overby.me
this repo has no description
1
2/*
3 * erf.h
4 *
5 * by Ian Ollmann, based on original work by Ali Sazegari, Steve Peters and Jeff Kidder
6 *
7 * Copyright (c) 2007 Apple Inc. All Rights Reserved.
8 *
9 */
10
11#include <math.h>
12
13static const double xxbig = 27.2e+0;
14static const double Maximum = 2.53e+307;
15static const double _HUGE = 6.71e+7;
16static const double InvSqrtPI = 5.6418958354775628695e-1;
17
18static const double a[5] = { 3.16112374387056560e+0,
19 1.13864154151050156e+2,
20 3.77485237685302021e+2,
21 3.20937758913846947e+3,
22 1.85777706184603153e-1 };
23
24static const double b[4] = { 2.36012909523441209e+1,
25 2.44024637934444173e+2,
26 1.28261652607737228e+3,
27 2.84423683343917062e+3 };
28
29static const double ccc[9] = { 5.64188496988670089e-1,
30 8.88314979438837594e+0,
31 6.61191906371416295e+1,
32 2.98635138197400131e+2,
33 8.81952221241769090e+2,
34 1.71204761263407058e+3,
35 2.05107837782607147e+3,
36 1.23033935479799725e+3,
37 2.15311535474403846e-8 };
38
39static const double d[8] = { 1.57449261107098347e+1,
40 1.17693950891312499e+2,
41 5.37181101862009858e+2,
42 1.62138957456669019e+3,
43 3.29079923573345963e+3,
44 4.36261909014324716e+3,
45 3.43936767414372164e+3,
46 1.23033935480374942e+3 };
47
48static const double pp[6] = { 3.05326634961232344e-1,
49 3.60344899949804439e-1,
50 1.25781726111229246e-1,
51 1.60837851487422766e-2,
52 6.58749161529837803e-4,
53 1.63153871373020978e-2 };
54
55static const double qq[5] = { 2.56852019228982242e+0,
56 1.87295284992346047e+0,
57 5.27905102951428412e-1,
58 6.05183413124413191e-2,
59 2.33520497626869185e-3 };
60
61
62static inline double ErrFunApprox ( double arg, double result, int which ) __attribute__ ((always_inline));
63static inline double ErrFunApprox ( double arg, double result, int which )
64{
65 register int i;
66 register double x, y, ysquared, numerator, denominator, del;
67
68 x = arg;
69 y = __builtin_fabs( x );
70
71/*******************************************************************************
72* Evaluate erfc for |x| <= 0.46875. *
73*******************************************************************************/
74
75 if ( y <= 0.46875e+0 )
76 {
77 if ( y > 1.11e-16 )
78 {
79 ysquared = y * y;
80 numerator = a[4] * ysquared;
81 denominator = ysquared;
82 for ( i = 0; i < 3; i++ )
83 {
84 numerator = ( numerator + a[i] ) * ysquared;
85 denominator = ( denominator + b[i] ) * ysquared;
86 }
87 result = y * ( numerator + a[3] ) / ( denominator + b[3] );
88 }
89 else
90 result = y * a[3] / b[3];
91
92 if ( which )
93 result = 1.0 - result;
94
95 return result;
96 }
97
98/*******************************************************************************
99* Evaluate erfc for 0.46875 < |x| <= 4.0 *
100*******************************************************************************/
101
102 else if ( y <= 4.0 )
103 {
104 numerator = ccc[8] * y;
105 denominator = y;
106 for ( i = 0; i < 7; i++ )
107 {
108 numerator = ( numerator + ccc[i] ) * y;
109 denominator = ( denominator + d[i] ) * y;
110 }
111 result = ( numerator + ccc[7] ) / ( denominator + d[7] );
112 ysquared = trunc ( y * 16.0 ) / 16.0;
113 del = ( y - ysquared ) * ( y + ysquared );
114 result = exp ( - ysquared * ysquared ) * exp ( - del ) * result;
115 }
116
117/*******************************************************************************
118* Evaluate erfc for |x| > 4.0 *
119*******************************************************************************/
120
121 else
122 {
123 if ( y >= xxbig )
124 {
125 if ( ( which != 2 ) || ( y >= Maximum ) )
126 {
127 if ( which == 1 )
128 {
129 double oldresult = result;
130 result *= 0x1.0000000000001p-1022;
131 result *= 0x1.0000000000001p-1022;
132 result *= 0x1.0000000000001p-1022; //set underflow
133 result += oldresult;
134 }
135
136 return result;
137 }
138 if ( y >= _HUGE )
139 {
140 result = InvSqrtPI / y;
141 return result;
142 }
143 }
144 ysquared = 1.0 / ( y * y );
145 numerator = pp[5] * ysquared;
146 denominator = ysquared;
147 for ( i = 0; i < 4; i++ )
148 {
149 numerator = ( numerator + pp[i] ) * ysquared;
150 denominator = ( denominator + qq[i] ) * ysquared;
151 }
152 result = ysquared * ( numerator + pp[4] ) / ( denominator + qq[4] );
153 result = ( InvSqrtPI - result ) / y;
154 ysquared = trunc ( y * 16.0 ) / 16.0;
155 del = ( y - ysquared ) * ( y + ysquared );
156 result = exp ( - ysquared * ysquared ) * exp ( - del ) * result;
157 }
158
159/*******************************************************************************
160* if the calling function is erf instead of erfc, take care of the *
161* underserved underflow. otherwise, the computation will produce the *
162* exception for erfc. *
163*******************************************************************************/
164
165
166 return ( which ) ? result : ( 0.5 - result ) + 0.5;
167}