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

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/C/erf.c
          +++ new/usr/src/lib/libm/common/C/erf.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  #pragma weak __erf = erf
  31   32  #pragma weak __erfc = erfc
  32   33  
  33      -/* INDENT OFF */
       34 +
  34   35  /*
  35   36   * double erf(double x)
  36   37   * double erfc(double x)
  37   38   *                           x
  38   39   *                    2      |\
  39   40   *     erf(x)  =  ---------  | exp(-t*t)dt
  40   41   *                 sqrt(pi) \|
  41   42   *                           0
  42   43   *
  43   44   *     erfc(x) =  1-erf(x)
↓ open down ↓ 16 lines elided ↑ open up ↑
  60   61   *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
  61   62   *         and that
  62   63   *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
  63   64   *         is close to one. The interval is chosen because the fix
  64   65   *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
  65   66   *         near 0.6174), and by some experiment, 0.84375 is chosen to
  66   67   *         guarantee the error is less than one ulp for erf.
  67   68   *
  68   69   *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
  69   70   *         c = 0.84506291151 rounded to single (24 bits)
  70      - *              erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
  71      - *              erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
       71 + *              erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
       72 + *              erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
  72   73   *                        1+(c+P1(s)/Q1(s))    if x < 0
  73      - *              |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
       74 + *              |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
  74   75   *         Remark: here we use the taylor series expansion at x=1.
  75   76   *              erf(1+s) = erf(1) + s*Poly(s)
  76   77   *                       = 0.845.. + P1(s)/Q1(s)
  77   78   *         That is, we use rational approximation to approximate
  78   79   *                      erf(1+s) - (c = (single)0.84506291151)
  79   80   *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
  80   81   *         where
  81   82   *              P1(s) = degree 6 poly in s
  82   83   *              Q1(s) = degree 6 poly in s
  83   84   *
  84   85   *      3. For x in [1.25,1/0.35(~2.857143)],
  85      - *              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
  86      - *              erf(x)  = 1 - erfc(x)
       86 + *              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
       87 + *              erf(x)  = 1 - erfc(x)
  87   88   *         where
  88   89   *              R1(z) = degree 7 poly in z, (z=1/x^2)
  89   90   *              S1(z) = degree 8 poly in z
  90   91   *
  91   92   *      4. For x in [1/0.35,28]
  92      - *              erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
       93 + *              erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
  93   94   *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
  94   95   *                      = 2.0 - tiny            (if x <= -6)
  95      - *              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
  96      - *              erf(x)  = sign(x)*(1.0 - tiny)
       96 + *              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
       97 + *              erf(x)  = sign(x)*(1.0 - tiny)
  97   98   *         where
  98   99   *              R2(z) = degree 6 poly in z, (z=1/x^2)
  99  100   *              S2(z) = degree 7 poly in z
 100  101   *
 101  102   *      Note1:
 102  103   *         To compute exp(-x*x-0.5625+R/S), let s be a single
 103  104   *         precision number and s := x; then
 104  105   *              -x*x = -s*s + (s-x)*(s+x)
 105  106   *              exp(-x*x-0.5626+R/S) =
 106  107   *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
 107  108   *      Note2:
 108  109   *         Here 4 and 5 make use of the asymptotic series
 109  110   *                        exp(-x*x)
 110  111   *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
 111  112   *                        x*sqrt(pi)
 112  113   *         We use rational approximation to approximate
 113      - *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
      114 + *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
 114  115   *         Here is the error bound for R1/S1 and R2/S2
 115      - *              |R1/S1 - f(x)|  < 2**(-62.57)
 116      - *              |R2/S2 - f(x)|  < 2**(-61.52)
      116 + *              |R1/S1 - f(x)|  < 2**(-62.57)
      117 + *              |R2/S2 - f(x)|  < 2**(-61.52)
 117  118   *
 118  119   *      5. For inf > x >= 28
 119      - *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
 120      - *              erfc(x) = tiny*tiny (raise underflow) if x > 0
      120 + *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
      121 + *              erfc(x) = tiny*tiny (raise underflow) if x > 0
 121  122   *                      = 2 - tiny if x<0
 122  123   *
 123  124   *      7. Special case:
 124      - *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
 125      - *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
 126      - *      erfc/erf(NaN) is NaN
      125 + *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
      126 + *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
      127 + *      erfc/erf(NaN) is NaN
 127  128   */
 128      -/* INDENT ON */
 129  129  
 130  130  #include "libm_macros.h"
 131  131  #include <math.h>
 132  132  
 133  133  static const double xxx[] = {
 134      -/* tiny */      1e-300,
 135      -/* half */      5.00000000000000000000e-01,     /* 3FE00000, 00000000 */
 136      -/* one */       1.00000000000000000000e+00,     /* 3FF00000, 00000000 */
 137      -/* two */       2.00000000000000000000e+00,     /* 40000000, 00000000 */
 138      -/* erx */       8.45062911510467529297e-01,     /* 3FEB0AC1, 60000000 */
      134 +/* tiny */
      135 +        1e-300,
      136 +/* half */ 5.00000000000000000000e-01,  /* 3FE00000, 00000000 */
      137 +/* one */ 1.00000000000000000000e+00,   /* 3FF00000, 00000000 */
      138 +/* two */ 2.00000000000000000000e+00,   /* 40000000, 00000000 */
      139 +/* erx */ 8.45062911510467529297e-01,   /* 3FEB0AC1, 60000000 */
      140 +
 139  141  /*
 140  142   * Coefficients for approximation to  erf on [0,0.84375]
 141  143   */
 142      -/* efx */        1.28379167095512586316e-01,    /* 3FC06EBA, 8214DB69 */
 143      -/* efx8 */       1.02703333676410069053e+00,    /* 3FF06EBA, 8214DB69 */
 144      -/* pp0 */        1.28379167095512558561e-01,    /* 3FC06EBA, 8214DB68 */
 145      -/* pp1 */       -3.25042107247001499370e-01,    /* BFD4CD7D, 691CB913 */
 146      -/* pp2 */       -2.84817495755985104766e-02,    /* BF9D2A51, DBD7194F */
 147      -/* pp3 */       -5.77027029648944159157e-03,    /* BF77A291, 236668E4 */
 148      -/* pp4 */       -2.37630166566501626084e-05,    /* BEF8EAD6, 120016AC */
 149      -/* qq1 */        3.97917223959155352819e-01,    /* 3FD97779, CDDADC09 */
 150      -/* qq2 */        6.50222499887672944485e-02,    /* 3FB0A54C, 5536CEBA */
 151      -/* qq3 */        5.08130628187576562776e-03,    /* 3F74D022, C4D36B0F */
 152      -/* qq4 */        1.32494738004321644526e-04,    /* 3F215DC9, 221C1A10 */
 153      -/* qq5 */       -3.96022827877536812320e-06,    /* BED09C43, 42A26120 */
      144 +/* efx */ 1.28379167095512586316e-01,   /* 3FC06EBA, 8214DB69 */
      145 +/* efx8 */ 1.02703333676410069053e+00,  /* 3FF06EBA, 8214DB69 */
      146 +/* pp0 */ 1.28379167095512558561e-01,   /* 3FC06EBA, 8214DB68 */
      147 +/* pp1 */ -3.25042107247001499370e-01,  /* BFD4CD7D, 691CB913 */
      148 +/* pp2 */ -2.84817495755985104766e-02,  /* BF9D2A51, DBD7194F */
      149 +/* pp3 */ -5.77027029648944159157e-03,  /* BF77A291, 236668E4 */
      150 +/* pp4 */ -2.37630166566501626084e-05,  /* BEF8EAD6, 120016AC */
      151 +/* qq1 */ 3.97917223959155352819e-01,   /* 3FD97779, CDDADC09 */
      152 +/* qq2 */ 6.50222499887672944485e-02,   /* 3FB0A54C, 5536CEBA */
      153 +/* qq3 */ 5.08130628187576562776e-03,   /* 3F74D022, C4D36B0F */
      154 +/* qq4 */ 1.32494738004321644526e-04,   /* 3F215DC9, 221C1A10 */
      155 +/* qq5 */ -3.96022827877536812320e-06,  /* BED09C43, 42A26120 */
      156 +
 154  157  /*
 155  158   * Coefficients for approximation to  erf  in [0.84375,1.25]
 156  159   */
 157      -/* pa0 */       -2.36211856075265944077e-03,    /* BF6359B8, BEF77538 */
 158      -/* pa1 */        4.14856118683748331666e-01,    /* 3FDA8D00, AD92B34D */
 159      -/* pa2 */       -3.72207876035701323847e-01,    /* BFD7D240, FBB8C3F1 */
 160      -/* pa3 */        3.18346619901161753674e-01,    /* 3FD45FCA, 805120E4 */
 161      -/* pa4 */       -1.10894694282396677476e-01,    /* BFBC6398, 3D3E28EC */
 162      -/* pa5 */        3.54783043256182359371e-02,    /* 3FA22A36, 599795EB */
 163      -/* pa6 */       -2.16637559486879084300e-03,    /* BF61BF38, 0A96073F */
 164      -/* qa1 */        1.06420880400844228286e-01,    /* 3FBB3E66, 18EEE323 */
 165      -/* qa2 */        5.40397917702171048937e-01,    /* 3FE14AF0, 92EB6F33 */
 166      -/* qa3 */        7.18286544141962662868e-02,    /* 3FB2635C, D99FE9A7 */
 167      -/* qa4 */        1.26171219808761642112e-01,    /* 3FC02660, E763351F */
 168      -/* qa5 */        1.36370839120290507362e-02,    /* 3F8BEDC2, 6B51DD1C */
 169      -/* qa6 */        1.19844998467991074170e-02,    /* 3F888B54, 5735151D */
      160 +/* pa0 */ -2.36211856075265944077e-03,  /* BF6359B8, BEF77538 */
      161 +/* pa1 */ 4.14856118683748331666e-01,   /* 3FDA8D00, AD92B34D */
      162 +/* pa2 */ -3.72207876035701323847e-01,  /* BFD7D240, FBB8C3F1 */
      163 +/* pa3 */ 3.18346619901161753674e-01,   /* 3FD45FCA, 805120E4 */
      164 +/* pa4 */ -1.10894694282396677476e-01,  /* BFBC6398, 3D3E28EC */
      165 +/* pa5 */ 3.54783043256182359371e-02,   /* 3FA22A36, 599795EB */
      166 +/* pa6 */ -2.16637559486879084300e-03,  /* BF61BF38, 0A96073F */
      167 +/* qa1 */ 1.06420880400844228286e-01,   /* 3FBB3E66, 18EEE323 */
      168 +/* qa2 */ 5.40397917702171048937e-01,   /* 3FE14AF0, 92EB6F33 */
      169 +/* qa3 */ 7.18286544141962662868e-02,   /* 3FB2635C, D99FE9A7 */
      170 +/* qa4 */ 1.26171219808761642112e-01,   /* 3FC02660, E763351F */
      171 +/* qa5 */ 1.36370839120290507362e-02,   /* 3F8BEDC2, 6B51DD1C */
      172 +/* qa6 */ 1.19844998467991074170e-02,   /* 3F888B54, 5735151D */
      173 +
 170  174  /*
 171  175   * Coefficients for approximation to  erfc in [1.25,1/0.35]
 172  176   */
 173      -/* ra0 */       -9.86494403484714822705e-03,    /* BF843412, 600D6435 */
 174      -/* ra1 */       -6.93858572707181764372e-01,    /* BFE63416, E4BA7360 */
 175      -/* ra2 */       -1.05586262253232909814e+01,    /* C0251E04, 41B0E726 */
 176      -/* ra3 */       -6.23753324503260060396e+01,    /* C04F300A, E4CBA38D */
 177      -/* ra4 */       -1.62396669462573470355e+02,    /* C0644CB1, 84282266 */
 178      -/* ra5 */       -1.84605092906711035994e+02,    /* C067135C, EBCCABB2 */
 179      -/* ra6 */       -8.12874355063065934246e+01,    /* C0545265, 57E4D2F2 */
 180      -/* ra7 */       -9.81432934416914548592e+00,    /* C023A0EF, C69AC25C */
 181      -/* sa1 */        1.96512716674392571292e+01,    /* 4033A6B9, BD707687 */
 182      -/* sa2 */        1.37657754143519042600e+02,    /* 4061350C, 526AE721 */
 183      -/* sa3 */        4.34565877475229228821e+02,    /* 407B290D, D58A1A71 */
 184      -/* sa4 */        6.45387271733267880336e+02,    /* 40842B19, 21EC2868 */
 185      -/* sa5 */        4.29008140027567833386e+02,    /* 407AD021, 57700314 */
 186      -/* sa6 */        1.08635005541779435134e+02,    /* 405B28A3, EE48AE2C */
 187      -/* sa7 */        6.57024977031928170135e+00,    /* 401A47EF, 8E484A93 */
 188      -/* sa8 */       -6.04244152148580987438e-02,    /* BFAEEFF2, EE749A62 */
      177 +/* ra0 */ -9.86494403484714822705e-03,  /* BF843412, 600D6435 */
      178 +/* ra1 */ -6.93858572707181764372e-01,  /* BFE63416, E4BA7360 */
      179 +/* ra2 */ -1.05586262253232909814e+01,  /* C0251E04, 41B0E726 */
      180 +/* ra3 */ -6.23753324503260060396e+01,  /* C04F300A, E4CBA38D */
      181 +/* ra4 */ -1.62396669462573470355e+02,  /* C0644CB1, 84282266 */
      182 +/* ra5 */ -1.84605092906711035994e+02,  /* C067135C, EBCCABB2 */
      183 +/* ra6 */ -8.12874355063065934246e+01,  /* C0545265, 57E4D2F2 */
      184 +/* ra7 */ -9.81432934416914548592e+00,  /* C023A0EF, C69AC25C */
      185 +/* sa1 */ 1.96512716674392571292e+01,   /* 4033A6B9, BD707687 */
      186 +/* sa2 */ 1.37657754143519042600e+02,   /* 4061350C, 526AE721 */
      187 +/* sa3 */ 4.34565877475229228821e+02,   /* 407B290D, D58A1A71 */
      188 +/* sa4 */ 6.45387271733267880336e+02,   /* 40842B19, 21EC2868 */
      189 +/* sa5 */ 4.29008140027567833386e+02,   /* 407AD021, 57700314 */
      190 +/* sa6 */ 1.08635005541779435134e+02,   /* 405B28A3, EE48AE2C */
      191 +/* sa7 */ 6.57024977031928170135e+00,   /* 401A47EF, 8E484A93 */
      192 +/* sa8 */ -6.04244152148580987438e-02,  /* BFAEEFF2, EE749A62 */
      193 +
 189  194  /*
 190  195   * Coefficients for approximation to  erfc in [1/.35,28]
 191  196   */
 192      -/* rb0 */       -9.86494292470009928597e-03,    /* BF843412, 39E86F4A */
 193      -/* rb1 */       -7.99283237680523006574e-01,    /* BFE993BA, 70C285DE */
 194      -/* rb2 */       -1.77579549177547519889e+01,    /* C031C209, 555F995A */
 195      -/* rb3 */       -1.60636384855821916062e+02,    /* C064145D, 43C5ED98 */
 196      -/* rb4 */       -6.37566443368389627722e+02,    /* C083EC88, 1375F228 */
 197      -/* rb5 */       -1.02509513161107724954e+03,    /* C0900461, 6A2E5992 */
 198      -/* rb6 */       -4.83519191608651397019e+02,    /* C07E384E, 9BDC383F */
 199      -/* sb1 */        3.03380607434824582924e+01,    /* 403E568B, 261D5190 */
 200      -/* sb2 */        3.25792512996573918826e+02,    /* 40745CAE, 221B9F0A */
 201      -/* sb3 */        1.53672958608443695994e+03,    /* 409802EB, 189D5118 */
 202      -/* sb4 */        3.19985821950859553908e+03,    /* 40A8FFB7, 688C246A */
 203      -/* sb5 */        2.55305040643316442583e+03,    /* 40A3F219, CEDF3BE6 */
 204      -/* sb6 */        4.74528541206955367215e+02,    /* 407DA874, E79FE763 */
 205      -/* sb7 */       -2.24409524465858183362e+01     /* C03670E2, 42712D62 */
      197 +/* rb0 */ -9.86494292470009928597e-03,  /* BF843412, 39E86F4A */
      198 +/* rb1 */ -7.99283237680523006574e-01,  /* BFE993BA, 70C285DE */
      199 +/* rb2 */ -1.77579549177547519889e+01,  /* C031C209, 555F995A */
      200 +/* rb3 */ -1.60636384855821916062e+02,  /* C064145D, 43C5ED98 */
      201 +/* rb4 */ -6.37566443368389627722e+02,  /* C083EC88, 1375F228 */
      202 +/* rb5 */ -1.02509513161107724954e+03,  /* C0900461, 6A2E5992 */
      203 +/* rb6 */ -4.83519191608651397019e+02,  /* C07E384E, 9BDC383F */
      204 +/* sb1 */ 3.03380607434824582924e+01,   /* 403E568B, 261D5190 */
      205 +/* sb2 */ 3.25792512996573918826e+02,   /* 40745CAE, 221B9F0A */
      206 +/* sb3 */ 1.53672958608443695994e+03,   /* 409802EB, 189D5118 */
      207 +/* sb4 */ 3.19985821950859553908e+03,   /* 40A8FFB7, 688C246A */
      208 +/* sb5 */ 2.55305040643316442583e+03,   /* 40A3F219, CEDF3BE6 */
      209 +/* sb6 */ 4.74528541206955367215e+02,   /* 407DA874, E79FE763 */
      210 +/* sb7 */ -2.24409524465858183362e+01   /* C03670E2, 42712D62 */
 206  211  };
 207  212  
 208      -#define tiny    xxx[0]
 209      -#define half    xxx[1]
 210      -#define one     xxx[2]
 211      -#define two     xxx[3]
 212      -#define erx     xxx[4]
      213 +#define tiny            xxx[0]
      214 +#define half            xxx[1]
      215 +#define one             xxx[2]
      216 +#define two             xxx[3]
      217 +#define erx             xxx[4]
      218 +
 213  219  /*
 214  220   * Coefficients for approximation to  erf on [0,0.84375]
 215  221   */
 216      -#define efx     xxx[5]
 217      -#define efx8    xxx[6]
 218      -#define pp0     xxx[7]
 219      -#define pp1     xxx[8]
 220      -#define pp2     xxx[9]
 221      -#define pp3     xxx[10]
 222      -#define pp4     xxx[11]
 223      -#define qq1     xxx[12]
 224      -#define qq2     xxx[13]
 225      -#define qq3     xxx[14]
 226      -#define qq4     xxx[15]
 227      -#define qq5     xxx[16]
      222 +#define efx             xxx[5]
      223 +#define efx8            xxx[6]
      224 +#define pp0             xxx[7]
      225 +#define pp1             xxx[8]
      226 +#define pp2             xxx[9]
      227 +#define pp3             xxx[10]
      228 +#define pp4             xxx[11]
      229 +#define qq1             xxx[12]
      230 +#define qq2             xxx[13]
      231 +#define qq3             xxx[14]
      232 +#define qq4             xxx[15]
      233 +#define qq5             xxx[16]
      234 +
 228  235  /*
 229  236   * Coefficients for approximation to  erf  in [0.84375,1.25]
 230  237   */
 231      -#define pa0     xxx[17]
 232      -#define pa1     xxx[18]
 233      -#define pa2     xxx[19]
 234      -#define pa3     xxx[20]
 235      -#define pa4     xxx[21]
 236      -#define pa5     xxx[22]
 237      -#define pa6     xxx[23]
 238      -#define qa1     xxx[24]
 239      -#define qa2     xxx[25]
 240      -#define qa3     xxx[26]
 241      -#define qa4     xxx[27]
 242      -#define qa5     xxx[28]
 243      -#define qa6     xxx[29]
      238 +#define pa0             xxx[17]
      239 +#define pa1             xxx[18]
      240 +#define pa2             xxx[19]
      241 +#define pa3             xxx[20]
      242 +#define pa4             xxx[21]
      243 +#define pa5             xxx[22]
      244 +#define pa6             xxx[23]
      245 +#define qa1             xxx[24]
      246 +#define qa2             xxx[25]
      247 +#define qa3             xxx[26]
      248 +#define qa4             xxx[27]
      249 +#define qa5             xxx[28]
      250 +#define qa6             xxx[29]
      251 +
 244  252  /*
 245  253   * Coefficients for approximation to  erfc in [1.25,1/0.35]
 246  254   */
 247      -#define ra0     xxx[30]
 248      -#define ra1     xxx[31]
 249      -#define ra2     xxx[32]
 250      -#define ra3     xxx[33]
 251      -#define ra4     xxx[34]
 252      -#define ra5     xxx[35]
 253      -#define ra6     xxx[36]
 254      -#define ra7     xxx[37]
 255      -#define sa1     xxx[38]
 256      -#define sa2     xxx[39]
 257      -#define sa3     xxx[40]
 258      -#define sa4     xxx[41]
 259      -#define sa5     xxx[42]
 260      -#define sa6     xxx[43]
 261      -#define sa7     xxx[44]
 262      -#define sa8     xxx[45]
      255 +#define ra0             xxx[30]
      256 +#define ra1             xxx[31]
      257 +#define ra2             xxx[32]
      258 +#define ra3             xxx[33]
      259 +#define ra4             xxx[34]
      260 +#define ra5             xxx[35]
      261 +#define ra6             xxx[36]
      262 +#define ra7             xxx[37]
      263 +#define sa1             xxx[38]
      264 +#define sa2             xxx[39]
      265 +#define sa3             xxx[40]
      266 +#define sa4             xxx[41]
      267 +#define sa5             xxx[42]
      268 +#define sa6             xxx[43]
      269 +#define sa7             xxx[44]
      270 +#define sa8             xxx[45]
      271 +
 263  272  /*
 264  273   * Coefficients for approximation to  erfc in [1/.35,28]
 265  274   */
 266      -#define rb0     xxx[46]
 267      -#define rb1     xxx[47]
 268      -#define rb2     xxx[48]
 269      -#define rb3     xxx[49]
 270      -#define rb4     xxx[50]
 271      -#define rb5     xxx[51]
 272      -#define rb6     xxx[52]
 273      -#define sb1     xxx[53]
 274      -#define sb2     xxx[54]
 275      -#define sb3     xxx[55]
 276      -#define sb4     xxx[56]
 277      -#define sb5     xxx[57]
 278      -#define sb6     xxx[58]
 279      -#define sb7     xxx[59]
      275 +#define rb0             xxx[46]
      276 +#define rb1             xxx[47]
      277 +#define rb2             xxx[48]
      278 +#define rb3             xxx[49]
      279 +#define rb4             xxx[50]
      280 +#define rb5             xxx[51]
      281 +#define rb6             xxx[52]
      282 +#define sb1             xxx[53]
      283 +#define sb2             xxx[54]
      284 +#define sb3             xxx[55]
      285 +#define sb4             xxx[56]
      286 +#define sb5             xxx[57]
      287 +#define sb6             xxx[58]
      288 +#define sb7             xxx[59]
 280  289  
 281  290  double
 282      -erf(double x) {
      291 +erf(double x)
      292 +{
 283  293          int hx, ix, i;
 284  294          double R, S, P, Q, s, y, z, r;
 285  295  
 286      -        hx = ((int *) &x)[HIWORD];
      296 +        hx = ((int *)&x)[HIWORD];
 287  297          ix = hx & 0x7fffffff;
 288      -        if (ix >= 0x7ff00000) { /* erf(nan)=nan */
      298 +
      299 +        if (ix >= 0x7ff00000) {         /* erf(nan)=nan */
 289  300  #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
 290      -                if (ix >= 0x7ff80000)           /* assumes sparc-like QNaN */
      301 +                if (ix >= 0x7ff80000)   /* assumes sparc-like QNaN */
 291  302                          return (x);
 292  303  #endif
 293      -                i = ((unsigned) hx >> 31) << 1;
 294      -                return ((double) (1 - i) + one / x);    /* erf(+-inf)=+-1 */
      304 +                i = ((unsigned)hx >> 31) << 1;
      305 +                return ((double)(1 - i) + one / x);     /* erf(+-inf)=+-1 */
 295  306          }
 296  307  
 297      -        if (ix < 0x3feb0000) {  /* |x|<0.84375 */
 298      -                if (ix < 0x3e300000) {  /* |x|<2**-28 */
 299      -                        if (ix < 0x00800000)    /* avoid underflow */
      308 +        if (ix < 0x3feb0000) {                          /* |x|<0.84375 */
      309 +                if (ix < 0x3e300000) {                  /* |x|<2**-28 */
      310 +                        if (ix < 0x00800000)            /* avoid underflow */
 300  311                                  return (0.125 * (8.0 * x + efx8 * x));
      312 +
 301  313                          return (x + efx * x);
 302  314                  }
      315 +
 303  316                  z = x * x;
 304  317                  r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
 305      -                s = one +
 306      -                        z *(qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
      318 +                s = one + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z *
      319 +                    qq5))));
 307  320                  y = r / s;
 308  321                  return (x + x * y);
 309  322          }
 310      -        if (ix < 0x3ff40000) {  /* 0.84375 <= |x| < 1.25 */
      323 +
      324 +        if (ix < 0x3ff40000) {          /* 0.84375 <= |x| < 1.25 */
 311  325                  s = fabs(x) - one;
 312      -                P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 +
 313      -                        s * (pa5 + s * pa6)))));
 314      -                Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 +
 315      -                        s * (qa5 + s * qa6)))));
      326 +                P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 +
      327 +                    s * pa6)))));
      328 +                Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 +
      329 +                    s * qa6)))));
      330 +
 316  331                  if (hx >= 0)
 317  332                          return (erx + P / Q);
 318  333                  else
 319  334                          return (-erx - P / Q);
 320  335          }
 321      -        if (ix >= 0x40180000) { /* inf > |x| >= 6 */
      336 +
      337 +        if (ix >= 0x40180000) {         /* inf > |x| >= 6 */
 322  338                  if (hx >= 0)
 323  339                          return (one - tiny);
 324  340                  else
 325  341                          return (tiny - one);
 326  342          }
      343 +
 327  344          x = fabs(x);
 328  345          s = one / (x * x);
 329      -        if (ix < 0x4006DB6E) {  /* |x| < 1/0.35 */
 330      -                R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 +
 331      -                        s * (ra5 + s * (ra6 + s * ra7))))));
 332      -                S = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 +
 333      -                        s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
      346 +
      347 +        if (ix < 0x4006DB6E) {          /* |x| < 1/0.35 */
      348 +                R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 + s * (ra5 +
      349 +                    s * (ra6 + s * ra7))))));
      350 +                S = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 + s * (sa5 +
      351 +                    s * (sa6 + s * (sa7 + s * sa8)))))));
 334  352          } else {                        /* |x| >= 1/0.35 */
 335      -                R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 +
 336      -                        s * (rb5 + s * rb6)))));
 337      -                S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 +
 338      -                        s * (sb5 + s * (sb6 + s * sb7))))));
      353 +                R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * (rb5 +
      354 +                    s * rb6)))));
      355 +                S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 + s * (sb5 +
      356 +                    s * (sb6 + s * sb7))))));
 339  357          }
      358 +
 340  359          z = x;
 341      -        ((int *) &z)[LOWORD] = 0;
      360 +        ((int *)&z)[LOWORD] = 0;
 342  361          r = exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S);
      362 +
 343  363          if (hx >= 0)
 344  364                  return (one - r / x);
 345  365          else
 346  366                  return (r / x - one);
 347  367  }
 348  368  
 349  369  double
 350      -erfc(double x) {
      370 +erfc(double x)
      371 +{
 351  372          int hx, ix;
 352  373          double R, S, P, Q, s, y, z, r;
 353  374  
 354      -        hx = ((int *) &x)[HIWORD];
      375 +        hx = ((int *)&x)[HIWORD];
 355  376          ix = hx & 0x7fffffff;
 356      -        if (ix >= 0x7ff00000) { /* erfc(nan)=nan */
      377 +
      378 +        if (ix >= 0x7ff00000) {         /* erfc(nan)=nan */
 357  379  #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
 358      -                if (ix >= 0x7ff80000)           /* assumes sparc-like QNaN */
      380 +                if (ix >= 0x7ff80000)   /* assumes sparc-like QNaN */
 359  381                          return (x);
 360  382  #endif
 361  383                  /* erfc(+-inf)=0,2 */
 362      -                return ((double) (((unsigned) hx >> 31) << 1) + one / x);
      384 +                return ((double)(((unsigned)hx >> 31) << 1) + one / x);
 363  385          }
 364  386  
 365      -        if (ix < 0x3feb0000) {  /* |x| < 0.84375 */
      387 +        if (ix < 0x3feb0000) {          /* |x| < 0.84375 */
 366  388                  if (ix < 0x3c700000)    /* |x| < 2**-56 */
 367  389                          return (one - x);
      390 +
 368  391                  z = x * x;
 369  392                  r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4)));
 370      -                s = one +
 371      -                        z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5))));
      393 +                s = one + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z *
      394 +                    qq5))));
 372  395                  y = r / s;
      396 +
 373  397                  if (hx < 0x3fd00000) {  /* x < 1/4 */
 374  398                          return (one - (x + x * y));
 375  399                  } else {
 376  400                          r = x * y;
 377  401                          r += (x - half);
 378  402                          return (half - r);
 379  403                  }
 380  404          }
 381      -        if (ix < 0x3ff40000) {  /* 0.84375 <= |x| < 1.25 */
      405 +
      406 +        if (ix < 0x3ff40000) {          /* 0.84375 <= |x| < 1.25 */
 382  407                  s = fabs(x) - one;
 383      -                P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 +
 384      -                        s * (pa5 + s * pa6)))));
 385      -                Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 +
 386      -                        s * (qa5 + s * qa6)))));
      408 +                P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 +
      409 +                    s * pa6)))));
      410 +                Q = one + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 +
      411 +                    s * qa6)))));
      412 +
 387  413                  if (hx >= 0) {
 388  414                          z = one - erx;
 389  415                          return (z - P / Q);
 390  416                  } else {
 391  417                          z = erx + P / Q;
 392  418                          return (one + z);
 393  419                  }
 394  420          }
 395      -        if (ix < 0x403c0000) {  /* |x|<28 */
      421 +
      422 +        if (ix < 0x403c0000) {          /* |x|<28 */
 396  423                  x = fabs(x);
 397  424                  s = one / (x * x);
      425 +
 398  426                  if (ix < 0x4006DB6D) {  /* |x| < 1/.35 ~ 2.857143 */
 399  427                          R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 +
 400      -                                s * (ra5 + s * (ra6 + s * ra7))))));
      428 +                            s * (ra5 + s * (ra6 + s * ra7))))));
 401  429                          S = one + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 +
 402      -                                s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
      430 +                            s * (sa5 + s * (sa6 + s * (sa7 + s * sa8)))))));
 403  431                  } else {
 404  432                          /* |x| >= 1/.35 ~ 2.857143 */
 405  433                          if (hx < 0 && ix >= 0x40180000)
 406  434                                  return (two - tiny);    /* x < -6 */
 407  435  
 408  436                          R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 +
 409      -                                s * (rb5 + s * rb6)))));
      437 +                            s * (rb5 + s * rb6)))));
 410  438                          S = one + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 +
 411      -                                s * (sb5 + s * (sb6 + s * sb7))))));
      439 +                            s * (sb5 + s * (sb6 + s * sb7))))));
 412  440                  }
      441 +
 413  442                  z = x;
 414      -                ((int *) &z)[LOWORD] = 0;
      443 +                ((int *)&z)[LOWORD] = 0;
 415  444                  r = exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S);
      445 +
 416  446                  if (hx > 0)
 417  447                          return (r / x);
 418  448                  else
 419  449                          return (two - r / x);
 420  450          } else {
 421  451                  if (hx > 0)
 422  452                          return (tiny * tiny);
 423  453                  else
 424  454                          return (two - tiny);
 425  455          }
 426  456  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX