package riso.numerical;
/** The special function math library.
* This is copied from JSci version 0.876, and redistributed under GPL by permission of Mark Hale.
* See <a href="http://fourier.dur.ac.uk:8000/~dma3mjh/jsci/">the JSci project home page</a>.
* See also <a href="http://fourier.dur.ac.uk:8000/~dma3mjh/">Mark Hale's home page</a>.
*
* <p> In response to my question,
* <pre>
* > Do I have your permission to redistribute the
* > SpecialMath.java file under the GPL?
* </pre>
* Mark Hale wrote:
* <pre>
* Date: Fri, 04 Jan 1980 12:44:25 +0000
* From: "Mark" <kk8mark@netscape.net>
* Reply-to: mark.hale@physics.org
* To: "Robert Dodier" <robert_dodier@yahoo.com>
* Subject: Re: Permission to redistribute SpecialMath code?
*
* Yes, sure. I will probably unbundle xml4j at some stage so I can
* explicitly include a GPL license with JSci.
*
* Mark
* </pre>
*
* @author Mark Hale
* @author Jaco van Kooten
*/
public final class SpecialMath {
private SpecialMath() {}
/** Square root of 2 pi. Copied from the <tt>NumericalConstants</tt> interface.
*/
final static double SQRT2PI=2.5066282746310005024157652848110452530069867406099;
// Some IEEE machine constants
/**
* Relative machine precision.
*/
private final static double EPS=2.22e-16;
/**
* The smallest positive floating-point number such that 1/xminin is machine representable.
*/
private final static double XMININ=2.23e-308;
// CHEBYSHEV SERIES
// series for ai0 on the interval 1.25000d-01 to 3.33333d-01
// with weighted error 7.87e-17
// log weighted error 16.10
// significant figures required 14.69
// decimal places required 16.76
private final static double ai0cs[]={
0.07575994494023796,
0.00759138081082334,
0.00041531313389237,
0.00001070076463439,
-0.00000790117997921,
-0.00000078261435014,
0.00000027838499429,
0.00000000825247260,
-0.00000001204463945,
0.00000000155964859,
0.00000000022925563,
-0.00000000011916228,
0.00000000001757854,
0.00000000000112822,
-0.00000000000114684,
0.00000000000027155,
-0.00000000000002415,
-0.00000000000000608,
0.00000000000000314,
-0.00000000000000071,
0.00000000000000007};
// series for ai02 on the interval 0. to 1.25000d-01
// with weighted error 3.79e-17
// log weighted error 16.42
// significant figures required 14.86
// decimal places required 17.09
private final static double ai02cs[]={
0.05449041101410882,
0.00336911647825569,
0.00006889758346918,
0.00000289137052082,
0.00000020489185893,
0.00000002266668991,
0.00000000339623203,
0.00000000049406022,
0.00000000001188914,
-0.00000000003149915,
-0.00000000001321580,
-0.00000000000179419,
0.00000000000071801,
0.00000000000038529,
0.00000000000001539,
-0.00000000000004151,
-0.00000000000000954,
0.00000000000000382,
0.00000000000000176,
-0.00000000000000034,
-0.00000000000000027,
0.00000000000000003};
// series for ai1 on the interval 1.25000d-01 to 3.33333d-01
// with weighted error 6.98e-17
// log weighted error 16.16
// significant figures required 14.53
// decimal places required 16.82
private final static double ai1cs[]={
-0.02846744181881479,
-0.01922953231443221,
-0.00061151858579437,
-0.00002069971253350,
0.00000858561914581,
0.00000104949824671,
-0.00000029183389184,
-0.00000001559378146,
0.00000001318012367,
-0.00000000144842341,
-0.00000000029085122,
0.00000000012663889,
-0.00000000001664947,
-0.00000000000166665,
0.00000000000124260,
-0.00000000000027315,
0.00000000000002023,
0.00000000000000730,
-0.00000000000000333,
0.00000000000000071,
-0.00000000000000006};
// series for ai12 on the interval 0. to 1.25000d-01
// with weighted error 3.55e-17
// log weighted error 16.45
// significant figures required 14.69
// decimal places required 17.12
private final static double ai12cs[]={
0.02857623501828014,
-0.00976109749136147,
-0.00011058893876263,
-0.00000388256480887,
-0.00000025122362377,
-0.00000002631468847,
-0.00000000383538039,
-0.00000000055897433,
-0.00000000001897495,
0.00000000003252602,
0.00000000001412580,
0.00000000000203564,
-0.00000000000071985,
-0.00000000000040836,
-0.00000000000002101,
0.00000000000004273,
0.00000000000001041,
-0.00000000000000382,
-0.00000000000000186,
0.00000000000000033,
0.00000000000000028,
-0.00000000000000003};
// series for aif on the interval -1.00000d+00 to 1.00000d+00
// with weighted error 1.09e-19
// log weighted error 18.96
// significant figures required 17.76
// decimal places required 19.44
private final static double aifcs[]={
-0.03797135849666999750,
0.05919188853726363857,
0.00098629280577279975,
0.00000684884381907656,
0.00000002594202596219,
0.00000000006176612774,
0.00000000000010092454,
0.00000000000000012014,
0.00000000000000000010};
// series for aig on the interval -1.00000d+00 to 1.00000d+00
// with weighted error 1.51e-17
// log weighted error 16.82
// significant figures required 15.19
// decimal places required 17.27
private final static double aigcs[]={
0.01815236558116127,
0.02157256316601076,
0.00025678356987483,
0.00000142652141197,
0.00000000457211492,
0.00000000000952517,
0.00000000000001392,
0.00000000000000001};
// series for aip on the interval 0. to 1.00000d+00
// with weighted error 5.10e-17
// log weighted error 16.29
// significant figures required 14.41
// decimal places required 17.06
private final static double aipcs[]={
-0.0187519297793868,
-0.0091443848250055,
0.0009010457337825,
-0.0001394184127221,
0.0000273815815785,
-0.0000062750421119,
0.0000016064844184,
-0.0000004476392158,
0.0000001334635874,
-0.0000000420735334,
0.0000000139021990,
-0.0000000047831848,
0.0000000017047897,
-0.0000000006268389,
0.0000000002369824,
-0.0000000000918641,
0.0000000000364278,
-0.0000000000147475,
0.0000000000060851,
-0.0000000000025552,
0.0000000000010906,
-0.0000000000004725,
0.0000000000002076,
-0.0000000000000924,
0.0000000000000417,
-0.0000000000000190,
0.0000000000000087,
-0.0000000000000040,
0.0000000000000019,
-0.0000000000000009,
0.0000000000000004,
-0.0000000000000002,
0.0000000000000001,
-0.0000000000000000};
// series for am21 on the interval -1.25000d-01 to 0.
// with weighted error 2.89e-17
// log weighted error 16.54
// significant figures required 14.15
// decimal places required 17.34
private final static double am21cs[]={
0.0065809191761485,
0.0023675984685722,
0.0001324741670371,
0.0000157600904043,
0.0000027529702663,
0.0000006102679017,
0.0000001595088468,
0.0000000471033947,
0.0000000152933871,
0.0000000053590722,
0.0000000020000910,
0.0000000007872292,
0.0000000003243103,
0.0000000001390106,
0.0000000000617011,
0.0000000000282491,
0.0000000000132979,
0.0000000000064188,
0.0000000000031697,
0.0000000000015981,
0.0000000000008213,
0.0000000000004296,
0.0000000000002284,
0.0000000000001232,
0.0000000000000675,
0.0000000000000374,
0.0000000000000210,
0.0000000000000119,
0.0000000000000068,
0.0000000000000039,
0.0000000000000023,
0.0000000000000013,
0.0000000000000008,
0.0000000000000005,
0.0000000000000003,
0.0000000000000001,
0.0000000000000001,
0.0000000000000000,
0.0000000000000000,
0.0000000000000000};
// series for ath1 on the interval -1.25000d-01 to 0.
// with weighted error 2.53e-17
// log weighted error 16.60
// significant figures required 15.15
// decimal places required 17.38
private final static double ath1cs[]={
-0.07125837815669365,
-0.00590471979831451,
-0.00012114544069499,
-0.00000988608542270,
-0.00000138084097352,
-0.00000026142640172,
-0.00000006050432589,
-0.00000001618436223,
-0.00000000483464911,
-0.00000000157655272,
-0.00000000055231518,
-0.00000000020545441,
-0.00000000008043412,
-0.00000000003291252,
-0.00000000001399875,
-0.00000000000616151,
-0.00000000000279614,
-0.00000000000130428,
-0.00000000000062373,
-0.00000000000030512,
-0.00000000000015239,
-0.00000000000007758,
-0.00000000000004020,
-0.00000000000002117,
-0.00000000000001132,
-0.00000000000000614,
-0.00000000000000337,
-0.00000000000000188,
-0.00000000000000105,
-0.00000000000000060,
-0.00000000000000034,
-0.00000000000000020,
-0.00000000000000011,
-0.00000000000000007,
-0.00000000000000004,
-0.00000000000000002};
// series for am22 on the interval -1.00000d+00 to -1.25000d-01
// with weighted error 2.99e-17
// log weighted error 16.52
// significant figures required 14.57
// decimal places required 17.28
private final static double am22cs[]={
-0.01562844480625341,
0.00778336445239681,
0.00086705777047718,
0.00015696627315611,
0.00003563962571432,
0.00000924598335425,
0.00000262110161850,
0.00000079188221651,
0.00000025104152792,
0.00000008265223206,
0.00000002805711662,
0.00000000976821090,
0.00000000347407923,
0.00000000125828132,
0.00000000046298826,
0.00000000017272825,
0.00000000006523192,
0.00000000002490471,
0.00000000000960156,
0.00000000000373448,
0.00000000000146417,
0.00000000000057826,
0.00000000000022991,
0.00000000000009197,
0.00000000000003700,
0.00000000000001496,
0.00000000000000608,
0.00000000000000248,
0.00000000000000101,
0.00000000000000041,
0.00000000000000017,
0.00000000000000007,
0.00000000000000002};
// series for ath2 on the interval -1.00000d+00 to -1.25000d-01
// with weighted error 2.57e-17
// log weighted error 16.59
// significant figures required 15.07
// decimal places required 17.34
private final static double ath2cs[]={
0.00440527345871877,
-0.03042919452318455,
-0.00138565328377179,
-0.00018044439089549,
-0.00003380847108327,
-0.00000767818353522,
-0.00000196783944371,
-0.00000054837271158,
-0.00000016254615505,
-0.00000005053049981,
-0.00000001631580701,
-0.00000000543420411,
-0.00000000185739855,
-0.00000000064895120,
-0.00000000023105948,
-0.00000000008363282,
-0.00000000003071196,
-0.00000000001142367,
-0.00000000000429811,
-0.00000000000163389,
-0.00000000000062693,
-0.00000000000024260,
-0.00000000000009461,
-0.00000000000003716,
-0.00000000000001469,
-0.00000000000000584,
-0.00000000000000233,
-0.00000000000000093,
-0.00000000000000037,
-0.00000000000000015,
-0.00000000000000006,
-0.00000000000000002};
// series for bi0 on the interval 0. to 9.00000d+00
// with weighted error 2.46e-18
// log weighted error 17.61
// significant figures required 17.90
// decimal places required 18.15
private final static double bi0cs[]={
-0.07660547252839144951,
1.927337953993808270,
0.2282644586920301339,
0.01304891466707290428,
0.00043442709008164874,
0.00000942265768600193,
0.00000014340062895106,
0.00000000161384906966,
0.00000000001396650044,
0.00000000000009579451,
0.00000000000000053339,
0.00000000000000000245};
// series for bj0 on the interval 0. to 1.60000d+01
// with weighted error 7.47e-18
// log weighted error 17.13
// significant figures required 16.98
// decimal places required 17.68
private final static double bj0cs[]={
0.100254161968939137,
-0.665223007764405132,
0.248983703498281314,
-0.0332527231700357697,
0.0023114179304694015,
-0.0000991127741995080,
0.0000028916708643998,
-0.0000000612108586630,
0.0000000009838650793,
-0.0000000000124235515,
0.0000000000001265433,
-0.0000000000000010619,
0.0000000000000000074};
// series for bm0 on the interval 0. to 6.25000d-02
// with weighted error 4.98e-17
// log weighted error 16.30
// significant figures required 14.97
// decimal places required 16.96
private final static double bm0cs[]={
0.09284961637381644,
-0.00142987707403484,
0.00002830579271257,
-0.00000143300611424,
0.00000012028628046,
-0.00000001397113013,
0.00000000204076188,
-0.00000000035399669,
0.00000000007024759,
-0.00000000001554107,
0.00000000000376226,
-0.00000000000098282,
0.00000000000027408,
-0.00000000000008091,
0.00000000000002511,
-0.00000000000000814,
0.00000000000000275,
-0.00000000000000096,
0.00000000000000034,
-0.00000000000000012,
0.00000000000000004};
// series for bth0 on the interval 0. to 6.25000d-02
// with weighted error 3.67e-17
// log weighted error 16.44
// significant figures required 15.53
// decimal places required 17.13
private final static double bth0cs[]={
-0.24639163774300119,
0.001737098307508963,
-0.000062183633402968,
0.000004368050165742,
-0.000000456093019869,
0.000000062197400101,
-0.000000010300442889,
0.000000001979526776,
-0.000000000428198396,
0.000000000102035840,
-0.000000000026363898,
0.000000000007297935,
-0.000000000002144188,
0.000000000000663693,
-0.000000000000215126,
0.000000000000072659,
-0.000000000000025465,
0.000000000000009229,
-0.000000000000003448,
0.000000000000001325,
-0.000000000000000522,
0.000000000000000210,
-0.000000000000000087,
0.000000000000000036};
// series for by0 on the interval 0. to 1.60000d+01
// with weighted error 1.20e-17
// log weighted error 16.92
// significant figures required 16.15
// decimal places required 17.48
private final static double by0cs[]={
-0.011277839392865573,
-0.12834523756042035,
-0.10437884799794249,
0.023662749183969695,
-0.002090391647700486,
0.000103975453939057,
-0.000003369747162423,
0.000000077293842676,
-0.000000001324976772,
0.000000000017648232,
-0.000000000000188105,
0.000000000000001641,
-0.000000000000000011};
// series for bi1 on the interval 0. to 9.00000d+00
// with weighted error 2.40e-17
// log weighted error 16.62
// significant figures required 16.23
// decimal places required 17.14
private final static double bi1cs[]={
-0.001971713261099859,
0.40734887667546481,
0.034838994299959456,
0.001545394556300123,
0.000041888521098377,
0.000000764902676483,
0.000000010042493924,
0.000000000099322077,
0.000000000000766380,
0.000000000000004741,
0.000000000000000024};
// series for bj1 on the interval 0. to 1.60000d+01
// with weighted error 4.48e-17
// log weighted error 16.35
// significant figures required 15.77
// decimal places required 16.89
private final static double bj1cs[]={
-0.11726141513332787,
-0.25361521830790640,
0.050127080984469569,
-0.004631514809625081,
0.000247996229415914,
-0.000008678948686278,
0.000000214293917143,
-0.000000003936093079,
0.000000000055911823,
-0.000000000000632761,
0.000000000000005840,
-0.000000000000000044};
// series for bm1 on the interval 0. to 6.25000d-02
// with weighted error 5.61e-17
// log weighted error 16.25
// significant figures required 14.97
// decimal places required 16.91
private final static double bm1cs[]={
0.1047362510931285,
0.00442443893702345,
-0.00005661639504035,
0.00000231349417339,
-0.00000017377182007,
0.00000001893209930,
-0.00000000265416023,
0.00000000044740209,
-0.00000000008691795,
0.00000000001891492,
-0.00000000000451884,
0.00000000000116765,
-0.00000000000032265,
0.00000000000009450,
-0.00000000000002913,
0.00000000000000939,
-0.00000000000000315,
0.00000000000000109,
-0.00000000000000039,
0.00000000000000014,
-0.00000000000000005};
// series for bth1 on the interval 0. to 6.25000d-02
// with weighted error 4.10e-17
// log weighted error 16.39
// significant figures required 15.96
// decimal places required 17.08
private final static double bth1cs[]={
0.74060141026313850,
-0.004571755659637690,
0.000119818510964326,
-0.000006964561891648,
0.000000655495621447,
-0.000000084066228945,
0.000000013376886564,
-0.000000002499565654,
0.000000000529495100,
-0.000000000124135944,
0.000000000031656485,
-0.000000000008668640,
0.000000000002523758,
-0.000000000000775085,
0.000000000000249527,
-0.000000000000083773,
0.000000000000029205,
-0.000000000000010534,
0.000000000000003919,
-0.000000000000001500,
0.000000000000000589,
-0.000000000000000237,
0.000000000000000097,
-0.000000000000000040};
// series for by1 on the interval 0. to 1.60000d+01
// with weighted error 1.87e-18
// log weighted error 17.73
// significant figures required 17.83
// decimal places required 18.30
private final static double by1cs[]={
0.03208047100611908629,
1.262707897433500450,
0.00649996189992317500,
-0.08936164528860504117,
0.01325088122175709545,
-0.00089790591196483523,
0.00003647361487958306,
-0.00000100137438166600,
0.00000001994539657390,
-0.00000000030230656018,
0.00000000000360987815,
-0.00000000000003487488,
0.00000000000000027838,
-0.00000000000000000186};
/**
* Evaluates a Chebyshev series.
* @param x value at which to evaluate series
* @param series the coefficients of the series
*/
public static double chebyshev(double x, double series[]) {
double twox,b0=0.0,b1=0.0,b2=0.0;
twox=2*x;
for(int i=series.length-1;i>-1;i--) {
b2=b1;
b1=b0;
b0=twox*b1-b2+series[i];
}
return 0.5*(b0-b2);
}
/**
* Airy function.
* Based on the NETLIB Fortran function ai written by W. Fullerton.
*/
public static double airy(double x) {
if(x<-1.0) {
double mp[]=airyModPhase(x);
return mp[0]*Math.cos(mp[1]);
} else if(x>1.0)
return expAiry(x)*Math.exp(-2.0*x*Math.sqrt(x)/3.0);
else {
final double z=x*x*x;
return 0.375+(chebyshev(z,aifcs)-x*(0.25+chebyshev(z,aigcs)));
}
}
/**
* Airy modulus and phase.
* Based on the NETLIB Fortran subroutine r9aimp written by W. Fullerton.
* @return an array with [0] containing the modulus and [1] containing the phase.
*/
private static double[] airyModPhase(double x) {
double z,mp[]=new double[2];
if(x<-2.0) {
z=16.0/(x*x*x)+1.0;
mp[0]=0.3125+chebyshev(z,am21cs);
mp[1] =-0.625+chebyshev(z,ath1cs);
} else {
z=(16.0/(x*x*x)+9.0)/7.0;
mp[0]=0.3125+chebyshev(z,am22cs);
mp[1]=-0.625+chebyshev(z,ath2cs);
}
final double sqrtx=Math.sqrt(-x);
mp[0]=Math.sqrt(mp[0]/sqrtx);
mp[1]=Math.PI/4.0-x*sqrtx*mp[1];
return mp;
}
/**
* Exponential scaled Airy function.
* Based on the NETLIB Fortran function aie written by W. Fullerton.
*/
private static double expAiry(double x) {
if(x<-1.0) {
double mp[]=airyModPhase(x);
return mp[0]*Math.cos(mp[1]);
} else if(x<=1.0) {
final double z=x*x*x;
return 0.375+(chebyshev(z,aifcs)-x*(0.25+chebyshev(z,aigcs)))*Math.exp(2.0*x*Math.sqrt(x)/3.0);
} else {
final double sqrtx=Math.sqrt(x);
final double z=2.0/(x*sqrtx)-1.0;
return (0.28125+chebyshev(z,aipcs))/Math.sqrt(sqrtx);
}
}
/**
* Bessel function of first kind, order zero.
* Based on the NETLIB Fortran function besj0 written by W. Fullerton.
*/
public static double besselFirstZero(double x) {
double y=Math.abs(x);
if(y>4.0) {
final double z=32/(y*y)-1;
final double amplitude=(0.75+chebyshev(z,bm0cs))/Math.sqrt(y);
final double theta=y-Math.PI/4.0+chebyshev(z,bth0cs)/y;
return amplitude*Math.cos(theta);
} else if(y==0.0)
return 1.0;
else
return chebyshev(0.125*y*y-1,bj0cs);
}
/**
* Modified Bessel function of first kind, order zero.
* Based on the NETLIB Fortran function besi0 written by W. Fullerton.
*/
public static double modBesselFirstZero(double x) {
double y=Math.abs(x);
if(y>3.0)
return Math.exp(y)*expModBesselFirstZero(x);
else
return 2.75+chebyshev(y*y/4.5-1.0,bi0cs);
}
/**
* Exponential scaled modified Bessel function of first kind, order zero.
* Based on the NETLIB Fortran function besi0e written by W. Fullerton.
*/
private static double expModBesselFirstZero(double x) {
final double y=Math.abs(x);
if(y>3.0) {
if(y>8.0)
return (0.375+chebyshev(16.0/y-1.0,ai02cs))/Math.sqrt(y);
else
return (0.375+chebyshev((48.0/y-11.0)/5.0,ai0cs))/Math.sqrt(y);
} else
return Math.exp(-y)*(2.75+chebyshev(y*y/4.5-1.0,bi0cs));
}
/**
* Bessel function of first kind, order one.
* Based on the NETLIB Fortran function besj1 written by W. Fullerton.
*/
public static double besselFirstOne(double x) {
double y=Math.abs(x);
if(y>4.0) {
final double z=32.0/(y*y)-1.0;
final double amplitude=(0.75+chebyshev(z,bm1cs))/Math.sqrt(y);
final double theta=y-3.0*Math.PI/4.0+chebyshev(z,bth1cs)/y;
return Math.abs(amplitude)*x*Math.cos(theta)/Math.abs(x);
} else if(y==0.0)
return 0.0;
else
return x*(0.25+chebyshev(0.125*y*y-1.0,bj1cs));
}
/**
* Modified Bessel function of first kind, order one.
* Based on the NETLIB Fortran function besi1 written by W. Fullerton.
*/
public static double modBesselFirstOne(double x) {
final double y=Math.abs(x);
if(y>3.0)
return Math.exp(y)*expModBesselFirstOne(x);
else if(y==0.0)
return 0.0;
else
return x*(0.875+chebyshev(y*y/4.5-1.0,bi1cs));
}
/**
* Exponential scaled modified Bessel function of first kind, order one.
* Based on the NETLIB Fortran function besi1e written by W. Fullerton.
*/
private static double expModBesselFirstOne(double x) {
final double y=Math.abs(x);
if(y>3.0) {
if(y>8.0)
return x/y*(0.375+chebyshev(16.0/y-1.0,ai12cs))/Math.sqrt(y);
else
return x/y*(0.375+chebyshev((48.0/y-11.0)/5.0,ai1cs))/Math.sqrt(y);
} else if(y==0.0)
return 0.0;
else
return Math.exp(-y)*x*(0.875+chebyshev(y*y/4.5-1.0,bi1cs));
}
/**
* Bessel function of second kind, order zero.
* Based on the NETLIB Fortran function besy0 written by W. Fullerton.
*/
public static double besselSecondZero(double x) {
if(x>4.0) {
final double z=32.0/(x*x)-1.0;
final double amplitude=(0.75+chebyshev(z,bm0cs))/Math.sqrt(x);
final double theta=x-Math.PI/4+chebyshev(z,bth0cs)/x;
return amplitude*Math.sin(theta);
} else
return (Math.log(0.5)+Math.log(x))*besselFirstZero(x)+0.375+chebyshev(0.125*x*x-1.0,by0cs)*2.0/Math.PI;
}
/**
* Bessel function of second kind, order one.
* Based on the NETLIB Fortran function besy1 written by W. Fullerton.
*/
public static double besselSecondOne(double x) {
if(x>4.0) {
final double z=32.0/(x*x)-1.0;
final double amplitude=(0.75+chebyshev(z,bm1cs))/Math.sqrt(x);
final double theta=x-3.0*Math.PI/4.0+chebyshev(z,bth1cs)/x;
return amplitude*Math.sin(theta);
} else
return 2.0*Math.log(0.5*x)*besselFirstOne(x)/Math.PI+(0.5+chebyshev(0.125*x*x-1.0,by1cs))/x;
}
private final static double LOGSQRT2PI=Math.log(SQRT2PI);
/**
* Gamma function.
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
* Applied Mathematics Division<BR>
* Argonne National Laboratory<BR>
* Argonne, IL 60439<BR>
* <P>
* References:
* <OL>
* <LI>"An Overview of Software Development for Special Functions", W. J. Cody, Lecture Notes in Mathematics, 506, Numerical Analysis Dundee, 1975, G. A. Watson (ed.), Springer Verlag, Berlin, 1976.
* <LI>Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
* </OL></P><P>
* From the original documentation:
* </P><P>
* This routine calculates the GAMMA function for a real argument X.
* Computation is based on an algorithm outlined in reference 1.
* The program uses rational functions that approximate the GAMMA
* function to at least 20 significant decimal digits. Coefficients
* for the approximation over the interval (1,2) are unpublished.
* Those for the approximation for X .GE. 12 are from reference 2.
* The accuracy achieved depends on the arithmetic system, the
* compiler, the intrinsic functions, and proper selection of the
* machine-dependent constants.
* </P><P>
* Error returns:<BR>
* The program returns the value XINF for singularities or when overflow would occur.
* The computation is believed to be free of underflow and overflow.
* </P>
* @return Double.MAX_VALUE if overflow would occur, i.e. if abs(x) > 171.624
* @author Jaco van Kooten
*/
public static double gamma(double x) {
// Gamma function related constants
final double g_p[] = { -1.71618513886549492533811,
24.7656508055759199108314,-379.804256470945635097577,
629.331155312818442661052,866.966202790413211295064,
-31451.2729688483675254357,-36144.4134186911729807069,
66456.1438202405440627855 };
final double g_q[] = { -30.8402300119738975254353,
315.350626979604161529144,-1015.15636749021914166146,
-3107.77167157231109440444,22538.1184209801510330112,
4755.84627752788110767815,-134659.959864969306392456,
-115132.259675553483497211 };
final double g_c[] = { -.001910444077728,8.4171387781295e-4,
-5.952379913043012e-4,7.93650793500350248e-4,
-.002777777777777681622553,.08333333333333333331554247,
.0057083835261 };
final double g_xbig = 171.624;
double fact=1.0, xden, xnum;
int i, n=0;
double y=x, z, y1;
boolean parity=false;
double res, sum, ysq;
if (y <= 0.0) {
// ----------------------------------------------------------------------
// Argument is negative
// ----------------------------------------------------------------------
y = -(x);
y1 = (int)y;
res = y - y1;
if (res != 0.0) {
if (y1 != (((int)(y1*0.5)) * 2.0))
parity = true;
fact = -Math.PI/ Math.sin(Math.PI * res);
y++;
} else
return Double.MAX_VALUE;
}
// ----------------------------------------------------------------------
// Argument is positive
// ----------------------------------------------------------------------
if (y < EPS) {
// ----------------------------------------------------------------------
// Argument .LT. EPS
// ----------------------------------------------------------------------
if (y >= XMININ)
res = 1.0 / y;
else
return Double.MAX_VALUE;
} else if (y < 12.0) {
y1 = y;
if (y < 1.0) {
// ----------------------------------------------------------------------
// 0.0 .LT. argument .LT. 1.0
// ----------------------------------------------------------------------
z = y;
y++;
} else {
// ----------------------------------------------------------------------
// 1.0 .LT. argument .LT. 12.0, reduce argument if necessary
// ----------------------------------------------------------------------
n = (int)y - 1;
y -= (double) n;
z = y - 1.0;
}
// ----------------------------------------------------------------------
// Evaluate approximation for 1.0 .LT. argument .LT. 2.0
// ----------------------------------------------------------------------
xnum = 0.0;
xden = 1.0;
for (i = 0; i < 8; ++i) {
xnum = (xnum + g_p[i]) * z;
xden = xden * z + g_q[i];
}
res = xnum / xden + 1.0;
if (y1 < y)
// ----------------------------------------------------------------------
// Adjust result for case 0.0 .LT. argument .LT. 1.0
// ----------------------------------------------------------------------
res /= y1;
else if (y1 > y) {
// ----------------------------------------------------------------------
// Adjust result for case 2.0 .LT. argument .LT. 12.0
// ----------------------------------------------------------------------
for (i = 0; i < n; ++i) {
res *= y;
y++;
}
}
} else {
// ----------------------------------------------------------------------
// Evaluate for argument .GE. 12.0
// ----------------------------------------------------------------------
if (y <= g_xbig) {
ysq = y * y;
sum = g_c[6];
for (i = 0; i < 6; ++i)
sum = sum / ysq + g_c[i];
sum = sum / y - y + LOGSQRT2PI;
sum += (y - 0.5) * Math.log(y);
res = Math.exp(sum);
} else
return Double.MAX_VALUE;
}
// ----------------------------------------------------------------------
// Final adjustments and return
// ----------------------------------------------------------------------
if (parity)
res = -res;
if (fact != 1.0)
res = fact / res;
return res;
}
/**
* The largest argument for which logGamma(x) is representable in the machine.
*/
private final static double logGamma_xBig=2.55e305;
// Function cache for logGamma
private static double logGammaCache_res=0.0;
private static double logGammaCache_x=0.0;
/**
* The natural logarithm of the gamma function.
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
* Applied Mathematics Division<BR>
* Argonne National Laboratory<BR>
* Argonne, IL 60439<BR>
* <P>
* References:
* <OL>
* <LI>W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.
* <LI>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.
* <LI>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.
* </OL></P><P>
* From the original documentation:
* </P><P>
* This routine calculates the LOG(GAMMA) function for a positive real argument X.
* Computation is based on an algorithm outlined in references 1 and 2.
* The program uses rational functions that theoretically approximate LOG(GAMMA)
* to at least 18 significant decimal digits. The approximation for X > 12 is from reference 3,
* while approximations for X < 12.0 are similar to those in reference 1, but are unpublished.
* The accuracy achieved depends on the arithmetic system, the compiler, the intrinsic functions,
* and proper selection of the machine-dependent constants.
* </P><P>
* Error returns:<BR>
* The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
* The computation is believed to be free of underflow and overflow.
* </P>
* @return Double.MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
* @author Jaco van Kooten
*/
public static double logGamma(double x) {
// Log Gamma related constants
final double lg_d1 = -.5772156649015328605195174,
lg_d2 = .4227843350984671393993777,
lg_d4 = 1.791759469228055000094023;
final double lg_p1[] = { 4.945235359296727046734888,
201.8112620856775083915565,2290.838373831346393026739,
11319.67205903380828685045,28557.24635671635335736389,
38484.96228443793359990269,26377.48787624195437963534,
7225.813979700288197698961 };
final double lg_q1[] = { 67.48212550303777196073036,
1113.332393857199323513008,7738.757056935398733233834,
27639.87074403340708898585,54993.10206226157329794414,
61611.22180066002127833352,36351.27591501940507276287,
8785.536302431013170870835 };
final double lg_p2[] = { 4.974607845568932035012064,
542.4138599891070494101986,15506.93864978364947665077,
184793.2904445632425417223,1088204.76946882876749847,
3338152.967987029735917223,5106661.678927352456275255,
3074109.054850539556250927 };
final double lg_q2[] = { 183.0328399370592604055942,
7765.049321445005871323047,133190.3827966074194402448,
1136705.821321969608938755,5267964.117437946917577538,
13467014.54311101692290052,17827365.30353274213975932,
9533095.591844353613395747 };
final double lg_p4[] = { 14745.02166059939948905062,
2426813.369486704502836312,121475557.4045093227939592,
2663432449.630976949898078,29403789566.34553899906876,
170266573776.5398868392998,492612579337.743088758812,
560625185622.3951465078242 };
final double lg_q4[] = { 2690.530175870899333379843,
639388.5654300092398984238,41355999.30241388052042842,
1120872109.61614794137657,14886137286.78813811542398,
101680358627.2438228077304,341747634550.7377132798597,
446315818741.9713286462081 };
final double lg_c[] = { -0.001910444077728,8.4171387781295e-4,
-5.952379913043012e-4,7.93650793500350248e-4,
-0.002777777777777681622553,0.08333333333333333331554247,
0.0057083835261 };
// Rough estimate of the fourth root of logGamma_xBig
final double lg_frtbig = 2.25e76;
final double pnt68 = 0.6796875;
double xden, corr, xnum;
int i;
double y, xm1, xm2, xm4, res, ysq;
if (x==logGammaCache_x)
return logGammaCache_res;
y = x;
if (y > 0.0 && y <= logGamma_xBig) {
if (y <= EPS) {
res = -Math.log(y);
} else if (y <= 1.5) {
// ----------------------------------------------------------------------
// EPS .LT. X .LE. 1.5
// ----------------------------------------------------------------------
if (y < pnt68) {
corr = -Math.log(y);
xm1 = y;
} else {
corr = 0.0;
xm1 = y - 1.0;
}
if (y <= 0.5 || y >= pnt68) {
xden = 1.0;
xnum = 0.0;
for (i = 0; i < 8; i++) {
xnum = xnum * xm1 + lg_p1[i];
xden = xden * xm1 + lg_q1[i];
}
res = corr + xm1 * (lg_d1 + xm1 * (xnum / xden));
} else {
xm2 = y - 1.0;
xden = 1.0;
xnum = 0.0;
for (i = 0; i < 8; i++) {
xnum = xnum * xm2 + lg_p2[i];
xden = xden * xm2 + lg_q2[i];
}
res = corr + xm2 * (lg_d2 + xm2 * (xnum / xden));
}
} else if (y <= 4.0) {
// ----------------------------------------------------------------------
// 1.5 .LT. X .LE. 4.0
// ----------------------------------------------------------------------
xm2 = y - 2.0;
xden = 1.0;
xnum = 0.0;
for (i = 0; i < 8; i++) {
xnum = xnum * xm2 + lg_p2[i];
xden = xden * xm2 + lg_q2[i];
}
res = xm2 * (lg_d2 + xm2 * (xnum / xden));
} else if (y <= 12.0) {
// ----------------------------------------------------------------------
// 4.0 .LT. X .LE. 12.0
// ----------------------------------------------------------------------
xm4 = y - 4.0;
xden = -1.0;
xnum = 0.0;
for (i = 0; i < 8; i++) {
xnum = xnum * xm4 + lg_p4[i];
xden = xden * xm4 + lg_q4[i];
}
res = lg_d4 + xm4 * (xnum / xden);
} else {
// ----------------------------------------------------------------------
// Evaluate for argument .GE. 12.0
// ----------------------------------------------------------------------
res = 0.0;
if (y <= lg_frtbig) {
res = lg_c[6];
ysq = y * y;
for (i = 0; i < 6; i++)
res = res / ysq + lg_c[i];
}
res /= y;
corr = Math.log(y);
res = res + LOGSQRT2PI - 0.5 * corr;
res += y * (corr - 1.0);
}
} else {
// ----------------------------------------------------------------------
// Return for bad arguments
// ----------------------------------------------------------------------
res = Double.MAX_VALUE;
}
// ----------------------------------------------------------------------
// Final adjustments and return
// ----------------------------------------------------------------------
logGammaCache_x = x;
logGammaCache_res = res;
return res;
}
private final static int MAX_ITERATIONS = 150;
private final static double PRECISION = 4.0*EPS;
/**
* Incomplete gamma function.
* The computation is based on approximations presented in Numerical Recipes, Chapter 6.2 (W.H. Press et al, 1992).
* @param a require a>=0
* @param x require x>=0
* @return 0 if x<0, a<=0 or a>2.55E305 to avoid errors and over/underflow
* @author Jaco van Kooten
*/
public static double incompleteGamma(double a, double x) {
if (x <= 0.0 || a <= 0.0 || a > logGamma_xBig)
return 0.0;
if (x < (a+1.0))
return gammaSeriesExpansion(a,x);
else
return 1.0-gammaFraction(a,x);
}
/**
* @author Jaco van Kooten
*/
private static double gammaSeriesExpansion(double a, double x) {
double ap=a;
double del=1.0/a;
double sum = del;
for (int n=1; n < MAX_ITERATIONS; n++) {
++ap;
del *= x/ap;
sum += del;
if (del < sum*PRECISION)
return sum*Math.exp(-x + a*Math.log(x) - logGamma(a));
}
return 0.0;
}
/**
* @author Jaco van Kooten
*/
private static double gammaFraction(double a, double x) {
double b=x+1.0-a;
double c=1.0/XMININ;
double d=1.0/b;
double h=d;
double del=0.0;
double an;
for (int i=1; i<MAX_ITERATIONS && Math.abs(del-1.0)>PRECISION; i++) {
an = -i*(i-a);
b += 2.0;
d=an*d+b;
c=b+an/c;
if (Math.abs(c) < XMININ)
c=XMININ;
if (Math.abs(d) < XMININ)
c=XMININ;
d=1.0/d;
del=d*c;
h *= del;
}
return Math.exp(-x + a*Math.log(x) - logGamma(a))*h;
}
/**
* Beta function.
* @param p require p>0
* @param q require q>0
* @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
* @author Jaco van Kooten
*/
public static double beta(double p, double q) {
if (p <= 0.0 || q <= 0.0 || (p+q) > logGamma_xBig)
return 0.0;
else
return Math.exp(logBeta(p,q));
}
// Function cache for logBeta
private static double logBetaCache_res=0.0;
private static double logBetaCache_p=0.0;
private static double logBetaCache_q=0.0;
/**
* The natural logarithm of the beta function.
* @param p require p>0
* @param q require q>0
* @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
* @author Jaco van Kooten
*/
public static double logBeta(double p, double q) {
if (p != logBetaCache_p || q != logBetaCache_q) {
logBetaCache_p = p;
logBetaCache_q = q;
if (p <= 0.0 || q <= 0.0 || (p+q) > logGamma_xBig)
logBetaCache_res = 0.0;
else
logBetaCache_res = logGamma(p)+logGamma(q)-logGamma(p+q);
}
return logBetaCache_res;
}
/**
* Incomplete beta function.
* The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
* @param x require 0<=x<=1
* @param p require p>0
* @param q require q>0
* @return 0 if x<0, p<=0, q<=0 or p+q>2.55E305 and 1 if x>1 to avoid errors and over/underflow
* @author Jaco van Kooten
*/
public static double incompleteBeta(double x, double p, double q) {
if (x <= 0.0)
return 0.0;
else if (x >= 1.0)
return 1.0;
else if (p <= 0.0 || q <= 0.0 || (p+q) > logGamma_xBig)
return 0.0;
else {
final double beta_gam=Math.exp(-logBeta(p,q) + p*Math.log(x) + q*Math.log(1.0-x));
if (x < (p+1.0)/(p+q+2.0))
return beta_gam*betaFraction(x,p,q)/p;
else
return 1.0-(beta_gam*betaFraction(1.0-x,q,p)/q);
}
}
/**
* Evaluates of continued fraction part of incomplete beta function.
* Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
* @author Jaco van Kooten
*/
private static double betaFraction(double x, double p, double q) {
int m, m2;
double sum_pq, p_plus, p_minus, c =1.0 , d, delta, h, frac;
sum_pq = p + q;
p_plus = p + 1.0;
p_minus = p - 1.0;
h=1.0-sum_pq*x/p_plus;
if (Math.abs(h) < XMININ)
h=XMININ;
h=1.0/h;
frac = h;
m=1;
delta = 0.0;
while (m <= MAX_ITERATIONS && Math.abs(delta-1.0) > PRECISION ) {
m2=2*m;
// even index for d
d=m*(q-m)*x/((p_minus+m2)*(p+m2));
h=1.0+d*h;
if (Math.abs(h) < XMININ)
h=XMININ;
h=1.0/h;
c=1.0+d/c;
if (Math.abs(c) < XMININ)
c=XMININ;
frac *= h*c;
// odd index for d
d = -(p+m)*(sum_pq+m)*x/((p+m2)*(p_plus+m2));
h=1.0+d*h;
if (Math.abs(h) < XMININ)
h=XMININ;
h=1.0/h;
c=1.0+d/c;
if (Math.abs(c) < XMININ)
c=XMININ;
delta=h*c;
frac *= delta;
m++;
}
return frac;
}
// ====================================================
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
//
// Developed at SunSoft, a Sun Microsystems, Inc. business.
// Permission to use, copy, modify, and distribute this
// software is freely granted, provided that this notice
// is preserved.
// ====================================================
//
// x
// 2 |\
// erf(x) = --------- | exp(-t*t)dt
// sqrt(pi) \|
// 0
//
// erfc(x) = 1-erf(x)
// Note that
// erf(-x) = -erf(x)
// erfc(-x) = 2 - erfc(x)
//
// Method:
// 1. For |x| in [0, 0.84375]
// erf(x) = x + x*R(x^2)
// erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
// = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
// where R = P/Q where P is an odd poly of degree 8 and
// Q is an odd poly of degree 10.
// -57.90
// | R - (erf(x)-x)/x | <= 2
//
//
// Remark. The formula is derived by noting
// erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
// and that
// 2/sqrt(pi) = 1.128379167095512573896158903121545171688
// is close to one. The interval is chosen because the fix
// point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
// near 0.6174), and by some experiment, 0.84375 is chosen to
// guarantee the error is less than one ulp for erf.
//
// 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
// c = 0.84506291151 rounded to single (24 bits)
// erf(x) = sign(x) * (c + P1(s)/Q1(s))
// erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
// 1+(c+P1(s)/Q1(s)) if x < 0
// |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
// Remark: here we use the taylor series expansion at x=1.
// erf(1+s) = erf(1) + s*Poly(s)
// = 0.845.. + P1(s)/Q1(s)
// That is, we use rational approximation to approximate
// erf(1+s) - (c = (single)0.84506291151)
// Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
// where
// P1(s) = degree 6 poly in s
// Q1(s) = degree 6 poly in s
//
// 3. For x in [1.25,1/0.35(~2.857143)],
// erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
// erf(x) = 1 - erfc(x)
// where
// R1(z) = degree 7 poly in z, (z=1/x^2)
// S1(z) = degree 8 poly in z
//
// 4. For x in [1/0.35,28]
// erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
// = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
// = 2.0 - tiny (if x <= -6)
// erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
// erf(x) = sign(x)*(1.0 - tiny)
// where
// R2(z) = degree 6 poly in z, (z=1/x^2)
// S2(z) = degree 7 poly in z
//
// Note1:
// To compute exp(-x*x-0.5625+R/S), let s be a single
// precision number and s := x; then
// -x*x = -s*s + (s-x)*(s+x)
// exp(-x*x-0.5626+R/S) =
// exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
// Note2:
// Here 4 and 5 make use of the asymptotic series
// exp(-x*x)
// erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
// x*sqrt(pi)
// We use rational approximation to approximate
// g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
// Here is the error bound for R1/S1 and R2/S2
// |R1/S1 - f(x)| < 2**(-62.57)
// |R2/S2 - f(x)| < 2**(-61.52)
//
// 5. For inf > x >= 28
// erf(x) = sign(x) *(1 - tiny) (raise inexact)
// erfc(x) = tiny*tiny (raise underflow) if x > 0
// = 2 - tiny if x<0
//
// 7. Special case:
// erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
// erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
// erfc/erf(NaN) is NaN
//
/**
* Error function.
* Based on C-code for the error function developed at Sun Microsystems.
* @author Jaco van Kooten
*/
public static double error(double x) {
// Coefficients for approximation to erf on [0,0.84375]
final double e_efx=1.28379167095512586316e-01;
//double efx8=1.02703333676410069053e00;
final double ePp[]={
1.28379167095512558561e-01,
-3.25042107247001499370e-01,
-2.84817495755985104766e-02,
-5.77027029648944159157e-03,
-2.37630166566501626084e-05};
final double eQq[]={
3.97917223959155352819e-01,
6.50222499887672944485e-02,
5.08130628187576562776e-03,
1.32494738004321644526e-04,
-3.96022827877536812320e-06};
// Coefficients for approximation to erf in [0.84375,1.25]
final double ePa[]={
-2.36211856075265944077e-03,
4.14856118683748331666e-01,
-3.72207876035701323847e-01,
3.18346619901161753674e-01,
-1.10894694282396677476e-01,
3.54783043256182359371e-02,
-2.16637559486879084300e-03};
final double eQa[]={
1.06420880400844228286e-01,
5.40397917702171048937e-01,
7.18286544141962662868e-02,
1.26171219808761642112e-01,
1.36370839120290507362e-02,
1.19844998467991074170e-02};
final double e_erx=8.45062911510467529297e-01;
double P,Q,s,retval;
final double abs_x = (x >= 0.0 ? x : -x);
if ( abs_x < 0.84375 ) { // 0 < |x| < 0.84375
if (abs_x < 3.7252902984619141e-9 ) // |x| < 2**-28
retval = abs_x + abs_x*e_efx;
else {
s = x*x;
P = ePp[0]+s*(ePp[1]+s*(ePp[2]+s*(ePp[3]+s*ePp[4])));
Q = 1.0+s*(eQq[0]+s*(eQq[1]+s*(eQq[2]+s*(eQq[3]+s*eQq[4]))));
retval = abs_x + abs_x*(P/Q);
}
} else if (abs_x < 1.25) { // 0.84375 < |x| < 1.25
s = abs_x-1.0;
P = ePa[0]+s*(ePa[1]+s*(ePa[2]+s*(ePa[3]+s*(ePa[4]+s*(ePa[5]+s*ePa[6])))));
Q = 1.0+s*(eQa[0]+s*(eQa[1]+s*(eQa[2]+s*(eQa[3]+s*(eQa[4]+s*eQa[5])))));
retval = e_erx + P/Q;
} else if (abs_x >= 6.0)
retval = 1.0;
else // 1.25 < |x| < 6.0
retval = 1.0-complementaryError(abs_x);
return (x >= 0.0) ? retval : -retval;
}
/**
* Complementary error function.
* Based on C-code for the error function developed at Sun Microsystems.
* @author Jaco van Kooten
*/
public static double complementaryError(double x) {
// Coefficients for approximation to erfc in [1.25,1/.35]
final double eRa[]={
-9.86494403484714822705e-03,
-6.93858572707181764372e-01,
-1.05586262253232909814e01,
-6.23753324503260060396e01,
-1.62396669462573470355e02,
-1.84605092906711035994e02,
-8.12874355063065934246e01,
-9.81432934416914548592e00};
final double eSa[]={
1.96512716674392571292e01,
1.37657754143519042600e02,
4.34565877475229228821e02,
6.45387271733267880336e02,
4.29008140027567833386e02,
1.08635005541779435134e02,
6.57024977031928170135e00,
-6.04244152148580987438e-02};
// Coefficients for approximation to erfc in [1/.35,28]
final double eRb[]={
-9.86494292470009928597e-03,
-7.99283237680523006574e-01,
-1.77579549177547519889e01,
-1.60636384855821916062e02,
-6.37566443368389627722e02,
-1.02509513161107724954e03,
-4.83519191608651397019e02};
final double eSb[]={
3.03380607434824582924e01,
3.25792512996573918826e02,
1.53672958608443695994e03,
3.19985821950859553908e03,
2.55305040643316442583e03,
4.74528541206955367215e02,
-2.24409524465858183362e01};
double s,retval,R,S;
final double abs_x =(x>=0.0 ? x : -x);
if (abs_x < 1.25)
retval = 1.0-error(abs_x);
else if (abs_x > 28.0)
retval=0.0;
else { // 1.25 < |x| < 28
s = 1.0/(abs_x*abs_x);
if (abs_x < 2.8571428) { // ( |x| < 1/0.35 )
R=eRa[0]+s*(eRa[1]+s*(eRa[2]+s*(eRa[3]+s*(eRa[4]+s*(eRa[5]+s*(eRa[6]+s*eRa[7]))))));
S=1.0+s*(eSa[0]+s*(eSa[1]+s*(eSa[2]+s*(eSa[3]+s*(eSa[4]+s*(eSa[5]+s*(eSa[6]+s*eSa[7])))))));
} else { // ( |x| > 1/0.35 )
R=eRb[0]+s*(eRb[1]+s*(eRb[2]+s*(eRb[3]+s*(eRb[4]+s*(eRb[5]+s*eRb[6])))));
S=1.0+s*(eSb[0]+s*(eSb[1]+s*(eSb[2]+s*(eSb[3]+s*(eSb[4]+s*(eSb[5]+s*eSb[6]))))));
}
retval = Math.exp(-x*x - 0.5625 + R/S)/abs_x;
}
return (x >= 0.0) ? retval : 2.0-retval;
}
}