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