API tools faq
paste
Advertisement
NaZaRa

Untitled

NaZaRa
Jan 20th, 2019
206
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 13.03 KB | None | 0 0
raw download clone embed print report
  1. /* @(#)s_erf.c 5.1 93/09/24 */
  2. /*
  3.  * ====================================================
  4.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  5.  *
  6.  * Developed at SunPro, a Sun Microsystems, Inc. business.
  7.  * Permission to use, copy, modify, and distribute this
  8.  * software is freely granted, provided that this notice
  9.  * is preserved.
  10.  * ====================================================
  11.  */
  12.  
  13. /* double erf(double x)
  14.  * double erfc(double x)
  15.  *                             x
  16.  *                      2      |\
  17.  *     erf(x)  =  ---------  | exp(-t*t)dt
  18.  *                    sqrt(pi) \|
  19.  *                             0
  20.  *
  21.  *     erfc(x) =  1-erf(x)
  22.  *  Note that
  23.  *                erf(-x) = -erf(x)
  24.  *                erfc(-x) = 2 - erfc(x)
  25.  *
  26.  * Method:
  27.  *        1. For |x| in [0, 0.84375]
  28.  *            erf(x)  = x + x*R(x^2)
  29.  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
  30.  *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
  31.  *           where R = P/Q where P is an odd poly of degree 8 and
  32.  *           Q is an odd poly of degree 10.
  33.  *                                                 -57.90
  34.  *                        | R - (erf(x)-x)/x | <= 2
  35.  *        
  36.  *
  37.  *           Remark. The formula is derived by noting
  38.  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
  39.  *           and that
  40.  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
  41.  *           is close to one. The interval is chosen because the fix
  42.  *           point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
  43.  *           near 0.6174), and by some experiment, 0.84375 is chosen to
  44.  *            guarantee the error is less than one ulp for erf.
  45.  *
  46.  *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
  47.  *         c = 0.84506291151 rounded to single (24 bits)
  48.  *                 erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
  49.  *                 erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
  50.  *                          1+(c+P1(s)/Q1(s))    if x < 0
  51.  *                 |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
  52.  *           Remark: here we use the taylor series expansion at x=1.
  53.  *                erf(1+s) = erf(1) + s*Poly(s)
  54.  *                         = 0.845.. + P1(s)/Q1(s)
  55.  *           That is, we use rational approximation to approximate
  56.  *                        erf(1+s) - (c = (single)0.84506291151)
  57.  *           Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
  58.  *           where
  59.  *                P1(s) = degree 6 poly in s
  60.  *                Q1(s) = degree 6 poly in s
  61.  *
  62.  *      3. For x in [1.25,1/0.35(~2.857143)],
  63.  *                 erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
  64.  *                 erf(x)  = 1 - erfc(x)
  65.  *           where
  66.  *                R1(z) = degree 7 poly in z, (z=1/x^2)
  67.  *                S1(z) = degree 8 poly in z
  68.  *
  69.  *      4. For x in [1/0.35,28]
  70.  *                 erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
  71.  *                        = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
  72.  *                        = 2.0 - tiny                (if x <= -6)
  73.  *                 erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
  74.  *                 erf(x)  = sign(x)*(1.0 - tiny)
  75.  *           where
  76.  *                R2(z) = degree 6 poly in z, (z=1/x^2)
  77.  *                S2(z) = degree 7 poly in z
  78.  *
  79.  *      Note1:
  80.  *           To compute exp(-x*x-0.5625+R/S), let s be a single
  81.  *           precision number and s := x; then
  82.  *                -x*x = -s*s + (s-x)*(s+x)
  83.  *                exp(-x*x-0.5626+R/S) =
  84.  *                        exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
  85.  *      Note2:
  86.  *           Here 4 and 5 make use of the asymptotic series
  87.  *                          exp(-x*x)
  88.  *                erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
  89.  *                          x*sqrt(pi)
  90.  *           We use rational approximation to approximate
  91.  *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
  92.  *           Here is the error bound for R1/S1 and R2/S2
  93.  *              |R1/S1 - f(x)|  < 2**(-62.57)
  94.  *              |R2/S2 - f(x)|  < 2**(-61.52)
  95.  *
  96.  *      5. For inf > x >= 28
  97.  *                 erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
  98.  *                 erfc(x) = tiny*tiny (raise underflow) if x > 0
  99.  *                        = 2 - tiny if x<0
  100.  *
  101.  *      7. Special case:
  102.  *                 erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
  103.  *                 erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
  104.  *                   erfc/erf(NaN) is NaN
  105.  */
  106.  
  107. #include <float.h>
  108. #include <math.h>
  109.  
  110. #include "math_private.h"
  111.  
  112. static const double
  113. tiny            = 1e-300,
  114. half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
  115. one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
  116. two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
  117.         /* c = (float)0.84506291151 */
  118. erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
  119. /*
  120.  * Coefficients for approximation to  erf on [0,0.84375]
  121.  */
  122. efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
  123. efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
  124. pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
  125. pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
  126. pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
  127. pp3  = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
  128. pp4  = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
  129. qq1  =  3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
  130. qq2  =  6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
  131. qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
  132. qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
  133. qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
  134. /*
  135.  * Coefficients for approximation to  erf  in [0.84375,1.25]
  136.  */
  137. pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
  138. pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
  139. pa2  = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
  140. pa3  =  3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
  141. pa4  = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
  142. pa5  =  3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
  143. pa6  = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
  144. qa1  =  1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
  145. qa2  =  5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
  146. qa3  =  7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
  147. qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
  148. qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
  149. qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
  150. /*
  151.  * Coefficients for approximation to  erfc in [1.25,1/0.35]
  152.  */
  153. ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
  154. ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
  155. ra2  = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
  156. ra3  = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
  157. ra4  = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
  158. ra5  = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
  159. ra6  = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
  160. ra7  = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
  161. sa1  =  1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
  162. sa2  =  1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
  163. sa3  =  4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
  164. sa4  =  6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
  165. sa5  =  4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
  166. sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
  167. sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
  168. sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
  169. /*
  170.  * Coefficients for approximation to  erfc in [1/.35,28]
  171.  */
  172. rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
  173. rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
  174. rb2  = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
  175. rb3  = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
  176. rb4  = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
  177. rb5  = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
  178. rb6  = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
  179. sb1  =  3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
  180. sb2  =  3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
  181. sb3  =  1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
  182. sb4  =  3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
  183. sb5  =  2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
  184. sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
  185. sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
  186.  
  187. double
  188. erf(double x)
  189. {
  190.         int32_t hx,ix,i;
  191.         double R,S,P,Q,s,y,z,r;
  192.         GET_HIGH_WORD(hx,x);
  193.         ix = hx&0x7fffffff;
  194.         if(ix>=0x7ff00000) {                /* erf(nan)=nan */
  195.             i = ((u_int32_t)hx>>31)<<1;
  196.             return (double)(1-i)+one/x;        /* erf(+-inf)=+-1 */
  197.         }
  198.  
  199.         if(ix < 0x3feb0000) {                /* |x|<0.84375 */
  200.             if(ix < 0x3e300000) {         /* |x|<2**-28 */
  201.                 if (ix < 0x00800000)
  202.                     return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
  203.                 return x + efx*x;
  204.             }
  205.             z = x*x;
  206.             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
  207.             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
  208.             y = r/s;
  209.             return x + x*y;
  210.         }
  211.         if(ix < 0x3ff40000) {                /* 0.84375 <= |x| < 1.25 */
  212.             s = fabs(x)-one;
  213.             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
  214.             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
  215.             if(hx>=0) return erx + P/Q; else return -erx - P/Q;
  216.         }
  217.         if (ix >= 0x40180000) {                /* inf>|x|>=6 */
  218.             if(hx>=0) return one-tiny; else return tiny-one;
  219.         }
  220.         x = fabs(x);
  221.          s = one/(x*x);
  222.         if(ix< 0x4006DB6E) {        /* |x| < 1/0.35 */
  223.             R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
  224.                                 ra5+s*(ra6+s*ra7))))));
  225.             S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
  226.                                 sa5+s*(sa6+s*(sa7+s*sa8)))))));
  227.         } else {        /* |x| >= 1/0.35 */
  228.             R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
  229.                                 rb5+s*rb6)))));
  230.             S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
  231.                                 sb5+s*(sb6+s*sb7))))));
  232.         }
  233.         z  = x;  
  234.         SET_LOW_WORD(z,0);
  235.         r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
  236.         if(hx>=0) return one-r/x; else return  r/x-one;
  237. }
  238. DEF_STD(erf);
  239. LDBL_MAYBE_CLONE(erf);
  240.  
  241. double
  242. erfc(double x)
  243. {
  244.         int32_t hx,ix;
  245.         double R,S,P,Q,s,y,z,r;
  246.         GET_HIGH_WORD(hx,x);
  247.         ix = hx&0x7fffffff;
  248.         if(ix>=0x7ff00000) {                        /* erfc(nan)=nan */
  249.                                                 /* erfc(+-inf)=0,2 */
  250.             return (double)(((u_int32_t)hx>>31)<<1)+one/x;
  251.         }
  252.  
  253.         if(ix < 0x3feb0000) {                /* |x|<0.84375 */
  254.             if(ix < 0x3c700000)          /* |x|<2**-56 */
  255.                 return one-x;
  256.             z = x*x;
  257.             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
  258.             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
  259.             y = r/s;
  260.             if(hx < 0x3fd00000) {          /* x<1/4 */
  261.                 return one-(x+x*y);
  262.             } else {
  263.                 r = x*y;
  264.                 r += (x-half);
  265.                 return half - r ;
  266.             }
  267.         }
  268.         if(ix < 0x3ff40000) {                /* 0.84375 <= |x| < 1.25 */
  269.             s = fabs(x)-one;
  270.             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
  271.             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
  272.             if(hx>=0) {
  273.                 z  = one-erx; return z - P/Q;
  274.             } else {
  275.                 z = erx+P/Q; return one+z;
  276.             }
  277.         }
  278.         if (ix < 0x403c0000) {                /* |x|<28 */
  279.             x = fabs(x);
  280.              s = one/(x*x);
  281.             if(ix< 0x4006DB6D) {        /* |x| < 1/.35 ~ 2.857143*/
  282.                 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
  283.                                 ra5+s*(ra6+s*ra7))))));
  284.                 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
  285.                                 sa5+s*(sa6+s*(sa7+s*sa8)))))));
  286.             } else {                        /* |x| >= 1/.35 ~ 2.857143 */
  287.                 if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
  288.                 R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
  289.                                 rb5+s*rb6)))));
  290.                 S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
  291.                                 sb5+s*(sb6+s*sb7))))));
  292.             }
  293.             z  = x;
  294.             SET_LOW_WORD(z,0);
  295.             r  =  exp(-z*z-0.5625) * exp((z-x)*(z+x)+R/S);
  296.             if(hx>0) return r/x; else return two-r/x;
  297.         } else {
  298.             if(hx>0) return tiny*tiny; else return two-tiny;
  299.         }
  300. }
  301. DEF_STD(erfc);
  302. LDBL_MAYBE_CLONE(erfc);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!