JSci.maths.SpecialMath Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of jsci Show documentation
Show all versions of jsci Show documentation
JSci is a set of open source Java packages. The aim is to encapsulate scientific methods/principles in the most natural way possible. As such they should greatly aid the development of scientific based software.
It offers: abstract math interfaces, linear algebra (support for various matrix and vector types), statistics (including probability distributions), wavelets, newtonian mechanics, chart/graph components (AWT and Swing), MathML DOM implementation, ...
Note: some packages, like javax.comm, for the astro and instruments package aren't listed as dependencies (not available).
The newest version!
package JSci.maths;
/**
* The special function math library.
* This class cannot be subclassed or instantiated because all methods are static.
* @version 1.2
* @author Mark Hale
*/
public final class SpecialMath extends AbstractMath implements NumericalConstants {
private SpecialMath() {}
// 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) {
return airyModPhase(x);
} 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 the real part, i.e. modulus*cos(phase).
*/
private static double airyModPhase(double x) {
double modulus, phase;
if(x < -2.0) {
double z = 16.0/(x*x*x)+1.0;
modulus = 0.3125+chebyshev(z, am21cs);
phase = -0.625+chebyshev(z, ath1cs);
} else {
double z = (16.0/(x*x*x)+9.0)/7.0;
modulus = 0.3125+chebyshev(z, am22cs);
phase = -0.625+chebyshev(z, ath2cs);
}
final double sqrtx = Math.sqrt(-x);
modulus = Math.sqrt(modulus/sqrtx);
phase = Math.PI/4.0-x*sqrtx*phase;
return modulus*Math.cos(phase);
}
/**
* 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) {
return airyModPhase(x);
} 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 related constants
private final static double g_p[] = { -1.71618513886549492533811,
24.7656508055759199108314, -379.804256470945635097577,
629.331155312818442661052, 866.966202790413211295064,
-31451.2729688483675254357, -36144.4134186911729807069,
66456.1438202405440627855 };
private final static double g_q[] = { -30.8402300119738975254353,
315.350626979604161529144, -1015.15636749021914166146,
-3107.77167157231109440444, 22538.1184209801510330112,
4755.84627752788110767815, -134659.959864969306392456,
-115132.259675553483497211 };
private final static double g_c[] = { -0.001910444077728,8.4171387781295e-4,
-5.952379913043012e-4, 7.93650793500350248e-4,
-0.002777777777777681622553, 0.08333333333333333331554247,
0.0057083835261 };
/**
* The largest argument for which gamma(x) is representable in the machine.
*/
public final static double GAMMA_X_MAX_VALUE = 171.624;
/**
* Gamma function.
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz
* Applied Mathematics Division
* Argonne National Laboratory
* Argonne, IL 60439
*
* References:
*
* - "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.
*
- Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
*
* From the original documentation:
*
* 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.
*
* Error returns:
* The program returns the value XINF for singularities or when overflow would occur.
* The computation is believed to be free of underflow and overflow.
*
* @return Double.MAX_VALUE if overflow would occur, i.e. if abs(x) > 171.624
* @jsci.planetmath GammaFunction
* @author Jaco van Kooten
*/
public static double gamma(double x) {
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 <= GAMMA_X_MAX_VALUE) {
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.
*/
public final static double LOG_GAMMA_X_MAX_VALUE = 2.55e305;
// Log Gamma related constants
private final static double lg_d1 = -0.5772156649015328605195174;
private final static double lg_d2 = 0.4227843350984671393993777;
private final static double lg_d4 = 1.791759469228055000094023;
private final static double lg_p1[] = { 4.945235359296727046734888,
201.8112620856775083915565, 2290.838373831346393026739,
11319.67205903380828685045, 28557.24635671635335736389,
38484.96228443793359990269, 26377.48787624195437963534,
7225.813979700288197698961 };
private final static double lg_q1[] = { 67.48212550303777196073036,
1113.332393857199323513008, 7738.757056935398733233834,
27639.87074403340708898585, 54993.10206226157329794414,
61611.22180066002127833352, 36351.27591501940507276287,
8785.536302431013170870835 };
private final static double lg_p2[] = { 4.974607845568932035012064,
542.4138599891070494101986, 15506.93864978364947665077,
184793.2904445632425417223, 1088204.76946882876749847,
3338152.967987029735917223, 5106661.678927352456275255,
3074109.054850539556250927 };
private final static double lg_q2[] = { 183.0328399370592604055942,
7765.049321445005871323047, 133190.3827966074194402448,
1136705.821321969608938755, 5267964.117437946917577538,
13467014.54311101692290052, 17827365.30353274213975932,
9533095.591844353613395747 };
private final static double lg_p4[] = { 14745.02166059939948905062,
2426813.369486704502836312, 121475557.4045093227939592,
2663432449.630976949898078, 29403789566.34553899906876,
170266573776.5398868392998, 492612579337.743088758812,
560625185622.3951465078242 };
private final static double lg_q4[] = { 2690.530175870899333379843,
639388.5654300092398984238, 41355999.30241388052042842,
1120872109.61614794137657, 14886137286.78813811542398,
101680358627.2438228077304, 341747634550.7377132798597,
446315818741.9713286462081 };
private final static 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
private final static double lg_frtbig = 2.25e76;
private final static double pnt68 = 0.6796875;
// Function cache for logGamma
private static final ThreadLocal logGammaCache_res=new ThreadLocal() {
protected Object initialValue() {
return new Double(0.0);
}
};
private static final ThreadLocal logGammaCache_x=new ThreadLocal() {
protected Object initialValue() {
return new Double(0.0);
}
};
/**
* The natural logarithm of the gamma function.
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz
* Applied Mathematics Division
* Argonne National Laboratory
* Argonne, IL 60439
*
* References:
*
* - W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.
*
- K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.
*
- Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.
*
* From the original documentation:
*
* 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.
*
* Error returns:
* 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.
*
* @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) {
double xden, corr, xnum;
int i;
double y, xm1, xm2, xm4, res, ysq;
if (x == ((Double) logGammaCache_x.get()).doubleValue())
return ((Double) logGammaCache_res.get()).doubleValue();
y = x;
if (y > 0.0 && y <= LOG_GAMMA_X_MAX_VALUE) {
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.set(new Double(x));
logGammaCache_res.set(new Double(res));
return res;
}
private final static int MAX_ITERATIONS = 1000000000;
// lower value = higher precision
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 > LOG_GAMMA_X_MAX_VALUE)
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));
}
throw new RuntimeException("Maximum iterations exceeded: please file a bug report.");
}
/**
* @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; iPRECISION; 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) > LOG_GAMMA_X_MAX_VALUE)
return 0.0;
else
return Math.exp(logBeta(p,q));
}
// Function cache for logBeta
private static final ThreadLocal logBetaCache_res=new ThreadLocal() {
protected Object initialValue() {
return new Double(0.0);
}
};
private static final ThreadLocal logBetaCache_p=new ThreadLocal() {
protected Object initialValue() {
return new Double(0.0);
}
};
private static final ThreadLocal logBetaCache_q=new ThreadLocal() {
protected Object initialValue() {
return new Double(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 != ((Double) logBetaCache_p.get()).doubleValue()
|| q != ((Double) logBetaCache_q.get()).doubleValue()) {
logBetaCache_p.set(new Double(p));
logBetaCache_q.set(new Double(q));
double res;
if (p <= 0.0 || q <= 0.0 || (p+q) > LOG_GAMMA_X_MAX_VALUE)
res = 0.0;
else
res = logGamma(p)+logGamma(q)-logGamma(p+q);
logBetaCache_res.set(new Double(res));
return res;
} else {
return ((Double) logBetaCache_res.get()).doubleValue();
}
}
/**
* 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) > LOG_GAMMA_X_MAX_VALUE)
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 >= 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
//
// Coefficients for approximation to erf on [0,0.84375]
private final static double e_efx=1.28379167095512586316e-01;
// private final static double efx8=1.02703333676410069053e00;
private final static double ePp[]={
1.28379167095512558561e-01,
-3.25042107247001499370e-01,
-2.84817495755985104766e-02,
-5.77027029648944159157e-03,
-2.37630166566501626084e-05};
private final static 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]
private final static double ePa[]={
-2.36211856075265944077e-03,
4.14856118683748331666e-01,
-3.72207876035701323847e-01,
3.18346619901161753674e-01,
-1.10894694282396677476e-01,
3.54783043256182359371e-02,
-2.16637559486879084300e-03};
private final static double eQa[]={
1.06420880400844228286e-01,
5.40397917702171048937e-01,
7.18286544141962662868e-02,
1.26171219808761642112e-01,
1.36370839120290507362e-02,
1.19844998467991074170e-02};
private final static double e_erx=8.45062911510467529297e-01;
/**
* Error function.
* Based on C-code for the error function developed at Sun Microsystems.
* @author Jaco van Kooten
*/
public static double error(double x) {
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;
}
// Coefficients for approximation to erfc in [1.25,1/.35]
private final static double eRa[]={
-9.86494403484714822705e-03,
-6.93858572707181764372e-01,
-1.05586262253232909814e01,
-6.23753324503260060396e01,
-1.62396669462573470355e02,
-1.84605092906711035994e02,
-8.12874355063065934246e01,
-9.81432934416914548592e00};
private final static 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]
private final static double eRb[]={
-9.86494292470009928597e-03,
-7.99283237680523006574e-01,
-1.77579549177547519889e01,
-1.60636384855821916062e02,
-6.37566443368389627722e02,
-1.02509513161107724954e03,
-4.83519191608651397019e02};
private final static double eSb[]={
3.03380607434824582924e01,
3.25792512996573918826e02,
1.53672958608443695994e03,
3.19985821950859553908e03,
2.55305040643316442583e03,
4.74528541206955367215e02,
-2.24409524465858183362e01};
/**
* 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) {
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;
}
}