WavesDefBU/ 0040755 0057434 0000145 00000000000 07634547076 0012014 5 ustar 00 mit 0000000 0011064 WavesDefBU/JSci/ 0040700 0057434 0000145 00000000000 07634547074 0012630 5 ustar 00 mit 0000000 0011064 WavesDefBU/JSci/maths/ 0040700 0057434 0000145 00000000000 07634547074 0013744 5 ustar 00 mit 0000000 0011064 WavesDefBU/JSci/maths/SpecialMath.java 0100600 0057434 0000145 00000205343 07634547074 0017007 0 ustar 00 mit 0000000 0011064 package JSci.maths;
/**
* The special function math library.
* This class cannot be subclassed or instantiated because all methods are static.
* @version 1.0
* @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) {
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
* Applied Mathematics Division
* Argonne National Laboratory
* Argonne, IL 60439
*
* References: *
* 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.
*
* References: *
* 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.
*