KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > JSci > maths > SpecialMath


1 package JSci.maths;
2
3 /**
4 * The special function math library.
5 * This class cannot be subclassed or instantiated because all methods are static.
6 * @version 1.1
7 * @author Mark Hale
8 */
9 public final class SpecialMath extends AbstractMath implements NumericalConstants {
10         private SpecialMath() {}
11
12 // Some IEEE machine constants
13 /**
14         * Relative machine precision.
15         */
16         private final static double EPS=2.22e-16;
17         /**
18         * The smallest positive floating-point number such that 1/xminin is machine representable.
19         */
20         private final static double XMININ=2.23e-308;
21
22
23 // CHEBYSHEV SERIES
24
25 // series for ai0 on the interval 1.25000d-01 to 3.33333d-01
26 // with weighted error 7.87e-17
27 // log weighted error 16.10
28 // significant figures required 14.69
29 // decimal places required 16.76
30
31         private final static double ai0cs[]={
32                 0.07575994494023796,
33                 0.00759138081082334,
34                 0.00041531313389237,
35                 0.00001070076463439,
36                 -0.00000790117997921,
37                 -0.00000078261435014,
38                 0.00000027838499429,
39                 0.00000000825247260,
40                 -0.00000001204463945,
41                 0.00000000155964859,
42                 0.00000000022925563,
43                 -0.00000000011916228,
44                 0.00000000001757854,
45                 0.00000000000112822,
46                 -0.00000000000114684,
47                 0.00000000000027155,
48                 -0.00000000000002415,
49                 -0.00000000000000608,
50                 0.00000000000000314,
51                 -0.00000000000000071,
52                 0.00000000000000007};
53
54 // series for ai02 on the interval 0. to 1.25000d-01
55 // with weighted error 3.79e-17
56 // log weighted error 16.42
57 // significant figures required 14.86
58 // decimal places required 17.09
59 private final static double ai02cs[]={
60                 0.05449041101410882,
61                 0.00336911647825569,
62                 0.00006889758346918,
63                 0.00000289137052082,
64                 0.00000020489185893,
65                 0.00000002266668991,
66                 0.00000000339623203,
67                 0.00000000049406022,
68                 0.00000000001188914,
69                 -0.00000000003149915,
70                 -0.00000000001321580,
71                 -0.00000000000179419,
72                 0.00000000000071801,
73                 0.00000000000038529,
74                 0.00000000000001539,
75                 -0.00000000000004151,
76                 -0.00000000000000954,
77                 0.00000000000000382,
78                 0.00000000000000176,
79                 -0.00000000000000034,
80                 -0.00000000000000027,
81                 0.00000000000000003};
82
83 // series for ai1 on the interval 1.25000d-01 to 3.33333d-01
84 // with weighted error 6.98e-17
85 // log weighted error 16.16
86 // significant figures required 14.53
87 // decimal places required 16.82
88
89         private final static double ai1cs[]={
90                 -0.02846744181881479,
91                 -0.01922953231443221,
92                 -0.00061151858579437,
93                 -0.00002069971253350,
94                 0.00000858561914581,
95                 0.00000104949824671,
96                 -0.00000029183389184,
97                 -0.00000001559378146,
98                 0.00000001318012367,
99                 -0.00000000144842341,
100                 -0.00000000029085122,
101                 0.00000000012663889,
102                 -0.00000000001664947,
103                 -0.00000000000166665,
104                 0.00000000000124260,
105                 -0.00000000000027315,
106                 0.00000000000002023,
107                 0.00000000000000730,
108                 -0.00000000000000333,
109                 0.00000000000000071,
110                 -0.00000000000000006};
111
112 // series for ai12 on the interval 0. to 1.25000d-01
113 // with weighted error 3.55e-17
114 // log weighted error 16.45
115 // significant figures required 14.69
116 // decimal places required 17.12
117
118         private final static double ai12cs[]={
119                 0.02857623501828014,
120                 -0.00976109749136147,
121                 -0.00011058893876263,
122                 -0.00000388256480887,
123                 -0.00000025122362377,
124                 -0.00000002631468847,
125                 -0.00000000383538039,
126                 -0.00000000055897433,
127                 -0.00000000001897495,
128                 0.00000000003252602,
129                 0.00000000001412580,
130                 0.00000000000203564,
131                 -0.00000000000071985,
132                 -0.00000000000040836,
133                 -0.00000000000002101,
134                 0.00000000000004273,
135                 0.00000000000001041,
136                 -0.00000000000000382,
137                 -0.00000000000000186,
138                 0.00000000000000033,
139                 0.00000000000000028,
140                 -0.00000000000000003};
141
142 // series for aif on the interval -1.00000d+00 to 1.00000d+00
143 // with weighted error 1.09e-19
144 // log weighted error 18.96
145 // significant figures required 17.76
146 // decimal places required 19.44
147
148         private final static double aifcs[]={
149                 -0.03797135849666999750,
150                 0.05919188853726363857,
151                 0.00098629280577279975,
152                 0.00000684884381907656,
153                 0.00000002594202596219,
154                 0.00000000006176612774,
155                 0.00000000000010092454,
156                 0.00000000000000012014,
157                 0.00000000000000000010};
158
159 // series for aig on the interval -1.00000d+00 to 1.00000d+00
160 // with weighted error 1.51e-17
161 // log weighted error 16.82
162 // significant figures required 15.19
163 // decimal places required 17.27
164
165         private final static double aigcs[]={
166                 0.01815236558116127,
167                 0.02157256316601076,
168                 0.00025678356987483,
169                 0.00000142652141197,
170                 0.00000000457211492,
171                 0.00000000000952517,
172                 0.00000000000001392,
173                 0.00000000000000001};
174
175 // series for aip on the interval 0. to 1.00000d+00
176 // with weighted error 5.10e-17
177 // log weighted error 16.29
178 // significant figures required 14.41
179 // decimal places required 17.06
180
181         private final static double aipcs[]={
182                 -0.0187519297793868,
183                 -0.0091443848250055,
184                 0.0009010457337825,
185                 -0.0001394184127221,
186                 0.0000273815815785,
187                 -0.0000062750421119,
188                 0.0000016064844184,
189                 -0.0000004476392158,
190                 0.0000001334635874,
191                 -0.0000000420735334,
192                 0.0000000139021990,
193                 -0.0000000047831848,
194                 0.0000000017047897,
195                 -0.0000000006268389,
196                 0.0000000002369824,
197                 -0.0000000000918641,
198                 0.0000000000364278,
199                 -0.0000000000147475,
200                 0.0000000000060851,
201                 -0.0000000000025552,
202                 0.0000000000010906,
203                 -0.0000000000004725,
204                 0.0000000000002076,
205                 -0.0000000000000924,
206                 0.0000000000000417,
207                 -0.0000000000000190,
208                 0.0000000000000087,
209                 -0.0000000000000040,
210                 0.0000000000000019,
211                 -0.0000000000000009,
212                 0.0000000000000004,
213                 -0.0000000000000002,
214                 0.0000000000000001,
215                 -0.0000000000000000};
216
217 // series for am21 on the interval -1.25000d-01 to 0.
218 // with weighted error 2.89e-17
219 // log weighted error 16.54
220 // significant figures required 14.15
221 // decimal places required 17.34
222
223         private final static double am21cs[]={
224                 0.0065809191761485,
225                 0.0023675984685722,
226                 0.0001324741670371,
227                 0.0000157600904043,
228                 0.0000027529702663,
229                 0.0000006102679017,
230                 0.0000001595088468,
231                 0.0000000471033947,
232                 0.0000000152933871,
233                 0.0000000053590722,
234                 0.0000000020000910,
235                 0.0000000007872292,
236                 0.0000000003243103,
237                 0.0000000001390106,
238                 0.0000000000617011,
239                 0.0000000000282491,
240                 0.0000000000132979,
241                 0.0000000000064188,
242                 0.0000000000031697,
243                 0.0000000000015981,
244                 0.0000000000008213,
245                 0.0000000000004296,
246                 0.0000000000002284,
247                 0.0000000000001232,
248                 0.0000000000000675,
249                 0.0000000000000374,
250                 0.0000000000000210,
251                 0.0000000000000119,
252                 0.0000000000000068,
253                 0.0000000000000039,
254                 0.0000000000000023,
255                 0.0000000000000013,
256                 0.0000000000000008,
257                 0.0000000000000005,
258                 0.0000000000000003,
259                 0.0000000000000001,
260                 0.0000000000000001,
261                 0.0000000000000000,
262                 0.0000000000000000,
263                 0.0000000000000000};
264
265 // series for ath1 on the interval -1.25000d-01 to 0.
266 // with weighted error 2.53e-17
267 // log weighted error 16.60
268 // significant figures required 15.15
269 // decimal places required 17.38
270
271         private final static double ath1cs[]={
272                 -0.07125837815669365,
273                 -0.00590471979831451,
274                 -0.00012114544069499,
275                 -0.00000988608542270,
276                 -0.00000138084097352,
277                 -0.00000026142640172,
278                 -0.00000006050432589,
279                 -0.00000001618436223,
280                 -0.00000000483464911,
281                 -0.00000000157655272,
282                 -0.00000000055231518,
283                 -0.00000000020545441,
284                 -0.00000000008043412,
285                 -0.00000000003291252,
286                 -0.00000000001399875,
287                 -0.00000000000616151,
288                 -0.00000000000279614,
289                 -0.00000000000130428,
290                 -0.00000000000062373,
291                 -0.00000000000030512,
292                 -0.00000000000015239,
293                 -0.00000000000007758,
294                 -0.00000000000004020,
295                 -0.00000000000002117,
296                 -0.00000000000001132,
297                 -0.00000000000000614,
298                 -0.00000000000000337,
299                 -0.00000000000000188,
300                 -0.00000000000000105,
301                 -0.00000000000000060,
302                 -0.00000000000000034,
303                 -0.00000000000000020,
304                 -0.00000000000000011,
305                 -0.00000000000000007,
306                 -0.00000000000000004,
307                 -0.00000000000000002};
308
309 // series for am22 on the interval -1.00000d+00 to -1.25000d-01
310 // with weighted error 2.99e-17
311 // log weighted error 16.52
312 // significant figures required 14.57
313 // decimal places required 17.28
314
315         private final static double am22cs[]={
316                 -0.01562844480625341,
317                 0.00778336445239681,
318                 0.00086705777047718,
319                 0.00015696627315611,
320                 0.00003563962571432,
321                 0.00000924598335425,
322                 0.00000262110161850,
323                 0.00000079188221651,
324                 0.00000025104152792,
325                 0.00000008265223206,
326                 0.00000002805711662,
327                 0.00000000976821090,
328                 0.00000000347407923,
329                 0.00000000125828132,
330                 0.00000000046298826,
331                 0.00000000017272825,
332                 0.00000000006523192,
333                 0.00000000002490471,
334                 0.00000000000960156,
335                 0.00000000000373448,
336                 0.00000000000146417,
337                 0.00000000000057826,
338                 0.00000000000022991,
339                 0.00000000000009197,
340                 0.00000000000003700,
341                 0.00000000000001496,
342                 0.00000000000000608,
343                 0.00000000000000248,
344                 0.00000000000000101,
345                 0.00000000000000041,
346                 0.00000000000000017,
347                 0.00000000000000007,
348                 0.00000000000000002};
349
350 // series for ath2 on the interval -1.00000d+00 to -1.25000d-01
351 // with weighted error 2.57e-17
352 // log weighted error 16.59
353 // significant figures required 15.07
354 // decimal places required 17.34
355
356         private final static double ath2cs[]={
357                 0.00440527345871877,
358                 -0.03042919452318455,
359                 -0.00138565328377179,
360                 -0.00018044439089549,
361                 -0.00003380847108327,
362                 -0.00000767818353522,
363                 -0.00000196783944371,
364                 -0.00000054837271158,
365                 -0.00000016254615505,
366                 -0.00000005053049981,
367                 -0.00000001631580701,
368                 -0.00000000543420411,
369                 -0.00000000185739855,
370                 -0.00000000064895120,
371                 -0.00000000023105948,
372                 -0.00000000008363282,
373                 -0.00000000003071196,
374                 -0.00000000001142367,
375                 -0.00000000000429811,
376                 -0.00000000000163389,
377                 -0.00000000000062693,
378                 -0.00000000000024260,
379                 -0.00000000000009461,
380                 -0.00000000000003716,
381                 -0.00000000000001469,
382                 -0.00000000000000584,
383                 -0.00000000000000233,
384                 -0.00000000000000093,
385                 -0.00000000000000037,
386                 -0.00000000000000015,
387                 -0.00000000000000006,
388                 -0.00000000000000002};
389
390 // series for bi0 on the interval 0. to 9.00000d+00
391 // with weighted error 2.46e-18
392 // log weighted error 17.61
393 // significant figures required 17.90
394 // decimal places required 18.15
395
396         private final static double bi0cs[]={
397                 -0.07660547252839144951,
398                 1.927337953993808270,
399                 0.2282644586920301339,
400                 0.01304891466707290428,
401                 0.00043442709008164874,
402                 0.00000942265768600193,
403                 0.00000014340062895106,
404                 0.00000000161384906966,
405                 0.00000000001396650044,
406                 0.00000000000009579451,
407                 0.00000000000000053339,
408                 0.00000000000000000245};
409
410 // series for bj0 on the interval 0. to 1.60000d+01
411 // with weighted error 7.47e-18
412 // log weighted error 17.13
413 // significant figures required 16.98
414 // decimal places required 17.68
415
416         private final static double bj0cs[]={
417                 0.100254161968939137,
418                 -0.665223007764405132,
419                 0.248983703498281314,
420                 -0.0332527231700357697,
421                 0.0023114179304694015,
422                 -0.0000991127741995080,
423                 0.0000028916708643998,
424                 -0.0000000612108586630,
425                 0.0000000009838650793,
426                 -0.0000000000124235515,
427                 0.0000000000001265433,
428                 -0.0000000000000010619,
429                 0.0000000000000000074};
430
431 // series for bm0 on the interval 0. to 6.25000d-02
432 // with weighted error 4.98e-17
433 // log weighted error 16.30
434 // significant figures required 14.97
435 // decimal places required 16.96
436
437         private final static double bm0cs[]={
438                 0.09284961637381644,
439                 -0.00142987707403484,
440                 0.00002830579271257,
441                 -0.00000143300611424,
442                 0.00000012028628046,
443                 -0.00000001397113013,
444                 0.00000000204076188,
445                 -0.00000000035399669,
446                 0.00000000007024759,
447                 -0.00000000001554107,
448                 0.00000000000376226,
449                 -0.00000000000098282,
450                 0.00000000000027408,
451                 -0.00000000000008091,
452                 0.00000000000002511,
453                 -0.00000000000000814,
454                 0.00000000000000275,
455                 -0.00000000000000096,
456                 0.00000000000000034,
457                 -0.00000000000000012,
458                 0.00000000000000004};
459
460 // series for bth0 on the interval 0. to 6.25000d-02
461 // with weighted error 3.67e-17
462 // log weighted error 16.44
463 // significant figures required 15.53
464 // decimal places required 17.13
465
466         private final static double bth0cs[]={
467                 -0.24639163774300119,
468                 0.001737098307508963,
469                 -0.000062183633402968,
470                 0.000004368050165742,
471                 -0.000000456093019869,
472                 0.000000062197400101,
473                 -0.000000010300442889,
474                 0.000000001979526776,
475                 -0.000000000428198396,
476                 0.000000000102035840,
477                 -0.000000000026363898,
478                 0.000000000007297935,
479                 -0.000000000002144188,
480                 0.000000000000663693,
481                 -0.000000000000215126,
482                 0.000000000000072659,
483                 -0.000000000000025465,
484                 0.000000000000009229,
485                 -0.000000000000003448,
486                 0.000000000000001325,
487                 -0.000000000000000522,
488                 0.000000000000000210,
489                 -0.000000000000000087,
490                 0.000000000000000036};
491
492 // series for by0 on the interval 0. to 1.60000d+01
493 // with weighted error 1.20e-17
494 // log weighted error 16.92
495 // significant figures required 16.15
496 // decimal places required 17.48
497
498         private final static double by0cs[]={
499                 -0.011277839392865573,
500                 -0.12834523756042035,
501                 -0.10437884799794249,
502                 0.023662749183969695,
503                 -0.002090391647700486,
504                 0.000103975453939057,
505                 -0.000003369747162423,
506                 0.000000077293842676,
507                 -0.000000001324976772,
508                 0.000000000017648232,
509                 -0.000000000000188105,
510                 0.000000000000001641,
511                 -0.000000000000000011};
512
513 // series for bi1 on the interval 0. to 9.00000d+00
514 // with weighted error 2.40e-17
515 // log weighted error 16.62
516 // significant figures required 16.23
517 // decimal places required 17.14
518
519         private final static double bi1cs[]={
520                 -0.001971713261099859,
521                 0.40734887667546481,
522                 0.034838994299959456,
523                 0.001545394556300123,
524                 0.000041888521098377,
525                 0.000000764902676483,
526                 0.000000010042493924,
527                 0.000000000099322077,
528                 0.000000000000766380,
529                 0.000000000000004741,
530                 0.000000000000000024};
531
532 // series for bj1 on the interval 0. to 1.60000d+01
533 // with weighted error 4.48e-17
534 // log weighted error 16.35
535 // significant figures required 15.77
536 // decimal places required 16.89
537
538         private final static double bj1cs[]={
539                 -0.11726141513332787,
540                 -0.25361521830790640,
541                 0.050127080984469569,
542                 -0.004631514809625081,
543                 0.000247996229415914,
544                 -0.000008678948686278,
545                 0.000000214293917143,
546                 -0.000000003936093079,
547                 0.000000000055911823,
548                 -0.000000000000632761,
549                 0.000000000000005840,
550                 -0.000000000000000044};
551
552 // series for bm1 on the interval 0. to 6.25000d-02
553 // with weighted error 5.61e-17
554 // log weighted error 16.25
555 // significant figures required 14.97
556 // decimal places required 16.91
557
558         private final static double bm1cs[]={
559                 0.1047362510931285,
560                 0.00442443893702345,
561                 -0.00005661639504035,
562                 0.00000231349417339,
563                 -0.00000017377182007,
564                 0.00000001893209930,
565                 -0.00000000265416023,
566                 0.00000000044740209,
567                 -0.00000000008691795,
568                 0.00000000001891492,
569                 -0.00000000000451884,
570                 0.00000000000116765,
571                 -0.00000000000032265,
572                 0.00000000000009450,
573                 -0.00000000000002913,
574                 0.00000000000000939,
575                 -0.00000000000000315,
576                 0.00000000000000109,
577                 -0.00000000000000039,
578                 0.00000000000000014,
579                 -0.00000000000000005};
580
581 // series for bth1 on the interval 0. to 6.25000d-02
582 // with weighted error 4.10e-17
583 // log weighted error 16.39
584 // significant figures required 15.96
585 // decimal places required 17.08
586
587         private final static double bth1cs[]={
588                 0.74060141026313850,
589                 -0.004571755659637690,
590                 0.000119818510964326,
591                 -0.000006964561891648,
592                 0.000000655495621447,
593                 -0.000000084066228945,
594                 0.000000013376886564,
595                 -0.000000002499565654,
596                 0.000000000529495100,
597                 -0.000000000124135944,
598                 0.000000000031656485,
599                 -0.000000000008668640,
600                 0.000000000002523758,
601                 -0.000000000000775085,
602                 0.000000000000249527,
603                 -0.000000000000083773,
604                 0.000000000000029205,
605                 -0.000000000000010534,
606                 0.000000000000003919,
607                 -0.000000000000001500,
608                 0.000000000000000589,
609                 -0.000000000000000237,
610                 0.000000000000000097,
611                 -0.000000000000000040};
612
613 // series for by1 on the interval 0. to 1.60000d+01
614 // with weighted error 1.87e-18
615 // log weighted error 17.73
616 // significant figures required 17.83
617 // decimal places required 18.30
618
619         private final static double by1cs[]={
620                 0.03208047100611908629,
621                 1.262707897433500450,
622                 0.00649996189992317500,
623                 -0.08936164528860504117,
624                 0.01325088122175709545,
625                 -0.00089790591196483523,
626                 0.00003647361487958306,
627                 -0.00000100137438166600,
628                 0.00000001994539657390,
629                 -0.00000000030230656018,
630                 0.00000000000360987815,
631                 -0.00000000000003487488,
632                 0.00000000000000027838,
633                 -0.00000000000000000186};
634
635
636
637         /**
638         * Evaluates a Chebyshev series.
639         * @param x value at which to evaluate series
640         * @param series the coefficients of the series
641         */
642         public static double chebyshev(double x, double series[]) {
643                 double twox,b0=0.0,b1=0.0,b2=0.0;
644                 twox=2*x;
645                 for(int i=series.length-1;i>-1;i--) {
646                         b2=b1;
647                         b1=b0;
648                         b0=twox*b1-b2+series[i];
649                 }
650                 return 0.5*(b0-b2);
651         }
652         /**
653         * Airy function.
654         * Based on the NETLIB Fortran function ai written by W. Fullerton.
655         */
656         public static double airy(double x) {
657                 if(x<-1.0) {
658                         return airyModPhase(x);
659                 } else if(x>1.0)
660                         return expAiry(x)*Math.exp(-2.0*x*Math.sqrt(x)/3.0);
661                 else {
662                         final double z=x*x*x;
663                         return 0.375+(chebyshev(z,aifcs)-x*(0.25+chebyshev(z,aigcs)));
664                 }
665         }
666         /**
667         * Airy modulus and phase.
668         * Based on the NETLIB Fortran subroutine r9aimp written by W. Fullerton.
669         * @return the real part, i.e. modulus*cos(phase).
670         */
671         private static double airyModPhase(double x) {
672                 double modulus, phase;
673                 if(x < -2.0) {
674                         double z = 16.0/(x*x*x)+1.0;
675                         modulus = 0.3125+chebyshev(z, am21cs);
676                         phase = -0.625+chebyshev(z, ath1cs);
677                 } else {
678                         double z = (16.0/(x*x*x)+9.0)/7.0;
679                         modulus = 0.3125+chebyshev(z, am22cs);
680                         phase = -0.625+chebyshev(z, ath2cs);
681                 }
682                 final double sqrtx = Math.sqrt(-x);
683                 modulus = Math.sqrt(modulus/sqrtx);
684                 phase = Math.PI/4.0-x*sqrtx*phase;
685                 return modulus*Math.cos(phase);
686         }
687         /**
688         * Exponential scaled Airy function.
689         * Based on the NETLIB Fortran function aie written by W. Fullerton.
690         */
691         private static double expAiry(double x) {
692                 if(x<-1.0) {
693                         return airyModPhase(x);
694                 } else if(x<=1.0) {
695                         final double z=x*x*x;
696                         return 0.375+(chebyshev(z,aifcs)-x*(0.25+chebyshev(z,aigcs)))*Math.exp(2.0*x*Math.sqrt(x)/3.0);
697                 } else {
698                         final double sqrtx=Math.sqrt(x);
699                         final double z=2.0/(x*sqrtx)-1.0;
700                         return (0.28125+chebyshev(z,aipcs))/Math.sqrt(sqrtx);
701                 }
702         }
703         /**
704         * Bessel function of first kind, order zero.
705         * Based on the NETLIB Fortran function besj0 written by W. Fullerton.
706         */
707         public static double besselFirstZero(double x) {
708                 double y=Math.abs(x);
709                 if(y>4.0) {
710                         final double z=32/(y*y)-1;
711                         final double amplitude=(0.75+chebyshev(z,bm0cs))/Math.sqrt(y);
712                         final double theta=y-Math.PI/4.0+chebyshev(z,bth0cs)/y;
713                         return amplitude*Math.cos(theta);
714                 } else if(y==0.0)
715                         return 1.0;
716                 else
717                         return chebyshev(0.125*y*y-1,bj0cs);
718         }
719         /**
720         * Modified Bessel function of first kind, order zero.
721         * Based on the NETLIB Fortran function besi0 written by W. Fullerton.
722         */
723         public static double modBesselFirstZero(double x) {
724                 double y=Math.abs(x);
725                 if(y>3.0)
726                         return Math.exp(y)*expModBesselFirstZero(x);
727                 else
728                         return 2.75+chebyshev(y*y/4.5-1.0, bi0cs);
729         }
730         /**
731         * Exponential scaled modified Bessel function of first kind, order zero.
732         * Based on the NETLIB Fortran function besi0e written by W. Fullerton.
733         */
734         private static double expModBesselFirstZero(double x) {
735                 final double y=Math.abs(x);
736                 if(y>3.0) {
737                         if(y>8.0)
738                                 return (0.375+chebyshev(16.0/y-1.0, ai02cs))/Math.sqrt(y);
739                         else
740                                 return (0.375+chebyshev((48.0/y-11.0)/5.0, ai0cs))/Math.sqrt(y);
741                 } else
742                         return Math.exp(-y)*(2.75+chebyshev(y*y/4.5-1.0, bi0cs));
743         }
744         /**
745         * Bessel function of first kind, order one.
746         * Based on the NETLIB Fortran function besj1 written by W. Fullerton.
747         */
748         public static double besselFirstOne(double x) {
749                 double y=Math.abs(x);
750                 if(y>4.0) {
751                         final double z=32.0/(y*y)-1.0;
752                         final double amplitude=(0.75+chebyshev(z, bm1cs))/Math.sqrt(y);
753                         final double theta=y-3.0*Math.PI/4.0+chebyshev(z, bth1cs)/y;
754                         return Math.abs(amplitude)*x*Math.cos(theta)/Math.abs(x);
755                 } else if(y==0.0)
756                         return 0.0;
757                 else
758                         return x*(0.25+chebyshev(0.125*y*y-1.0, bj1cs));
759         }
760         /**
761         * Modified Bessel function of first kind, order one.
762         * Based on the NETLIB Fortran function besi1 written by W. Fullerton.
763         */
764         public static double modBesselFirstOne(double x) {
765                 final double y=Math.abs(x);
766                 if(y>3.0)
767                         return Math.exp(y)*expModBesselFirstOne(x);
768                 else if(y==0.0)
769                         return 0.0;
770                 else
771                         return x*(0.875+chebyshev(y*y/4.5-1.0, bi1cs));
772         }
773         /**
774         * Exponential scaled modified Bessel function of first kind, order one.
775         * Based on the NETLIB Fortran function besi1e written by W. Fullerton.
776         */
777         private static double expModBesselFirstOne(double x) {
778                 final double y=Math.abs(x);
779                 if(y>3.0) {
780                         if(y>8.0)
781                                 return x/y*(0.375+chebyshev(16.0/y-1.0, ai12cs))/Math.sqrt(y);
782                         else
783                                 return x/y*(0.375+chebyshev((48.0/y-11.0)/5.0, ai1cs))/Math.sqrt(y);
784                 } else if(y==0.0)
785                         return 0.0;
786                 else
787                         return Math.exp(-y)*x*(0.875+chebyshev(y*y/4.5-1.0, bi1cs));
788         }
789         /**
790         * Bessel function of second kind, order zero.
791         * Based on the NETLIB Fortran function besy0 written by W. Fullerton.
792         */
793         public static double besselSecondZero(double x) {
794                 if(x>4.0) {
795                         final double z=32.0/(x*x)-1.0;
796                         final double amplitude=(0.75+chebyshev(z, bm0cs))/Math.sqrt(x);
797                         final double theta=x-Math.PI/4+chebyshev(z, bth0cs)/x;
798                         return amplitude*Math.sin(theta);
799                 } else
800                         return (Math.log(0.5)+Math.log(x))*besselFirstZero(x)+0.375+chebyshev(0.125*x*x-1.0,by0cs)*2.0/Math.PI;
801         }
802         /**
803         * Bessel function of second kind, order one.
804         * Based on the NETLIB Fortran function besy1 written by W. Fullerton.
805         */
806         public static double besselSecondOne(double x) {
807                 if(x>4.0) {
808                         final double z=32.0/(x*x)-1.0;
809                         final double amplitude=(0.75+chebyshev(z, bm1cs))/Math.sqrt(x);
810                         final double theta=x-3.0*Math.PI/4.0+chebyshev(z, bth1cs)/x;
811                         return amplitude*Math.sin(theta);
812                 } else
813                         return 2.0*Math.log(0.5*x)*besselFirstOne(x)/Math.PI+(0.5+chebyshev(0.125*x*x-1.0, by1cs))/x;
814         }
815
816         private final static double LOGSQRT2PI=Math.log(SQRT2PI);
817
818 // Gamma function related constants
819 private final static double g_p[] = { -1.71618513886549492533811,
820                 24.7656508055759199108314, -379.804256470945635097577,
821                 629.331155312818442661052, 866.966202790413211295064,
822                 -31451.2729688483675254357, -36144.4134186911729807069,
823                 66456.1438202405440627855 };
824         private final static double g_q[] = { -30.8402300119738975254353,
825                 315.350626979604161529144, -1015.15636749021914166146,
826                 -3107.77167157231109440444, 22538.1184209801510330112,
827                 4755.84627752788110767815, -134659.959864969306392456,
828                 -115132.259675553483497211 };
829         private final static double g_c[] = { -0.001910444077728,8.4171387781295e-4,
830                 -5.952379913043012e-4, 7.93650793500350248e-4,
831                 -0.002777777777777681622553, 0.08333333333333333331554247,
832                 0.0057083835261 };
833         /**
834         * The largest argument for which <code>gamma(x)</code> is representable in the machine.
835         */
836         public final static double GAMMA_X_MAX_VALUE = 171.624;
837
838         /**
839         * Gamma function.
840         * Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
841         * Applied Mathematics Division<BR>
842         * Argonne National Laboratory<BR>
843         * Argonne, IL 60439<BR>
844         * <P>
845         * References:
846         * <OL>
847         * <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.
848         * <LI>Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
849         * </OL></P><P>
850         * From the original documentation:
851         * </P><P>
852         * This routine calculates the GAMMA function for a real argument X.
853     * Computation is based on an algorithm outlined in reference 1.
854     * The program uses rational functions that approximate the GAMMA
855     * function to at least 20 significant decimal digits. Coefficients
856     * for the approximation over the interval (1,2) are unpublished.
857     * Those for the approximation for X .GE. 12 are from reference 2.
858     * The accuracy achieved depends on the arithmetic system, the
859     * compiler, the intrinsic functions, and proper selection of the
860     * machine-dependent constants.
861         * </P><P>
862     * Error returns:<BR>
863     * The program returns the value XINF for singularities or when overflow would occur.
864         * The computation is believed to be free of underflow and overflow.
865     * </P>
866         * @return Double.MAX_VALUE if overflow would occur, i.e. if abs(x) > 171.624
867         * @jsci.planetmath GammaFunction
868         * @author Jaco van Kooten
869     */
870         public static double gamma(double x) {
871                 double fact=1.0, xden, xnum;
872                 int i, n=0;
873                 double y=x, z, y1;
874                 boolean parity=false;
875                 double res, sum, ysq;
876
877                 if (y <= 0.0) {
878 // ----------------------------------------------------------------------
879 // Argument is negative
880 // ----------------------------------------------------------------------
881 y = -(x);
882                         y1 = (int)y;
883                         res = y - y1;
884                         if (res != 0.0) {
885                                 if (y1 != (((int)(y1*0.5)) * 2.0))
886                                         parity = true;
887                                 fact = -Math.PI/ Math.sin(Math.PI * res);
888                                 y++;
889                         } else
890                                 return Double.MAX_VALUE;
891                 }
892 // ----------------------------------------------------------------------
893 // Argument is positive
894 // ----------------------------------------------------------------------
895 if (y < EPS) {
896 // ----------------------------------------------------------------------
897 // Argument .LT. EPS
898 // ----------------------------------------------------------------------
899 if (y >= XMININ)
900                                 res = 1.0 / y;
901                         else
902                                 return Double.MAX_VALUE;
903                 } else if (y < 12.0) {
904                         y1 = y;
905                         if (y < 1.0) {
906 // ----------------------------------------------------------------------
907 // 0.0 .LT. argument .LT. 1.0
908 // ----------------------------------------------------------------------
909 z = y;
910                                 y++;
911                         } else {
912 // ----------------------------------------------------------------------
913 // 1.0 .LT. argument .LT. 12.0, reduce argument if necessary
914 // ----------------------------------------------------------------------
915 n = (int)y - 1;
916                                 y -= (double) n;
917                                 z = y - 1.0;
918                         }
919 // ----------------------------------------------------------------------
920 // Evaluate approximation for 1.0 .LT. argument .LT. 2.0
921 // ----------------------------------------------------------------------
922 xnum = 0.0;
923                         xden = 1.0;
924                         for (i = 0; i < 8; ++i) {
925                                 xnum = (xnum + g_p[i]) * z;
926                                 xden = xden * z + g_q[i];
927                         }
928                         res = xnum / xden + 1.0;
929                         if (y1 < y)
930 // ----------------------------------------------------------------------
931 // Adjust result for case 0.0 .LT. argument .LT. 1.0
932 // ----------------------------------------------------------------------
933 res /= y1;
934                         else if (y1 > y) {
935 // ----------------------------------------------------------------------
936 // Adjust result for case 2.0 .LT. argument .LT. 12.0
937 // ----------------------------------------------------------------------
938 for (i = 0; i < n; ++i) {
939                                         res *= y;
940                                         y++;
941                                 }
942                         }
943                 } else {
944 // ----------------------------------------------------------------------
945 // Evaluate for argument .GE. 12.0
946 // ----------------------------------------------------------------------
947 if (y <= GAMMA_X_MAX_VALUE) {
948                                 ysq = y * y;
949                                 sum = g_c[6];
950                                 for (i = 0; i < 6; ++i)
951                                         sum = sum / ysq + g_c[i];
952                                 sum = sum / y - y + LOGSQRT2PI;
953                                 sum += (y - 0.5) * Math.log(y);
954                                 res = Math.exp(sum);
955                         } else
956                                 return Double.MAX_VALUE;
957                 }
958 // ----------------------------------------------------------------------
959 // Final adjustments and return
960 // ----------------------------------------------------------------------
961 if (parity)
962                         res = -res;
963                 if (fact != 1.0)
964                         res = fact / res;
965                 return res;
966         }
967
968         /**
969         * The largest argument for which <code>logGamma(x)</code> is representable in the machine.
970         */
971         public final static double LOG_GAMMA_X_MAX_VALUE = 2.55e305;
972
973 // Log Gamma related constants
974 private final static double lg_d1 = -0.5772156649015328605195174;
975         private final static double lg_d2 = 0.4227843350984671393993777;
976         private final static double lg_d4 = 1.791759469228055000094023;
977         private final static double lg_p1[] = { 4.945235359296727046734888,
978                 201.8112620856775083915565, 2290.838373831346393026739,
979                 11319.67205903380828685045, 28557.24635671635335736389,
980                 38484.96228443793359990269, 26377.48787624195437963534,
981                 7225.813979700288197698961 };
982         private final static double lg_q1[] = { 67.48212550303777196073036,
983                 1113.332393857199323513008, 7738.757056935398733233834,
984                 27639.87074403340708898585, 54993.10206226157329794414,
985                 61611.22180066002127833352, 36351.27591501940507276287,
986                 8785.536302431013170870835 };
987         private final static double lg_p2[] = { 4.974607845568932035012064,
988                 542.4138599891070494101986, 15506.93864978364947665077,
989                 184793.2904445632425417223, 1088204.76946882876749847,
990                 3338152.967987029735917223, 5106661.678927352456275255,
991                 3074109.054850539556250927 };
992         private final static double lg_q2[] = { 183.0328399370592604055942,
993                 7765.049321445005871323047, 133190.3827966074194402448,
994                 1136705.821321969608938755, 5267964.117437946917577538,
995                 13467014.54311101692290052, 17827365.30353274213975932,
996                 9533095.591844353613395747 };
997         private final static double lg_p4[] = { 14745.02166059939948905062,
998                 2426813.369486704502836312, 121475557.4045093227939592,
999                 2663432449.630976949898078, 29403789566.34553899906876,
1000                170266573776.5398868392998, 492612579337.743088758812,
1001                560625185622.3951465078242 };
1002        private final static double lg_q4[] = { 2690.530175870899333379843,
1003                639388.5654300092398984238, 41355999.30241388052042842,
1004                1120872109.61614794137657, 14886137286.78813811542398,
1005                101680358627.2438228077304, 341747634550.7377132798597,
1006                446315818741.9713286462081 };
1007        private final static double lg_c[] = { -0.001910444077728,8.4171387781295e-4,
1008                -5.952379913043012e-4, 7.93650793500350248e-4,
1009                -0.002777777777777681622553, 0.08333333333333333331554247,
1010                0.0057083835261 };
1011// Rough estimate of the fourth root of logGamma_xBig
1012 private final static double lg_frtbig = 2.25e76;
1013        private final static double pnt68 = 0.6796875;
1014
1015// Function cache for logGamma
1016 private static double logGammaCache_res=0.0;
1017        private static double logGammaCache_x=0.0;
1018
1019    /**
1020        * The natural logarithm of the gamma function.
1021        * Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR>
1022        * Applied Mathematics Division<BR>
1023        * Argonne National Laboratory<BR>
1024        * Argonne, IL 60439<BR>
1025        * <P>
1026        * References:
1027        * <OL>
1028        * <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.
1029        * <LI>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.
1030        * <LI>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.
1031        * </OL></P><P>
1032        * From the original documentation:
1033        * </P><P>
1034        * This routine calculates the LOG(GAMMA) function for a positive real argument X.
1035        * Computation is based on an algorithm outlined in references 1 and 2.
1036        * The program uses rational functions that theoretically approximate LOG(GAMMA)
1037        * to at least 18 significant decimal digits. The approximation for X > 12 is from reference 3,
1038        * while approximations for X < 12.0 are similar to those in reference 1, but are unpublished.
1039        * The accuracy achieved depends on the arithmetic system, the compiler, the intrinsic functions,
1040        * and proper selection of the machine-dependent constants.
1041        * </P><P>
1042        * Error returns:<BR>
1043        * The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
1044        * The computation is believed to be free of underflow and overflow.
1045        * </P>
1046        * @return Double.MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
1047        * @author Jaco van Kooten
1048    */
1049        public static double logGamma(double x) {
1050                double xden, corr, xnum;
1051                int i;
1052                double y, xm1, xm2, xm4, res, ysq;
1053
1054                if (x==logGammaCache_x)
1055                        return logGammaCache_res;
1056                y = x;
1057                if (y > 0.0 && y <= LOG_GAMMA_X_MAX_VALUE) {
1058                        if (y <= EPS) {
1059                                res = -Math.log(y);
1060                        } else if (y <= 1.5) {
1061// ----------------------------------------------------------------------
1062// EPS .LT. X .LE. 1.5
1063// ----------------------------------------------------------------------
1064 if (y < pnt68) {
1065                                        corr = -Math.log(y);
1066                                        xm1 = y;
1067                                } else {
1068                                        corr = 0.0;
1069                                        xm1 = y - 1.0;
1070                                }
1071                                if (y <= 0.5 || y >= pnt68) {
1072                                        xden = 1.0;
1073                                        xnum = 0.0;
1074                                        for (i = 0; i < 8; i++) {
1075                                                xnum = xnum * xm1 + lg_p1[i];
1076                                                xden = xden * xm1 + lg_q1[i];
1077                                        }
1078                                        res = corr + xm1 * (lg_d1 + xm1 * (xnum / xden));
1079                                } else {
1080                                        xm2 = y - 1.0;
1081                                        xden = 1.0;
1082                                        xnum = 0.0;
1083                                        for (i = 0; i < 8; i++) {
1084                                                xnum = xnum * xm2 + lg_p2[i];
1085                                                xden = xden * xm2 + lg_q2[i];
1086                                        }
1087                                        res = corr + xm2 * (lg_d2 + xm2 * (xnum / xden));
1088                                }
1089                        } else if (y <= 4.0) {
1090// ----------------------------------------------------------------------
1091// 1.5 .LT. X .LE. 4.0
1092// ----------------------------------------------------------------------
1093 xm2 = y - 2.0;
1094                                xden = 1.0;
1095                                xnum = 0.0;
1096                                for (i = 0; i < 8; i++) {
1097                                        xnum = xnum * xm2 + lg_p2[i];
1098                                        xden = xden * xm2 + lg_q2[i];
1099                                }
1100                                res = xm2 * (lg_d2 + xm2 * (xnum / xden));
1101                        } else if (y <= 12.0) {
1102// ----------------------------------------------------------------------
1103// 4.0 .LT. X .LE. 12.0
1104// ----------------------------------------------------------------------
1105 xm4 = y - 4.0;
1106                                xden = -1.0;
1107                                xnum = 0.0;
1108                                for (i = 0; i < 8; i++) {
1109                                        xnum = xnum * xm4 + lg_p4[i];
1110                                        xden = xden * xm4 + lg_q4[i];
1111                                }
1112                                res = lg_d4 + xm4 * (xnum / xden);
1113                        } else {
1114// ----------------------------------------------------------------------
1115// Evaluate for argument .GE. 12.0
1116// ----------------------------------------------------------------------
1117 res = 0.0;
1118                                if (y <= lg_frtbig) {
1119                            res = lg_c[6];
1120                            ysq = y * y;
1121                            for (i = 0; i < 6; i++)
1122                                        res = res / ysq + lg_c[i];
1123                                }
1124                                res /= y;
1125                                corr = Math.log(y);
1126                                res = res + LOGSQRT2PI - 0.5 * corr;
1127                                res += y * (corr - 1.0);
1128                        }
1129                } else {
1130// ----------------------------------------------------------------------
1131// Return for bad arguments
1132// ----------------------------------------------------------------------
1133 res = Double.MAX_VALUE;
1134                }
1135// ----------------------------------------------------------------------
1136// Final adjustments and return
1137// ----------------------------------------------------------------------
1138 logGammaCache_x = x;
1139                logGammaCache_res = res;
1140                return res;
1141        }
1142
1143        private final static int MAX_ITERATIONS = 1000;
1144        // lower value = higher precision
1145 private final static double PRECISION = 4.0*EPS;
1146
1147        /**
1148        * Incomplete gamma function.
1149        * The computation is based on approximations presented in Numerical Recipes, Chapter 6.2 (W.H. Press et al, 1992).
1150        * @param a require a>=0
1151        * @param x require x>=0
1152        * @return 0 if x<0, a<=0 or a>2.55E305 to avoid errors and over/underflow
1153        * @author Jaco van Kooten
1154        */
1155        public static double incompleteGamma(double a, double x) {
1156                if (x <= 0.0 || a <= 0.0 || a > LOG_GAMMA_X_MAX_VALUE)
1157                        return 0.0;
1158                if (x < (a+1.0))
1159                        return gammaSeriesExpansion(a,x);
1160                else
1161                        return 1.0-gammaFraction(a,x);
1162        }
1163        /**
1164        * @author Jaco van Kooten
1165        */
1166        private static double gammaSeriesExpansion(double a, double x) {
1167                double ap = a;
1168                double del = 1.0/a;
1169                double sum = del;
1170                for (int n=1; n < MAX_ITERATIONS; n++) {
1171                        ++ap;
1172                        del *= x/ap;
1173                        sum += del;
1174                        if (del < sum*PRECISION)
1175                                return sum*Math.exp(-x + a*Math.log(x) - logGamma(a));
1176                }
1177                throw new RuntimeException JavaDoc("Maximum iterations exceeded: please file a bug report.");
1178        }
1179        /**
1180        * @author Jaco van Kooten
1181        */
1182        private static double gammaFraction(double a, double x) {
1183                double b=x+1.0-a;
1184                double c=1.0/XMININ;
1185                double d=1.0/b;
1186                double h=d;
1187                double del=0.0;
1188                double an;
1189                for (int i=1; i<MAX_ITERATIONS && Math.abs(del-1.0)>PRECISION; i++) {
1190                        an = -i*(i-a);
1191                        b += 2.0;
1192                        d=an*d+b;
1193                        c=b+an/c;
1194                        if (Math.abs(c) < XMININ)
1195                                c=XMININ;
1196                        if (Math.abs(d) < XMININ)
1197                                c=XMININ;
1198                        d=1.0/d;
1199                        del=d*c;
1200                        h *= del;
1201                }
1202                return Math.exp(-x + a*Math.log(x) - logGamma(a))*h;
1203        }
1204        /**
1205        * Beta function.
1206        * @param p require p>0
1207        * @param q require q>0
1208        * @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
1209        * @author Jaco van Kooten
1210        */
1211        public static double beta(double p, double q) {
1212                if (p <= 0.0 || q <= 0.0 || (p+q) > LOG_GAMMA_X_MAX_VALUE)
1213                        return 0.0;
1214                else
1215                        return Math.exp(logBeta(p,q));
1216    }
1217
1218// Function cache for logBeta
1219 private static double logBetaCache_res=0.0;
1220        private static double logBetaCache_p=0.0;
1221        private static double logBetaCache_q=0.0;
1222
1223        /**
1224        * The natural logarithm of the beta function.
1225        * @param p require p>0
1226        * @param q require q>0
1227        * @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
1228        * @author Jaco van Kooten
1229        */
1230        public static double logBeta(double p, double q) {
1231                if (p != logBetaCache_p || q != logBetaCache_q) {
1232                        logBetaCache_p = p;
1233                        logBetaCache_q = q;
1234                        if (p <= 0.0 || q <= 0.0 || (p+q) > LOG_GAMMA_X_MAX_VALUE)
1235                                logBetaCache_res = 0.0;
1236                        else
1237                                logBetaCache_res = logGamma(p)+logGamma(q)-logGamma(p+q);
1238                }
1239                return logBetaCache_res;
1240        }
1241        /**
1242        * Incomplete beta function.
1243        * The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
1244        * @param x require 0<=x<=1
1245        * @param p require p>0
1246        * @param q require q>0
1247        * @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
1248        * @author Jaco van Kooten
1249        */
1250    public static double incompleteBeta(double x, double p, double q) {
1251                if (x <= 0.0)
1252                return 0.0;
1253                else if (x >= 1.0)
1254                        return 1.0;
1255                else if (p <= 0.0 || q <= 0.0 || (p+q) > LOG_GAMMA_X_MAX_VALUE)
1256                return 0.0;
1257                else {
1258                        final double beta_gam=Math.exp(-logBeta(p,q) + p*Math.log(x) + q*Math.log(1.0-x));
1259                        if (x < (p+1.0)/(p+q+2.0))
1260                                return beta_gam*betaFraction(x,p,q)/p;
1261                else
1262                    return 1.0-(beta_gam*betaFraction(1.0-x,q,p)/q);
1263                }
1264        }
1265        /**
1266    * Evaluates of continued fraction part of incomplete beta function.
1267        * Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
1268        * @author Jaco van Kooten
1269        */
1270    private static double betaFraction(double x, double p, double q) {
1271            int m, m2;
1272            double sum_pq, p_plus, p_minus, c =1.0 , d, delta, h, frac;
1273            sum_pq = p + q;
1274            p_plus = p + 1.0;
1275            p_minus = p - 1.0;
1276            h=1.0-sum_pq*x/p_plus;
1277            if (Math.abs(h) < XMININ)
1278                        h=XMININ;
1279            h=1.0/h;
1280                frac = h;
1281            m=1;
1282            delta = 0.0;
1283            while (m <= MAX_ITERATIONS && Math.abs(delta-1.0) > PRECISION ) {
1284                        m2=2*m;
1285                // even index for d
1286 d=m*(q-m)*x/((p_minus+m2)*(p+m2));
1287                h=1.0+d*h;
1288                if (Math.abs(h) < XMININ)
1289                                h=XMININ;
1290                h=1.0/h;
1291                c=1.0+d/c;
1292                if (Math.abs(c) < XMININ)
1293                                c=XMININ;
1294                frac *= h*c;
1295                // odd index for d
1296 d = -(p+m)*(sum_pq+m)*x/((p+m2)*(p_plus+m2));
1297                h=1.0+d*h;
1298                if (Math.abs(h) < XMININ)
1299                                h=XMININ;
1300                h=1.0/h;
1301                c=1.0+d/c;
1302                if (Math.abs(c) < XMININ)
1303                                c=XMININ;
1304                        delta=h*c;
1305                        frac *= delta;
1306                        m++;
1307                }
1308                return frac;
1309        }
1310
1311// ====================================================
1312// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1313//
1314// Developed at SunSoft, a Sun Microsystems, Inc. business.
1315// Permission to use, copy, modify, and distribute this
1316// software is freely granted, provided that this notice
1317// is preserved.
1318// ====================================================
1319//
1320// x
1321// 2 |\
1322// erf(x) = --------- | exp(-t*t)dt
1323// sqrt(pi) \|
1324// 0
1325//
1326// erfc(x) = 1-erf(x)
1327// Note that
1328// erf(-x) = -erf(x)
1329// erfc(-x) = 2 - erfc(x)
1330//
1331// Method:
1332// 1. For |x| in [0, 0.84375]
1333// erf(x) = x + x*R(x^2)
1334// erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
1335// = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
1336// where R = P/Q where P is an odd poly of degree 8 and
1337// Q is an odd poly of degree 10.
1338// -57.90
1339// | R - (erf(x)-x)/x | <= 2
1340//
1341//
1342// Remark. The formula is derived by noting
1343// erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
1344// and that
1345// 2/sqrt(pi) = 1.128379167095512573896158903121545171688
1346// is close to one. The interval is chosen because the fix
1347// point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
1348// near 0.6174), and by some experiment, 0.84375 is chosen to
1349// guarantee the error is less than one ulp for erf.
1350//
1351// 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
1352// c = 0.84506291151 rounded to single (24 bits)
1353// erf(x) = sign(x) * (c + P1(s)/Q1(s))
1354// erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
1355// 1+(c+P1(s)/Q1(s)) if x < 0
1356// |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
1357// Remark: here we use the taylor series expansion at x=1.
1358// erf(1+s) = erf(1) + s*Poly(s)
1359// = 0.845.. + P1(s)/Q1(s)
1360// That is, we use rational approximation to approximate
1361// erf(1+s) - (c = (single)0.84506291151)
1362// Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
1363// where
1364// P1(s) = degree 6 poly in s
1365// Q1(s) = degree 6 poly in s
1366//
1367// 3. For x in [1.25,1/0.35(~2.857143)],
1368// erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
1369// erf(x) = 1 - erfc(x)
1370// where
1371// R1(z) = degree 7 poly in z, (z=1/x^2)
1372// S1(z) = degree 8 poly in z
1373//
1374// 4. For x in [1/0.35,28]
1375// erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
1376// = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
1377// = 2.0 - tiny (if x <= -6)
1378// erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
1379// erf(x) = sign(x)*(1.0 - tiny)
1380// where
1381// R2(z) = degree 6 poly in z, (z=1/x^2)
1382// S2(z) = degree 7 poly in z
1383//
1384// Note1:
1385// To compute exp(-x*x-0.5625+R/S), let s be a single
1386// precision number and s := x; then
1387// -x*x = -s*s + (s-x)*(s+x)
1388// exp(-x*x-0.5626+R/S) =
1389// exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
1390// Note2:
1391// Here 4 and 5 make use of the asymptotic series
1392// exp(-x*x)
1393// erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
1394// x*sqrt(pi)
1395// We use rational approximation to approximate
1396// g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
1397// Here is the error bound for R1/S1 and R2/S2
1398// |R1/S1 - f(x)| < 2**(-62.57)
1399// |R2/S2 - f(x)| < 2**(-61.52)
1400//
1401// 5. For inf > x >= 28
1402// erf(x) = sign(x) *(1 - tiny) (raise inexact)
1403// erfc(x) = tiny*tiny (raise underflow) if x > 0
1404// = 2 - tiny if x<0
1405//
1406// 7. Special case:
1407// erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
1408// erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
1409// erfc/erf(NaN) is NaN
1410//
1411
1412// Coefficients for approximation to erf on [0,0.84375]
1413 private final static double e_efx=1.28379167095512586316e-01;
1414// private final static double efx8=1.02703333676410069053e00;
1415 private final static double ePp[]={
1416                1.28379167095512558561e-01,
1417                -3.25042107247001499370e-01,
1418                -2.84817495755985104766e-02,
1419                -5.77027029648944159157e-03,
1420                -2.37630166566501626084e-05};
1421        private final static double eQq[]={
1422                3.97917223959155352819e-01,
1423                6.50222499887672944485e-02,
1424                5.08130628187576562776e-03,
1425                1.32494738004321644526e-04,
1426                -3.96022827877536812320e-06};
1427// Coefficients for approximation to erf in [0.84375,1.25]
1428 private final static double ePa[]={
1429                -2.36211856075265944077e-03,
1430                4.14856118683748331666e-01,
1431                -3.72207876035701323847e-01,
1432                3.18346619901161753674e-01,
1433                -1.10894694282396677476e-01,
1434                3.54783043256182359371e-02,
1435                -2.16637559486879084300e-03};
1436        private final static double eQa[]={
1437                1.06420880400844228286e-01,
1438                5.40397917702171048937e-01,
1439                7.18286544141962662868e-02,
1440                1.26171219808761642112e-01,
1441                1.36370839120290507362e-02,
1442                1.19844998467991074170e-02};
1443        private final static double e_erx=8.45062911510467529297e-01;
1444
1445    /**
1446        * Error function.
1447        * Based on C-code for the error function developed at Sun Microsystems.
1448        * @author Jaco van Kooten
1449    */
1450        public static double error(double x) {
1451            double P,Q,s,retval;
1452            final double abs_x = (x >= 0.0 ? x : -x);
1453                if ( abs_x < 0.84375 ) { // 0 < |x| < 0.84375
1454 if (abs_x < 3.7252902984619141e-9 ) // |x| < 2**-28
1455 retval = abs_x + abs_x*e_efx;
1456                        else {
1457                    s = x*x;
1458                    P = ePp[0]+s*(ePp[1]+s*(ePp[2]+s*(ePp[3]+s*ePp[4])));
1459                    Q = 1.0+s*(eQq[0]+s*(eQq[1]+s*(eQq[2]+s*(eQq[3]+s*eQq[4]))));
1460                    retval = abs_x + abs_x*(P/Q);
1461                }
1462            } else if (abs_x < 1.25) { // 0.84375 < |x| < 1.25
1463 s = abs_x-1.0;
1464                        P = ePa[0]+s*(ePa[1]+s*(ePa[2]+s*(ePa[3]+s*(ePa[4]+s*(ePa[5]+s*ePa[6])))));
1465                        Q = 1.0+s*(eQa[0]+s*(eQa[1]+s*(eQa[2]+s*(eQa[3]+s*(eQa[4]+s*eQa[5])))));
1466                        retval = e_erx + P/Q;
1467           } else if (abs_x >= 6.0)
1468                        retval = 1.0;
1469                else // 1.25 < |x| < 6.0
1470 retval = 1.0-complementaryError(abs_x);
1471                return (x >= 0.0) ? retval : -retval;
1472        }
1473
1474// Coefficients for approximation to erfc in [1.25,1/.35]
1475 private final static double eRa[]={
1476                -9.86494403484714822705e-03,
1477                -6.93858572707181764372e-01,
1478                -1.05586262253232909814e01,
1479                -6.23753324503260060396e01,
1480                -1.62396669462573470355e02,
1481                -1.84605092906711035994e02,
1482                -8.12874355063065934246e01,
1483                -9.81432934416914548592e00};
1484        private final static double eSa[]={
1485                1.96512716674392571292e01,
1486                1.37657754143519042600e02,
1487                4.34565877475229228821e02,
1488                6.45387271733267880336e02,
1489                4.29008140027567833386e02,
1490                1.08635005541779435134e02,
1491                6.57024977031928170135e00,
1492                -6.04244152148580987438e-02};
1493// Coefficients for approximation to erfc in [1/.35,28]
1494 private final static double eRb[]={
1495                -9.86494292470009928597e-03,
1496                -7.99283237680523006574e-01,
1497                -1.77579549177547519889e01,
1498                -1.60636384855821916062e02,
1499                -6.37566443368389627722e02,
1500                -1.02509513161107724954e03,
1501                -4.83519191608651397019e02};
1502        private final static double eSb[]={
1503                3.03380607434824582924e01,
1504                3.25792512996573918826e02,
1505                1.53672958608443695994e03,
1506                3.19985821950859553908e03,
1507                2.55305040643316442583e03,
1508                4.74528541206955367215e02,
1509                -2.24409524465858183362e01};
1510
1511    /**
1512        * Complementary error function.
1513        * Based on C-code for the error function developed at Sun Microsystems.
1514        * @author Jaco van Kooten
1515    */
1516    public static double complementaryError(double x) {
1517            double s,retval,R,S;
1518                final double abs_x =(x>=0.0 ? x : -x);
1519                if (abs_x < 1.25)
1520                retval = 1.0-error(abs_x);
1521                else if (abs_x > 28.0)
1522                        retval=0.0;
1523                else { // 1.25 < |x| < 28
1524 s = 1.0/(abs_x*abs_x);
1525                        if (abs_x < 2.8571428) { // ( |x| < 1/0.35 )
1526 R=eRa[0]+s*(eRa[1]+s*(eRa[2]+s*(eRa[3]+s*(eRa[4]+s*(eRa[5]+s*(eRa[6]+s*eRa[7]))))));
1527                        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])))))));
1528                } else { // ( |x| > 1/0.35 )
1529 R=eRb[0]+s*(eRb[1]+s*(eRb[2]+s*(eRb[3]+s*(eRb[4]+s*(eRb[5]+s*eRb[6])))));
1530                        S=1.0+s*(eSb[0]+s*(eSb[1]+s*(eSb[2]+s*(eSb[3]+s*(eSb[4]+s*(eSb[5]+s*eSb[6]))))));
1531            }
1532                        retval = Math.exp(-x*x - 0.5625 + R/S)/abs_x;
1533                }
1534                return (x >= 0.0) ? retval : 2.0-retval;
1535    }
1536}
1537
1538
Popular Tags