erfl.c 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386
  1. /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_erfl.c */
  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. * Copyright (c) 2008 Stephen L. Moshier <[email protected]>
  14. *
  15. * Permission to use, copy, modify, and distribute this software for any
  16. * purpose with or without fee is hereby granted, provided that the above
  17. * copyright notice and this permission notice appear in all copies.
  18. *
  19. * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  20. * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  21. * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  22. * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  23. * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  24. * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  25. * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  26. */
  27. /* double erf(double x)
  28. * double erfc(double x)
  29. * x
  30. * 2 |\
  31. * erf(x) = --------- | exp(-t*t)dt
  32. * sqrt(pi) \|
  33. * 0
  34. *
  35. * erfc(x) = 1-erf(x)
  36. * Note that
  37. * erf(-x) = -erf(x)
  38. * erfc(-x) = 2 - erfc(x)
  39. *
  40. * Method:
  41. * 1. For |x| in [0, 0.84375]
  42. * erf(x) = x + x*R(x^2)
  43. * erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
  44. * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
  45. * Remark. The formula is derived by noting
  46. * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
  47. * and that
  48. * 2/sqrt(pi) = 1.128379167095512573896158903121545171688
  49. * is close to one. The interval is chosen because the fix
  50. * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
  51. * near 0.6174), and by some experiment, 0.84375 is chosen to
  52. * guarantee the error is less than one ulp for erf.
  53. *
  54. * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
  55. * c = 0.84506291151 rounded to single (24 bits)
  56. * erf(x) = sign(x) * (c + P1(s)/Q1(s))
  57. * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
  58. * 1+(c+P1(s)/Q1(s)) if x < 0
  59. * Remark: here we use the taylor series expansion at x=1.
  60. * erf(1+s) = erf(1) + s*Poly(s)
  61. * = 0.845.. + P1(s)/Q1(s)
  62. * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
  63. *
  64. * 3. For x in [1.25,1/0.35(~2.857143)],
  65. * erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
  66. * z=1/x^2
  67. * erf(x) = 1 - erfc(x)
  68. *
  69. * 4. For x in [1/0.35,107]
  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(z)/S2(z))
  72. * if -6.666<x<0
  73. * = 2.0 - tiny (if x <= -6.666)
  74. * z=1/x^2
  75. * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
  76. * erf(x) = sign(x)*(1.0 - tiny)
  77. * Note1:
  78. * To compute exp(-x*x-0.5625+R/S), let s be a single
  79. * precision number and s := x; then
  80. * -x*x = -s*s + (s-x)*(s+x)
  81. * exp(-x*x-0.5626+R/S) =
  82. * exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
  83. * Note2:
  84. * Here 4 and 5 make use of the asymptotic series
  85. * exp(-x*x)
  86. * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
  87. * x*sqrt(pi)
  88. *
  89. * 5. For inf > x >= 107
  90. * erf(x) = sign(x) *(1 - tiny) (raise inexact)
  91. * erfc(x) = tiny*tiny (raise underflow) if x > 0
  92. * = 2 - tiny if x<0
  93. *
  94. * 7. Special case:
  95. * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
  96. * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
  97. * erfc/erf(NaN) is NaN
  98. */
  99. #include "libm.h"
  100. #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  101. long double erfl(long double x)
  102. {
  103. return erf(x);
  104. }
  105. #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
  106. static const long double
  107. tiny = 1e-4931L,
  108. /* c = (float)0.84506291151 */
  109. erx = 0.845062911510467529296875L,
  110. /*
  111. * Coefficients for approximation to erf on [0,0.84375]
  112. */
  113. /* 2/sqrt(pi) - 1 */
  114. efx = 1.2837916709551257389615890312154517168810E-1L,
  115. /* 8 * (2/sqrt(pi) - 1) */
  116. efx8 = 1.0270333367641005911692712249723613735048E0L,
  117. pp[6] = {
  118. 1.122751350964552113068262337278335028553E6L,
  119. -2.808533301997696164408397079650699163276E6L,
  120. -3.314325479115357458197119660818768924100E5L,
  121. -6.848684465326256109712135497895525446398E4L,
  122. -2.657817695110739185591505062971929859314E3L,
  123. -1.655310302737837556654146291646499062882E2L,
  124. },
  125. qq[6] = {
  126. 8.745588372054466262548908189000448124232E6L,
  127. 3.746038264792471129367533128637019611485E6L,
  128. 7.066358783162407559861156173539693900031E5L,
  129. 7.448928604824620999413120955705448117056E4L,
  130. 4.511583986730994111992253980546131408924E3L,
  131. 1.368902937933296323345610240009071254014E2L,
  132. /* 1.000000000000000000000000000000000000000E0 */
  133. },
  134. /*
  135. * Coefficients for approximation to erf in [0.84375,1.25]
  136. */
  137. /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
  138. -0.15625 <= x <= +.25
  139. Peak relative error 8.5e-22 */
  140. pa[8] = {
  141. -1.076952146179812072156734957705102256059E0L,
  142. 1.884814957770385593365179835059971587220E2L,
  143. -5.339153975012804282890066622962070115606E1L,
  144. 4.435910679869176625928504532109635632618E1L,
  145. 1.683219516032328828278557309642929135179E1L,
  146. -2.360236618396952560064259585299045804293E0L,
  147. 1.852230047861891953244413872297940938041E0L,
  148. 9.394994446747752308256773044667843200719E-2L,
  149. },
  150. qa[7] = {
  151. 4.559263722294508998149925774781887811255E2L,
  152. 3.289248982200800575749795055149780689738E2L,
  153. 2.846070965875643009598627918383314457912E2L,
  154. 1.398715859064535039433275722017479994465E2L,
  155. 6.060190733759793706299079050985358190726E1L,
  156. 2.078695677795422351040502569964299664233E1L,
  157. 4.641271134150895940966798357442234498546E0L,
  158. /* 1.000000000000000000000000000000000000000E0 */
  159. },
  160. /*
  161. * Coefficients for approximation to erfc in [1.25,1/0.35]
  162. */
  163. /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
  164. 1/2.85711669921875 < 1/x < 1/1.25
  165. Peak relative error 3.1e-21 */
  166. ra[] = {
  167. 1.363566591833846324191000679620738857234E-1L,
  168. 1.018203167219873573808450274314658434507E1L,
  169. 1.862359362334248675526472871224778045594E2L,
  170. 1.411622588180721285284945138667933330348E3L,
  171. 5.088538459741511988784440103218342840478E3L,
  172. 8.928251553922176506858267311750789273656E3L,
  173. 7.264436000148052545243018622742770549982E3L,
  174. 2.387492459664548651671894725748959751119E3L,
  175. 2.220916652813908085449221282808458466556E2L,
  176. },
  177. sa[] = {
  178. -1.382234625202480685182526402169222331847E1L,
  179. -3.315638835627950255832519203687435946482E2L,
  180. -2.949124863912936259747237164260785326692E3L,
  181. -1.246622099070875940506391433635999693661E4L,
  182. -2.673079795851665428695842853070996219632E4L,
  183. -2.880269786660559337358397106518918220991E4L,
  184. -1.450600228493968044773354186390390823713E4L,
  185. -2.874539731125893533960680525192064277816E3L,
  186. -1.402241261419067750237395034116942296027E2L,
  187. /* 1.000000000000000000000000000000000000000E0 */
  188. },
  189. /*
  190. * Coefficients for approximation to erfc in [1/.35,107]
  191. */
  192. /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
  193. 1/6.6666259765625 < 1/x < 1/2.85711669921875
  194. Peak relative error 4.2e-22 */
  195. rb[] = {
  196. -4.869587348270494309550558460786501252369E-5L,
  197. -4.030199390527997378549161722412466959403E-3L,
  198. -9.434425866377037610206443566288917589122E-2L,
  199. -9.319032754357658601200655161585539404155E-1L,
  200. -4.273788174307459947350256581445442062291E0L,
  201. -8.842289940696150508373541814064198259278E0L,
  202. -7.069215249419887403187988144752613025255E0L,
  203. -1.401228723639514787920274427443330704764E0L,
  204. },
  205. sb[] = {
  206. 4.936254964107175160157544545879293019085E-3L,
  207. 1.583457624037795744377163924895349412015E-1L,
  208. 1.850647991850328356622940552450636420484E0L,
  209. 9.927611557279019463768050710008450625415E0L,
  210. 2.531667257649436709617165336779212114570E1L,
  211. 2.869752886406743386458304052862814690045E1L,
  212. 1.182059497870819562441683560749192539345E1L,
  213. /* 1.000000000000000000000000000000000000000E0 */
  214. },
  215. /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
  216. 1/107 <= 1/x <= 1/6.6666259765625
  217. Peak relative error 1.1e-21 */
  218. rc[] = {
  219. -8.299617545269701963973537248996670806850E-5L,
  220. -6.243845685115818513578933902532056244108E-3L,
  221. -1.141667210620380223113693474478394397230E-1L,
  222. -7.521343797212024245375240432734425789409E-1L,
  223. -1.765321928311155824664963633786967602934E0L,
  224. -1.029403473103215800456761180695263439188E0L,
  225. },
  226. sc[] = {
  227. 8.413244363014929493035952542677768808601E-3L,
  228. 2.065114333816877479753334599639158060979E-1L,
  229. 1.639064941530797583766364412782135680148E0L,
  230. 4.936788463787115555582319302981666347450E0L,
  231. 5.005177727208955487404729933261347679090E0L,
  232. /* 1.000000000000000000000000000000000000000E0 */
  233. };
  234. long double erfl(long double x)
  235. {
  236. long double R, S, P, Q, s, y, z, r;
  237. int32_t ix, i;
  238. uint32_t se, i0, i1;
  239. GET_LDOUBLE_WORDS(se, i0, i1, x);
  240. ix = se & 0x7fff;
  241. if (ix >= 0x7fff) { /* erf(nan)=nan */
  242. i = ((se & 0xffff) >> 15) << 1;
  243. return (long double)(1 - i) + 1.0 / x; /* erf(+-inf)=+-1 */
  244. }
  245. ix = (ix << 16) | (i0 >> 16);
  246. if (ix < 0x3ffed800) { /* |x| < 0.84375 */
  247. if (ix < 0x3fde8000) { /* |x| < 2**-33 */
  248. if (ix < 0x00080000)
  249. return 0.125 * (8.0 * x + efx8 * x); /* avoid underflow */
  250. return x + efx * x;
  251. }
  252. z = x * x;
  253. r = pp[0] + z * (pp[1] +
  254. z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
  255. s = qq[0] + z * (qq[1] +
  256. z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
  257. y = r / s;
  258. return x + x * y;
  259. }
  260. if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
  261. s = fabsl(x) - 1.0;
  262. P = pa[0] + s * (pa[1] + s * (pa[2] +
  263. s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
  264. Q = qa[0] + s * (qa[1] + s * (qa[2] +
  265. s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
  266. if ((se & 0x8000) == 0)
  267. return erx + P / Q;
  268. return -erx - P / Q;
  269. }
  270. if (ix >= 0x4001d555) { /* inf > |x| >= 6.6666259765625 */
  271. if ((se & 0x8000) == 0)
  272. return 1.0 - tiny;
  273. return tiny - 1.0;
  274. }
  275. x = fabsl (x);
  276. s = 1.0 / (x * x);
  277. if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
  278. R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
  279. s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
  280. S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
  281. s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
  282. } else { /* 2.857 <= |x| < 6.667 */
  283. R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
  284. s * (rb[5] + s * (rb[6] + s * rb[7]))))));
  285. S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
  286. s * (sb[5] + s * (sb[6] + s))))));
  287. }
  288. z = x;
  289. GET_LDOUBLE_WORDS(i, i0, i1, z);
  290. i1 = 0;
  291. SET_LDOUBLE_WORDS(z, i, i0, i1);
  292. r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
  293. if ((se & 0x8000) == 0)
  294. return 1.0 - r / x;
  295. return r / x - 1.0;
  296. }
  297. long double erfcl(long double x)
  298. {
  299. int32_t hx, ix;
  300. long double R, S, P, Q, s, y, z, r;
  301. uint32_t se, i0, i1;
  302. GET_LDOUBLE_WORDS(se, i0, i1, x);
  303. ix = se & 0x7fff;
  304. if (ix >= 0x7fff) { /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
  305. return (long double)(((se & 0xffff) >> 15) << 1) + 1.0 / x;
  306. }
  307. ix = (ix << 16) | (i0 >> 16);
  308. if (ix < 0x3ffed800) { /* |x| < 0.84375 */
  309. if (ix < 0x3fbe0000) /* |x| < 2**-65 */
  310. return 1.0 - x;
  311. z = x * x;
  312. r = pp[0] + z * (pp[1] +
  313. z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
  314. s = qq[0] + z * (qq[1] +
  315. z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
  316. y = r / s;
  317. if (ix < 0x3ffd8000) /* x < 1/4 */
  318. return 1.0 - (x + x * y);
  319. r = x * y;
  320. r += x - 0.5L;
  321. return 0.5L - r;
  322. }
  323. if (ix < 0x3fffa000) { /* 0.84375 <= |x| < 1.25 */
  324. s = fabsl(x) - 1.0;
  325. P = pa[0] + s * (pa[1] + s * (pa[2] +
  326. s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
  327. Q = qa[0] + s * (qa[1] + s * (qa[2] +
  328. s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
  329. if ((se & 0x8000) == 0) {
  330. z = 1.0 - erx;
  331. return z - P / Q;
  332. }
  333. z = erx + P / Q;
  334. return 1.0 + z;
  335. }
  336. if (ix < 0x4005d600) { /* |x| < 107 */
  337. x = fabsl(x);
  338. s = 1.0 / (x * x);
  339. if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
  340. R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
  341. s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
  342. S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
  343. s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
  344. } else if (ix < 0x4001d555) { /* 6.6666259765625 > |x| >= 1/.35 ~ 2.857143 */
  345. R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
  346. s * (rb[5] + s * (rb[6] + s * rb[7]))))));
  347. S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
  348. s * (sb[5] + s * (sb[6] + s))))));
  349. } else { /* 107 > |x| >= 6.666 */
  350. if (se & 0x8000)
  351. return 2.0 - tiny;/* x < -6.666 */
  352. R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
  353. s * (rc[4] + s * rc[5]))));
  354. S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
  355. s * (sc[4] + s))));
  356. }
  357. z = x;
  358. GET_LDOUBLE_WORDS(hx, i0, i1, z);
  359. i1 = 0;
  360. i0 &= 0xffffff00;
  361. SET_LDOUBLE_WORDS(z, hx, i0, i1);
  362. r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
  363. if ((se & 0x8000) == 0)
  364. return r / x;
  365. return 2.0 - r / x;
  366. }
  367. if ((se & 0x8000) == 0)
  368. return tiny * tiny;
  369. return 2.0 - tiny;
  370. }
  371. #endif