C中的'erf'功能 [英] 'erf' function in C
问题描述
是否有任何具有erf功能的C编译器（来自
概率）作为数学库的一部分？或者是否有任何数学可供下载实现此功能的
库？
谢谢，
James Harlacher
Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher
推荐答案
" James Harlacher" < JP **** @ hotmail.com> écritdansle message de
news：82 ************************** @ posting.google.c om ...
"James Harlacher" <jp****@hotmail.com> a écrit dans le message de
news:82**************************@posting.google.c om...
是否有任何具有erf功能（来自
概率）的C编译器作为其数学库的一部分？或者是否有可供下载的数学库实现此功能？
谢谢，
James Harlacher
Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
Thanks,
James Harlacher
是的。 lcc-win32使用Sun的实现。
/ * @（＃）s_erf.c 5.1 93/09/24 * /
/ *
* ========================================== ======== ==
* Sun Microsystems，Inc。版权所有（C）1993。保留所有权利。
*
*在Sun Microsystems，Inc。业务的SunPro开发。
*允许使用，复制，修改和分发此
*软件是免费授予的，前提是这个通知
*被保留。
* =========================== ======================= ==
* /
/ * double erf（双x）
*双erfc（双x）
* x
* 2 | \
* erf（x）= --------- | exp（-t * t）dt
* sqrt（pi）\ |
* 0
*
* erfc（x）= 1-erf（x）
*注意
* erf（-x）= -erf（x）
* erfc（-x）= 2 - erfc（x）
*
*方法：
* 1.对于| x | in [0,0.84375]
* erf（x）= x + x * R（x ^ 2）
* erfc（x）= 1 - erf（x）如果x在[-.84375,0.25]
* = 0.5 +（（0.5-x）-x * R），如果x在[0.25,0.84375]
*其中R = P / Q，其中P是8度的奇数多边形，
* Q是10度的奇数多边形。
* -57.90
* | R - （erf（x）-x）/ x | < = 2
*
*
*备注。通过注意来得出公式
* erf（x）=（2 / sqrt（pi））*（x - x ^ 3/3 + x ^ 5/10-x ^ 7/42 + ....）
*和
* 2 / sqrt（pi）= 1.128379167095512573896158903121545171688
*接近1。选择间隔是因为erf（x）的修正点接近0.6174（即，当x为
*接近0.6174时，erf（x）= x），通过一些实验，选择0.84375作为
*保证erf的误差小于1 ulp。
*
* 2.对于| X |在[0.84375,1.25]中，设s = | x | - 1，
* c = 0.84506291151四舍五入到单（24位）
* erf（x）= sign（x）*（c + P1（s）/ Q1（s））
* erfc（x）=（1-c）-P1（s）/ Q1（s）如果x> 0
* 1+（c + P1（s）/ Q1（s））如果x <0。 0
* | P1 / Q1 - （erf（| x |） - c）| < = 2 ** - 59.06
*备注：这里我们在x = 1时使用泰勒系列扩展。
* erf（1 + s）= erf（1 ）+ s * Poly（s）
* = 0.845 .. + P1（s）/ Q1（s）
*也就是说，我们使用有理逼近近似
* erf（1 + s） - （c =（单）0.84506291151）
*注意| P1 / Q1 |< 0.078 for x in [0.84375,1.25]
* where
* P1（s）= degree 6 poly in s
* Q1（s ）= 6度聚合物s
*
* 3.对于[1.25,1 / 0.35（~2.857143）]中的x，
* erfc（x）=（1 / x）* exp（-x * x-0.5625 + R1 / S1）
* erf（x）= 1 - erfc（x）
*其中
* R1（z）=度数7在z中，（z = 1 / x ^ 2）
* S1（z）= 8级poly in z
*
* 4.对于[1 / 0.35,28]中的x
* erfc（x）=（1 / x）* exp（-x * x-0.5625 + R2 / S2）如果x> 0
* = 2.0 - （1 / x）* exp（-x * x-0.5625 + R2 / S2）如果-6< x< 0
* = 2.0 - 微小（如果x <= -6）
* erf（x）= sign（x）*（1.0 - erfc（x））如果x< 6，否则
* erf（x）= sign（x）*（1.0 - 微小）
*其中
* R2（z） =度数6在z中，（z = 1 / x ^ 2）
* S2（z）=度数7在z中
*
*注1：
*要计算exp（-x * x-0.5625 + R / S），设s为单个
*精度数，s：= X;然后
* -x * x = -s * s +（sx）*（s + x）
* exp（-x * x-0.5626 + R / S ）=
* exp（-s * s-0.5625）* exp（（sx）*（s + x）+ R / S）;
*注2：
*这里4和5使用渐近系列
* exp（-x * x）
* erfc（x）〜 - -------- *（1 + Poly（1 / x ^ 2））
* x * sqrt（pi）
*我们使用有理逼近近似值
* g（s）= f（1 / x ^ 2）= log（erfc（x）* x） - x * x + 0.5625
*这是R1 / S1和R2 / S2的误差范围
* | R1 / S1 - f（x）| < 2 **（ - 62.57）
* | R2 / S2 - f（x）| < 2 **（ - 61.52）
*
* 5.对于inf> x> = 28
* erf（x）= sign（x）*（1 - 微小）（提高不精确度）
* erfc（x）= tiny * tiny （提高下溢）如果x> 0
* = 2 - 微小如果x <0
*
* 7.特殊情况：
* erf（0）= 0，erf（inf）= 1，erf（-inf）= -1，
* erfc（0）= 1，erfc（inf）= 0，erfc（-inf ）= 2，
* erfc / erf（NaN）是NaN
* /
静态const double tiny = 1e-300，
half = 5.00000000000000000000e-01，/ * 0x3FE00000,0x00000000 * /
one = 1.00000000000000000000e + 00，/ * 0x3FF00000,0x00000000 * /
2 = 2.00000000000000000000e + 00，/ * 0x40000000,0x00000000 * /
/ * c =（浮动）0.84506291151 * /
erx = 8.45062911510467529297e-01，/ * 0x3FEB0AC1 ，0x60000000 * /
/ *
*近似于[0,0.84375]的erf的系数
* /
efx = 1.28379167095512586316e-01，/ * 0x3FC06EBA，0x8214DB69 * /
efx8 = 1.02703333676410069053e + 00，/ * 0x3FF06EBA，0x8214DB69 * /
pp0 = 1.28379167095512558561e -01，/ * 0x3FC06EBA，0x8214DB68 * /
pp1 = -3.25042107247001499370e-01，/ * 0xBFD4CD7D，0x691CB913 * /
pp2 = -2.84817495755985104766e-02，/ * 0xBF9D2A51,0xDBD7194F * /
pp3 = -5.770270296489441515157e- 03，/ * 0xBF77A291,0x236668E4 * /
pp4 = -2.37630166566501626084e-05，/ * 0xBEF8EAD6,0x120016AC * /
qq1 = 3.97917223959155352819e-01，/ * 0x3FD97779 ，0xCDDADC09 * /
qq2 = 6.50222499887672944485e-02，/ * 0x3FB0A54C，0x5536CEBA * /
qq3 = 5.08130628187576562776e-03，/ * 0x3F74D022,0xC4D36B0F * /
qq4 = 1.32494738004321644526e-04，/ * 0x3F215DC9,0x221C1A10 * /
qq5 = -3.96022827877536812320e-06，/ * 0xBED09C43,0x42A26120 * /
/ *
*近似于erf的系数[0.84375,1.25]
* /
pa0 = -2.36211856075265944077e-03，/ * 0xBF6359B8,0xBEF77538 * /
pa1 = 4.14856118683748331666e-01，/ * 0x3FDA8D00,0xAD92B34D * /
pa2 = -3.72207876035701323847e-01，/ * 0xBFD7D240,0xFBB8C3F1 * /
pa3 = 3.18346619901161753674e-01，/ * 0x3FD45FCA，0x805120E4 * /
pa4 = -1.10894694282396677476e-01，/ * 0xBFBC6398,0x3D3E28EC * /
pa5 = 3.54783043256182359371e-02，/ * 0x3FA22A36,0x599795EB * /
pa6 = -2.16637559486879084300e-03，/ * 0xBF61BF38,0x0A96073F * /
qa1 = 1.06420880400844228286e-01 ，/ * 0x3FBB3E66,0x18EEE323 * /
qa2 = 5.40397917702171048937e-01，/ * 0x3FE14AF0,0x92EB6F33 * /
qa3 = 7.18286544141962662868e-02，/ * 0x3FB2635C，0xD99FE9A7 * /
qa4 = 1.26171219808761642112e-01，/ * 0x3FC02660,0xE763351F * /
qa5 = 1.36370839120290507362e-02，/ * 0x3F8BEDC2,0x6B51DD1C * /
qa6 = 1.19844998467991074170e-02，/ * 0x3F888B54,0x5735151D * /
/ *
*在[1.25,1 / 0.35]中逼近erfc的系数
* /
ra0 = -9.86494403484714822705e-03，/ * 0xBF843412,0x600D6435 * /
ra1 = -6.93858572707181764372e-01，/ * 0xBFE63416,0xE4BA7360 * /
ra2 = -1.05586262253232909814e + 01，/ * 0xC0251E04,0x41B0E726 * /
ra3 = -6.23753324503260060396e + 01，/ * 0xC04F300A，0xE4CBA38D * /
ra4 = -1.62396669462573470355e + 02，/ * 0xC0644CB1,0x84282266 * /
ra5 = -1.84605092906711035994e + 02，/ * 0xC067135C，0xEBCCABB2 * /
ra6 = -8.12874355063065934246e + 01，/ * 0xC0545265,0x57E4D2F2 * /
ra7 = -9.81432934416914548592e + 00，/ * 0xC023A0EF，0xC69AC25C * /
sa1 = 1.96512716674392571292e + 01，/ * 0x4033A6B9,0xBD707687 * /
sa2 = 1.37657754143519042600e + 02，/ * 0x4061350C，0x526AE721 * /
sa3 = 4.34565877475229228821e + 02，/ * 0x407B290D，0xD58A1A71 * /
sa4 = 6.45387271733267880336e + 02，/ * 0x40842B19,0x21EC2868 * /
sa5 = 4.29008140027567833386e + 02，/ * 0x407AD021,0x57700314 * /
sa6 = 1.08635005541779435134e + 02，/ * 0x405B28A3,0xEE48AE2C * /
sa7 = 6.57024977031928170135e + 00，/ * 0x401A47EF，0x8E484A93 * /
sa8 = -6.04244152148580987438e-02，/ * 0xBFAEEFF2,0xEE749A62 * /
/ *
*在[1 / .35,28]中逼近erfc的系数
* /
rb0 = -9.86494292470009928597e-03，/ * 0xBF843412,0x39E86F4A * /
rb1 = -7.99283237680523006574e-01，/ * 0xBFE993BA，0x70C285DE * /
rb2 = -1.77579549177547519889e + 01，/ * 0xC031C209,0x555F995A * /
rb3 = -1.60636384855821916062e + 02，/ * 0xC064145D，0x43C5ED98 * /
rb4 = -6.37566443368389627722e + 02，/ * 0xC083EC88,0x1375F228 * /
rb5 = -1.02509513161107724954e + 03，/ * 0xC0900461,0x6A2E5992 * /
rb6 = -4.83519191608651397019e + 02，/ * 0xC07E384E，0x9BDC383F * /
sb1 = 3.03380607434824582924e + 01，/ * 0x403E568B，0x261D5190 * /
sb2 = 3.25792512996573918826e + 02，/ * 0x40745CAE，0x221B9F0A * /
sb3 = 1.53672958608443695994e + 03，/ * 0x409802EB，0x189D5118 * /
sb4 = 3.19985821950859553908e + 03，/ * 0x40A8FFB7,0x688C246A * /
sb5 = 2.55305040643316442583e + 03，/ * 0x40A3F219,0xCEDF3BE6 * /
sb6 = 4.74528541206955367215e + 02，/ * 0x407DA874,0xE79FE763 * /
sb7 = -2.24409524465858183362e + 01; / * 0xC03670E2,0x42712D62 * /
extern double exp（double）;
extern double fabs（double）;
double erf（双x）
{
int n0，hx，ix，i;
双R，S，P，Q，s， y，z，r;
n0 =（（*（int *）& one）>> 29）^ 1;
hx = *（n0 +（ int *）& x）;
ix = hx& 0x7fffffff;
if（ix> = 0x7ff00000）{/ * erf（nan）= nan * /
i =（（unsigned）hx>> 31）<< 1;
return（double）（1-i）+ one / x; / * erf（+ - inf）= + - 1 * /
}
if（ix< 0x3feb0000）{/ * | x |< 0.84375 * /
if（ix< 0x3e300000）{/ * | x |< 2 ** - 28 * /
if（ix< 0x00800000）
返回0.125 *（8.0 * x + efx8 * x）; / *避免下溢* /
返回x + efx * x;
}
z = x * x;
r = pp0 + z *（pp1 + z *（pp2 + z *（pp3 + z * pp4）））;
s = 1 + z *（qq1 + z *（qq2 + z *（ qq3 + z *（qq4 + z * qq5））））;
y = r / s;
返回x + x * y;
}
if（ix< 0x3ff40000）{/ * 0.84375< = | x | < 1.25 * /
s = fabs（x）-one;
P = pa0 + s *（pa1 + s *（pa2 + s *（pa3 + s *（pa4） + s *（pa5 + s * pa6）））））;
Q = 1 + s *（qa1 + s *（qa2 + s *（qa3 + s *（qa4 + s *（ qa5 + s * qa6）））））;
if（hx> = 0）返回erx + P / Q;否则返回-erx - P / Q;
}
if（ix> = 0x40180000）{/ * inf> | x |> = 6 * /
if（hx> = 0）返回one-tiny;否则返回tiny-one;
}
x = fabs（x）;
s = 1 /（x * x）;
if（ix< 0x4006DB6E）{/ * | x | < 1 / 0.35 * /
R = ra0 + s *（ra1 + s *（ra2 + s *（ra3 + s *（ra4 + s *）（
ra5 + s *（ra6 + s * ra7）））））;
S = 1 + s *（sa1 + s *（sa2 + s *（sa3 + s *（sa4 + s *（
sa5 + s *（sa6 + s *（sa7 + s * sa8）））））））;
} else {/ * | x | > = 1 / 0.35 * /
R = rb0 + s *（rb1 + s *（rb2 + s *（rb3 + s *）（rb4 + s *（
$ b $） b rb5 + s * rb6））））;
S = 1 + s *（sb1 + s *（sb2 + s *）（sb3 + s *（sb4 + s *（
sb5 + s *（sb6 + s * sb7））））））;
}
z = x;
*（1-n0 +（int *）& z）= 0;
r = exp（-z * z-0.5625）* exp（（zx）*（z + x）+ R / S ）;
if（hx> = 0）返回1-r / x;否则返回r / x-one;
}
双倍erfc（双x）
{
int n0，hx，ix;
双R，S，P，Q，s，y，z，r;
n0 =（（*（int * ）& one）>> 29）^ 1;
hx = *（n0 +（int *）& x）;
ix = hx& 0x7fffffff;
if（ix> = 0x7ff00000）{/ * erfc（nan）= nan * /
/ * erfc（+ - inf）= 0,2 * /
return（double）（（（unsigned）hx>> 31）<< 1）+ one / x;
}
if（ix< 0x3feb0000）{/ * | x |< 0.84375 * /
if（ix< 0x3c700000）/ * | x |< 2 ** - 56 * /
返回1-x;
z = x * x;
r = pp0 + z *（pp1 + z *（pp2 + z *（ pp3 + z * pp4）））;
s = 1 + z *（qq1 + z *（qq2 + z *（qq3 + z *（qq4 + z * qq5））））;
y = r / s;
if（hx< 0x3fd00000）{/ * x< 1/4 * /
返回一个 - （x + x * y）;
}否则{
r = x * y;
r + =（x-half）;
返回半 - r;
}
}
if（ix< 0x3ff40000）{/ * 0.84375< = | x | < 1.25 * /
s = fabs（x）-one;
P = pa0 + s *（pa1 + s *（pa2 + s *（pa3 + s *（pa4） + s *（pa5 + s * pa6）））））;
Q = 1 + s *（qa1 + s *（qa2 + s *（qa3 + s *（qa4 + s *（ qa5 + s * qa6）））））;
if（hx> = 0）{
z = one-erx;返回z - P / Q;
}否则{
z = erx + P / Q;如果（ix< 0x403c0000）{/ * | x |< 28 * /
x = fabs（x）;
s = 1 /（x * x）;
if（ix< 0x4006DB6D）{/ * | x | < 1 / .35~2.857143 * /
R = ra0 + s *（ra1 + s *（ra2 + s *（ra3 + s *（ra4 + s *（
$ b $） b ra5 + s *（ra6 + s * ra7）））））;
S = 1 + s *（sa1 + s *（sa2 + s *（sa3 + s *（sa4 +） s *（
sa5 + s *（sa6 + s *（sa7 + s * sa8）））））））;
} else {/ * | x | > = 1 / .35~2.857143 * /
if（hx< 0&& ix> = 0x40180000）返回两个微小; / * x< -6 * /
R = rb0 + s *（rb1 + s *（rb2 + s *（rb3 + s *）（rb4 + s *（
rb5 + s） * rb6））））;
S = 1 + s *（sb1 + s *（sb2 + s *（sb3 + s *）（sb4 + s *（
$ b $） b sb5 + s *（sb6 + s * sb7））））））;
}
z = x;
*（1- n0 +（int *）& z）= 0;
r = exp（-z * z-0.5625）*
exp（（zx）*（z + x） + R / S）;
如果（hx> 0）返回r / x;否则返回2-r / x;
}否则{
如果（hx> 0）返回tiny * tiny;否则返回两个小;
}
}
Yes. lcc-win32 uses Sun''s implementation.
/* @(#)s_erf.c 5.1 93/09/24 */
/*
* ================================================== ==
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ================================================== ==
*/
/* double erf(double x)
* double erfc(double x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
* sqrt(pi) \|
* 0
*
* erfc(x) = 1-erf(x)
* Note that
* erf(-x) = -erf(x)
* erfc(-x) = 2 - erfc(x)
*
* Method:
* 1. For |x| in [0, 0.84375]
* erf(x) = x + x*R(x^2)
* erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
* = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
* where R = P/Q where P is an odd poly of degree 8 and
* Q is an odd poly of degree 10.
* -57.90
* | R - (erf(x)-x)/x | <= 2
*
*
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
* and that
* 2/sqrt(pi) = 1.128379167095512573896158903121545171688
* is close to one. The interval is chosen because the fix
* point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
* near 0.6174), and by some experiment, 0.84375 is chosen to
* guarantee the error is less than one ulp for erf.
*
* 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
* erf(x) = sign(x) * (c + P1(s)/Q1(s))
* erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
* 1+(c+P1(s)/Q1(s)) if x < 0
* |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
* That is, we use rational approximation to approximate
* erf(1+s) - (c = (single)0.84506291151)
* Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
* where
* P1(s) = degree 6 poly in s
* Q1(s) = degree 6 poly in s
*
* 3. For x in [1.25,1/0.35(~2.857143)],
* erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
* erf(x) = 1 - erfc(x)
* where
* R1(z) = degree 7 poly in z, (z=1/x^2)
* S1(z) = degree 8 poly in z
*
* 4. For x in [1/0.35,28]
* erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
* = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
* = 2.0 - tiny (if x <= -6)
* erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
* erf(x) = sign(x)*(1.0 - tiny)
* where
* R2(z) = degree 6 poly in z, (z=1/x^2)
* S2(z) = degree 7 poly in z
*
* Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
* exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
* Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
* We use rational approximation to approximate
* g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
* Here is the error bound for R1/S1 and R2/S2
* |R1/S1 - f(x)| < 2**(-62.57)
* |R2/S2 - f(x)| < 2**(-61.52)
*
* 5. For inf > x >= 28
* erf(x) = sign(x) *(1 - tiny) (raise inexact)
* erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
* 7. Special case:
* erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
* erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
* erfc/erf(NaN) is NaN
*/
static const double tiny = 1e-300,
half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
/* c = (float)0.84506291151 */
erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
/*
* Coefficients for approximation to erf on [0,0.84375]
*/
efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
/*
* Coefficients for approximation to erf in [0.84375,1.25]
*/
pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
/*
* Coefficients for approximation to erfc in [1.25,1/0.35]
*/
ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
/*
* Coefficients for approximation to erfc in [1/.35,28]
*/
rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
extern double exp(double);
extern double fabs(double);
double erf(double x)
{
int n0,hx,ix,i;
double R,S,P,Q,s,y,z,r;
n0 = ((*(int*)&one)>>29)^1;
hx = *(n0+(int*)&x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erf(nan)=nan */
i = ((unsigned)hx>>31)<<1;
return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
}
if(ix < 0x3feb0000) { /* |x|<0.84375 */
if(ix < 0x3e300000) { /* |x|<2**-28 */
if (ix < 0x00800000)
return 0.125*(8.0*x+efx8*x); /*avoid underflow */
return x + efx*x;
}
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
return x + x*y;
}
if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx>=0) return erx + P/Q; else return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf>|x|>=6 */
if(hx>=0) return one-tiny; else return tiny-one;
}
x = fabs(x);
s = one/(x*x);
if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */
R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
*(1-n0+(int*)&z) = 0;
r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
if(hx>=0) return one-r/x; else return r/x-one;
}
double erfc(double x)
{
int n0,hx,ix;
double R,S,P,Q,s,y,z,r;
n0 = ((*(int*)&one)>>29)^1;
hx = *(n0+(int*)&x);
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (double)(((unsigned)hx>>31)<<1)+one/x;
}
if(ix < 0x3feb0000) { /* |x|<0.84375 */
if(ix < 0x3c700000) /* |x|<2**-56 */
return one-x;
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
if(hx < 0x3fd00000) { /* x<1/4 */
return one-(x+x*y);
} else {
r = x*y;
r += (x-half);
return half - r ;
}
}
if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
if(hx>=0) {
z = one-erx; return z - P/Q;
} else {
z = erx+P/Q; return one+z;
}
}
if (ix < 0x403c0000) { /* |x|<28 */
x = fabs(x);
s = one/(x*x);
if(ix< 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
*(1-n0+(int*)&z) = 0;
r = exp(-z*z-0.5625)*
exp((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x;
} else {
if(hx>0) return tiny*tiny; else return two-tiny;
}
}
" James Harlacher" < JP **** @ hotmail.com>在消息中写道
news：82 ************************** @ posting.google.c om ...
"James Harlacher" <jp****@hotmail.com> wrote in message
news:82**************************@posting.google.c om...
是否有任何具有erf功能（来自
概率）的C编译器作为其数学库的一部分？或者是否有可供下载实现此功能的数学库？
Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
您的网络浏览器发生了什么变化？ gcc实现的大多数库（包括glibc和newlib）都有这个功能。
What happened to your web browser? Most libraries on which gcc is
implemented, including glibc and newlib, have this function.
James Harlacher写道：
James Harlacher wrote:
是否有任何具有erf功能（来自
概率）的C编译器作为其数学库的一部分？或者是否有可供下载的数学库实现此功能？
Hi,
Are there any C compilers which have the erf function (from
probability) as part of their math libraries? Or are there any math
libraries available to download which implement this function?
任何C99编译器都有它; gcc有它。
祝你好运，
Sidney
Any C99 compiler has it; gcc has it.
Best regards,
Sidney
这篇关于C中的'erf'功能的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
