Mercurial > projects > ldc

comparison tango/tango/math/ErrorFunction.d @ 132:1700239cab2e trunk

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