Source code: graph/SpecialFunction.java
1 package graph;
2
3 import java.lang.Math;
4 import java.lang.ArithmeticException;
5
6 /*
7 **************************************************************************
8 **
9 ** Class SpecialFunction
10 **
11 **************************************************************************
12 ** Copyright (C) 1996 Leigh Brookshaw
13 **
14 ** This program is free software; you can redistribute it and/or modify
15 ** it under the terms of the GNU General Public License as published by
16 ** the Free Software Foundation; either version 2 of the License, or
17 ** (at your option) any later version.
18 **
19 ** This program is distributed in the hope that it will be useful,
20 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
21 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 ** GNU General Public License for more details.
23 **
24 ** You should have received a copy of the GNU General Public License
25 ** along with this program; if not, write to the Free Software
26 ** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
27 **************************************************************************
28 ** Modified by Jim Cochrane, last changed in July, 2000
29 **************************************************************************
30 **
31 ** This class is an extension of java.lang.Math. It includes a number
32 ** of special functions not found in the Math class.
33 **
34 *************************************************************************/
35
36
37 /**
38 * This class contains physical constants and special functions not found
39 * in the java.lang.Math class.
40 * Like the java.lang.Math class this class is final and cannot be
41 * subclassed.
42 * All physical constants are in cgs units.
43 * <P>
44 * <B>NOTE:</B> These special functions do not necessarily use the fastest
45 * or most accurate algorithms.
46 *
47 * @version $Revision: 2.4 $, $Date: 2001/05/29 09:36:38 $
48 * @author Leigh Brookshaw
49 */
50
51
52 public final class SpecialFunction extends Object {
53
54 /*
55 ** machine constants
56 */
57 private static final double MACHEP = 1.11022302462515654042E-16;
58 private static final double MAXLOG = 7.09782712893383996732E2;
59 private static final double MINLOG = -7.451332191019412076235E2;
60 private static final double MAXGAM = 171.624376956302725;
61 private static final double SQTPI = 2.50662827463100050242E0;
62 private static final double SQRTH = 7.07106781186547524401E-1;
63 private static final double LOGPI = 1.14472988584940017414;
64
65
66 /*
67 ** Physical Constants in cgs Units
68 */
69
70
71 /**
72 * Boltzman Constant. Units erg/deg(K)
73 */
74 public static final double BOLTZMAN = 1.3807e-16;
75 /**
76 * Elementary Charge. Units statcoulomb
77 */
78 public static final double ECHARGE = 4.8032e-10;
79 /**
80 * Electron Mass. Units g
81 */
82 public static final double EMASS = 9.1095e-28;
83 /**
84 * Proton Mass. Units g
85 */
86 public static final double PMASS = 1.6726e-24;
87 /**
88 * Gravitational Constant. Units dyne-cm^2/g^2
89 */
90 public static final double GRAV = 6.6720e-08;
91 /**
92 * Planck constant. Units erg-sec
93 */
94 public static final double PLANCK = 6.6262e-27;
95 /**
96 * Speed of Light in a Vacuum. Units cm/sec
97 */
98 public static final double LIGHTSPEED = 2.9979e10;
99 /**
100 * Stefan-Boltzman Constant. Units erg/cm^2-sec-deg^4
101 */
102 public static final double STEFANBOLTZ = 5.6703e-5;
103 /**
104 * Avogadro Number. Units 1/mol
105 */
106 public static final double AVOGADRO = 6.0220e23;
107 /**
108 * Gas Constant. Units erg/deg-mol
109 */
110 public static final double GASCONSTANT = 8.3144e07;
111 /**
112 * Gravitational Acceleration at the Earths surface. Units cm/sec^2
113 */
114 public static final double GRAVACC = 980.67;
115
116 /**
117 * Solar Mass. Units g
118 */
119 public static final double SOLARMASS = 1.99e33;
120 /**
121 * Solar Radius. Units cm
122 */
123 public static final double SOLARRADIUS = 6.96e10;
124 /**
125 * Solar Luminosity. Units erg/sec
126 */
127 public static final double SOLARLUM = 3.90e33;
128 /**
129 * Solar Flux. Units erg/cm^2-sec
130 */
131 public static final double SOLARFLUX = 6.41e10;
132 /**
133 * Astronomical Unit (radius of the Earth's orbit). Units cm
134 */
135 public static final double AU = 1.50e13;
136
137
138 /**
139 * Don't let anyone instantiate this class.
140 */
141 private SpecialFunction() {}
142
143 /*
144 ** Function Methods
145 */
146
147 /**
148 * @param x a double value
149 * @return The log<sub>10</sub>
150 */
151 static public double log10(double x) throws ArithmeticException {
152 if( x <= 0.0 ) throw new ArithmeticException("range exception");
153 return Math.log(x)/2.30258509299404568401;
154 }
155
156
157 /**
158 * @param x a double value
159 * @return the hyperbolic cosine of the argument
160 */
161
162 static public double cosh(double x) throws ArithmeticException {
163 double a;
164 a = x;
165 if( a < 0.0 ) a = Math.abs(x);
166 a = Math.exp(a);
167 return 0.5*(a+1/a);
168 }
169
170 /**
171 * @param x a double value
172 * @return the hyperbolic sine of the argument
173 */
174 static public double sinh(double x) throws ArithmeticException {
175 double a;
176 if(x == 0.0) return x;
177 a = x;
178 if( a < 0.0 ) a = Math.abs(x);
179 a = Math.exp(a);
180 if( x < 0.0 ) return -0.5*(a-1/a);
181 else return 0.5*(a-1/a);
182 }
183
184 /**
185 * @param x a double value
186 * @return the hyperbolic tangent of the argument
187 */
188 static public double tanh(double x) throws ArithmeticException {
189 double a;
190 if( x == 0.0 ) return x;
191 a = x;
192 if( a < 0.0 ) a = Math.abs(x);
193 a = Math.exp(2.0*a);
194 if(x < 0.0 ) return -( 1.0-2.0/(a+1.0) );
195 else return ( 1.0-2.0/(a+1.0) );
196 }
197
198 /**
199 * @param x a double value
200 * @return the hyperbolic arc cosine of the argument
201 */
202
203 static public double acosh(double x) throws ArithmeticException {
204 if( x < 1.0 ) throw new ArithmeticException("range exception");
205 return Math.log( x + Math.sqrt(x*x-1));
206 }
207
208 /**
209 * @param x a double value
210 * @return the hyperbolic arc sine of the argument
211 */
212 static public double asinh(double xx) throws ArithmeticException {
213 double x;
214 int sign;
215 if(xx == 0.0) return xx;
216 if( xx < 0.0 ) {
217 sign = -1;
218 x = -xx;
219 } else {
220 sign = 1;
221 x = xx;
222 }
223 return sign*Math.log( x + Math.sqrt(x*x+1));
224 }
225
226 /**
227 * @param x a double value
228 * @return the hyperbolic arc tangent of the argument
229 */
230 static public double atanh(double x) throws ArithmeticException {
231 if( x > 1.0 || x < -1.0 ) throw
232 new ArithmeticException("range exception");
233 return 0.5 * Math.log( (1.0+x)/(1.0-x) );
234 }
235
236 /**
237 * @param x a double value
238 * @return the Bessel function of order 0 of the argument.
239 */
240
241 static public double j0(double x) throws ArithmeticException {
242 double ax;
243
244 if( (ax=Math.abs(x)) < 8.0 ) {
245 double y=x*x;
246 double ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
247 +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
248 double ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
249 +y*(59272.64853+y*(267.8532712+y*1.0))));
250
251 return ans1/ans2;
252
253 } else {
254 double z=8.0/ax;
255 double y=z*z;
256 double xx=ax-0.785398164;
257 double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
258 +y*(-0.2073370639e-5+y*0.2093887211e-6)));
259 double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
260 +y*(-0.6911147651e-5+y*(0.7621095161e-6
261 -y*0.934935152e-7)));
262
263 return Math.sqrt(0.636619772/ax)*
264 (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
265 }
266 }
267 /**
268 * @param x a double value
269 * @return the Bessel function of order 1 of the argument.
270 */
271
272 static public double j1(double x) throws ArithmeticException {
273
274 double ax;
275 double y;
276 double ans1, ans2;
277
278 if ((ax=Math.abs(x)) < 8.0) {
279 y=x*x;
280 ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
281 +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
282 ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
283 +y*(99447.43394+y*(376.9991397+y*1.0))));
284 return ans1/ans2;
285 } else {
286 double z=8.0/ax;
287 double xx=ax-2.356194491;
288 y=z*z;
289
290 ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
291 +y*(0.2457520174e-5+y*(-0.240337019e-6))));
292 ans2=0.04687499995+y*(-0.2002690873e-3
293 +y*(0.8449199096e-5+y*(-0.88228987e-6
294 +y*0.105787412e-6)));
295 double ans=Math.sqrt(0.636619772/ax)*
296 (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
297 if (x < 0.0) ans = -ans;
298 return ans;
299 }
300 }
301
302 /**
303 * @param n integer order
304 * @param x a double value
305 * @return the Bessel function of order n of the argument.
306 */
307 static public double jn(int n, double x) throws ArithmeticException {
308 int j,m;
309 double ax,bj,bjm,bjp,sum,tox,ans;
310 boolean jsum;
311
312 double ACC = 40.0;
313 double BIGNO = 1.0e+10;
314 double BIGNI = 1.0e-10;
315
316 if(n == 0) return j0(x);
317 if(n == 1) return j1(x);
318
319 ax=Math.abs(x);
320 if(ax == 0.0) return 0.0;
321 else
322 if (ax > (double)n) {
323 tox=2.0/ax;
324 bjm=j0(ax);
325 bj=j1(ax);
326 for (j=1;j<n;j++) {
327 bjp=j*tox*bj-bjm;
328 bjm=bj;
329 bj=bjp;
330 }
331 ans=bj;
332 } else {
333 tox=2.0/ax;
334 m=2*((n+(int)Math.sqrt(ACC*n))/2);
335 jsum=false;
336 bjp=ans=sum=0.0;
337 bj=1.0;
338 for (j=m;j>0;j--) {
339 bjm=j*tox*bj-bjp;
340 bjp=bj;
341 bj=bjm;
342 if (Math.abs(bj) > BIGNO) {
343 bj *= BIGNI;
344 bjp *= BIGNI;
345 ans *= BIGNI;
346 sum *= BIGNI;
347 }
348 if (jsum) sum += bj;
349 jsum=!jsum;
350 if (j == n) ans=bjp;
351 }
352 sum=2.0*sum-bj;
353 ans /= sum;
354 }
355 return x < 0.0 && n%2 == 1 ? -ans : ans;
356 }
357 /**
358 * @param x a double value
359 * @return the Bessel function of the second kind,
360 * of order 0 of the argument.
361 */
362
363 static public double y0(double x) throws ArithmeticException {
364
365 if (x < 8.0) {
366 double y=x*x;
367
368 double ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
369 +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
370 double ans2=40076544269.0+y*(745249964.8+y*(7189466.438
371 +y*(47447.26470+y*(226.1030244+y*1.0))));
372
373 return (ans1/ans2)+0.636619772*j0(x)*Math.log(x);
374 } else {
375 double z=8.0/x;
376 double y=z*z;
377 double xx=x-0.785398164;
378
379 double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
380 +y*(-0.2073370639e-5+y*0.2093887211e-6)));
381 double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
382 +y*(-0.6911147651e-5+y*(0.7621095161e-6
383 +y*(-0.934945152e-7))));
384 return Math.sqrt(0.636619772/x)*
385 (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
386 }
387 }
388
389 /**
390 * @param x a double value
391 * @return the Bessel function of the second kind,
392 * of order 1 of the argument.
393 */
394 static public double y1(double x) throws ArithmeticException {
395
396 if (x < 8.0) {
397 double y=x*x;
398 double ans1=x*(-0.4900604943e13+y*(0.1275274390e13
399 +y*(-0.5153438139e11+y*(0.7349264551e9
400 +y*(-0.4237922726e7+y*0.8511937935e4)))));
401 double ans2=0.2499580570e14+y*(0.4244419664e12
402 +y*(0.3733650367e10+y*(0.2245904002e8
403 +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
404 return (ans1/ans2)+0.636619772*(j1(x)*Math.log(x)-1.0/x);
405 } else {
406 double z=8.0/x;
407 double y=z*z;
408 double xx=x-2.356194491;
409 double ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
410 +y*(0.2457520174e-5+y*(-0.240337019e-6))));
411 double ans2=0.04687499995+y*(-0.2002690873e-3
412 +y*(0.8449199096e-5+y*(-0.88228987e-6
413 +y*0.105787412e-6)));
414 return Math.sqrt(0.636619772/x)*
415 (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
416 }
417 }
418 /**
419 * @param n integer order
420 * @param x a double value
421 * @return the Bessel function of the second kind,
422 * of order n of the argument.
423 */
424 static public double yn(int n, double x) throws ArithmeticException {
425 double by,bym,byp,tox;
426
427 if(n == 0) return y0(x);
428 if(n == 1) return y1(x);
429
430 tox=2.0/x;
431 by=y1(x);
432 bym=y0(x);
433 for (int j=1;j<n;j++) {
434 byp=j*tox*by-bym;
435 bym=by;
436 by=byp;
437 }
438 return by;
439 }
440
441
442
443 /**
444 * @param x a double value
445 * @return the factorial of the argument
446 */
447 static public double fac(double x) throws ArithmeticException {
448 double d = Math.abs(x);
449 if(Math.floor(d) == d) return (double)fac( (int)x );
450 else return gamma(x+1.0);
451 }
452
453 /**
454 * @param x an integer value
455 * @return the factorial of the argument
456 */
457 static public int fac(int j) throws ArithmeticException {
458 int i = j;
459 int d = 1;
460 if(j < 0) i = Math.abs(j);
461 while( i > 1) { d *= i--; }
462 if(j < 0) return -d;
463 else return d;
464 }
465
466
467
468
469
470 /**
471 * @param x a double value
472 * @return the Gamma function of the value.
473 * <P>
474 * <FONT size=2>
475 * Converted to Java from<BR>
476 * Cephes Math Library Release 2.2: July, 1992<BR>
477 * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>
478 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
479 **/
480 static public double gamma(double x) throws ArithmeticException {
481
482 double P[] = {
483 1.60119522476751861407E-4,
484 1.19135147006586384913E-3,
485 1.04213797561761569935E-2,
486 4.76367800457137231464E-2,
487 2.07448227648435975150E-1,
488 4.94214826801497100753E-1,
489 9.99999999999999996796E-1
490 };
491 double Q[] = {
492 -2.31581873324120129819E-5,
493 5.39605580493303397842E-4,
494 -4.45641913851797240494E-3,
495 1.18139785222060435552E-2,
496 3.58236398605498653373E-2,
497 -2.34591795718243348568E-1,
498 7.14304917030273074085E-2,
499 1.00000000000000000320E0
500 };
501 double MAXGAM = 171.624376956302725;
502 double LOGPI = 1.14472988584940017414;
503
504 double p, z;
505 int i;
506
507 double q = Math.abs(x);
508
509 if( q > 33.0 ) {
510 if( x < 0.0 ) {
511 p = Math.floor(q);
512 if( p == q ) throw new ArithmeticException("gamma: overflow");
513 i = (int)p;
514 z = q - p;
515 if( z > 0.5 ) {
516 p += 1.0;
517 z = q - p;
518 }
519 z = q * Math.sin( Math.PI * z );
520 if( z == 0.0 ) throw new ArithmeticException("gamma: overflow");
521 z = Math.abs(z);
522 z = Math.PI/(z * stirf(q) );
523
524 return -z;
525 } else {
526 return stirf(x);
527 }
528 }
529
530 z = 1.0;
531 while( x >= 3.0 ) {
532 x -= 1.0;
533 z *= x;
534 }
535
536 while( x < 0.0 ) {
537 if( x == 0.0 ) {
538 throw new ArithmeticException("gamma: singular");
539 } else
540 if( x > -1.E-9 ) {
541 return( z/((1.0 + 0.5772156649015329 * x) * x) );
542 }
543 z /= x;
544 x += 1.0;
545 }
546
547 while( x < 2.0 ) {
548 if( x == 0.0 ) {
549 throw new ArithmeticException("gamma: singular");
550 } else
551 if( x < 1.e-9 ) {
552 return( z/((1.0 + 0.5772156649015329 * x) * x) );
553 }
554 z /= x;
555 x += 1.0;
556 }
557
558 if( (x == 2.0) || (x == 3.0) ) return z;
559
560 x -= 2.0;
561 p = polevl( x, P, 6 );
562 q = polevl( x, Q, 7 );
563 return z * p / q;
564
565 }
566
567 /* Gamma function computed by Stirling's formula.
568 * The polynomial STIR is valid for 33 <= x <= 172.
569
570 Cephes Math Library Release 2.2: July, 1992
571 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
572 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
573 */
574 static private double stirf(double x) throws ArithmeticException {
575 double STIR[] = {
576 7.87311395793093628397E-4,
577 -2.29549961613378126380E-4,
578 -2.68132617805781232825E-3,
579 3.47222221605458667310E-3,
580 8.33333333333482257126E-2,
581 };
582 double MAXSTIR = 143.01608;
583
584 double w = 1.0/x;
585 double y = Math.exp(x);
586
587 w = 1.0 + w * polevl( w, STIR, 4 );
588
589 if( x > MAXSTIR ) {
590 /* Avoid overflow in Math.pow() */
591 double v = Math.pow( x, 0.5 * x - 0.25 );
592 y = v * (v / y);
593 } else {
594 y = Math.pow( x, x - 0.5 ) / y;
595 }
596 y = SQTPI * y * w;
597 return y;
598 }
599
600 /**
601 * @param a double value
602 * @param x double value
603 * @return the Complemented Incomplete Gamma function.
604 * <P>
605 * <FONT size=2>
606 * Converted to Java from<BR>
607 * Cephes Math Library Release 2.2: July, 1992<BR>
608 * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>
609 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
610 **/
611
612 static public double igamc( double a, double x )
613 throws ArithmeticException {
614 double big = 4.503599627370496e15;
615 double biginv = 2.22044604925031308085e-16;
616 double ans, ax, c, yc, r, t, y, z;
617 double pk, pkm1, pkm2, qk, qkm1, qkm2;
618
619 if( x <= 0 || a <= 0 ) return 1.0;
620
621 if( x < 1.0 || x < a ) return 1.0 - igam(a,x);
622
623 ax = a * Math.log(x) - x - lgamma(a);
624 if( ax < -MAXLOG ) return 0.0;
625
626 ax = Math.exp(ax);
627
628 /* continued fraction */
629 y = 1.0 - a;
630 z = x + y + 1.0;
631 c = 0.0;
632 pkm2 = 1.0;
633 qkm2 = x;
634 pkm1 = x + 1.0;
635 qkm1 = z * x;
636 ans = pkm1/qkm1;
637
638 do {
639 c += 1.0;
640 y += 1.0;
641 z += 2.0;
642 yc = y * c;
643 pk = pkm1 * z - pkm2 * yc;
644 qk = qkm1 * z - qkm2 * yc;
645 if( qk != 0 ) {
646 r = pk/qk;
647 t = Math.abs( (ans - r)/r );
648 ans = r;
649 } else
650 t = 1.0;
651
652 pkm2 = pkm1;
653 pkm1 = pk;
654 qkm2 = qkm1;
655 qkm1 = qk;
656 if( Math.abs(pk) > big ) {
657 pkm2 *= biginv;
658 pkm1 *= biginv;
659 qkm2 *= biginv;
660 qkm1 *= biginv;
661 }
662 } while( t > MACHEP );
663
664 return ans * ax;
665 }
666
667
668 /**
669 * @param a double value
670 * @param x double value
671 * @return the Incomplete Gamma function.
672 * <P>
673 * <FONT size=2>
674 * Converted to Java from<BR>
675 * Cephes Math Library Release 2.2: July, 1992<BR>
676 * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>
677 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
678 **/
679 static public double igam(double a, double x)
680 throws ArithmeticException {
681
682
683 double ans, ax, c, r;
684
685 if( x <= 0 || a <= 0 ) return 0.0;
686
687 if( x > 1.0 && x > a ) return 1.0 - igamc(a,x);
688
689 /* Compute x**a * exp(-x) / gamma(a) */
690 ax = a * Math.log(x) - x - lgamma(a);
691 if( ax < -MAXLOG ) return( 0.0 );
692
693 ax = Math.exp(ax);
694
695 /* power series */
696 r = a;
697 c = 1.0;
698 ans = 1.0;
699
700 do {
701 r += 1.0;
702 c *= x/r;
703 ans += c;
704 }
705 while( c/ans > MACHEP );
706
707 return( ans * ax/a );
708
709 }
710
711 /**
712 * Returns the area under the left hand tail (from 0 to x)
713 * of the Chi square probability density function with
714 * v degrees of freedom.
715 *
716 * @param df degrees of freedom
717 * @param x double value
718 * @return the Chi-Square function.
719 **/
720
721 static public double chisq(double df, double x)
722 throws ArithmeticException {
723
724 if( x < 0.0 || df < 1.0 ) return 0.0;
725
726 return igam( df/2.0, x/2.0 );
727
728 }
729
730 /**
731 * Returns the area under the right hand tail (from x to
732 * infinity) of the Chi square probability density function
733 * with v degrees of freedom:
734 *
735 * @param df degrees of freedom
736 * @param x double value
737 * @return the Chi-Square function.
738 * <P>
739 **/
740
741 static public double chisqc(double df, double x)
742 throws ArithmeticException {
743
744 if( x < 0.0 || df < 1.0 ) return 0.0;
745
746 return igamc( df/2.0, x/2.0 );
747
748 }
749
750 /**
751 * Returns the sum of the first k terms of the Poisson
752 * distribution.
753 * @param k number of terms
754 * @param x double value
755 */
756
757 static public double poisson(int k, double x)
758 throws ArithmeticException {
759
760
761 if( k < 0 || x < 0 ) return 0.0;
762
763 return igamc((double)(k+1) ,x);
764 }
765
766 /**
767 * Returns the sum of the terms k+1 to infinity of the Poisson
768 * distribution.
769 * @param k start
770 * @param x double value
771 */
772
773 static public double poissonc(int k, double x)
774 throws ArithmeticException {
775
776
777 if( k < 0 || x < 0 ) return 0.0;
778
779 return igam((double)(k+1),x);
780 }
781
782
783
784 /**
785 * @param a double value
786 * @return The area under the Gaussian probability density
787 * function, integrated from minus infinity to x:
788 */
789
790 static public double normal( double a)
791 throws ArithmeticException {
792 double x, y, z;
793
794 x = a * SQRTH;
795 z = Math.abs(x);
796
797 if( z < SQRTH ) y = 0.5 + 0.5 * erf(x);
798 else {
799 y = 0.5 * erfc(z);
800 if( x > 0 ) y = 1.0 - y;
801 }
802
803 return y;
804 }
805
806
807 /**
808 * @param a double value
809 * @return The complementary Error function
810 * <P>
811 * <FONT size=2>
812 * Converted to Java from<BR>
813 * Cephes Math Library Release 2.2: July, 1992<BR>
814 * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>
815 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
816 */
817
818 static public double erfc(double a)
819 throws ArithmeticException {
820 double x,y,z,p,q;
821
822 double P[] = {
823 2.46196981473530512524E-10,
824 5.64189564831068821977E-1,
825 7.46321056442269912687E0,
826 4.86371970985681366614E1,
827 1.96520832956077098242E2,
828 5.26445194995477358631E2,
829 9.34528527171957607540E2,
830 1.02755188689515710272E3,
831 5.57535335369399327526E2
832 };
833 double Q[] = {
834 //1.0
835 1.32281951154744992508E1,
836 8.67072140885989742329E1,
837 3.54937778887819891062E2,
838 9.75708501743205489753E2,
839 1.82390916687909736289E3,
840 2.24633760818710981792E3,
841 1.65666309194161350182E3,
842 5.57535340817727675546E2
843 };
844
845 double R[] = {
846 5.64189583547755073984E-1,
847 1.27536670759978104416E0,
848 5.01905042251180477414E0,
849 6.16021097993053585195E0,
850 7.40974269950448939160E0,
851 2.97886665372100240670E0
852 };
853 double S[] = {
854 //1.00000000000000000000E0,
855 2.26052863220117276590E0,
856 9.39603524938001434673E0,
857 1.20489539808096656605E1,
858 1.70814450747565897222E1,
859 9.60896809063285878198E0,
860 3.36907645100081516050E0
861 };
862
863 if( a < 0.0 ) x = -a;
864 else x = a;
865
866 if( x < 1.0 ) return 1.0 - erf(a);
867
868 z = -a * a;
869
870 if( z < -MAXLOG ) {
871 if( a < 0 ) return( 2.0 );
872 else return( 0.0 );
873 }
874
875 z = Math.exp(z);
876
877 if( x < 8.0 ) {
878 p = polevl( x, P, 8 );
879 q = p1evl( x, Q, 8 );
880 } else {
881 p = polevl( x, R, 5 );
882 q = p1evl( x, S, 6 );
883 }
884
885 y = (z * p)/q;
886
887 if( a < 0 ) y = 2.0 - y;
888
889 if( y == 0.0 ) {
890 if( a < 0 ) return 2.0;
891 else return( 0.0 );
892 }
893
894
895 return y;
896 }
897
898 /**
899 * @param a double value
900 * @return The Error function
901 * <P>
902 * <FONT size=2>
903 * Converted to Java from<BR>
904 * Cephes Math Library Release 2.2: July, 1992<BR>
905 * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>
906 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
907 */
908
909 static public double erf(double x)
910 throws ArithmeticException {
911 double y, z;
912 double T[] = {
913 9.60497373987051638749E0,
914 9.00260197203842689217E1,
915 2.23200534594684319226E3,
916 7.00332514112805075473E3,
917 5.55923013010394962768E4
918 };
919 double U[] = {
920 //1.00000000000000000000E0,
921 3.35617141647503099647E1,
922 5.21357949780152679795E2,
923 4.59432382970980127987E3,
924 2.26290000613890934246E4,
925 4.92673942608635921086E4
926 };
927
928 if( Math.abs(x) > 1.0 ) return( 1.0 - erfc(x) );
929 z = x * x;
930 y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
931 return y;
932 }
933
934
935
936
937
938 static private double polevl( double x, double coef[], int N )
939 throws ArithmeticException {
940
941 double ans;
942
943 ans = coef[0];
944
945 for(int i=1; i<=N; i++) { ans = ans*x+coef[i]; }
946
947 return ans;
948 }
949
950 static private double p1evl( double x, double coef[], int N )
951 throws ArithmeticException {
952
953 double ans;
954
955 ans = x + coef[0];
956
957 for(int i=1; i<N; i++) { ans = ans*x+coef[i]; }
958
959 return ans;
960 }
961 /*
962 *
963 * Natural logarithm of gamma function
964 *
965 */
966 /*
967 Cephes Math Library Release 2.2: July, 1992
968 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
969 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
970 */
971
972
973 static private double lgamma(double x)
974 throws ArithmeticException {
975 double p, q, w, z;
976
977 double A[] = {
978 8.11614167470508450300E-4,
979 -5.95061904284301438324E-4,
980 7.93650340457716943945E-4,
981 -2.77777777730099687205E-3,
982 8.33333333333331927722E-2
983 };
984 double B[] = {
985 -1.37825152569120859100E3,
986 -3.88016315134637840924E4,
987 -3.31612992738871184744E5,
988 -1.16237097492762307383E6,
989 -1.72173700820839662146E6,
990 -8.53555664245765465627E5
991 };
992 double C[] = {
993 /* 1.00000000000000000000E0, */
994 -3.51815701436523470549E2,
995 -1.70642106651881159223E4,
996 -2.20528590553854454839E5,
997 -1.13933444367982507207E6,
998 -2.53252307177582951285E6,
999 -2.01889141433532773231E6
1000 };
1001
1002 if( x < -34.0 ) {
1003 q = -x;
1004 w = lgamma(q);
1005 p = Math.floor(q);
1006 if( p == q ) throw new ArithmeticException("lgam: Overflow");
1007 z = q - p;
1008 if( z > 0.5 ) {
1009 p += 1.0;
1010 z = p - q;
1011 }
1012 z = q * Math.sin( Math.PI * z );
1013 if( z == 0.0 ) throw new
1014 ArithmeticException("lgamma: Overflow");
1015 z = LOGPI - Math.log( z ) - w;
1016 return z;
1017 }
1018
1019 if( x < 13.0 ) {
1020 z = 1.0;
1021 while( x >= 3.0 ) {
1022 x -= 1.0;
1023 z *= x;
1024 }
1025 while( x < 2.0 ) {
1026 if( x == 0.0 ) throw new
1027 ArithmeticException("lgamma: Overflow");
1028 z /= x;
1029 x += 1.0;
1030 }
1031 if( z < 0.0 ) z = -z;
1032 if( x == 2.0 ) return Math.log(z);
1033 x -= 2.0;
1034 p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
1035 return( Math.log(z) + p );
1036 }
1037
1038 if( x > 2.556348e305 ) throw new
1039 ArithmeticException("lgamma: Overflow");
1040
1041 q = ( x - 0.5 ) * Math.log(x) - x + 0.91893853320467274178;
1042 if( x > 1.0e8 ) return( q );
1043
1044 p = 1.0/(x*x);
1045 if( x >= 1000.0 )
1046 q += (( 7.9365079365079365079365e-4 * p
1047 - 2.7777777777777777777778e-3) *p
1048 + 0.0833333333333333333333) / x;
1049 else
1050 q += polevl( p, A, 4 ) / x;
1051 return q;
1052 }
1053
1054
1055
1056 /**
1057 * @param aa double value
1058 * @param bb double value
1059 * @param xx double value
1060 * @return The Incomplete Beta Function evaluated from zero to xx.
1061 * <P>
1062 * <FONT size=2>
1063 * Converted to Java from<BR>
1064 * Cephes Math Library Release 2.3: July, 1995<BR>
1065 * Copyright 1984, 1995 by Stephen L. Moshier<BR>
1066 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
1067 */
1068
1069 public static double ibeta( double aa, double bb, double xx )
1070 throws ArithmeticException {
1071 double a, b, t, x, xc, w, y;
1072 boolean flag;
1073
1074 if( aa <= 0.0 || bb <= 0.0 ) throw new
1075 ArithmeticException("ibeta: Domain error!");
1076
1077 if( (xx <= 0.0) || ( xx >= 1.0) ) {
1078 if( xx == 0.0 ) return 0.0;
1079 if( xx == 1.0 ) return 1.0;
1080 throw new ArithmeticException("ibeta: Domain error!");
1081 }
1082
1083 flag = false;
1084 if( (bb * xx) <= 1.0 && xx <= 0.95) {
1085 t = pseries(aa, bb, xx);
1086 return t;
1087 }
1088
1089 w = 1.0 - xx;
1090
1091 /* Reverse a and b if x is greater than the mean. */
1092 if( xx > (aa/(aa+bb)) ) {
1093 flag = true;
1094 a = bb;
1095 b = aa;
1096 xc = xx;
1097 x = w;
1098 } else {
1099 a = aa;
1100 b = bb;
1101 xc = w;
1102 x = xx;
1103 }
1104
1105 if( flag && (b * x) <= 1.0 && x <= 0.95) {
1106 t = pseries(a, b, x);
1107 if( t <= MACHEP ) t = 1.0 - MACHEP;
1108 else t = 1.0 - t;
1109 return t;
1110 }
1111
1112 /* Choose expansion for better convergence. */
1113 y = x * (a+b-2.0) - (a-1.0);
1114 if( y < 0.0 )
1115 w = incbcf( a, b, x );
1116 else
1117 w = incbd( a, b, x ) / xc;
1118
1119 /* Multiply w by the factor
1120 a b _ _ _
1121 x (1-x) | (a+b) / ( a | (a) | (b) ) . */
1122
1123 y = a * Math.log(x);
1124 t = b * Math.log(xc);
1125 if( (a+b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG ) {
1126 t = Math.pow(xc,b);
1127 t *= Math.pow(x,a);
1128 t /= a;
1129 t *= w;
1130 t *= gamma(a+b) / (gamma(a) * gamma(b));
1131 if( flag ) {
1132 if( t <= MACHEP ) t = 1.0 - MACHEP;
1133 else t = 1.0 - t;
1134 }
1135 return t;
1136 }
1137 /* Resort to logarithms. */
1138 y += t + lgamma(a+b) - lgamma(a) - lgamma(b);
1139 y += Math.log(w/a);
1140 if( y < MINLOG )
1141 t = 0.0;
1142 else
1143 t = Math.exp(y);
1144
1145 if( flag ) {
1146 if( t <= MACHEP ) t = 1.0 - MACHEP;
1147 else t = 1.0 - t;
1148 }
1149 return t;
1150 }
1151
1152/* Continued fraction expansion #1
1153 * for incomplete beta integral
1154 */
1155
1156 private static double incbcf( double a, double b, double x )
1157 throws ArithmeticException {
1158 double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
1159 double k1, k2, k3, k4, k5, k6, k7, k8;
1160 double r, t, ans, thresh;
1161 int n;
1162 double big = 4.503599627370496e15;
1163 double biginv = 2.22044604925031308085e-16;
1164
1165 k1 = a;
1166 k2 = a + b;
1167 k3 = a;
1168 k4 = a + 1.0;
1169 k5 = 1.0;
1170 k6 = b - 1.0;
1171 k7 = k4;
1172 k8 = a + 2.0;
1173
1174 pkm2 = 0.0;
1175 qkm2 = 1.0;
1176 pkm1 = 1.0;
1177 qkm1 = 1.0;
1178 ans = 1.0;
1179 r = 1.0;
1180 n = 0;
1181 thresh = 3.0 * MACHEP;
1182 do {
1183 xk = -( x * k1 * k2 )/( k3 * k4 );
1184 pk = pkm1 + pkm2 * xk;
1185 qk = qkm1 + qkm2 * xk;
1186 pkm2 = pkm1;
1187 pkm1 = pk;
1188 qkm2 = qkm1;
1189 qkm1 = qk;
1190
1191 xk = ( x * k5 * k6 )/( k7 * k8 );
1192 pk = pkm1 + pkm2 * xk;
1193 qk = qkm1 + qkm2 * xk;
1194 pkm2 = pkm1;
1195 pkm1 = pk;
1196 qkm2 = qkm1;
1197 qkm1 = qk;
1198
1199 if( qk != 0 ) r = pk/qk;
1200 if( r != 0 ) {
1201 t = Math.abs( (ans - r)/r );
1202 ans = r;
1203 } else
1204 t = 1.0;
1205
1206 if( t < thresh ) return ans;
1207
1208 k1 += 1.0;
1209 k2 += 1.0;
1210 k3 += 2.0;
1211 k4 += 2.0;
1212 k5 += 1.0;
1213 k6 -= 1.0;
1214 k7 += 2.0;
1215 k8 += 2.0;
1216
1217 if( (Math.abs(qk) + Math.abs(pk)) > big ) {
1218 pkm2 *= biginv;
1219 pkm1 *= biginv;
1220 qkm2 *= biginv;
1221 qkm1 *= biginv;
1222 }
1223 if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
1224 pkm2 *= big;
1225 pkm1 *= big;
1226 qkm2 *= big;
1227 qkm1 *= big;
1228 }
1229 } while( ++n < 300 );
1230
1231 return ans;
1232 }
1233/* Continued fraction expansion #2
1234 * for incomplete beta integral
1235 */
1236
1237 static private double incbd( double a, double b, double x )
1238 throws ArithmeticException {
1239 double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
1240 double k1, k2, k3, k4, k5, k6, k7, k8;
1241 double r, t, ans, z, thresh;
1242 int n;
1243 double big = 4.503599627370496e15;
1244 double biginv = 2.22044604925031308085e-16;
1245
1246 k1 = a;
1247 k2 = b - 1.0;
1248 k3 = a;
1249 k4 = a + 1.0;
1250 k5 = 1.0;
1251 k6 = a + b;
1252 k7 = a + 1.0;;
1253 k8 = a + 2.0;
1254
1255 pkm2 = 0.0;
1256 qkm2 = 1.0;
1257 pkm1 = 1.0;
1258 qkm1 = 1.0;
1259 z = x / (1.0-x);
1260 ans = 1.0;
1261 r = 1.0;
1262 n = 0;
1263 thresh = 3.0 * MACHEP;
1264 do {
1265 xk = -( z * k1 * k2 )/( k3 * k4 );
1266 pk = pkm1 + pkm2 * xk;
1267 qk = qkm1 + qkm2 * xk;
1268 pkm2 = pkm1;
1269 pkm1 = pk;
1270 qkm2 = qkm1;
1271 qkm1 = qk;
1272
1273 xk = ( z * k5 * k6 )/( k7 * k8 );
1274 pk = pkm1 + pkm2 * xk;
1275 qk = qkm1 + qkm2 * xk;
1276 pkm2 = pkm1;
1277 pkm1 = pk;
1278 qkm2 = qkm1;
1279 qkm1 = qk;
1280
1281 if( qk != 0 ) r = pk/qk;
1282 if( r != 0 ) {
1283 t = Math.abs( (ans - r)/r );
1284 ans = r;
1285 } else
1286 t = 1.0;
1287
1288 if( t < thresh ) return ans;
1289
1290 k1 += 1.0;
1291 k2 -= 1.0;
1292 k3 += 2.0;
1293 k4 += 2.0;
1294 k5 += 1.0;
1295 k6 += 1.0;
1296 k7 += 2.0;
1297 k8 += 2.0;
1298
1299 if( (Math.abs(qk) + Math.abs(pk)) > big ) {
1300 pkm2 *= biginv;
1301 pkm1 *= biginv;
1302 qkm2 *= biginv;
1303 qkm1 *= biginv;
1304 }
1305 if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
1306 pkm2 *= big;
1307 pkm1 *= big;
1308 qkm2 *= big;
1309 qkm1 *= big;
1310 }
1311 } while( ++n < 300 );
1312
1313 return ans;
1314 }
1315/* Power series for incomplete beta integral.
1316 Use when b*x is small and x not too close to 1. */
1317
1318 static private double pseries( double a, double b, double x )
1319 throws ArithmeticException {
1320 double s, t, u, v, n, t1, z, ai;
1321
1322 ai = 1.0 / a;
1323 u = (1.0 - b) * x;
1324 v = u / (a + 1.0);
1325 t1 = v;
1326 t = u;
1327 n = 2.0;
1328 s = 0.0;
1329 z = MACHEP * ai;
1330 while( Math.abs(v) > z ) {
1331 u = (n - b) * x / n;
1332 t *= u;
1333 v = t / (a + n);
1334 s += v;
1335 n += 1.0;
1336 }
1337 s += t1;
1338 s += ai;
1339
1340 u = a * Math.log(x);
1341 if( (a+b) < MAXGAM && Math.abs(u) < MAXLOG ) {
1342 t = gamma(a+b)/(gamma(a)*gamma(b));
1343 s = s * t * Math.pow(x,a);
1344 } else {
1345 t = lgamma(a+b) - lgamma(a) - lgamma(b) + u + Math.log(s);
1346 if( t < MINLOG ) s = 0.0;
1347 else s = Math.exp(t);
1348 }
1349 return s;
1350 }
1351
1352}