Print this page
11210 libm should be cstyle(1ONBLD) clean

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/Q/erfl.c
          +++ new/usr/src/lib/libm/common/Q/erfl.c
↓ open down ↓ 14 lines elided ↑ open up ↑
  15   15   * If applicable, add the following below this CDDL HEADER, with the
  16   16   * fields enclosed by brackets "[]" replaced with your own identifying
  17   17   * information: Portions Copyright [yyyy] [name of copyright owner]
  18   18   *
  19   19   * CDDL HEADER END
  20   20   */
  21   21  
  22   22  /*
  23   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24   24   */
       25 +
  25   26  /*
  26   27   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   28   * Use is subject to license terms.
  28   29   */
  29   30  
  30   31  /*
  31   32   * long double function erf,erfc (long double x)
  32   33   * K.C. Ng, September, 1989.
  33   34   *                           x
  34   35   *                    2      |\
  35   36   *     erf(x)  =  ---------  | exp(-t*t)dt
  36      - *                 sqrt(pi) \|
       37 + *                 sqrt(pi) \|
  37   38   *                           0
  38   39   *
  39   40   *     erfc(x) =  1-erf(x)
  40   41   *
  41   42   * method:
  42      - *      Since erf(-x) = -erf(x), we assume x>=0.
       43 + *      Since erf(-x) = -erf(x), we assume x>=0.
  43   44   *      For x near 0, we have the expansion
  44   45   *
  45      - *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....).
       46 + *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....).
  46   47   *
  47      - *      Since 2/sqrt(pi) = 1.128379167095512573896158903121545171688,
       48 + *      Since 2/sqrt(pi) = 1.128379167095512573896158903121545171688,
  48   49   *      we use x + x*P(x^2) to approximate erf(x). This formula will
  49   50   *      guarantee the error less than one ulp where x is not too far
  50   51   *      away from 0. We note that erf(x)=x at x = 0.6174...... After
  51   52   *      some experiment, we choose the following approximation on
  52   53   *      interval [0,0.84375].
  53   54   *
  54   55   *      For x in [0,0.84375]
  55   56   *                 2                2        4               40
  56      - *         P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p20 * x  )
       57 + *         P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p20 * x  )
  57   58   *
  58   59   *         erf(x)  = x + x*P
  59      - *         erfc(x) = 1 - erf(x)           if x<=0.25
       60 + *         erfc(x) = 1 - erf(x)           if x<=0.25
  60   61   *                 = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]
  61   62   *      precision: |P(x^2)-(erf(x)-x)/x| <= 2**-122.50
  62   63   *
  63   64   *      For x in [0.84375,1.25], let s = x - 1, and
  64   65   *      c = 0.84506291151 rounded to single (24 bits)
  65   66   *         erf(x)  = c  + P1(s)/Q1(s)
  66   67   *         erfc(x) = (1-c)  - P1(s)/Q1(s)
  67   68   *      precision: |P1/Q1 - (erf(x)-c)| <= 2**-118.41
  68   69   *
  69   70   *
↓ open down ↓ 22 lines elided ↑ open up ↑
  92   93   *         erf(inf)  = 1
  93   94   *         erfc(inf) = 0
  94   95   */
  95   96  
  96   97  #pragma weak __erfl = erfl
  97   98  #pragma weak __erfcl = erfcl
  98   99  
  99  100  #include "libm.h"
 100  101  #include "longdouble.h"
 101  102  
 102      -static const long double
 103      -        tiny        = 1e-40L,
 104      -        nearunfl    = 1e-4000L,
 105      -        half        = 0.5L,
 106      -        one         = 1.0L,
 107      -        onehalf     = 1.5L,
 108      -        L16_3       = 16.0L/3.0L;
      103 +static const long double tiny = 1e-40L,
      104 +        nearunfl = 1e-4000L,
      105 +        half = 0.5L,
      106 +        one = 1.0L,
      107 +        onehalf = 1.5L,
      108 +        L16_3 = 16.0L / 3.0L;
      109 +
 109  110  /*
 110  111   * Coefficients for even polynomial P for erf(x)=x+x*P(x^2) on [0,0.84375]
 111  112   */
 112      -static const long double P[] = {        /* 21 coeffs */
 113      -   1.283791670955125738961589031215451715556e-0001L,
 114      -  -3.761263890318375246320529677071815594603e-0001L,
 115      -   1.128379167095512573896158903121205899135e-0001L,
 116      -  -2.686617064513125175943235483344625046092e-0002L,
 117      -   5.223977625442187842111846652980454568389e-0003L,
 118      -  -8.548327023450852832546626271083862724358e-0004L,
 119      -   1.205533298178966425102164715902231976672e-0004L,
 120      -  -1.492565035840625097674944905027897838996e-0005L,
 121      -   1.646211436588924733604648849172936692024e-0006L,
 122      -  -1.636584469123491976815834704799733514987e-0007L,
 123      -   1.480719281587897445302529007144770739305e-0008L,
 124      -  -1.229055530170782843046467986464722047175e-0009L,
 125      -   9.422759064320307357553954945760654341633e-0011L,
 126      -  -6.711366846653439036162105104991433380926e-0012L,
 127      -   4.463224090341893165100275380693843116240e-0013L,
 128      -  -2.783513452582658245422635662559779162312e-0014L,
 129      -   1.634227412586960195251346878863754661546e-0015L,
 130      -  -9.060782672889577722765711455623117802795e-0017L,
 131      -   4.741341801266246873412159213893613602354e-0018L,
 132      -  -2.272417596497826188374846636534317381203e-0019L,
 133      -   8.069088733716068462496835658928566920933e-0021L,
      113 +static const long double P[] = {        /* 21 coeffs */
      114 +        1.283791670955125738961589031215451715556e-0001L,
      115 +        -3.761263890318375246320529677071815594603e-0001L,
      116 +        1.128379167095512573896158903121205899135e-0001L,
      117 +        -2.686617064513125175943235483344625046092e-0002L,
      118 +        5.223977625442187842111846652980454568389e-0003L,
      119 +        -8.548327023450852832546626271083862724358e-0004L,
      120 +        1.205533298178966425102164715902231976672e-0004L,
      121 +        -1.492565035840625097674944905027897838996e-0005L,
      122 +        1.646211436588924733604648849172936692024e-0006L,
      123 +        -1.636584469123491976815834704799733514987e-0007L,
      124 +        1.480719281587897445302529007144770739305e-0008L,
      125 +        -1.229055530170782843046467986464722047175e-0009L,
      126 +        9.422759064320307357553954945760654341633e-0011L,
      127 +        -6.711366846653439036162105104991433380926e-0012L,
      128 +        4.463224090341893165100275380693843116240e-0013L,
      129 +        -2.783513452582658245422635662559779162312e-0014L,
      130 +        1.634227412586960195251346878863754661546e-0015L,
      131 +        -9.060782672889577722765711455623117802795e-0017L,
      132 +        4.741341801266246873412159213893613602354e-0018L,
      133 +        -2.272417596497826188374846636534317381203e-0019L,
      134 +        8.069088733716068462496835658928566920933e-0021L,
 134  135  };
 135  136  
 136  137  /*
 137  138   * Rational erf(x) = ((float)0.84506291151) + P1(x-1)/Q1(x-1) on [0.84375,1.25]
 138  139   */
 139      -static const long double C1   = (long double)((float)0.84506291151);
 140      -static const long double P1[] = {       /*  12 top coeffs */
 141      -  -2.362118560752659955654364917390741930316e-0003L,
 142      -   4.129623379624420034078926610650759979146e-0001L,
 143      -  -3.973857505403547283109417923182669976904e-0002L,
 144      -   4.357503184084022439763567513078036755183e-0002L,
 145      -   8.015593623388421371247676683754171456950e-0002L,
 146      -  -1.034459310403352486685467221776778474602e-0002L,
 147      -   5.671850295381046679675355719017720821383e-0003L,
 148      -   1.219262563232763998351452194968781174318e-0003L,
 149      -   5.390833481581033423020320734201065475098e-0004L,
 150      -  -1.978853912815115495053119023517805528300e-0004L,
 151      -   6.184234513953600118335017885706420552487e-0005L,
 152      -  -5.331802711697810861017518515816271808286e-0006L,
      140 +static const long double C1 = (long double)((float)0.84506291151);
      141 +static const long double P1[] = {       /*  12 top coeffs */
      142 +        -2.362118560752659955654364917390741930316e-0003L,
      143 +        4.129623379624420034078926610650759979146e-0001L,
      144 +        -3.973857505403547283109417923182669976904e-0002L,
      145 +        4.357503184084022439763567513078036755183e-0002L,
      146 +        8.015593623388421371247676683754171456950e-0002L,
      147 +        -1.034459310403352486685467221776778474602e-0002L,
      148 +        5.671850295381046679675355719017720821383e-0003L,
      149 +        1.219262563232763998351452194968781174318e-0003L,
      150 +        5.390833481581033423020320734201065475098e-0004L,
      151 +        -1.978853912815115495053119023517805528300e-0004L,
      152 +        6.184234513953600118335017885706420552487e-0005L,
      153 +        -5.331802711697810861017518515816271808286e-0006L,
 153  154  };
 154      -static const long double Q1[] = {       /*  12 bottom coeffs with leading 1.0 hidden */
 155      -   9.081506296064882195280178373107623196655e-0001L,
 156      -   6.821049531968204097604392183650687642520e-0001L,
 157      -   4.067869178233539502315055970743271822838e-0001L,
 158      -   1.702332233546316765818144723063881095577e-0001L,
 159      -   7.498098377690553934266423088708614219356e-0002L,
 160      -   2.050154396918178697056927234366372760310e-0002L,
 161      -   7.012988534031999899054782333851905939379e-0003L,
 162      -   1.149904787014400354649843451234570731076e-0003L,
 163      -   3.185620255011299476196039491205159718620e-0004L,
 164      -   1.273405072153008775426376193374105840517e-0005L,
 165      -   4.753866999959432971956781228148402971454e-0006L,
 166      -  -1.002287602111660026053981728549540200683e-0006L,
      155 +
      156 +static const long double Q1[] = { /* 12 bottom coeffs with leading 1.0 hidden */
      157 +        9.081506296064882195280178373107623196655e-0001L,
      158 +        6.821049531968204097604392183650687642520e-0001L,
      159 +        4.067869178233539502315055970743271822838e-0001L,
      160 +        1.702332233546316765818144723063881095577e-0001L,
      161 +        7.498098377690553934266423088708614219356e-0002L,
      162 +        2.050154396918178697056927234366372760310e-0002L,
      163 +        7.012988534031999899054782333851905939379e-0003L,
      164 +        1.149904787014400354649843451234570731076e-0003L,
      165 +        3.185620255011299476196039491205159718620e-0004L,
      166 +        1.273405072153008775426376193374105840517e-0005L,
      167 +        4.753866999959432971956781228148402971454e-0006L,
      168 +        -1.002287602111660026053981728549540200683e-0006L,
 167  169  };
      170 +
 168  171  /*
 169  172   * Rational erf(x) = ((float)0.95478588343) + P2(x-1.5)/Q2(x-1.5)
 170  173   * on [1.25,1.75]
 171  174   */
 172      -static const long double C2   = (long double)((float)0.95478588343);
 173      -static const long double P2[] = {       /*  12 top coeffs */
 174      -   1.131926304864446730135126164594785863512e-0002L,
 175      -   1.273617996967754151544330055186210322832e-0001L,
 176      -  -8.169980734667512519897816907190281143423e-0002L,
 177      -   9.512267486090321197833634271787944271746e-0002L,
 178      -  -2.394251569804872160005274999735914368170e-0002L,
 179      -   1.108768660227528667525252333184520222905e-0002L,
 180      -   3.527435492933902414662043314373277494221e-0004L,
 181      -   4.946116273341953463584319006669474625971e-0004L,
 182      -  -4.289851942513144714600285769022420962418e-0005L,
 183      -   8.304719841341952705874781636002085119978e-0005L,
 184      -  -1.040460226177309338781902252282849903189e-0005L,
 185      -   2.122913331584921470381327583672044434087e-0006L,
      175 +static const long double C2 = (long double)((float)0.95478588343);
      176 +static const long double P2[] = {       /*  12 top coeffs */
      177 +        1.131926304864446730135126164594785863512e-0002L,
      178 +        1.273617996967754151544330055186210322832e-0001L,
      179 +        -8.169980734667512519897816907190281143423e-0002L,
      180 +        9.512267486090321197833634271787944271746e-0002L,
      181 +        -2.394251569804872160005274999735914368170e-0002L,
      182 +        1.108768660227528667525252333184520222905e-0002L,
      183 +        3.527435492933902414662043314373277494221e-0004L,
      184 +        4.946116273341953463584319006669474625971e-0004L,
      185 +        -4.289851942513144714600285769022420962418e-0005L,
      186 +        8.304719841341952705874781636002085119978e-0005L,
      187 +        -1.040460226177309338781902252282849903189e-0005L,
      188 +        2.122913331584921470381327583672044434087e-0006L,
 186  189  };
 187      -static const long double Q2[] = {       /*  13 bottom coeffs with leading 1.0 hidden */
 188      -   7.448815737306992749168727691042003832150e-0001L,
 189      -   7.161813850236008294484744312430122188043e-0001L,
 190      -   3.603134756584225766144922727405641236121e-0001L,
 191      -   1.955811609133766478080550795194535852653e-0001L,
 192      -   7.253059963716225972479693813787810711233e-0002L,
 193      -   2.752391253757421424212770221541238324978e-0002L,
 194      -   7.677654852085240257439050673446546828005e-0003L,
 195      -   2.141102244555509687346497060326630061069e-0003L,
 196      -   4.342123013830957093949563339130674364271e-0004L,
 197      -   8.664587895570043348530991997272212150316e-0005L,
 198      -   1.109201582511752087060167429397033701988e-0005L,
 199      -   1.357834375781831062713347000030984364311e-0006L,
 200      -   4.957746280594384997273090385060680016451e-0008L,
      190 +
      191 +static const long double Q2[] = { /* 13 bottom coeffs with leading 1.0 hidden */
      192 +        7.448815737306992749168727691042003832150e-0001L,
      193 +        7.161813850236008294484744312430122188043e-0001L,
      194 +        3.603134756584225766144922727405641236121e-0001L,
      195 +        1.955811609133766478080550795194535852653e-0001L,
      196 +        7.253059963716225972479693813787810711233e-0002L,
      197 +        2.752391253757421424212770221541238324978e-0002L,
      198 +        7.677654852085240257439050673446546828005e-0003L,
      199 +        2.141102244555509687346497060326630061069e-0003L,
      200 +        4.342123013830957093949563339130674364271e-0004L,
      201 +        8.664587895570043348530991997272212150316e-0005L,
      202 +        1.109201582511752087060167429397033701988e-0005L,
      203 +        1.357834375781831062713347000030984364311e-0006L,
      204 +        4.957746280594384997273090385060680016451e-0008L,
 201  205  };
      206 +
 202  207  /*
 203  208   * erfc(x) = exp(-x*x)/x * R1(1/x)/S1(1/x) on [1.75, 16/3]
 204  209   */
 205      -static const long double R1[] = {       /*  14 top coeffs */
 206      -   4.630195122654315016370705767621550602948e+0006L,
 207      -   1.257949521746494830700654204488675713628e+0007L,
 208      -   1.704153822720260272814743497376181625707e+0007L,
 209      -   1.502600568706061872381577539537315739943e+0007L,
 210      -   9.543710793431995284827024445387333922861e+0006L,
 211      -   4.589344808584091011652238164935949522427e+0006L,
 212      -   1.714660662941745791190907071920671844289e+0006L,
 213      -   5.034802147768798894307672256192466283867e+0005L,
 214      -   1.162286400443554670553152110447126850725e+0005L,
 215      -   2.086643834548901681362757308058660399137e+0004L,
 216      -   2.839793161868140305907004392890348777338e+0003L,
 217      -   2.786687241658423601778258694498655680778e+0002L,
 218      -   1.779177837102695602425897452623985786464e+0001L,
 219      -   5.641895835477470769043614623819144434731e-0001L,
      210 +static const long double R1[] = {       /*  14 top coeffs */
      211 +        4.630195122654315016370705767621550602948e+0006L,
      212 +        1.257949521746494830700654204488675713628e+0007L,
      213 +        1.704153822720260272814743497376181625707e+0007L,
      214 +        1.502600568706061872381577539537315739943e+0007L,
      215 +        9.543710793431995284827024445387333922861e+0006L,
      216 +        4.589344808584091011652238164935949522427e+0006L,
      217 +        1.714660662941745791190907071920671844289e+0006L,
      218 +        5.034802147768798894307672256192466283867e+0005L,
      219 +        1.162286400443554670553152110447126850725e+0005L,
      220 +        2.086643834548901681362757308058660399137e+0004L,
      221 +        2.839793161868140305907004392890348777338e+0003L,
      222 +        2.786687241658423601778258694498655680778e+0002L,
      223 +        1.779177837102695602425897452623985786464e+0001L,
      224 +        5.641895835477470769043614623819144434731e-0001L,
 220  225  };
 221      -static const long double S1[] = {       /* 15 bottom coeffs with leading 1.0 hidden */
 222      -   4.630195122654331529595606896287596843110e+0006L,
 223      -   1.780411093345512024324781084220509055058e+0007L,
 224      -   3.250113097051800703707108623715776848283e+0007L,
 225      -   3.737857099176755050912193712123489115755e+0007L,
 226      -   3.029787497516578821459174055870781168593e+0007L,
 227      -   1.833850619965384765005769632103205777227e+0007L,
 228      -   8.562719999736915722210391222639186586498e+0006L,
 229      -   3.139684562074658971315545539760008136973e+0006L,
 230      -   9.106421313731384880027703627454366930945e+0005L,
 231      -   2.085108342384266508613267136003194920001e+0005L,
 232      -   3.723126272693120340730491416449539290600e+0004L,
 233      -   5.049169878567344046145695360784436929802e+0003L,
 234      -   4.944274532748010767670150730035392093899e+0002L,
 235      -   3.153510608818213929982940249162268971412e+0001L,
 236      -   1.0e00L,
      226 +
      227 +static const long double S1[] = { /* 15 bottom coeffs with leading 1.0 hidden */
      228 +        4.630195122654331529595606896287596843110e+0006L,
      229 +        1.780411093345512024324781084220509055058e+0007L,
      230 +        3.250113097051800703707108623715776848283e+0007L,
      231 +        3.737857099176755050912193712123489115755e+0007L,
      232 +        3.029787497516578821459174055870781168593e+0007L,
      233 +        1.833850619965384765005769632103205777227e+0007L,
      234 +        8.562719999736915722210391222639186586498e+0006L,
      235 +        3.139684562074658971315545539760008136973e+0006L,
      236 +        9.106421313731384880027703627454366930945e+0005L,
      237 +        2.085108342384266508613267136003194920001e+0005L,
      238 +        3.723126272693120340730491416449539290600e+0004L,
      239 +        5.049169878567344046145695360784436929802e+0003L,
      240 +        4.944274532748010767670150730035392093899e+0002L,
      241 +        3.153510608818213929982940249162268971412e+0001L,
      242 +        1.0e00L,
 237  243  };
 238  244  
 239  245  /*
 240  246   * erfc(x) = exp(-x*x)/x * R2(1/x)/S2(1/x) on [16/3, 107]
 241  247   */
 242      -static const long double R2[] = {       /*  15 top coeffs in reverse order!!*/
 243      -   2.447288012254302966796326587537136931669e+0005L,
 244      -   8.768592567189861896653369912716538739016e+0005L,
 245      -   1.552293152581780065761497908005779524953e+0006L,
 246      -   1.792075924835942935864231657504259926729e+0006L,
 247      -   1.504001463155897344947500222052694835875e+0006L,
 248      -   9.699485556326891411801230186016013019935e+0005L,
 249      -   4.961449933661807969863435013364796037700e+0005L,
 250      -   2.048726544693474028061176764716228273791e+0005L,
 251      -   6.891532964330949722479061090551896886635e+0004L,
 252      -   1.888014709010307507771964047905823237985e+0004L,
 253      -   4.189692064988957745054734809642495644502e+0003L,
 254      -   7.362346487427048068212968889642741734621e+0002L,
 255      -   9.980359714211411423007641056580813116207e+0001L,
 256      -   9.426910895135379181107191962193485174159e+0000L,
 257      -   5.641895835477562869480794515623601280429e-0001L,
      248 +static const long double R2[] = { /*  15 top coeffs in reverse order!! */
      249 +        2.447288012254302966796326587537136931669e+0005L,
      250 +        8.768592567189861896653369912716538739016e+0005L,
      251 +        1.552293152581780065761497908005779524953e+0006L,
      252 +        1.792075924835942935864231657504259926729e+0006L,
      253 +        1.504001463155897344947500222052694835875e+0006L,
      254 +        9.699485556326891411801230186016013019935e+0005L,
      255 +        4.961449933661807969863435013364796037700e+0005L,
      256 +        2.048726544693474028061176764716228273791e+0005L,
      257 +        6.891532964330949722479061090551896886635e+0004L,
      258 +        1.888014709010307507771964047905823237985e+0004L,
      259 +        4.189692064988957745054734809642495644502e+0003L,
      260 +        7.362346487427048068212968889642741734621e+0002L,
      261 +        9.980359714211411423007641056580813116207e+0001L,
      262 +        9.426910895135379181107191962193485174159e+0000L,
      263 +        5.641895835477562869480794515623601280429e-0001L,
 258  264  };
 259      -static const long double S2[] = {       /* 16 coefficients */
 260      -   2.447282203601902971246004716790604686880e+0005L,
 261      -   1.153009852759385309367759460934808489833e+0006L,
 262      -   2.608580649612639131548966265078663384849e+0006L,
 263      -   3.766673917346623308850202792390569025740e+0006L,
 264      -   3.890566255138383910789924920541335370691e+0006L,
 265      -   3.052882073900746207613166259994150527732e+0006L,
 266      -   1.885574519970380988460241047248519418407e+0006L,
 267      -   9.369722034759943185851450846811445012922e+0005L,
 268      -   3.792278350536686111444869752624492443659e+0005L,
 269      -   1.257750606950115799965366001773094058720e+0005L,
 270      -   3.410830600242369370645608634643620355058e+0004L,
 271      -   7.513984469742343134851326863175067271240e+0003L,
 272      -   1.313296320593190002554779998138695507840e+0003L,
 273      -   1.773972700887629157006326333696896516769e+0002L,
 274      -   1.670876451822586800422009013880457094162e+0001L,
 275      -   1.000L,
      265 +
      266 +static const long double S2[] = {       /* 16 coefficients */
      267 +        2.447282203601902971246004716790604686880e+0005L,
      268 +        1.153009852759385309367759460934808489833e+0006L,
      269 +        2.608580649612639131548966265078663384849e+0006L,
      270 +        3.766673917346623308850202792390569025740e+0006L,
      271 +        3.890566255138383910789924920541335370691e+0006L,
      272 +        3.052882073900746207613166259994150527732e+0006L,
      273 +        1.885574519970380988460241047248519418407e+0006L,
      274 +        9.369722034759943185851450846811445012922e+0005L,
      275 +        3.792278350536686111444869752624492443659e+0005L,
      276 +        1.257750606950115799965366001773094058720e+0005L,
      277 +        3.410830600242369370645608634643620355058e+0004L,
      278 +        7.513984469742343134851326863175067271240e+0003L,
      279 +        1.313296320593190002554779998138695507840e+0003L,
      280 +        1.773972700887629157006326333696896516769e+0002L,
      281 +        1.670876451822586800422009013880457094162e+0001L,
      282 +        1.000L,
 276  283  };
 277  284  
 278      -long double erfl(x)
 279      -long double x;
      285 +long double
      286 +erfl(long double x)
 280  287  {
 281      -        long double s,y,t;
      288 +        long double s, y, t;
 282  289  
 283  290          if (!finitel(x)) {
 284      -            if (x != x) return x+x;     /* NaN */
 285      -            return copysignl(one,x);    /* return +-1.0 is x=Inf */
      291 +                if (x != x)
      292 +                        return (x + x);         /* NaN */
      293 +
      294 +                return (copysignl(one, x));     /* return +-1.0 is x=Inf */
 286  295          }
 287  296  
 288  297          y = fabsl(x);
      298 +
 289  299          if (y <= 0.84375L) {
 290      -            if (y<=tiny) return x+P[0]*x;
 291      -            s = y*y;
 292      -            t = __poly_libmq(s,21,P);
 293      -            return  x+x*t;
      300 +                if (y <= tiny)
      301 +                        return (x + P[0] * x);
      302 +
      303 +                s = y * y;
      304 +                t = __poly_libmq(s, 21, P);
      305 +                return (x + x * t);
 294  306          }
 295      -        if (y<=1.25L) {
 296      -            s = y-one;
 297      -            t = C1+__poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
 298      -            return (signbitl(x))? -t: t;
 299      -        } else if (y<=1.75L) {
 300      -            s = y-onehalf;
 301      -            t = C2+__poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
 302      -            return (signbitl(x))? -t: t;
      307 +
      308 +        if (y <= 1.25L) {
      309 +                s = y - one;
      310 +                t = C1 + __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s,
      311 +                    12, Q1));
      312 +                return ((signbitl(x)) ? -t : t);
      313 +        } else if (y <= 1.75L) {
      314 +                s = y - onehalf;
      315 +                t = C2 + __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s,
      316 +                    13, Q2));
      317 +                return ((signbitl(x)) ? -t : t);
 303  318          }
 304      -        if (y<=9.0L) t = erfcl(y); else t = tiny;
 305      -        return (signbitl(x))? t-one: one-t;
      319 +
      320 +        if (y <= 9.0L)
      321 +                t = erfcl(y);
      322 +        else
      323 +                t = tiny;
      324 +
      325 +        return ((signbitl(x)) ? t - one : one - t);
 306  326  }
 307  327  
 308      -long double erfcl(x)
 309      -long double x;
      328 +long double
      329 +erfcl(long double x)
 310  330  {
 311      -        long double s,y,t;
      331 +        long double s, y, t;
 312  332  
 313  333          if (!finitel(x)) {
 314      -            if (x != x) return x+x;     /* NaN */
 315      -            /* return 2.0 if x= -inf; 0.0 if x= +inf */
 316      -            if (x < 0.0L) return 2.0L; else return 0.0L;
      334 +                if (x != x)
      335 +                        return (x + x);         /* NaN */
      336 +
      337 +                /* return 2.0 if x= -inf; 0.0 if x= +inf */
      338 +                if (x < 0.0L)
      339 +                        return (2.0L);
      340 +                else
      341 +                        return (0.0L);
 317  342          }
 318  343  
 319  344          if (x <= 0.84375L) {
 320      -            if (x<=0.25) return one-erfl(x);
 321      -            s = x*x;
 322      -            t = half-x;
 323      -            t = t - x*__poly_libmq(s,21,P);
 324      -            return  half+t;
      345 +                if (x <= 0.25)
      346 +                        return (one - erfl(x));
      347 +
      348 +                s = x * x;
      349 +                t = half - x;
      350 +                t = t - x * __poly_libmq(s, 21, P);
      351 +                return (half + t);
 325  352          }
 326      -        if (x<=1.25L) {
 327      -            s = x-one;
 328      -            t = one-C1;
 329      -            return t - __poly_libmq(s,12,P1)/(one+s*__poly_libmq(s,12,Q1));
 330      -        } else if (x<=1.75L) {
 331      -            s = x-onehalf;
 332      -            t = one-C2;
 333      -            return t - __poly_libmq(s,12,P2)/(one+s*__poly_libmq(s,13,Q2));
      353 +
      354 +        if (x <= 1.25L) {
      355 +                s = x - one;
      356 +                t = one - C1;
      357 +                return (t - __poly_libmq(s, 12, P1) / (one + s * __poly_libmq(s,
      358 +                    12, Q1)));
      359 +        } else if (x <= 1.75L) {
      360 +                s = x - onehalf;
      361 +                t = one - C2;
      362 +                return (t - __poly_libmq(s, 12, P2) / (one + s * __poly_libmq(s,
      363 +                    13, Q2)));
 334  364          }
 335      -        if (x>=107.0L) return nearunfl*nearunfl;                /* underflow */
 336      -        else if (x >= L16_3) {
 337      -            y = __poly_libmq(x,15,R2);
 338      -            t = y/__poly_libmq(x,16,S2);
      365 +
      366 +        if (x >= 107.0L) {
      367 +                return (nearunfl * nearunfl);   /* underflow */
      368 +        } else if (x >= L16_3) {
      369 +                y = __poly_libmq(x, 15, R2);
      370 +                t = y / __poly_libmq(x, 16, S2);
 339  371          } else {
 340      -            y = __poly_libmq(x,14,R1);
 341      -            t = y/__poly_libmq(x,15,S1);
      372 +                y = __poly_libmq(x, 14, R1);
      373 +                t = y / __poly_libmq(x, 15, S1);
 342  374          }
      375 +
 343  376          /*
 344  377           * Note that exp(-x*x+d) = exp(-x*x)*exp(d), so to compute
 345  378           * exp(-x*x) with a small relative error, we need to compute
 346  379           * -x*x with a small absolute error.  To this end, we set y
 347  380           * equal to the leading part of x but with enough trailing
 348  381           * zeros that y*y can be computed exactly and we rewrite x*x
 349  382           * as y*y + (x-y)*(x+y), distributing the latter expression
 350  383           * across the exponential.
 351  384           *
 352  385           * We could construct y in a portable way by setting
 353  386           *
 354  387           *   int i = (int)(x * ptwo);
 355  388           *   y = (long double)i * 1/ptwo;
 356  389           *
 357  390           * where ptwo is some power of two large enough to make x-y
 358  391           * small but not so large that the conversion to int overflows.
 359  392           * When long double arithmetic is slow, however, the following
 360  393           * non-portable code is preferable.
 361  394           */
 362  395          y = x;
 363      -        *(2+(int*)&y) = *(3+(int*)&y) = 0;
 364      -        t *= expl(-y*y)*expl(-(x-y)*(x+y));
 365      -        return  t;
      396 +        *(2 + (int *)&y) = *(3 + (int *)&y) = 0;
      397 +        t *= expl(-y * y) * expl(-(x - y) * (x + y));
      398 +        return (t);
 366  399  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX