overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22/*******************************************************************************
23* *
24* File erfcerf.c, *
25* Functions erf(x) and cerf(x), *
26* Implementation of the error and complementary error function for the *
27* PowerPC. *
28* *
29* Copyright c 1993 Apple Computer, Inc. All rights reserved. *
30* *
31* Written by Ali Sazegari, started on February 1993, *
32* *
33* based on FORTRAN routines in SpecFun by J. W. Cody of Applied *
34* Mathematics Division of Argonne National Laboratory, Argonne, IL *
35* 60439. *
36* *
37* W A R N I N G: This routine expects a 64 bit double model. *
38* *
39* February 22 1993: First C implementation for PowerPC. *
40* July 14 1993: added #pragma fenv_access, declared variables *
41* register automatics, changed feholdexcept to the *
42* internal MathLib function _feprocentry. *
43* October 07 1993: corrected the return of signed zero in erf. *
44* September 19 1994: replaced the environmental enquiries by __setflm. *
45* October 31 1995: corrected an error in the computation of erfc in *
46* the interval of [-0.46875,0.0]. changed the value *
47* of xbig to disallow flush to zero. *
48* November 01 1995: corrected an undeserved underfow flag in erf. *
49* *
50********************************************************************************
51* *
52* E R F ( X ) & C E R F ( X ) *
53* *
54********************************************************************************
55* *
56* Explanation of machine-dependent constants: *
57* *
58* xmin = the smallest positive floating-point number. *
59* xneg = the largest negative argument acceptable to ERFCX; *
60* the negative of the solution to the equation *
61* 2 * exp ( x * x ) = HUGE_VAL. *
62* xsmall = argument below which erf(x) may be represented by *
63* 2 * x / sqrt ( pi ) and above which x * x will not *
64* underflow. A conservative value is the largest machine *
65* number x such that 1.0 + x = 1.0 to machine precision. *
66* xbig = largest argument acceptable to erfc; solution to *
67* the equation: W(x) * ( 1 - 0.5 / x ** 2 ) = xmin, where *
68* W(x) = exp ( - x * x ) / ( x * sqrt ( pi ) ). *
69* HUGE = argument above which 1 - 1 / ( 2 * x * x ) = 1 to *
70* machine precision. A conservative value is *
71* 1 / ( 2 * sqrt ( xsmall ) ) . *
72* Maximum = largest acceptable argument to erfcx; the minimum *
73* of HUGE_VAL and 1 / ( sqrt ( pi ) * xmin ). *
74* *
75* Approximate values for the macintosh and the powerpc are: *
76* *
77* xmin xinf xneg xsmall *
78* *
79* Macintosh (E.P.) *
80* PowerPC (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 *
81* *
82* xbig HUGE Maximum *
83* *
84* Macintosh (E.P.) *
85* PowerPC (D.P.) 26.543 6.71D+7 2.53D+307 *
86* if not flush to 0 27.2 *
87* *
88********************************************************************************
89* *
90* Functions required are: *
91* *
92* Transandental: exp; *
93* Auxiliary: trunc, fabs and __fpclassifyd; *
94* Environmental: feholdexcept and feupdateenv. *
95* *
96* Reference: *
97* *
98* The main computation evaluates near-minimax approximations *
99* from "Rational Chebyshev approximations for the error function" *
100* by W. J. Cody, Math. Comp., 1969, PP. 631-638. *
101* *
102* This program uses rational functions that theoretically approximate *
103* erf(x) and erfc(x) to at least 18 significant decimal digits. *
104* *
105*******************************************************************************/
106
107#include "math.h"
108#include "fenv.h"
109#include "fp_private.h"
110#include "fenv_private.h"
111
112#define EXCEPT_MASK 0x1ff80000
113
114/*******************************************************************************
115* Functions needed for the computation. *
116*******************************************************************************/
117
118/* the following fp.h functions are used: */
119/* __fpclassifyd, copysign, trunc, fabs and exp; */
120/* the following environmental functions are used: */
121/* __setflm. */
122
123/*******************************************************************************
124* Coefficients for approximation to erf in [ -0.46875, 0.46875] in *
125* decreasing order. *
126*******************************************************************************/
127
128static const double a[5] = { 3.16112374387056560e+0,
129 1.13864154151050156e+2,
130 3.77485237685302021e+2,
131 3.20937758913846947e+3,
132 1.85777706184603153e-1 };
133
134static const double b[4] = { 2.36012909523441209e+1,
135 2.44024637934444173e+2,
136 1.28261652607737228e+3,
137 2.84423683343917062e+3 };
138
139/*******************************************************************************
140* Coefficients for approximation to erfc in [-4.0,-0.46875)U(0.46875,4.0] *
141* in decreasing order. *
142*******************************************************************************/
143
144static const double c[9] = { 5.64188496988670089e-1,
145 8.88314979438837594e+0,
146 6.61191906371416295e+1,
147 2.98635138197400131e+2,
148 8.81952221241769090e+2,
149 1.71204761263407058e+3,
150 2.05107837782607147e+3,
151 1.23033935479799725e+3,
152 2.15311535474403846e-8 };
153
154static const double d[8] = { 1.57449261107098347e+1,
155 1.17693950891312499e+2,
156 5.37181101862009858e+2,
157 1.62138957456669019e+3,
158 3.29079923573345963e+3,
159 4.36261909014324716e+3,
160 3.43936767414372164e+3,
161 1.23033935480374942e+3 };
162
163/*******************************************************************************
164* Coefficients for approximation to erfc in [-inf,-4.0)U(4.0,inf] in *
165* decreasing order. *
166*******************************************************************************/
167
168static const double p[6] = { 3.05326634961232344e-1,
169 3.60344899949804439e-1,
170 1.25781726111229246e-1,
171 1.60837851487422766e-2,
172 6.58749161529837803e-4,
173 1.63153871373020978e-2 };
174
175static const double q[5] = { 2.56852019228982242e+0,
176 1.87295284992346047e+0,
177 5.27905102951428412e-1,
178 6.05183413124413191e-2,
179 2.33520497626869185e-3 };
180
181static const double InvSqrtPI = 5.6418958354775628695e-1;
182static const double xbig = 27.2e+0;
183static const double Maximum = 2.53e+307;
184static const double _HUGE = 6.71e+7;
185
186#pragma STDC FENV_ACCESS ON
187
188static double ErrFunApprox ( double arg, double result, int which );
189
190/*******************************************************************************
191* E R R O R F U N C T I O N *
192*******************************************************************************/
193
194double erf ( double x )
195 {
196 register int which = 0;
197 register double result = 0.0;
198 hexdouble OldEnvironment, NewEnvironment;
199
200/*******************************************************************************
201* The next switch will decipher what sort of argument we have. If argument *
202* is SNaN then a QNaN has to be returned and the invalid flag signaled. *
203*******************************************************************************/
204
205 switch ( __fpclassifyd ( x ) )
206 {
207 case FP_NAN:
208 x *= 2.0; /* Sets invalid & quiets NaN */
209 return x;
210 case FP_ZERO:
211 return x;
212 case FP_INFINITE:
213 return (x > 0.0 ? 1.0 : - 1.0);
214 default: /* NORMALNUM and DENORMALNUM */
215 break;
216 }
217
218 FEGETENVD( OldEnvironment.d ); // save environment, set default
219 FESETENVD( 0.0 );
220
221 result = 1.0;
222 result = ErrFunApprox ( x, result, which );
223
224/*******************************************************************************
225* Take care of the negative argument. *
226*******************************************************************************/
227
228 result = copysign ( result, x);
229
230 FEGETENVD( NewEnvironment.d );
231 OldEnvironment.i.lo |= ( NewEnvironment.i.lo & EXCEPT_MASK ); // Merge new exceptions into old environment
232 FESETENVD( OldEnvironment.d ); // restore caller's environment
233
234 return ( result );
235 }
236
237/*******************************************************************************
238* C O M P L E M E N T A R Y E R R O R F U N C T I O N *
239*******************************************************************************/
240
241double erfc ( double x )
242 {
243 register int which = 1;
244 register double result = 0.0;
245 hexdouble OldEnvironment, NewEnvironment;
246
247
248/*******************************************************************************
249* The next switch will decipher what sort of argument we have. If argument *
250* is SNaN then a QNaN has to be returned and the invalid flag signaled. *
251*******************************************************************************/
252
253 switch ( __fpclassifyd ( x ) )
254 {
255 case FP_NAN:
256 x *= 2.0; /* Sets invalid & quiets NaN */
257 return x;
258 case FP_ZERO:
259 return 1.0;
260 case FP_INFINITE:
261 return (x > 0.0 ? 0.0 : 2.0);
262 default: /* NORMALNUM and DENORMALNUM */
263 break;
264 }
265
266 FEGETENVD( OldEnvironment.d ); // save environment, set default
267 FESETENVD( 0.0 );
268
269 result = 0.0;
270 result = ErrFunApprox ( x, result, which );
271
272/*******************************************************************************
273* Take care of the negative argument. *
274*******************************************************************************/
275
276 if ( x < 0.0 )
277 result = 2.0 - result;
278
279 FEGETENVD( NewEnvironment.d );
280 OldEnvironment.i.lo |= ( NewEnvironment.i.lo & EXCEPT_MASK ); // Merge new exceptions into old environment
281 FESETENVD( OldEnvironment.d ); // restore caller's environment
282
283 return ( result );
284 }
285
286
287/*******************************************************************************
288* C O R E A P P R O X I M A T I O N *
289*******************************************************************************/
290
291static double ErrFunApprox ( double arg, double result, int which )
292 {
293 register int i;
294 register double x, y, ysquared, numerator, denominator, del;
295
296 x = arg;
297 y = fabs ( x );
298
299/*******************************************************************************
300* Evaluate erfc for |x| <= 0.46875. *
301*******************************************************************************/
302
303 if ( y <= 0.46875e+0 )
304 {
305 ysquared = 0.0;
306 if ( y > 1.11e-16 )
307 ysquared = y * y;
308 numerator = a[4] * ysquared;
309 denominator = ysquared;
310 for ( i = 0; i < 3; i++ )
311 {
312 numerator = ( numerator + a[i] ) * ysquared;
313 denominator = ( denominator + b[i] ) * ysquared;
314 }
315 result = y * ( numerator + a[3] ) / ( denominator + b[3] );
316 if ( which ) result = 1.0 - result;
317 return result;
318 }
319
320/*******************************************************************************
321* Evaluate erfc for 0.46875 < |x| <= 4.0 *
322*******************************************************************************/
323
324 else if ( y <= 4.0 )
325 {
326 numerator = c[8] * y;
327 denominator = y;
328 for ( i = 0; i < 7; i++ )
329 {
330 numerator = ( numerator + c[i] ) * y;
331 denominator = ( denominator + d[i] ) * y;
332 }
333 result = ( numerator + c[7] ) / ( denominator + d[7] );
334 ysquared = trunc ( y * 16.0 ) / 16.0;
335 del = ( y - ysquared ) * ( y + ysquared );
336 result = exp ( - ysquared * ysquared ) * exp ( - del ) * result;
337 }
338
339/*******************************************************************************
340* Evaluate erfc for |x| > 4.0 *
341*******************************************************************************/
342
343 else
344 {
345 if ( y >= xbig )
346 {
347 if ( ( which != 2 ) || ( y >= Maximum ) )
348 {
349 if ( which == 1 )
350 {
351 hexdouble OldEnvironment;
352 FEGETENVD( OldEnvironment.d );
353 OldEnvironment.i.lo |= FE_UNDERFLOW;
354 FESETENVD( OldEnvironment.d );
355 }
356 return result;
357 }
358 if ( y >= _HUGE )
359 {
360 result = InvSqrtPI / y;
361 return result;
362 }
363 }
364 ysquared = 1.0 / ( y * y );
365 numerator = p[5] * ysquared;
366 denominator = ysquared;
367 for ( i = 0; i < 4; i++ )
368 {
369 numerator = ( numerator + p[i] ) * ysquared;
370 denominator = ( denominator + q[i] ) * ysquared;
371 }
372 result = ysquared * ( numerator + p[4] ) / ( denominator + q[4] );
373 result = ( InvSqrtPI - result ) / y;
374 ysquared = trunc ( y * 16.0 ) / 16.0;
375 del = ( y - ysquared ) * ( y + ysquared );
376 result = exp ( - ysquared * ysquared ) * exp ( - del ) * result;
377 }
378
379/*******************************************************************************
380* if the calling function is erf instead of erfc, take care of the *
381* underserved underflow. otherwise, the computation will produce the *
382* exception for erfc. *
383*******************************************************************************/
384
385 if ( which == 0 )
386 {
387 hexdouble OldEnvironment;
388 FEGETENVD( OldEnvironment.d );
389 OldEnvironment.i.lo &= ~FE_UNDERFLOW;
390 FESETENVD( OldEnvironment.d );
391 }
392
393 return ( which ) ? result : ( 0.5 - result ) + 0.5;
394 }