ocean.math.ErrorFunction source code

1 /**
2  * Error Functions and Normal Distribution.
3  *
4  * Copyright:
5  *     Copyright (C) 1984, 1995, 2000 Stephen L. Moshier
6  *     Code taken from the Cephes Math Library Release 2.3:  January, 1995
7  *     Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH.
8  *     All rights reserved.
9  *
10  * License:
11  *     Tango Dual License: 3-Clause BSD License / Academic Free License v3.0.
12  *     See LICENSE_TANGO.txt for details.
13  *
14  * Authors: Stephen L. Moshier, ported to D by Don Clugston
15  *
16  */
17 /**
18  * Macros:
19  *  NAN = $(RED NAN)
20  *  SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
21  *  GAMMA =  &#915;
22  *  INTEGRAL = &#8747;
23  *  INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
24  *  POWER = $1<sup>$2</sup>
25  *  BIGSUM = $(BIG &Sigma; <sup>$2</sup><sub>$(SMALL $1)</sub>)
26  *  CHOOSE = $(BIG &#40;) <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG &#41;)
27  *  TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
28  *      <caption>Special Values</caption>
29  *      $0</table>
30  *  SVH = $(TR $(TH $1) $(TH $2))
31  *  SV  = $(TR $(TD $1) $(TD $2))
32  */
33 module ocean.math.ErrorFunction;
34 
35 import ocean.math.Math;
36 import ocean.math.IEEE;
37 import ocean.core.Verify;
38 
39 version (unittest)
40 {
41     import ocean.core.Test;
42 }
43 
44 
45 static immutable real SQRT2PI = 0x1.40d931ff62705966p+1L;    // 2.5066282746310005024
46 static immutable real EXP_2  = 0.13533528323661269189L; /* exp(-2) */
47 
48 private {
49 
50 /* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
51    1/8 <= 1/x <= 1
52    Peak relative error 5.8e-21  */
53 static immutable real [] P = [ -0x1.30dfa809b3cc6676p-17, 0x1.38637cd0913c0288p+18,
54    0x1.2f015e047b4476bp+22, 0x1.24726f46aa9ab08p+25, 0x1.64b13c6395dc9c26p+27,
55    0x1.294c93046ad55b5p+29, 0x1.5962a82f92576dap+30, 0x1.11a709299faba04ap+31,
56    0x1.11028065b087be46p+31, 0x1.0d8ef40735b097ep+30
57 ];
58 
59 static immutable real [] Q = [ 0x1.14d8e2a72dec49f4p+19, 0x1.0c880ff467626e1p+23,
60    0x1.04417ef060b58996p+26, 0x1.404e61ba86df4ebap+28, 0x1.0f81887bc82b873ap+30,
61    0x1.4552a5e39fb49322p+31, 0x1.11779a0ceb2a01cep+32, 0x1.3544dd691b5b1d5cp+32,
62    0x1.a91781f12251f02ep+31, 0x1.0d8ef3da605a1c86p+30, 1.0
63 ];
64 
65 
66 /* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
67    1/128 <= 1/x < 1/8
68    Peak relative error 1.9e-21  */
69 static immutable real [] R = [ 0x1.b9f6d8b78e22459ep-6, 0x1.1b84686b0a4ea43ap-1,
70    0x1.b8f6aebe96000c2ap+1, 0x1.cb1dbedac27c8ec2p+2, 0x1.cf885f8f572a4c14p+1
71 ];
72 
73 static immutable real [] S = [
74    0x1.87ae3cae5f65eb5ep-5, 0x1.01616f266f306d08p+0, 0x1.a4abe0411eed6c22p+2,
75    0x1.eac9ce3da600abaap+3, 0x1.5752a9ac2faebbccp+3, 1.0
76 ];
77 
78 /* erf(x)  = x P(x^2)/Q(x^2)
79    0 <= x <= 1
80    Peak relative error 7.6e-23  */
81 static immutable real [] T = [ 0x1.0da01654d757888cp+20, 0x1.2eb7497bc8b4f4acp+17,
82    0x1.79078c19530f72a8p+15, 0x1.4eaf2126c0b2c23p+11, 0x1.1f2ea81c9d272a2ep+8,
83    0x1.59ca6e2d866e625p+2, 0x1.c188e0b67435faf4p-4
84 ];
85 
86 static immutable real [] U = [ 0x1.dde6025c395ae34ep+19, 0x1.c4bc8b6235df35aap+18,
87    0x1.8465900e88b6903ap+16, 0x1.855877093959ffdp+13, 0x1.e5c44395625ee358p+9,
88    0x1.6a0fed103f1c68a6p+5, 1.0
89 ];
90 
91 }
92 
93 /**
94  *  Complementary error function
95  *
96  * erfc(x) = 1 - erf(x), and has high relative accuracy for
97  * values of x far from zero. (For values near zero, use erf(x)).
98  *
99  *  1 - erf(x) =  2/ $(SQRT)(&pi;)
100  *     $(INTEGRAL x, $(INFINITY)) exp( - $(POWER t, 2)) dt
101  *
102  *
103  * For small x, erfc(x) = 1 - erf(x); otherwise rational
104  * approximations are computed.
105  *
106  * A special function expx2(x) is used to suppress error amplification
107  * in computing exp(-x^2).
108  */
109 real erfc(real a)
110 {
111     if (a == real.infinity)
112         return 0.0;
113     if (a == -real.infinity)
114         return 2.0;
115 
116     real x;
117 
118     if (a < 0.0L )
119         x = -a;
120     else
121         x = a;
122     if (x < 1.0)
123         return 1.0 - erf(a);
124 
125     real z = -a * a;
126 
127     if (z < -MAXLOG){
128 //    mtherr( "erfcl", UNDERFLOW );
129         if (a < 0) return 2.0;
130         else return 0.0;
131     }
132 
133     /* Compute z = exp(z).  */
134     z = expx2(a, -1);
135     real y = 1.0/x;
136 
137     real p, q;
138 
139     if( x < 8.0 ) y = z * rationalPoly(y, P, Q);
140     else          y = z * y * rationalPoly(y * y, R, S);
141 
142     if (a < 0.0L)
143         y = 2.0L - y;
144 
145     if (y == 0.0) {
146 //    mtherr( "erfcl", UNDERFLOW );
147         if (a < 0) return 2.0;
148         else return 0.0;
149     }
150 
151     return y;
152 }
153 
154 
155 private {
156 /* Exponentially scaled erfc function
157    exp(x^2) erfc(x)
158    valid for x > 1.
159    Use with normalDistribution and expx2.  */
160 
161 real erfce(real x)
162 {
163     real p, q;
164 
165     real y = 1.0/x;
166 
167     if (x < 8.0) {
168         return rationalPoly( y, P, Q);
169     } else {
170         return y * rationalPoly(y*y, R, S);
171     }
172 }
173 
174 }
175 
176 /**
177  *  Error function
178  *
179  * The integral is
180  *
181  *  erf(x) =  2/ $(SQRT)(&pi;)
182  *     $(INTEGRAL 0, x) exp( - $(POWER t, 2)) dt
183  *
184  * The magnitude of x is limited to about 106.56 for IEEE 80-bit
185  * arithmetic; 1 or -1 is returned outside this range.
186  *
187  * For 0 <= |x| < 1, a rational polynomials are used; otherwise
188  * erf(x) = 1 - erfc(x).
189  *
190  * ACCURACY:
191  *                      Relative error:
192  * arithmetic   domain     # trials      peak         rms
193  *    IEEE      0,1         50000       2.0e-19     5.7e-20
194  */
195 real erf(real x)
196 {
197     if (x == 0.0)
198         return x; // deal with negative zero
199     if (x == -real.infinity)
200         return -1.0;
201     if (x == real.infinity)
202         return 1.0;
203     if (abs(x) > 1.0L)
204         return 1.0L - erfc(x);
205 
206     real z = x * x;
207     return x * rationalPoly(z, T, U);
208 }
209 
210 unittest {
211    // High resolution test points.
212     static immutable real erfc0_250 = 0.723663330078125 + 1.0279753638067014931732235184287934646022E-5;
213     static immutable real erfc0_375 = 0.5958709716796875 + 1.2118885490201676174914080878232469565953E-5;
214     static immutable real erfc0_500 = 0.4794921875 + 7.9346869534623172533461080354712635484242E-6;
215     static immutable real erfc0_625 = 0.3767547607421875 + 4.3570693945275513594941232097252997287766E-6;
216     static immutable real erfc0_750 = 0.2888336181640625 + 1.0748182422368401062165408589222625794046E-5;
217     static immutable real erfc0_875 = 0.215911865234375 + 1.3073705765341685464282101150637224028267E-5;
218     static immutable real erfc1_000 = 0.15728759765625 + 1.1609394035130658779364917390740703933002E-5;
219     static immutable real erfc1_125 = 0.111602783203125 + 8.9850951672359304215530728365232161564636E-6;
220 
221     static immutable real erf0_875  = (1-0.215911865234375) - 1.3073705765341685464282101150637224028267E-5;
222 
223 
224     test(feqrel(erfc(0.250L), erfc0_250 )>=real.mant_dig-1);
225     test(feqrel(erfc(0.375L), erfc0_375 )>=real.mant_dig-0);
226     test(feqrel(erfc(0.500L), erfc0_500 )>=real.mant_dig-1);
227     test(feqrel(erfc(0.625L), erfc0_625 )>=real.mant_dig-1);
228     test(feqrel(erfc(0.750L), erfc0_750 )>=real.mant_dig-1);
229     test(feqrel(erfc(0.875L), erfc0_875 )>=real.mant_dig-4);
230     version(FailsOnLinux) test(feqrel(erfc(1.000L), erfc1_000 )>=real.mant_dig-0);
231     test(feqrel(erfc(1.125L), erfc1_125 )>=real.mant_dig-2);
232     test(feqrel(erf(0.875L), erf0_875 )>=real.mant_dig-1);
233     // The DMC implementation of erfc() fails this next test (just)
234     test(feqrel(erfc(4.1L),0.67000276540848983727e-8L)>=real.mant_dig-4);
235 
236     test(isIdentical(erf(0.0),0.0));
237     test(isIdentical(erf(-0.0),-0.0));
238     test(erf(real.infinity) == 1.0);
239     test(erf(-real.infinity) == -1.0);
240     test(isIdentical(erf(NaN(0xDEF)),NaN(0xDEF)));
241     test(isIdentical(erfc(NaN(0xDEF)),NaN(0xDEF)));
242     test(isIdentical(erfc(real.infinity),0.0));
243     test(erfc(-real.infinity) == 2.0);
244     test(erfc(0) == 1.0);
245 }
246 
247 /*
248  *  Exponential of squared argument
249  *
250  * Computes y = exp(x*x) while suppressing error amplification
251  * that would ordinarily arise from the inexactness of the
252  * exponential argument x*x.
253  *
254  * If sign < 0, the result is inverted; i.e., y = exp(-x*x) .
255  *
256  * ACCURACY:
257  *                      Relative error:
258  * arithmetic      domain        # trials      peak         rms
259  *   IEEE     -106.566, 106.566    10^5       1.6e-19     4.4e-20
260  */
261 
262 real expx2(real x, int sign)
263 {
264     /*
265     Cephes Math Library Release 2.9:  June, 2000
266     Copyright 2000 by Stephen L. Moshier
267     */
268     static immutable real M = 32768.0;
269     static immutable real MINV = 3.0517578125e-5L;
270 
271     x = abs(x);
272     if (sign < 0)
273         x = -x;
274 
275   /* Represent x as an exact multiple of M plus a residual.
276      M is a power of 2 chosen so that exp(m * m) does not overflow
277      or underflow and so that |x - m| is small.  */
278     real m = MINV * floor(M * x + 0.5L);
279     real f = x - m;
280 
281     /* x^2 = m^2 + 2mf + f^2 */
282     real u = m * m;
283     real u1 = 2 * m * f  +  f * f;
284 
285     if (sign < 0) {
286         u = -u;
287         u1 = -u1;
288     }
289 
290     if ((u+u1) > MAXLOG)
291         return real.infinity;
292 
293     /* u is exact, u1 is small.  */
294     return exp(u) * exp(u1);
295 }
296 
297 
298 package {
299 /*
300 Computes the normal distribution function.
301 
302 The normal (or Gaussian, or bell-shaped) distribution is
303 defined as:
304 
305 normalDist(x) = 1/$(SQRT) &pi; $(INTEGRAL -$(INFINITY), x) exp( - $(POWER t, 2)/2) dt
306     = 0.5 + 0.5 * erf(x/sqrt(2))
307     = 0.5 * erfc(- x/sqrt(2))
308 
309 To maintain accuracy at high values of x, use
310 normalDistribution(x) = 1 - normalDistribution(-x).
311 
312 Accuracy:
313 Within a few bits of machine resolution over the entire
314 range.
315 
316 References:
317 $(LINK http://www.netlib.org/cephes/ldoubdoc.html),
318 G. Marsaglia, "Evaluating the Normal Distribution",
319 Journal of Statistical Software <b>11</b>, (July 2004).
320 */
321 real normalDistributionImpl(real a)
322 {
323     real x = a * SQRT1_2;
324     real z = abs(x);
325 
326     if( z < 1.0 )
327         return 0.5L + 0.5L * erf(x);
328     else {
329         /* See below for erfce. */
330         real y = 0.5L * erfce(z);
331         /* Multiply by exp(-x^2 / 2)  */
332         z = expx2(a, -1);
333         y = y * sqrt(z);
334         if( x > 0.0L )
335             y = 1.0L - y;
336         return y;
337     }
338 }
339 
340 }
341 
342 unittest {
343     test(fabs(normalDistributionImpl(1L) - (0.841344746068543))< 0.0000000000000005);
344     test(isIdentical(normalDistributionImpl(NaN(0x325)), NaN(0x325)));
345 }
346 
347 package {
348 /*
349  * Inverse of Normal distribution function
350  *
351  * Returns the argument, x, for which the area under the
352  * Normal probability density function (integrated from
353  * minus infinity to x) is equal to p.
354  *
355  * For small arguments 0 < p < exp(-2), the program computes
356  * z = sqrt( -2 log(p) );  then the approximation is
357  * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z) .
358  * For larger arguments,  x/sqrt(2 pi) = w + w^3 R(w^2)/S(w^2)) ,
359  * where w = p - 0.5 .
360  */
361 real normalDistributionInvImpl(real p)
362 {
363 verify(p>=0.0L && p<=1.0L, "Domain error");
364 
365 static immutable real[] P0 = [ -0x1.758f4d969484bfdcp-7, 0x1.53cee17a59259dd2p-3,
366    -0x1.ea01e4400a9427a2p-1,  0x1.61f7504a0105341ap+1, -0x1.09475a594d0399f6p+2,
367     0x1.7c59e7a0df99e3e2p+1, -0x1.87a81da52edcdf14p-1,  0x1.1fb149fd3f83600cp-7
368 ];
369 
370 static immutable real[] Q0 = [ -0x1.64b92ae791e64bb2p-7, 0x1.7585c7d597298286p-3,
371    -0x1.40011be4f7591ce6p+0, 0x1.1fc067d8430a425ep+2, -0x1.21008ffb1e7ccdf2p+3,
372    0x1.3d1581cf9bc12fccp+3, -0x1.53723a89fd8f083cp+2, 1.0
373 ];
374 
375 static immutable real[] P1 = [ 0x1.20ceea49ea142f12p-13, 0x1.cbe8a7267aea80bp-7,
376    0x1.79fea765aa787c48p-2, 0x1.d1f59faa1f4c4864p+1, 0x1.1c22e426a013bb96p+4,
377    0x1.a8675a0c51ef3202p+5, 0x1.75782c4f83614164p+6, 0x1.7a2f3d90948f1666p+6,
378    0x1.5cd116ee4c088c3ap+5, 0x1.1361e3eb6e3cc20ap+2
379 ];
380 
381 static immutable real[] Q1 = [ 0x1.3a4ce1406cea98fap-13, 0x1.f45332623335cda2p-7,
382    0x1.98f28bbd4b98db1p-2, 0x1.ec3b24f9c698091cp+1, 0x1.1cc56ecda7cf58e4p+4,
383    0x1.92c6f7376bf8c058p+5, 0x1.4154c25aa47519b4p+6, 0x1.1b321d3b927849eap+6,
384    0x1.403a5f5a4ce7b202p+4, 1.0
385 ];
386 
387 static immutable real[] P2 = [ 0x1.8c124a850116a6d8p-21, 0x1.534abda3c2fb90bap-13,
388    0x1.29a055ec93a4718cp-7, 0x1.6468e98aad6dd474p-3, 0x1.3dab2ef4c67a601cp+0,
389    0x1.e1fb3a1e70c67464p+1, 0x1.b6cce8035ff57b02p+2, 0x1.9f4c9e749ff35f62p+1
390 ];
391 
392 static immutable real[] Q2 = [ 0x1.af03f4fc0655e006p-21, 0x1.713192048d11fb2p-13,
393    0x1.4357e5bbf5fef536p-7, 0x1.7fdac8749985d43cp-3, 0x1.4a080c813a2d8e84p+0,
394    0x1.c3a4b423cdb41bdap+1, 0x1.8160694e24b5557ap+2, 1.0
395 ];
396 
397 static immutable real[] P3 = [ -0x1.55da447ae3806168p-34, -0x1.145635641f8778a6p-24,
398  -0x1.abf46d6b48040128p-17, -0x1.7da550945da790fcp-11, -0x1.aa0b2a31157775fap-8,
399    0x1.b11d97522eed26bcp-3, 0x1.1106d22f9ae89238p+1, 0x1.029a358e1e630f64p+1
400 ];
401 
402 static immutable real[] Q3 = [ -0x1.74022dd5523e6f84p-34, -0x1.2cb60d61e29ee836p-24,
403    -0x1.d19e6ec03a85e556p-17, -0x1.9ea2a7b4422f6502p-11, -0x1.c54b1e852f107162p-8,
404    0x1.e05268dd3c07989ep-3, 0x1.239c6aff14afbf82p+1, 1.0
405 ];
406 
407   if(p<=0.0L || p>=1.0L) {
408         if (p == 0.0L) {
409             return -real.infinity;
410         }
411         if( p == 1.0L ) {
412             return real.infinity;
413         }
414         return NaN(TANGO_NAN.NORMALDISTRIBUTION_INV_DOMAIN);
415     }
416     int code = 1;
417     real y = p;
418     if( y > (1.0L - EXP_2) ) {
419         y = 1.0L - y;
420         code = 0;
421     }
422 
423     real x, z, y2, x0, x1;
424 
425     if ( y > EXP_2 ) {
426         y = y - 0.5L;
427         y2 = y * y;
428         x = y + y * (y2 * rationalPoly( y2, P0, Q0));
429         return x * SQRT2PI;
430     }
431 
432     x = sqrt( -2.0L * log(y) );
433     x0 = x - log(x)/x;
434     z = 1.0L/x;
435     if ( x < 8.0L ) {
436         x1 = z * rationalPoly( z, P1, Q1);
437     } else if( x < 32.0L ) {
438         x1 = z * rationalPoly( z, P2, Q2);
439     } else {
440         x1 = z * rationalPoly( z, P3, Q3);
441     }
442     x = x0 - x1;
443     if ( code != 0 ) {
444         x = -x;
445     }
446     return x;
447 }
448 
449 }
450 
451 
452 unittest {
453     // TODO: Use verified test points.
454     // The values below are from Excel 2003.
455     test(fabs(normalDistributionInvImpl(0.001) - (-3.09023230616779))< 0.00000000000005);
456     test(fabs(normalDistributionInvImpl(1e-50) - (-14.9333375347885))< 0.00000000000005);
457     test(feqrel(normalDistributionInvImpl(0.999), -normalDistributionInvImpl(0.001))>real.mant_dig-6);
458 
459     // Excel 2003 gets all the following values wrong!
460     test(normalDistributionInvImpl(0.0)==-real.infinity);
461     test(normalDistributionInvImpl(1.0)==real.infinity);
462     test(normalDistributionInvImpl(0.5)==0);
463     // (Excel 2003 returns norminv(p) = -30 for all p < 1e-200).
464     // The value tested here is the one the function returned in Jan 2006.
465     real unknown1 = normalDistributionInvImpl(1e-250L);
466     test( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005);
467 }