0001: package JSci.maths;
0002:
0003: /**
0004: * The special function math library.
0005: * This class cannot be subclassed or instantiated because all methods are static.
0006: * @version 1.1
0007: * @author Mark Hale
0008: */
0009: public final class SpecialMath extends AbstractMath implements
0010: NumericalConstants {
0011: private SpecialMath() {
0012: }
0013:
0014: // Some IEEE machine constants
0015: /**
0016: * Relative machine precision.
0017: */
0018: private final static double EPS = 2.22e-16;
0019: /**
0020: * The smallest positive floating-point number such that 1/xminin is machine representable.
0021: */
0022: private final static double XMININ = 2.23e-308;
0023:
0024: // CHEBYSHEV SERIES
0025:
0026: // series for ai0 on the interval 1.25000d-01 to 3.33333d-01
0027: // with weighted error 7.87e-17
0028: // log weighted error 16.10
0029: // significant figures required 14.69
0030: // decimal places required 16.76
0031:
0032: private final static double ai0cs[] = { 0.07575994494023796,
0033: 0.00759138081082334, 0.00041531313389237,
0034: 0.00001070076463439, -0.00000790117997921,
0035: -0.00000078261435014, 0.00000027838499429,
0036: 0.00000000825247260, -0.00000001204463945,
0037: 0.00000000155964859, 0.00000000022925563,
0038: -0.00000000011916228, 0.00000000001757854,
0039: 0.00000000000112822, -0.00000000000114684,
0040: 0.00000000000027155, -0.00000000000002415,
0041: -0.00000000000000608, 0.00000000000000314,
0042: -0.00000000000000071, 0.00000000000000007 };
0043:
0044: // series for ai02 on the interval 0. to 1.25000d-01
0045: // with weighted error 3.79e-17
0046: // log weighted error 16.42
0047: // significant figures required 14.86
0048: // decimal places required 17.09
0049: private final static double ai02cs[] = { 0.05449041101410882,
0050: 0.00336911647825569, 0.00006889758346918,
0051: 0.00000289137052082, 0.00000020489185893,
0052: 0.00000002266668991, 0.00000000339623203,
0053: 0.00000000049406022, 0.00000000001188914,
0054: -0.00000000003149915, -0.00000000001321580,
0055: -0.00000000000179419, 0.00000000000071801,
0056: 0.00000000000038529, 0.00000000000001539,
0057: -0.00000000000004151, -0.00000000000000954,
0058: 0.00000000000000382, 0.00000000000000176,
0059: -0.00000000000000034, -0.00000000000000027,
0060: 0.00000000000000003 };
0061:
0062: // series for ai1 on the interval 1.25000d-01 to 3.33333d-01
0063: // with weighted error 6.98e-17
0064: // log weighted error 16.16
0065: // significant figures required 14.53
0066: // decimal places required 16.82
0067:
0068: private final static double ai1cs[] = { -0.02846744181881479,
0069: -0.01922953231443221, -0.00061151858579437,
0070: -0.00002069971253350, 0.00000858561914581,
0071: 0.00000104949824671, -0.00000029183389184,
0072: -0.00000001559378146, 0.00000001318012367,
0073: -0.00000000144842341, -0.00000000029085122,
0074: 0.00000000012663889, -0.00000000001664947,
0075: -0.00000000000166665, 0.00000000000124260,
0076: -0.00000000000027315, 0.00000000000002023,
0077: 0.00000000000000730, -0.00000000000000333,
0078: 0.00000000000000071, -0.00000000000000006 };
0079:
0080: // series for ai12 on the interval 0. to 1.25000d-01
0081: // with weighted error 3.55e-17
0082: // log weighted error 16.45
0083: // significant figures required 14.69
0084: // decimal places required 17.12
0085:
0086: private final static double ai12cs[] = { 0.02857623501828014,
0087: -0.00976109749136147, -0.00011058893876263,
0088: -0.00000388256480887, -0.00000025122362377,
0089: -0.00000002631468847, -0.00000000383538039,
0090: -0.00000000055897433, -0.00000000001897495,
0091: 0.00000000003252602, 0.00000000001412580,
0092: 0.00000000000203564, -0.00000000000071985,
0093: -0.00000000000040836, -0.00000000000002101,
0094: 0.00000000000004273, 0.00000000000001041,
0095: -0.00000000000000382, -0.00000000000000186,
0096: 0.00000000000000033, 0.00000000000000028,
0097: -0.00000000000000003 };
0098:
0099: // series for aif on the interval -1.00000d+00 to 1.00000d+00
0100: // with weighted error 1.09e-19
0101: // log weighted error 18.96
0102: // significant figures required 17.76
0103: // decimal places required 19.44
0104:
0105: private final static double aifcs[] = { -0.03797135849666999750,
0106: 0.05919188853726363857, 0.00098629280577279975,
0107: 0.00000684884381907656, 0.00000002594202596219,
0108: 0.00000000006176612774, 0.00000000000010092454,
0109: 0.00000000000000012014, 0.00000000000000000010 };
0110:
0111: // series for aig on the interval -1.00000d+00 to 1.00000d+00
0112: // with weighted error 1.51e-17
0113: // log weighted error 16.82
0114: // significant figures required 15.19
0115: // decimal places required 17.27
0116:
0117: private final static double aigcs[] = { 0.01815236558116127,
0118: 0.02157256316601076, 0.00025678356987483,
0119: 0.00000142652141197, 0.00000000457211492,
0120: 0.00000000000952517, 0.00000000000001392,
0121: 0.00000000000000001 };
0122:
0123: // series for aip on the interval 0. to 1.00000d+00
0124: // with weighted error 5.10e-17
0125: // log weighted error 16.29
0126: // significant figures required 14.41
0127: // decimal places required 17.06
0128:
0129: private final static double aipcs[] = { -0.0187519297793868,
0130: -0.0091443848250055, 0.0009010457337825,
0131: -0.0001394184127221, 0.0000273815815785,
0132: -0.0000062750421119, 0.0000016064844184,
0133: -0.0000004476392158, 0.0000001334635874,
0134: -0.0000000420735334, 0.0000000139021990,
0135: -0.0000000047831848, 0.0000000017047897,
0136: -0.0000000006268389, 0.0000000002369824,
0137: -0.0000000000918641, 0.0000000000364278,
0138: -0.0000000000147475, 0.0000000000060851,
0139: -0.0000000000025552, 0.0000000000010906,
0140: -0.0000000000004725, 0.0000000000002076,
0141: -0.0000000000000924, 0.0000000000000417,
0142: -0.0000000000000190, 0.0000000000000087,
0143: -0.0000000000000040, 0.0000000000000019,
0144: -0.0000000000000009, 0.0000000000000004,
0145: -0.0000000000000002, 0.0000000000000001,
0146: -0.0000000000000000 };
0147:
0148: // series for am21 on the interval -1.25000d-01 to 0.
0149: // with weighted error 2.89e-17
0150: // log weighted error 16.54
0151: // significant figures required 14.15
0152: // decimal places required 17.34
0153:
0154: private final static double am21cs[] = { 0.0065809191761485,
0155: 0.0023675984685722, 0.0001324741670371, 0.0000157600904043,
0156: 0.0000027529702663, 0.0000006102679017, 0.0000001595088468,
0157: 0.0000000471033947, 0.0000000152933871, 0.0000000053590722,
0158: 0.0000000020000910, 0.0000000007872292, 0.0000000003243103,
0159: 0.0000000001390106, 0.0000000000617011, 0.0000000000282491,
0160: 0.0000000000132979, 0.0000000000064188, 0.0000000000031697,
0161: 0.0000000000015981, 0.0000000000008213, 0.0000000000004296,
0162: 0.0000000000002284, 0.0000000000001232, 0.0000000000000675,
0163: 0.0000000000000374, 0.0000000000000210, 0.0000000000000119,
0164: 0.0000000000000068, 0.0000000000000039, 0.0000000000000023,
0165: 0.0000000000000013, 0.0000000000000008, 0.0000000000000005,
0166: 0.0000000000000003, 0.0000000000000001, 0.0000000000000001,
0167: 0.0000000000000000, 0.0000000000000000, 0.0000000000000000 };
0168:
0169: // series for ath1 on the interval -1.25000d-01 to 0.
0170: // with weighted error 2.53e-17
0171: // log weighted error 16.60
0172: // significant figures required 15.15
0173: // decimal places required 17.38
0174:
0175: private final static double ath1cs[] = { -0.07125837815669365,
0176: -0.00590471979831451, -0.00012114544069499,
0177: -0.00000988608542270, -0.00000138084097352,
0178: -0.00000026142640172, -0.00000006050432589,
0179: -0.00000001618436223, -0.00000000483464911,
0180: -0.00000000157655272, -0.00000000055231518,
0181: -0.00000000020545441, -0.00000000008043412,
0182: -0.00000000003291252, -0.00000000001399875,
0183: -0.00000000000616151, -0.00000000000279614,
0184: -0.00000000000130428, -0.00000000000062373,
0185: -0.00000000000030512, -0.00000000000015239,
0186: -0.00000000000007758, -0.00000000000004020,
0187: -0.00000000000002117, -0.00000000000001132,
0188: -0.00000000000000614, -0.00000000000000337,
0189: -0.00000000000000188, -0.00000000000000105,
0190: -0.00000000000000060, -0.00000000000000034,
0191: -0.00000000000000020, -0.00000000000000011,
0192: -0.00000000000000007, -0.00000000000000004,
0193: -0.00000000000000002 };
0194:
0195: // series for am22 on the interval -1.00000d+00 to -1.25000d-01
0196: // with weighted error 2.99e-17
0197: // log weighted error 16.52
0198: // significant figures required 14.57
0199: // decimal places required 17.28
0200:
0201: private final static double am22cs[] = { -0.01562844480625341,
0202: 0.00778336445239681, 0.00086705777047718,
0203: 0.00015696627315611, 0.00003563962571432,
0204: 0.00000924598335425, 0.00000262110161850,
0205: 0.00000079188221651, 0.00000025104152792,
0206: 0.00000008265223206, 0.00000002805711662,
0207: 0.00000000976821090, 0.00000000347407923,
0208: 0.00000000125828132, 0.00000000046298826,
0209: 0.00000000017272825, 0.00000000006523192,
0210: 0.00000000002490471, 0.00000000000960156,
0211: 0.00000000000373448, 0.00000000000146417,
0212: 0.00000000000057826, 0.00000000000022991,
0213: 0.00000000000009197, 0.00000000000003700,
0214: 0.00000000000001496, 0.00000000000000608,
0215: 0.00000000000000248, 0.00000000000000101,
0216: 0.00000000000000041, 0.00000000000000017,
0217: 0.00000000000000007, 0.00000000000000002 };
0218:
0219: // series for ath2 on the interval -1.00000d+00 to -1.25000d-01
0220: // with weighted error 2.57e-17
0221: // log weighted error 16.59
0222: // significant figures required 15.07
0223: // decimal places required 17.34
0224:
0225: private final static double ath2cs[] = { 0.00440527345871877,
0226: -0.03042919452318455, -0.00138565328377179,
0227: -0.00018044439089549, -0.00003380847108327,
0228: -0.00000767818353522, -0.00000196783944371,
0229: -0.00000054837271158, -0.00000016254615505,
0230: -0.00000005053049981, -0.00000001631580701,
0231: -0.00000000543420411, -0.00000000185739855,
0232: -0.00000000064895120, -0.00000000023105948,
0233: -0.00000000008363282, -0.00000000003071196,
0234: -0.00000000001142367, -0.00000000000429811,
0235: -0.00000000000163389, -0.00000000000062693,
0236: -0.00000000000024260, -0.00000000000009461,
0237: -0.00000000000003716, -0.00000000000001469,
0238: -0.00000000000000584, -0.00000000000000233,
0239: -0.00000000000000093, -0.00000000000000037,
0240: -0.00000000000000015, -0.00000000000000006,
0241: -0.00000000000000002 };
0242:
0243: // series for bi0 on the interval 0. to 9.00000d+00
0244: // with weighted error 2.46e-18
0245: // log weighted error 17.61
0246: // significant figures required 17.90
0247: // decimal places required 18.15
0248:
0249: private final static double bi0cs[] = { -0.07660547252839144951,
0250: 1.927337953993808270, 0.2282644586920301339,
0251: 0.01304891466707290428, 0.00043442709008164874,
0252: 0.00000942265768600193, 0.00000014340062895106,
0253: 0.00000000161384906966, 0.00000000001396650044,
0254: 0.00000000000009579451, 0.00000000000000053339,
0255: 0.00000000000000000245 };
0256:
0257: // series for bj0 on the interval 0. to 1.60000d+01
0258: // with weighted error 7.47e-18
0259: // log weighted error 17.13
0260: // significant figures required 16.98
0261: // decimal places required 17.68
0262:
0263: private final static double bj0cs[] = { 0.100254161968939137,
0264: -0.665223007764405132, 0.248983703498281314,
0265: -0.0332527231700357697, 0.0023114179304694015,
0266: -0.0000991127741995080, 0.0000028916708643998,
0267: -0.0000000612108586630, 0.0000000009838650793,
0268: -0.0000000000124235515, 0.0000000000001265433,
0269: -0.0000000000000010619, 0.0000000000000000074 };
0270:
0271: // series for bm0 on the interval 0. to 6.25000d-02
0272: // with weighted error 4.98e-17
0273: // log weighted error 16.30
0274: // significant figures required 14.97
0275: // decimal places required 16.96
0276:
0277: private final static double bm0cs[] = { 0.09284961637381644,
0278: -0.00142987707403484, 0.00002830579271257,
0279: -0.00000143300611424, 0.00000012028628046,
0280: -0.00000001397113013, 0.00000000204076188,
0281: -0.00000000035399669, 0.00000000007024759,
0282: -0.00000000001554107, 0.00000000000376226,
0283: -0.00000000000098282, 0.00000000000027408,
0284: -0.00000000000008091, 0.00000000000002511,
0285: -0.00000000000000814, 0.00000000000000275,
0286: -0.00000000000000096, 0.00000000000000034,
0287: -0.00000000000000012, 0.00000000000000004 };
0288:
0289: // series for bth0 on the interval 0. to 6.25000d-02
0290: // with weighted error 3.67e-17
0291: // log weighted error 16.44
0292: // significant figures required 15.53
0293: // decimal places required 17.13
0294:
0295: private final static double bth0cs[] = { -0.24639163774300119,
0296: 0.001737098307508963, -0.000062183633402968,
0297: 0.000004368050165742, -0.000000456093019869,
0298: 0.000000062197400101, -0.000000010300442889,
0299: 0.000000001979526776, -0.000000000428198396,
0300: 0.000000000102035840, -0.000000000026363898,
0301: 0.000000000007297935, -0.000000000002144188,
0302: 0.000000000000663693, -0.000000000000215126,
0303: 0.000000000000072659, -0.000000000000025465,
0304: 0.000000000000009229, -0.000000000000003448,
0305: 0.000000000000001325, -0.000000000000000522,
0306: 0.000000000000000210, -0.000000000000000087,
0307: 0.000000000000000036 };
0308:
0309: // series for by0 on the interval 0. to 1.60000d+01
0310: // with weighted error 1.20e-17
0311: // log weighted error 16.92
0312: // significant figures required 16.15
0313: // decimal places required 17.48
0314:
0315: private final static double by0cs[] = { -0.011277839392865573,
0316: -0.12834523756042035, -0.10437884799794249,
0317: 0.023662749183969695, -0.002090391647700486,
0318: 0.000103975453939057, -0.000003369747162423,
0319: 0.000000077293842676, -0.000000001324976772,
0320: 0.000000000017648232, -0.000000000000188105,
0321: 0.000000000000001641, -0.000000000000000011 };
0322:
0323: // series for bi1 on the interval 0. to 9.00000d+00
0324: // with weighted error 2.40e-17
0325: // log weighted error 16.62
0326: // significant figures required 16.23
0327: // decimal places required 17.14
0328:
0329: private final static double bi1cs[] = { -0.001971713261099859,
0330: 0.40734887667546481, 0.034838994299959456,
0331: 0.001545394556300123, 0.000041888521098377,
0332: 0.000000764902676483, 0.000000010042493924,
0333: 0.000000000099322077, 0.000000000000766380,
0334: 0.000000000000004741, 0.000000000000000024 };
0335:
0336: // series for bj1 on the interval 0. to 1.60000d+01
0337: // with weighted error 4.48e-17
0338: // log weighted error 16.35
0339: // significant figures required 15.77
0340: // decimal places required 16.89
0341:
0342: private final static double bj1cs[] = { -0.11726141513332787,
0343: -0.25361521830790640, 0.050127080984469569,
0344: -0.004631514809625081, 0.000247996229415914,
0345: -0.000008678948686278, 0.000000214293917143,
0346: -0.000000003936093079, 0.000000000055911823,
0347: -0.000000000000632761, 0.000000000000005840,
0348: -0.000000000000000044 };
0349:
0350: // series for bm1 on the interval 0. to 6.25000d-02
0351: // with weighted error 5.61e-17
0352: // log weighted error 16.25
0353: // significant figures required 14.97
0354: // decimal places required 16.91
0355:
0356: private final static double bm1cs[] = { 0.1047362510931285,
0357: 0.00442443893702345, -0.00005661639504035,
0358: 0.00000231349417339, -0.00000017377182007,
0359: 0.00000001893209930, -0.00000000265416023,
0360: 0.00000000044740209, -0.00000000008691795,
0361: 0.00000000001891492, -0.00000000000451884,
0362: 0.00000000000116765, -0.00000000000032265,
0363: 0.00000000000009450, -0.00000000000002913,
0364: 0.00000000000000939, -0.00000000000000315,
0365: 0.00000000000000109, -0.00000000000000039,
0366: 0.00000000000000014, -0.00000000000000005 };
0367:
0368: // series for bth1 on the interval 0. to 6.25000d-02
0369: // with weighted error 4.10e-17
0370: // log weighted error 16.39
0371: // significant figures required 15.96
0372: // decimal places required 17.08
0373:
0374: private final static double bth1cs[] = { 0.74060141026313850,
0375: -0.004571755659637690, 0.000119818510964326,
0376: -0.000006964561891648, 0.000000655495621447,
0377: -0.000000084066228945, 0.000000013376886564,
0378: -0.000000002499565654, 0.000000000529495100,
0379: -0.000000000124135944, 0.000000000031656485,
0380: -0.000000000008668640, 0.000000000002523758,
0381: -0.000000000000775085, 0.000000000000249527,
0382: -0.000000000000083773, 0.000000000000029205,
0383: -0.000000000000010534, 0.000000000000003919,
0384: -0.000000000000001500, 0.000000000000000589,
0385: -0.000000000000000237, 0.000000000000000097,
0386: -0.000000000000000040 };
0387:
0388: // series for by1 on the interval 0. to 1.60000d+01
0389: // with weighted error 1.87e-18
0390: // log weighted error 17.73
0391: // significant figures required 17.83
0392: // decimal places required 18.30
0393:
0394: private final static double by1cs[] = { 0.03208047100611908629,
0395: 1.262707897433500450, 0.00649996189992317500,
0396: -0.08936164528860504117, 0.01325088122175709545,
0397: -0.00089790591196483523, 0.00003647361487958306,
0398: -0.00000100137438166600, 0.00000001994539657390,
0399: -0.00000000030230656018, 0.00000000000360987815,
0400: -0.00000000000003487488, 0.00000000000000027838,
0401: -0.00000000000000000186 };
0402:
0403: /**
0404: * Evaluates a Chebyshev series.
0405: * @param x value at which to evaluate series
0406: * @param series the coefficients of the series
0407: */
0408: public static double chebyshev(double x, double series[]) {
0409: double twox, b0 = 0.0, b1 = 0.0, b2 = 0.0;
0410: twox = 2 * x;
0411: for (int i = series.length - 1; i > -1; i--) {
0412: b2 = b1;
0413: b1 = b0;
0414: b0 = twox * b1 - b2 + series[i];
0415: }
0416: return 0.5 * (b0 - b2);
0417: }
0418:
0419: /**
0420: * Airy function.
0421: * Based on the NETLIB Fortran function ai written by W. Fullerton.
0422: */
0423: public static double airy(double x) {
0424: if (x < -1.0) {
0425: return airyModPhase(x);
0426: } else if (x > 1.0)
0427: return expAiry(x) * Math.exp(-2.0 * x * Math.sqrt(x) / 3.0);
0428: else {
0429: final double z = x * x * x;
0430: return 0.375 + (chebyshev(z, aifcs) - x
0431: * (0.25 + chebyshev(z, aigcs)));
0432: }
0433: }
0434:
0435: /**
0436: * Airy modulus and phase.
0437: * Based on the NETLIB Fortran subroutine r9aimp written by W. Fullerton.
0438: * @return the real part, i.e. modulus*cos(phase).
0439: */
0440: private static double airyModPhase(double x) {
0441: double modulus, phase;
0442: if (x < -2.0) {
0443: double z = 16.0 / (x * x * x) + 1.0;
0444: modulus = 0.3125 + chebyshev(z, am21cs);
0445: phase = -0.625 + chebyshev(z, ath1cs);
0446: } else {
0447: double z = (16.0 / (x * x * x) + 9.0) / 7.0;
0448: modulus = 0.3125 + chebyshev(z, am22cs);
0449: phase = -0.625 + chebyshev(z, ath2cs);
0450: }
0451: final double sqrtx = Math.sqrt(-x);
0452: modulus = Math.sqrt(modulus / sqrtx);
0453: phase = Math.PI / 4.0 - x * sqrtx * phase;
0454: return modulus * Math.cos(phase);
0455: }
0456:
0457: /**
0458: * Exponential scaled Airy function.
0459: * Based on the NETLIB Fortran function aie written by W. Fullerton.
0460: */
0461: private static double expAiry(double x) {
0462: if (x < -1.0) {
0463: return airyModPhase(x);
0464: } else if (x <= 1.0) {
0465: final double z = x * x * x;
0466: return 0.375
0467: + (chebyshev(z, aifcs) - x
0468: * (0.25 + chebyshev(z, aigcs)))
0469: * Math.exp(2.0 * x * Math.sqrt(x) / 3.0);
0470: } else {
0471: final double sqrtx = Math.sqrt(x);
0472: final double z = 2.0 / (x * sqrtx) - 1.0;
0473: return (0.28125 + chebyshev(z, aipcs)) / Math.sqrt(sqrtx);
0474: }
0475: }
0476:
0477: /**
0478: * Bessel function of first kind, order zero.
0479: * Based on the NETLIB Fortran function besj0 written by W. Fullerton.
0480: */
0481: public static double besselFirstZero(double x) {
0482: double y = Math.abs(x);
0483: if (y > 4.0) {
0484: final double z = 32 / (y * y) - 1;
0485: final double amplitude = (0.75 + chebyshev(z, bm0cs))
0486: / Math.sqrt(y);
0487: final double theta = y - Math.PI / 4.0
0488: + chebyshev(z, bth0cs) / y;
0489: return amplitude * Math.cos(theta);
0490: } else if (y == 0.0)
0491: return 1.0;
0492: else
0493: return chebyshev(0.125 * y * y - 1, bj0cs);
0494: }
0495:
0496: /**
0497: * Modified Bessel function of first kind, order zero.
0498: * Based on the NETLIB Fortran function besi0 written by W. Fullerton.
0499: */
0500: public static double modBesselFirstZero(double x) {
0501: double y = Math.abs(x);
0502: if (y > 3.0)
0503: return Math.exp(y) * expModBesselFirstZero(x);
0504: else
0505: return 2.75 + chebyshev(y * y / 4.5 - 1.0, bi0cs);
0506: }
0507:
0508: /**
0509: * Exponential scaled modified Bessel function of first kind, order zero.
0510: * Based on the NETLIB Fortran function besi0e written by W. Fullerton.
0511: */
0512: private static double expModBesselFirstZero(double x) {
0513: final double y = Math.abs(x);
0514: if (y > 3.0) {
0515: if (y > 8.0)
0516: return (0.375 + chebyshev(16.0 / y - 1.0, ai02cs))
0517: / Math.sqrt(y);
0518: else
0519: return (0.375 + chebyshev((48.0 / y - 11.0) / 5.0,
0520: ai0cs))
0521: / Math.sqrt(y);
0522: } else
0523: return Math.exp(-y)
0524: * (2.75 + chebyshev(y * y / 4.5 - 1.0, bi0cs));
0525: }
0526:
0527: /**
0528: * Bessel function of first kind, order one.
0529: * Based on the NETLIB Fortran function besj1 written by W. Fullerton.
0530: */
0531: public static double besselFirstOne(double x) {
0532: double y = Math.abs(x);
0533: if (y > 4.0) {
0534: final double z = 32.0 / (y * y) - 1.0;
0535: final double amplitude = (0.75 + chebyshev(z, bm1cs))
0536: / Math.sqrt(y);
0537: final double theta = y - 3.0 * Math.PI / 4.0
0538: + chebyshev(z, bth1cs) / y;
0539: return Math.abs(amplitude) * x * Math.cos(theta)
0540: / Math.abs(x);
0541: } else if (y == 0.0)
0542: return 0.0;
0543: else
0544: return x * (0.25 + chebyshev(0.125 * y * y - 1.0, bj1cs));
0545: }
0546:
0547: /**
0548: * Modified Bessel function of first kind, order one.
0549: * Based on the NETLIB Fortran function besi1 written by W. Fullerton.
0550: */
0551: public static double modBesselFirstOne(double x) {
0552: final double y = Math.abs(x);
0553: if (y > 3.0)
0554: return Math.exp(y) * expModBesselFirstOne(x);
0555: else if (y == 0.0)
0556: return 0.0;
0557: else
0558: return x * (0.875 + chebyshev(y * y / 4.5 - 1.0, bi1cs));
0559: }
0560:
0561: /**
0562: * Exponential scaled modified Bessel function of first kind, order one.
0563: * Based on the NETLIB Fortran function besi1e written by W. Fullerton.
0564: */
0565: private static double expModBesselFirstOne(double x) {
0566: final double y = Math.abs(x);
0567: if (y > 3.0) {
0568: if (y > 8.0)
0569: return x / y
0570: * (0.375 + chebyshev(16.0 / y - 1.0, ai12cs))
0571: / Math.sqrt(y);
0572: else
0573: return x
0574: / y
0575: * (0.375 + chebyshev((48.0 / y - 11.0) / 5.0,
0576: ai1cs)) / Math.sqrt(y);
0577: } else if (y == 0.0)
0578: return 0.0;
0579: else
0580: return Math.exp(-y) * x
0581: * (0.875 + chebyshev(y * y / 4.5 - 1.0, bi1cs));
0582: }
0583:
0584: /**
0585: * Bessel function of second kind, order zero.
0586: * Based on the NETLIB Fortran function besy0 written by W. Fullerton.
0587: */
0588: public static double besselSecondZero(double x) {
0589: if (x > 4.0) {
0590: final double z = 32.0 / (x * x) - 1.0;
0591: final double amplitude = (0.75 + chebyshev(z, bm0cs))
0592: / Math.sqrt(x);
0593: final double theta = x - Math.PI / 4 + chebyshev(z, bth0cs)
0594: / x;
0595: return amplitude * Math.sin(theta);
0596: } else
0597: return (Math.log(0.5) + Math.log(x)) * besselFirstZero(x)
0598: + 0.375 + chebyshev(0.125 * x * x - 1.0, by0cs)
0599: * 2.0 / Math.PI;
0600: }
0601:
0602: /**
0603: * Bessel function of second kind, order one.
0604: * Based on the NETLIB Fortran function besy1 written by W. Fullerton.
0605: */
0606: public static double besselSecondOne(double x) {
0607: if (x > 4.0) {
0608: final double z = 32.0 / (x * x) - 1.0;
0609: final double amplitude = (0.75 + chebyshev(z, bm1cs))
0610: / Math.sqrt(x);
0611: final double theta = x - 3.0 * Math.PI / 4.0
0612: + chebyshev(z, bth1cs) / x;
0613: return amplitude * Math.sin(theta);
0614: } else
0615: return 2.0 * Math.log(0.5 * x) * besselFirstOne(x)
0616: / Math.PI
0617: + (0.5 + chebyshev(0.125 * x * x - 1.0, by1cs)) / x;
0618: }
0619:
0620: private final static double LOGSQRT2PI = Math.log(SQRT2PI);
0621:
0622: // Gamma function related constants
0623: private final static double g_p[] = { -1.71618513886549492533811,
0624: 24.7656508055759199108314, -379.804256470945635097577,
0625: 629.331155312818442661052, 866.966202790413211295064,
0626: -31451.2729688483675254357, -36144.4134186911729807069,
0627: 66456.1438202405440627855 };
0628: private final static double g_q[] = { -30.8402300119738975254353,
0629: 315.350626979604161529144, -1015.15636749021914166146,
0630: -3107.77167157231109440444, 22538.1184209801510330112,
0631: 4755.84627752788110767815, -134659.959864969306392456,
0632: -115132.259675553483497211 };
0633: private final static double g_c[] = { -0.001910444077728,
0634: 8.4171387781295e-4, -5.952379913043012e-4,
0635: 7.93650793500350248e-4, -0.002777777777777681622553,
0636: 0.08333333333333333331554247, 0.0057083835261 };
0637: /**
0638: * The largest argument for which <code>gamma(x)</code> is representable in the machine.
0639: */
0640: public final static double GAMMA_X_MAX_VALUE = 171.624;
0641:
0642: /**
0643: * Gamma function.
0644: * Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
0645: * Applied Mathematics Division<BR>
0646: * Argonne National Laboratory<BR>
0647: * Argonne, IL 60439<BR>
0648: * <P>
0649: * References:
0650: * <OL>
0651: * <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.
0652: * <LI>Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
0653: * </OL></P><P>
0654: * From the original documentation:
0655: * </P><P>
0656: * This routine calculates the GAMMA function for a real argument X.
0657: * Computation is based on an algorithm outlined in reference 1.
0658: * The program uses rational functions that approximate the GAMMA
0659: * function to at least 20 significant decimal digits. Coefficients
0660: * for the approximation over the interval (1,2) are unpublished.
0661: * Those for the approximation for X .GE. 12 are from reference 2.
0662: * The accuracy achieved depends on the arithmetic system, the
0663: * compiler, the intrinsic functions, and proper selection of the
0664: * machine-dependent constants.
0665: * </P><P>
0666: * Error returns:<BR>
0667: * The program returns the value XINF for singularities or when overflow would occur.
0668: * The computation is believed to be free of underflow and overflow.
0669: * </P>
0670: * @return Double.MAX_VALUE if overflow would occur, i.e. if abs(x) > 171.624
0671: * @jsci.planetmath GammaFunction
0672: * @author Jaco van Kooten
0673: */
0674: public static double gamma(double x) {
0675: double fact = 1.0, xden, xnum;
0676: int i, n = 0;
0677: double y = x, z, y1;
0678: boolean parity = false;
0679: double res, sum, ysq;
0680:
0681: if (y <= 0.0) {
0682: // ----------------------------------------------------------------------
0683: // Argument is negative
0684: // ----------------------------------------------------------------------
0685: y = -(x);
0686: y1 = (int) y;
0687: res = y - y1;
0688: if (res != 0.0) {
0689: if (y1 != (((int) (y1 * 0.5)) * 2.0))
0690: parity = true;
0691: fact = -Math.PI / Math.sin(Math.PI * res);
0692: y++;
0693: } else
0694: return Double.MAX_VALUE;
0695: }
0696: // ----------------------------------------------------------------------
0697: // Argument is positive
0698: // ----------------------------------------------------------------------
0699: if (y < EPS) {
0700: // ----------------------------------------------------------------------
0701: // Argument .LT. EPS
0702: // ----------------------------------------------------------------------
0703: if (y >= XMININ)
0704: res = 1.0 / y;
0705: else
0706: return Double.MAX_VALUE;
0707: } else if (y < 12.0) {
0708: y1 = y;
0709: if (y < 1.0) {
0710: // ----------------------------------------------------------------------
0711: // 0.0 .LT. argument .LT. 1.0
0712: // ----------------------------------------------------------------------
0713: z = y;
0714: y++;
0715: } else {
0716: // ----------------------------------------------------------------------
0717: // 1.0 .LT. argument .LT. 12.0, reduce argument if necessary
0718: // ----------------------------------------------------------------------
0719: n = (int) y - 1;
0720: y -= (double) n;
0721: z = y - 1.0;
0722: }
0723: // ----------------------------------------------------------------------
0724: // Evaluate approximation for 1.0 .LT. argument .LT. 2.0
0725: // ----------------------------------------------------------------------
0726: xnum = 0.0;
0727: xden = 1.0;
0728: for (i = 0; i < 8; ++i) {
0729: xnum = (xnum + g_p[i]) * z;
0730: xden = xden * z + g_q[i];
0731: }
0732: res = xnum / xden + 1.0;
0733: if (y1 < y)
0734: // ----------------------------------------------------------------------
0735: // Adjust result for case 0.0 .LT. argument .LT. 1.0
0736: // ----------------------------------------------------------------------
0737: res /= y1;
0738: else if (y1 > y) {
0739: // ----------------------------------------------------------------------
0740: // Adjust result for case 2.0 .LT. argument .LT. 12.0
0741: // ----------------------------------------------------------------------
0742: for (i = 0; i < n; ++i) {
0743: res *= y;
0744: y++;
0745: }
0746: }
0747: } else {
0748: // ----------------------------------------------------------------------
0749: // Evaluate for argument .GE. 12.0
0750: // ----------------------------------------------------------------------
0751: if (y <= GAMMA_X_MAX_VALUE) {
0752: ysq = y * y;
0753: sum = g_c[6];
0754: for (i = 0; i < 6; ++i)
0755: sum = sum / ysq + g_c[i];
0756: sum = sum / y - y + LOGSQRT2PI;
0757: sum += (y - 0.5) * Math.log(y);
0758: res = Math.exp(sum);
0759: } else
0760: return Double.MAX_VALUE;
0761: }
0762: // ----------------------------------------------------------------------
0763: // Final adjustments and return
0764: // ----------------------------------------------------------------------
0765: if (parity)
0766: res = -res;
0767: if (fact != 1.0)
0768: res = fact / res;
0769: return res;
0770: }
0771:
0772: /**
0773: * The largest argument for which <code>logGamma(x)</code> is representable in the machine.
0774: */
0775: public final static double LOG_GAMMA_X_MAX_VALUE = 2.55e305;
0776:
0777: // Log Gamma related constants
0778: private final static double lg_d1 = -0.5772156649015328605195174;
0779: private final static double lg_d2 = 0.4227843350984671393993777;
0780: private final static double lg_d4 = 1.791759469228055000094023;
0781: private final static double lg_p1[] = { 4.945235359296727046734888,
0782: 201.8112620856775083915565, 2290.838373831346393026739,
0783: 11319.67205903380828685045, 28557.24635671635335736389,
0784: 38484.96228443793359990269, 26377.48787624195437963534,
0785: 7225.813979700288197698961 };
0786: private final static double lg_q1[] = { 67.48212550303777196073036,
0787: 1113.332393857199323513008, 7738.757056935398733233834,
0788: 27639.87074403340708898585, 54993.10206226157329794414,
0789: 61611.22180066002127833352, 36351.27591501940507276287,
0790: 8785.536302431013170870835 };
0791: private final static double lg_p2[] = { 4.974607845568932035012064,
0792: 542.4138599891070494101986, 15506.93864978364947665077,
0793: 184793.2904445632425417223, 1088204.76946882876749847,
0794: 3338152.967987029735917223, 5106661.678927352456275255,
0795: 3074109.054850539556250927 };
0796: private final static double lg_q2[] = { 183.0328399370592604055942,
0797: 7765.049321445005871323047, 133190.3827966074194402448,
0798: 1136705.821321969608938755, 5267964.117437946917577538,
0799: 13467014.54311101692290052, 17827365.30353274213975932,
0800: 9533095.591844353613395747 };
0801: private final static double lg_p4[] = { 14745.02166059939948905062,
0802: 2426813.369486704502836312, 121475557.4045093227939592,
0803: 2663432449.630976949898078, 29403789566.34553899906876,
0804: 170266573776.5398868392998, 492612579337.743088758812,
0805: 560625185622.3951465078242 };
0806: private final static double lg_q4[] = { 2690.530175870899333379843,
0807: 639388.5654300092398984238, 41355999.30241388052042842,
0808: 1120872109.61614794137657, 14886137286.78813811542398,
0809: 101680358627.2438228077304, 341747634550.7377132798597,
0810: 446315818741.9713286462081 };
0811: private final static double lg_c[] = { -0.001910444077728,
0812: 8.4171387781295e-4, -5.952379913043012e-4,
0813: 7.93650793500350248e-4, -0.002777777777777681622553,
0814: 0.08333333333333333331554247, 0.0057083835261 };
0815: // Rough estimate of the fourth root of logGamma_xBig
0816: private final static double lg_frtbig = 2.25e76;
0817: private final static double pnt68 = 0.6796875;
0818:
0819: // Function cache for logGamma
0820: private static final ThreadLocal logGammaCache_res = new ThreadLocal() {
0821: protected Object initialValue() {
0822: return new Double(0.0);
0823: }
0824:
0825: };
0826: private static final ThreadLocal logGammaCache_x = new ThreadLocal() {
0827: protected Object initialValue() {
0828: return new Double(0.0);
0829: }
0830: };
0831:
0832: /**
0833: * The natural logarithm of the gamma function.
0834: * Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
0835: * Applied Mathematics Division<BR>
0836: * Argonne National Laboratory<BR>
0837: * Argonne, IL 60439<BR>
0838: * <P>
0839: * References:
0840: * <OL>
0841: * <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.
0842: * <LI>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.
0843: * <LI>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.
0844: * </OL></P><P>
0845: * From the original documentation:
0846: * </P><P>
0847: * This routine calculates the LOG(GAMMA) function for a positive real argument X.
0848: * Computation is based on an algorithm outlined in references 1 and 2.
0849: * The program uses rational functions that theoretically approximate LOG(GAMMA)
0850: * to at least 18 significant decimal digits. The approximation for X > 12 is from reference 3,
0851: * while approximations for X < 12.0 are similar to those in reference 1, but are unpublished.
0852: * The accuracy achieved depends on the arithmetic system, the compiler, the intrinsic functions,
0853: * and proper selection of the machine-dependent constants.
0854: * </P><P>
0855: * Error returns:<BR>
0856: * The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
0857: * The computation is believed to be free of underflow and overflow.
0858: * </P>
0859: * @return Double.MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
0860: * @author Jaco van Kooten
0861: */
0862: public static double logGamma(double x) {
0863: double xden, corr, xnum;
0864: int i;
0865: double y, xm1, xm2, xm4, res, ysq;
0866:
0867: if (x == ((Double) logGammaCache_x.get()).doubleValue())
0868: return ((Double) logGammaCache_res.get()).doubleValue();
0869:
0870: y = x;
0871: if (y > 0.0 && y <= LOG_GAMMA_X_MAX_VALUE) {
0872: if (y <= EPS) {
0873: res = -Math.log(y);
0874: } else if (y <= 1.5) {
0875: // ----------------------------------------------------------------------
0876: // EPS .LT. X .LE. 1.5
0877: // ----------------------------------------------------------------------
0878: if (y < pnt68) {
0879: corr = -Math.log(y);
0880: xm1 = y;
0881: } else {
0882: corr = 0.0;
0883: xm1 = y - 1.0;
0884: }
0885: if (y <= 0.5 || y >= pnt68) {
0886: xden = 1.0;
0887: xnum = 0.0;
0888: for (i = 0; i < 8; i++) {
0889: xnum = xnum * xm1 + lg_p1[i];
0890: xden = xden * xm1 + lg_q1[i];
0891: }
0892: res = corr + xm1 * (lg_d1 + xm1 * (xnum / xden));
0893: } else {
0894: xm2 = y - 1.0;
0895: xden = 1.0;
0896: xnum = 0.0;
0897: for (i = 0; i < 8; i++) {
0898: xnum = xnum * xm2 + lg_p2[i];
0899: xden = xden * xm2 + lg_q2[i];
0900: }
0901: res = corr + xm2 * (lg_d2 + xm2 * (xnum / xden));
0902: }
0903: } else if (y <= 4.0) {
0904: // ----------------------------------------------------------------------
0905: // 1.5 .LT. X .LE. 4.0
0906: // ----------------------------------------------------------------------
0907: xm2 = y - 2.0;
0908: xden = 1.0;
0909: xnum = 0.0;
0910: for (i = 0; i < 8; i++) {
0911: xnum = xnum * xm2 + lg_p2[i];
0912: xden = xden * xm2 + lg_q2[i];
0913: }
0914: res = xm2 * (lg_d2 + xm2 * (xnum / xden));
0915: } else if (y <= 12.0) {
0916: // ----------------------------------------------------------------------
0917: // 4.0 .LT. X .LE. 12.0
0918: // ----------------------------------------------------------------------
0919: xm4 = y - 4.0;
0920: xden = -1.0;
0921: xnum = 0.0;
0922: for (i = 0; i < 8; i++) {
0923: xnum = xnum * xm4 + lg_p4[i];
0924: xden = xden * xm4 + lg_q4[i];
0925: }
0926: res = lg_d4 + xm4 * (xnum / xden);
0927: } else {
0928: // ----------------------------------------------------------------------
0929: // Evaluate for argument .GE. 12.0
0930: // ----------------------------------------------------------------------
0931: res = 0.0;
0932: if (y <= lg_frtbig) {
0933: res = lg_c[6];
0934: ysq = y * y;
0935: for (i = 0; i < 6; i++)
0936: res = res / ysq + lg_c[i];
0937: }
0938: res /= y;
0939: corr = Math.log(y);
0940: res = res + LOGSQRT2PI - 0.5 * corr;
0941: res += y * (corr - 1.0);
0942: }
0943: } else {
0944: // ----------------------------------------------------------------------
0945: // Return for bad arguments
0946: // ----------------------------------------------------------------------
0947: res = Double.MAX_VALUE;
0948: }
0949: // ----------------------------------------------------------------------
0950: // Final adjustments and return
0951: // ----------------------------------------------------------------------
0952: logGammaCache_x.set(new Double(x));
0953: logGammaCache_res.set(new Double(res));
0954: return res;
0955: }
0956:
0957: private final static int MAX_ITERATIONS = 1000;
0958: // lower value = higher precision
0959: private final static double PRECISION = 4.0 * EPS;
0960:
0961: /**
0962: * Incomplete gamma function.
0963: * The computation is based on approximations presented in Numerical Recipes, Chapter 6.2 (W.H. Press et al, 1992).
0964: * @param a require a>=0
0965: * @param x require x>=0
0966: * @return 0 if x<0, a<=0 or a>2.55E305 to avoid errors and over/underflow
0967: * @author Jaco van Kooten
0968: */
0969: public static double incompleteGamma(double a, double x) {
0970: if (x <= 0.0 || a <= 0.0 || a > LOG_GAMMA_X_MAX_VALUE)
0971: return 0.0;
0972: if (x < (a + 1.0))
0973: return gammaSeriesExpansion(a, x);
0974: else
0975: return 1.0 - gammaFraction(a, x);
0976: }
0977:
0978: /**
0979: * @author Jaco van Kooten
0980: */
0981: private static double gammaSeriesExpansion(double a, double x) {
0982: double ap = a;
0983: double del = 1.0 / a;
0984: double sum = del;
0985: for (int n = 1; n < MAX_ITERATIONS; n++) {
0986: ++ap;
0987: del *= x / ap;
0988: sum += del;
0989: if (del < sum * PRECISION)
0990: return sum
0991: * Math.exp(-x + a * Math.log(x) - logGamma(a));
0992: }
0993: throw new RuntimeException(
0994: "Maximum iterations exceeded: please file a bug report.");
0995: }
0996:
0997: /**
0998: * @author Jaco van Kooten
0999: */
1000: private static double gammaFraction(double a, double x) {
1001: double b = x + 1.0 - a;
1002: double c = 1.0 / XMININ;
1003: double d = 1.0 / b;
1004: double h = d;
1005: double del = 0.0;
1006: double an;
1007: for (int i = 1; i < MAX_ITERATIONS
1008: && Math.abs(del - 1.0) > PRECISION; i++) {
1009: an = -i * (i - a);
1010: b += 2.0;
1011: d = an * d + b;
1012: c = b + an / c;
1013: if (Math.abs(c) < XMININ)
1014: c = XMININ;
1015: if (Math.abs(d) < XMININ)
1016: c = XMININ;
1017: d = 1.0 / d;
1018: del = d * c;
1019: h *= del;
1020: }
1021: return Math.exp(-x + a * Math.log(x) - logGamma(a)) * h;
1022: }
1023:
1024: /**
1025: * Beta function.
1026: * @param p require p>0
1027: * @param q require q>0
1028: * @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
1029: * @author Jaco van Kooten
1030: */
1031: public static double beta(double p, double q) {
1032: if (p <= 0.0 || q <= 0.0 || (p + q) > LOG_GAMMA_X_MAX_VALUE)
1033: return 0.0;
1034: else
1035: return Math.exp(logBeta(p, q));
1036: }
1037:
1038: // Function cache for logBeta
1039: private static final ThreadLocal logBetaCache_res = new ThreadLocal() {
1040: protected Object initialValue() {
1041: return new Double(0.0);
1042: }
1043: };
1044: private static final ThreadLocal logBetaCache_p = new ThreadLocal() {
1045: protected Object initialValue() {
1046: return new Double(0.0);
1047: }
1048: };
1049: private static final ThreadLocal logBetaCache_q = new ThreadLocal() {
1050: protected Object initialValue() {
1051: return new Double(0.0);
1052: }
1053: };
1054:
1055: /**
1056: * The natural logarithm of the beta function.
1057: * @param p require p>0
1058: * @param q require q>0
1059: * @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
1060: * @author Jaco van Kooten
1061: */
1062: public static double logBeta(double p, double q) {
1063: if (p != ((Double) logBetaCache_p.get()).doubleValue()
1064: || q != ((Double) logBetaCache_q.get()).doubleValue()) {
1065: logBetaCache_p.set(new Double(p));
1066: logBetaCache_q.set(new Double(q));
1067: double res;
1068: if (p <= 0.0 || q <= 0.0 || (p + q) > LOG_GAMMA_X_MAX_VALUE)
1069: res = 0.0;
1070: else
1071: res = logGamma(p) + logGamma(q) - logGamma(p + q);
1072: logBetaCache_res.set(new Double(res));
1073: return res;
1074: } else {
1075: return ((Double) logBetaCache_res.get()).doubleValue();
1076: }
1077: }
1078:
1079: /**
1080: * Incomplete beta function.
1081: * The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
1082: * @param x require 0<=x<=1
1083: * @param p require p>0
1084: * @param q require q>0
1085: * @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
1086: * @author Jaco van Kooten
1087: */
1088: public static double incompleteBeta(double x, double p, double q) {
1089: if (x <= 0.0)
1090: return 0.0;
1091: else if (x >= 1.0)
1092: return 1.0;
1093: else if (p <= 0.0 || q <= 0.0
1094: || (p + q) > LOG_GAMMA_X_MAX_VALUE)
1095: return 0.0;
1096: else {
1097: final double beta_gam = Math.exp(-logBeta(p, q) + p
1098: * Math.log(x) + q * Math.log(1.0 - x));
1099: if (x < (p + 1.0) / (p + q + 2.0))
1100: return beta_gam * betaFraction(x, p, q) / p;
1101: else
1102: return 1.0 - (beta_gam * betaFraction(1.0 - x, q, p) / q);
1103: }
1104: }
1105:
1106: /**
1107: * Evaluates of continued fraction part of incomplete beta function.
1108: * Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
1109: * @author Jaco van Kooten
1110: */
1111: private static double betaFraction(double x, double p, double q) {
1112: int m, m2;
1113: double sum_pq, p_plus, p_minus, c = 1.0, d, delta, h, frac;
1114: sum_pq = p + q;
1115: p_plus = p + 1.0;
1116: p_minus = p - 1.0;
1117: h = 1.0 - sum_pq * x / p_plus;
1118: if (Math.abs(h) < XMININ)
1119: h = XMININ;
1120: h = 1.0 / h;
1121: frac = h;
1122: m = 1;
1123: delta = 0.0;
1124: while (m <= MAX_ITERATIONS && Math.abs(delta - 1.0) > PRECISION) {
1125: m2 = 2 * m;
1126: // even index for d
1127: d = m * (q - m) * x / ((p_minus + m2) * (p + m2));
1128: h = 1.0 + d * h;
1129: if (Math.abs(h) < XMININ)
1130: h = XMININ;
1131: h = 1.0 / h;
1132: c = 1.0 + d / c;
1133: if (Math.abs(c) < XMININ)
1134: c = XMININ;
1135: frac *= h * c;
1136: // odd index for d
1137: d = -(p + m) * (sum_pq + m) * x
1138: / ((p + m2) * (p_plus + m2));
1139: h = 1.0 + d * h;
1140: if (Math.abs(h) < XMININ)
1141: h = XMININ;
1142: h = 1.0 / h;
1143: c = 1.0 + d / c;
1144: if (Math.abs(c) < XMININ)
1145: c = XMININ;
1146: delta = h * c;
1147: frac *= delta;
1148: m++;
1149: }
1150: return frac;
1151: }
1152:
1153: // ====================================================
1154: // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1155: //
1156: // Developed at SunSoft, a Sun Microsystems, Inc. business.
1157: // Permission to use, copy, modify, and distribute this
1158: // software is freely granted, provided that this notice
1159: // is preserved.
1160: // ====================================================
1161: //
1162: // x
1163: // 2 |\
1164: // erf(x) = --------- | exp(-t*t)dt
1165: // sqrt(pi) \|
1166: // 0
1167: //
1168: // erfc(x) = 1-erf(x)
1169: // Note that
1170: // erf(-x) = -erf(x)
1171: // erfc(-x) = 2 - erfc(x)
1172: //
1173: // Method:
1174: // 1. For |x| in [0, 0.84375]
1175: // erf(x) = x + x*R(x^2)
1176: // erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
1177: // = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
1178: // where R = P/Q where P is an odd poly of degree 8 and
1179: // Q is an odd poly of degree 10.
1180: // -57.90
1181: // | R - (erf(x)-x)/x | <= 2
1182: //
1183: //
1184: // Remark. The formula is derived by noting
1185: // erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
1186: // and that
1187: // 2/sqrt(pi) = 1.128379167095512573896158903121545171688
1188: // is close to one. The interval is chosen because the fix
1189: // point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
1190: // near 0.6174), and by some experiment, 0.84375 is chosen to
1191: // guarantee the error is less than one ulp for erf.
1192: //
1193: // 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
1194: // c = 0.84506291151 rounded to single (24 bits)
1195: // erf(x) = sign(x) * (c + P1(s)/Q1(s))
1196: // erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
1197: // 1+(c+P1(s)/Q1(s)) if x < 0
1198: // |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
1199: // Remark: here we use the taylor series expansion at x=1.
1200: // erf(1+s) = erf(1) + s*Poly(s)
1201: // = 0.845.. + P1(s)/Q1(s)
1202: // That is, we use rational approximation to approximate
1203: // erf(1+s) - (c = (single)0.84506291151)
1204: // Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
1205: // where
1206: // P1(s) = degree 6 poly in s
1207: // Q1(s) = degree 6 poly in s
1208: //
1209: // 3. For x in [1.25,1/0.35(~2.857143)],
1210: // erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
1211: // erf(x) = 1 - erfc(x)
1212: // where
1213: // R1(z) = degree 7 poly in z, (z=1/x^2)
1214: // S1(z) = degree 8 poly in z
1215: //
1216: // 4. For x in [1/0.35,28]
1217: // erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
1218: // = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
1219: // = 2.0 - tiny (if x <= -6)
1220: // erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
1221: // erf(x) = sign(x)*(1.0 - tiny)
1222: // where
1223: // R2(z) = degree 6 poly in z, (z=1/x^2)
1224: // S2(z) = degree 7 poly in z
1225: //
1226: // Note1:
1227: // To compute exp(-x*x-0.5625+R/S), let s be a single
1228: // precision number and s := x; then
1229: // -x*x = -s*s + (s-x)*(s+x)
1230: // exp(-x*x-0.5626+R/S) =
1231: // exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
1232: // Note2:
1233: // Here 4 and 5 make use of the asymptotic series
1234: // exp(-x*x)
1235: // erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
1236: // x*sqrt(pi)
1237: // We use rational approximation to approximate
1238: // g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
1239: // Here is the error bound for R1/S1 and R2/S2
1240: // |R1/S1 - f(x)| < 2**(-62.57)
1241: // |R2/S2 - f(x)| < 2**(-61.52)
1242: //
1243: // 5. For inf > x >= 28
1244: // erf(x) = sign(x) *(1 - tiny) (raise inexact)
1245: // erfc(x) = tiny*tiny (raise underflow) if x > 0
1246: // = 2 - tiny if x<0
1247: //
1248: // 7. Special case:
1249: // erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
1250: // erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
1251: // erfc/erf(NaN) is NaN
1252: //
1253:
1254: // Coefficients for approximation to erf on [0,0.84375]
1255: private final static double e_efx = 1.28379167095512586316e-01;
1256: // private final static double efx8=1.02703333676410069053e00;
1257: private final static double ePp[] = { 1.28379167095512558561e-01,
1258: -3.25042107247001499370e-01, -2.84817495755985104766e-02,
1259: -5.77027029648944159157e-03, -2.37630166566501626084e-05 };
1260: private final static double eQq[] = { 3.97917223959155352819e-01,
1261: 6.50222499887672944485e-02, 5.08130628187576562776e-03,
1262: 1.32494738004321644526e-04, -3.96022827877536812320e-06 };
1263: // Coefficients for approximation to erf in [0.84375,1.25]
1264: private final static double ePa[] = { -2.36211856075265944077e-03,
1265: 4.14856118683748331666e-01, -3.72207876035701323847e-01,
1266: 3.18346619901161753674e-01, -1.10894694282396677476e-01,
1267: 3.54783043256182359371e-02, -2.16637559486879084300e-03 };
1268: private final static double eQa[] = { 1.06420880400844228286e-01,
1269: 5.40397917702171048937e-01, 7.18286544141962662868e-02,
1270: 1.26171219808761642112e-01, 1.36370839120290507362e-02,
1271: 1.19844998467991074170e-02 };
1272: private final static double e_erx = 8.45062911510467529297e-01;
1273:
1274: /**
1275: * Error function.
1276: * Based on C-code for the error function developed at Sun Microsystems.
1277: * @author Jaco van Kooten
1278: */
1279: public static double error(double x) {
1280: double P, Q, s, retval;
1281: final double abs_x = (x >= 0.0 ? x : -x);
1282: if (abs_x < 0.84375) { // 0 < |x| < 0.84375
1283: if (abs_x < 3.7252902984619141e-9) // |x| < 2**-28
1284: retval = abs_x + abs_x * e_efx;
1285: else {
1286: s = x * x;
1287: P = ePp[0]
1288: + s
1289: * (ePp[1] + s
1290: * (ePp[2] + s * (ePp[3] + s * ePp[4])));
1291: Q = 1.0
1292: + s
1293: * (eQq[0] + s
1294: * (eQq[1] + s
1295: * (eQq[2] + s
1296: * (eQq[3] + s * eQq[4]))));
1297: retval = abs_x + abs_x * (P / Q);
1298: }
1299: } else if (abs_x < 1.25) { // 0.84375 < |x| < 1.25
1300: s = abs_x - 1.0;
1301: P = ePa[0]
1302: + s
1303: * (ePa[1] + s
1304: * (ePa[2] + s
1305: * (ePa[3] + s
1306: * (ePa[4] + s
1307: * (ePa[5] + s
1308: * ePa[6])))));
1309: Q = 1.0
1310: + s
1311: * (eQa[0] + s
1312: * (eQa[1] + s
1313: * (eQa[2] + s
1314: * (eQa[3] + s
1315: * (eQa[4] + s
1316: * eQa[5])))));
1317: retval = e_erx + P / Q;
1318: } else if (abs_x >= 6.0)
1319: retval = 1.0;
1320: else
1321: // 1.25 < |x| < 6.0
1322: retval = 1.0 - complementaryError(abs_x);
1323: return (x >= 0.0) ? retval : -retval;
1324: }
1325:
1326: // Coefficients for approximation to erfc in [1.25,1/.35]
1327: private final static double eRa[] = { -9.86494403484714822705e-03,
1328: -6.93858572707181764372e-01, -1.05586262253232909814e01,
1329: -6.23753324503260060396e01, -1.62396669462573470355e02,
1330: -1.84605092906711035994e02, -8.12874355063065934246e01,
1331: -9.81432934416914548592e00 };
1332: private final static double eSa[] = { 1.96512716674392571292e01,
1333: 1.37657754143519042600e02, 4.34565877475229228821e02,
1334: 6.45387271733267880336e02, 4.29008140027567833386e02,
1335: 1.08635005541779435134e02, 6.57024977031928170135e00,
1336: -6.04244152148580987438e-02 };
1337: // Coefficients for approximation to erfc in [1/.35,28]
1338: private final static double eRb[] = { -9.86494292470009928597e-03,
1339: -7.99283237680523006574e-01, -1.77579549177547519889e01,
1340: -1.60636384855821916062e02, -6.37566443368389627722e02,
1341: -1.02509513161107724954e03, -4.83519191608651397019e02 };
1342: private final static double eSb[] = { 3.03380607434824582924e01,
1343: 3.25792512996573918826e02, 1.53672958608443695994e03,
1344: 3.19985821950859553908e03, 2.55305040643316442583e03,
1345: 4.74528541206955367215e02, -2.24409524465858183362e01 };
1346:
1347: /**
1348: * Complementary error function.
1349: * Based on C-code for the error function developed at Sun Microsystems.
1350: * @author Jaco van Kooten
1351: */
1352: public static double complementaryError(double x) {
1353: double s, retval, R, S;
1354: final double abs_x = (x >= 0.0 ? x : -x);
1355: if (abs_x < 1.25)
1356: retval = 1.0 - error(abs_x);
1357: else if (abs_x > 28.0)
1358: retval = 0.0;
1359: else { // 1.25 < |x| < 28
1360: s = 1.0 / (abs_x * abs_x);
1361: if (abs_x < 2.8571428) { // ( |x| < 1/0.35 )
1362: R = eRa[0]
1363: + s
1364: * (eRa[1] + s
1365: * (eRa[2] + s
1366: * (eRa[3] + s
1367: * (eRa[4] + s
1368: * (eRa[5] + s
1369: * (eRa[6] + s
1370: * eRa[7]))))));
1371: S = 1.0
1372: + s
1373: * (eSa[0] + s
1374: * (eSa[1] + s
1375: * (eSa[2] + s
1376: * (eSa[3] + s
1377: * (eSa[4] + s
1378: * (eSa[5] + s
1379: * (eSa[6] + s
1380: * eSa[7])))))));
1381: } else { // ( |x| > 1/0.35 )
1382: R = eRb[0]
1383: + s
1384: * (eRb[1] + s
1385: * (eRb[2] + s
1386: * (eRb[3] + s
1387: * (eRb[4] + s
1388: * (eRb[5] + s
1389: * eRb[6])))));
1390: S = 1.0
1391: + s
1392: * (eSb[0] + s
1393: * (eSb[1] + s
1394: * (eSb[2] + s
1395: * (eSb[3] + s
1396: * (eSb[4] + s
1397: * (eSb[5] + s
1398: * eSb[6]))))));
1399: }
1400: retval = Math.exp(-x * x - 0.5625 + R / S) / abs_x;
1401: }
1402: return (x >= 0.0) ? retval : 2.0 - retval;
1403: }
1404: }
|