View Javadoc

1   /*
2    * =========================================================================
3    * Copyright (C) 1997 - 1998 by Visual Numerics, Inc. All rights reserved.
4    *
5    * Permission to use, copy, modify, and distribute this software is freely
6    * granted by Visual Numerics, Inc., provided that the copyright notice
7    * above and the following warranty disclaimer are preserved in human
8    * readable form.
9    *
10   * Because this software is licenses free of charge, it is provided
11   * "AS IS", with NO WARRANTY.  TO THE EXTENT PERMITTED BY LAW, VNI
12   * DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
13   * TO ITS PERFORMANCE, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
14   * VNI WILL NOT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING OUT OF THE USE
15   * OF OR INABILITY TO USE THIS SOFTWARE, INCLUDING BUT NOT LIMITED TO DIRECT,
16   * INDIRECT, SPECIAL, CONSEQUENTIAL, PUNITIVE, AND EXEMPLARY DAMAGES, EVEN
17   * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
18   *
19   * =========================================================================
20   */
21  
22  //package com.visualnumerics.javagrande;
23  package edu.sc.seis.TauP;
24  
25  /***
26   *  Collection of special functions.
27   */
28  public class Sfun
29  {
30      /*** The smallest relative spacing for doubles.*/
31      public final static double EPSILON_SMALL = 1.1102230246252e-16;
32      
33      /*** The largest relative spacing for doubles. */
34      public final static double EPSILON_LARGE = 2.2204460492503e-16;
35  
36      
37      /***
38       *  Private contructor, so nobody can make an instance of this class.
39       */
40      private Sfun()
41      {
42      }
43  
44  
45      /*
46       *  Evaluate a Chebyschev series
47       */
48      static double csevl(double x, double coef[])
49      {
50          double  b0, b1, b2, twox;
51          int     i;
52          b1 = 0.0;
53          b0 = 0.0;
54          b2 = 0.0;
55          twox = 2.0*x;
56          for (i = coef.length-1;  i >= 0;  i--) {
57              b2 = b1;
58              b1 = b0;
59              b0 = twox*b1 - b2 + coef[i];
60          }
61          return 0.5*(b0-b2);
62      }
63  
64  
65      // Series on [0,0.0625]
66      private static final double COT_COEF[] = {
67          .240259160982956302509553617744970e+0,
68          -.165330316015002278454746025255758e-1,
69          -.429983919317240189356476228239895e-4,
70          -.159283223327541046023490851122445e-6,
71          -.619109313512934872588620579343187e-9,
72          -.243019741507264604331702590579575e-11,
73          -.956093675880008098427062083100000e-14,
74          -.376353798194580580416291539706666e-16,
75          -.148166574646746578852176794666666e-18
76      };
77  
78      /***
79       *  Returns the cotangent of a double.
80       *  @param  x   A double value.
81       *  @return  The cotangent of x.
82       *  If x is NaN, the result is NaN.
83       */
84      static public double cot(double x)
85      {
86          double ans, ainty, ainty2, prodbg, y, yrem;
87          double pi2rec = 0.011619772367581343075535053490057; //  2/PI - 0.625
88  
89          y = Math.abs(x);
90  
91          if (y > 4.5036e+15) {
92              // 4.5036e+15 = 1.0/EPSILON_LARGE
93              return Double.NaN;
94          }
95  
96          // Carefully compute
97          // Y * (2/PI) = (AINT(Y) + REM(Y)) * (.625 + PI2REC)
98          //      = AINT(.625*Y) + REM(.625*Y) + Y*PI2REC  =  AINT(.625*Y) + Z
99          //      = AINT(.625*Y) + AINT(Z) + REM(Z)
100         ainty  = (int)y;
101         yrem   = y - ainty;
102         prodbg = 0.625*ainty;
103         ainty  = (int)prodbg;
104         y      = (prodbg-ainty) + 0.625*yrem + y*pi2rec;
105         ainty2 = (int)y;
106         ainty  = ainty + ainty2;
107         y      = y - ainty2;
108 
109         int ifn = (int)(ainty%2.0);
110         if (ifn == 1) y = 1.0 - y;
111 
112         if (y == 0.0) {
113             ans = Double.POSITIVE_INFINITY;
114         } else if (y <= 1.82501e-08) {
115             // 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)
116             ans = 1.0/y;
117         } else if (y <= 0.25) {
118             ans = (0.5+csevl(32.0*y*y-1.0,COT_COEF))/y;
119         } else if (y <= 0.5) {
120             ans = (0.5+csevl(8.0*y*y-1.0,COT_COEF))/(0.5*y);
121             ans = (ans*ans-1.0)*0.5/ans;
122         } else {
123             ans = (0.5+csevl(2.0*y*y-1.0,COT_COEF))/(0.25*y);
124             ans = (ans*ans-1.0)*0.5/ans;
125             ans = (ans*ans-1.0)*0.5/ans;
126         }
127         if (x != 0.0) ans = sign(ans,x);
128         if (ifn == 1) ans = -ans;
129         return ans;
130     }
131 
132     /***
133      *  Returns the common (base 10) logarithm of a double.
134      *  @param  x   A double value.
135      *  @return  The common logarithm of x.
136      */
137     static public double log10(double x)
138     {
139         //if (Double.isNaN(x)) return Double.NaN;
140         return 0.43429448190325182765*Math.log(x);
141     }
142 
143     /*
144      *  Returns the value of x with the sign of y.
145      */
146     static private double sign(double x, double y)
147     {
148         double abs_x = ((x < 0) ? -x : x);
149         return (y < 0.0) ? -abs_x : abs_x;
150     }
151 
152     // Series on the interval [0,1]
153     private static final double SINH_COEF[] = {
154         0.1730421940471796,
155         0.08759422192276048,
156         0.00107947777456713,
157         0.00000637484926075,
158         0.00000002202366404,
159         0.00000000004987940,
160         0.00000000000007973,
161         0.00000000000000009};
162 
163     /***
164      *  Returns the inverse (arc) hyperbolic sine of a double.
165      *  @param  x   A double value.
166      *  @return  The arc hyperbolic sine of x.
167      *  If x is NaN or less than one, the result is NaN.
168      */
169     static public double sinh(double x)
170     {
171         double  ans;
172         double  y = Math.abs(x);
173         
174         if (Double.isNaN(x)) {
175             ans = Double.NaN;
176         } else if (Double.isInfinite(y)) {
177             return x;
178         } else if (y < 2.58096e-08) {
179             // 2.58096e-08 = Math.sqrt(6.0*EPSILON_SMALL)
180             ans = x;
181         } else if (y <= 1.0) {
182             ans = x*(1.0+csevl(2.0*x*x-1.0,SINH_COEF));
183         } else {
184             y = Math.exp(y);
185             if (y >= 94906265.62) {
186                 // 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)
187                 ans = sign(0.5*y,x);
188             } else {
189                 ans = sign(0.5*(y-1.0/y),x);
190             }
191         }
192         return ans;
193     }
194 
195     /***
196      *  Returns the hyperbolic cosine of a double.
197      *  @param  x   A double value.
198      *  @return  The hyperbolic cosine of x.
199      *  If x is NaN, the result is NaN.
200      */
201     static public double cosh(double x)
202     {
203         double  ans;
204         double  y = Math.exp(Math.abs(x));
205 
206         if (Double.isNaN(x)) {
207             ans = Double.NaN;
208         } else if (Double.isInfinite(x)) {
209             ans = x;
210         } else if (y < 94906265.62) {
211             // 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)
212             ans = 0.5*(y+1.0/y);
213         } else {
214             ans = 0.5*y;
215         }
216         return ans;
217     }
218     
219     // Series on [0,1]
220     private static final double TANH_COEF[] = {
221         -.25828756643634710,
222         -.11836106330053497,
223         .009869442648006398,
224         -.000835798662344582,
225         .000070904321198943,
226         -.000006016424318120,
227         .000000510524190800,
228         -.000000043320729077,
229         .000000003675999055,
230         -.000000000311928496,
231         .000000000026468828,
232         -.000000000002246023,
233         .000000000000190587,
234         -.000000000000016172,
235         .000000000000001372,
236         -.000000000000000116,
237         .000000000000000009};
238 
239     /***
240      *  Returns the hyperbolic tangent of a double.
241      *  @param  x   A double value.
242      *  @return  The hyperbolic tangent of x.
243      */
244     static public double tanh(double x)
245     {
246         double  ans, y;
247         y = Math.abs(x);
248         
249         if (Double.isNaN(x)) {
250             ans = Double.NaN;
251         } else if (y < 1.82501e-08) {
252             // 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)
253             ans = x;
254         } else if (y <= 1.0) {
255             ans = x*(1.0+csevl(2.0*x*x-1.0,TANH_COEF));
256         } else if (y < 7.977294885) {
257             // 7.977294885 = -0.5*Math.log(EPSILON_SMALL)
258             y = Math.exp(y);
259             ans = sign((y-1.0/y)/(y+1.0/y),x);
260         } else {
261             ans = sign(1.0,x);
262         }
263         return ans;
264     }
265     // Series on the interval [0,1]
266     private static final double ASINH_COEF[] = {
267          -.12820039911738186343372127359268e+0,
268         -.58811761189951767565211757138362e-1,
269         .47274654322124815640725249756029e-2,
270         -.49383631626536172101360174790273e-3,
271         .58506207058557412287494835259321e-4,
272         -.74669983289313681354755069217188e-5,
273         .10011693583558199265966192015812e-5,
274         -.13903543858708333608616472258886e-6,
275         .19823169483172793547317360237148e-7,
276         -.28847468417848843612747272800317e-8,
277         .42672965467159937953457514995907e-9,
278         -.63976084654366357868752632309681e-10,
279         .96991686089064704147878293131179e-11,
280         -.14844276972043770830246658365696e-11,
281         .22903737939027447988040184378983e-12,
282         -.35588395132732645159978942651310e-13,
283         .55639694080056789953374539088554e-14,
284         -.87462509599624678045666593520162e-15,
285         .13815248844526692155868802298129e-15,
286         -.21916688282900363984955142264149e-16,
287         .34904658524827565638313923706880e-17
288     };
289 
290     /***
291      *  Returns the inverse (arc) hyperbolic sine of a double.
292      *  @param  x   A double value.
293      *  @return  The arc hyperbolic sine of x.
294      *  If x is NaN, the result is NaN.
295      */
296     static public double asinh(double x)
297     {
298         double  ans;
299         double  y = Math.abs(x);
300     
301         if (Double.isNaN(x)) {
302             ans = Double.NaN;
303         } else if (y <= 1.05367e-08) {
304             // 1.05367e-08 = Math.sqrt(EPSILON_SMALL)
305             ans = x;
306         } else if (y <= 1.0) {
307             ans = x*(1.0+csevl(2.0*x*x-1.0,ASINH_COEF));
308         } else if (y < 94906265.62) {
309             // 94906265.62 = 1/Math.sqrt(EPSILON_SMALL)
310             ans = Math.log(y+Math.sqrt(y*y+1.0));
311         } else {
312             ans = 0.69314718055994530941723212145818 + Math.log(y);
313         }
314         if (x < 0.0) ans = -ans;
315         return ans;
316     }
317 
318     /***
319      *  Returns the inverse (arc) hyperbolic cosine of a double.
320      *  @param  x   A double value.
321      *  @return  The arc hyperbolic cosine of x.
322      *  If x is NaN or less than one, the result is NaN.
323      */
324     static public double acosh(double x)
325     {
326         double ans;
327         
328         if (Double.isNaN(x) || x < 1) {
329             ans = Double.NaN;
330         } else if (x < 94906265.62) {
331             // 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)
332             ans = Math.log(x+Math.sqrt(x*x-1.0));
333         } else {
334             ans = 0.69314718055994530941723212145818 + Math.log(x);
335         }
336         return ans;
337     }
338 
339 
340     // Series on the interval [0,0.25]
341     private static final double ATANH_COEF[] = {
342         .9439510239319549230842892218633e-1,
343         .4919843705578615947200034576668e-1,
344         .2102593522455432763479327331752e-2,
345         .1073554449776116584640731045276e-3,
346         .5978267249293031478642787517872e-5,
347         .3505062030889134845966834886200e-6,
348         .2126374343765340350896219314431e-7,
349         .1321694535715527192129801723055e-8,
350         .8365875501178070364623604052959e-10,
351         .5370503749311002163881434587772e-11,
352         .3486659470157107922971245784290e-12,
353         .2284549509603433015524024119722e-13,
354         .1508407105944793044874229067558e-14,
355         .1002418816804109126136995722837e-15,
356         .6698674738165069539715526882986e-17,
357         .4497954546494931083083327624533e-18
358     };
359 
360     /***
361      *  Returns the inverse (arc) hyperbolic tangent of a double.
362      *  @param  x   A double value.
363      *  @return  The arc hyperbolic tangent of x.
364      *  If x is NaN or |x|>1, the result is NaN.
365      */
366     static public double atanh(double x)
367     {
368         double  y = Math.abs(x);
369         double  ans;
370 
371         if (Double.isNaN(x)) {
372             ans = Double.NaN;
373         } else if (y < 1.82501e-08) {
374             // 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)
375             ans = x;
376         } else if (y <= 0.5) {
377             ans = x*(1.0+csevl(8.0*x*x-1.0,ATANH_COEF));
378         } else if (y < 1.0) {
379             ans = 0.5*Math.log((1.0+x)/(1.0-x));
380         } else if (y == 1.0) {
381             ans = x*Double.POSITIVE_INFINITY;
382         } else {
383             ans = Double.NaN;
384         }
385         return ans;
386     }
387 
388 
389     /***
390      *  Returns the factorial of an integer.
391      *  @param  n   An integer value.
392      *  @return  The factorial of n, n!.
393      *  If x is negative, the result is NaN.
394      */
395     static public double fact(int n)
396     {
397         double ans = 1;
398 
399         if (Double.isNaN(n) || n < 0) {
400             ans = Double.NaN;
401         } else if (n > 170) {
402             // The 171! is too large to fit in a double.
403             ans = Double.POSITIVE_INFINITY;
404         } else {
405             for (int k = 2;  k <= n;  k++)
406                 ans *= k;
407         }
408         return ans;
409     }
410 
411 
412     // Series on the interval [0,1]
413     private static final double GAMMA_COEF[] = {
414         .8571195590989331421920062399942e-2,
415         .4415381324841006757191315771652e-2,
416         .5685043681599363378632664588789e-1,
417         -.4219835396418560501012500186624e-2,
418         .1326808181212460220584006796352e-2,
419         -.1893024529798880432523947023886e-3,
420         .3606925327441245256578082217225e-4,
421         -.6056761904460864218485548290365e-5,
422         .1055829546302283344731823509093e-5,
423         -.1811967365542384048291855891166e-6,
424         .3117724964715322277790254593169e-7,
425         -.5354219639019687140874081024347e-8,
426         .9193275519859588946887786825940e-9,
427         -.1577941280288339761767423273953e-9,
428         .2707980622934954543266540433089e-10,
429         -.4646818653825730144081661058933e-11,
430         .7973350192007419656460767175359e-12,
431         -.1368078209830916025799499172309e-12,
432         .2347319486563800657233471771688e-13,
433         -.4027432614949066932766570534699e-14,
434         .6910051747372100912138336975257e-15,
435         -.1185584500221992907052387126192e-15,
436         .2034148542496373955201026051932e-16,
437         -.3490054341717405849274012949108e-17,
438         .5987993856485305567135051066026e-18,
439         -.1027378057872228074490069778431e-18
440     };
441 
442     /***
443      *  Returns the Gamma function of a double.
444      *  @param  x   A double value.
445      *  @return  The Gamma function of x.
446      *  If x is a negative integer, the result is NaN.
447      */
448     static public double gamma(double x)
449     {
450         double  ans;
451         double  y = Math.abs(x);
452 
453         if (y <= 10.0) {
454             /*
455              * Compute gamma(x) for |x|<=10.
456              * First reduce the interval and  find gamma(1+y) for 0 <= y < 1.
457              */
458             int n = (int)x;
459             if (x < 0.0)  n--;
460             y = x - n;
461             n--;
462             ans = 0.9375 + csevl(2.0*y-1.0, GAMMA_COEF);
463             if (n == 0) {
464             } else if (n < 0) {
465                 // Compute gamma(x) for x < 1
466                 n = -n;
467                 if (x == 0.0) {
468                     ans = Double.NaN;
469                 } else if (y < 1.0/Double.MAX_VALUE) {
470                     ans = Double.POSITIVE_INFINITY;
471                 } else {
472                     double xn = n - 2;
473                     if (x < 0.0 && x + xn == 0.0) {
474                         ans = Double.NaN;
475                     } else {
476                         for (int i = 0; i < n; i++) {
477                             ans /= x + i;
478                         }
479                     }
480                 }
481             } else {    // gamma(x) for x >= 2.0
482                 for (int i = 1; i <= n; i++) {
483                     ans *= y + i;
484                 }
485             }
486         } else {  // gamma(x) for |x| > 10
487             if (x > 171.614) {
488                 ans = Double.POSITIVE_INFINITY;
489             } else if (x < -170.56) {
490                 ans = 0.0; // underflows
491             } else {
492                 // 0.9189385332046727 = 0.5*log(2*PI)
493                 ans = Math.exp((y-0.5)*Math.log(y)-y+0.9189385332046727+r9lgmc(y));
494                 if (x < 0.0) {
495                     double sinpiy = Math.sin(Math.PI * y);
496                     if (sinpiy == 0 || Math.round(y) == y) {
497                         ans = Double.NaN;
498                     } else {
499                         ans = -Math.PI / (y * sinpiy * ans);
500                     }
501                 }
502             }
503         }
504         return ans;
505     }
506 
507     /***
508      *  Returns the logarithm of the Gamma function of a double.
509      *  @param  x   A double value.
510      *  @return  The natural logarithm of the Gamma function of x.
511      *  If x is a negative integer, the result is NaN.
512      */
513     static public double logGamma(double x)
514     {
515         double  ans, sinpiy, y;
516 
517         y = Math.abs(x);
518 
519         if (y <= 10) {
520             ans = Math.log(Math.abs(gamma(x)));
521         } else if (x > 0) {
522             // A&S 6.1.40
523             // 0.9189385332046727 = 0.5*log(2*PI)
524             ans = 0.9189385332046727 + (x-0.5)*Math.log(x) - x + r9lgmc(y);
525         } else {
526             sinpiy = Math.abs(Math.sin(Math.PI * y));
527             if (sinpiy == 0  || Math.round(y) == y) {
528                 // The argument for the function can not be a negative integer.
529                 ans = Double.NaN;
530             } else {
531                 ans = 0.22579135264472743236 + (x-0.5)*Math.log(y) - x - Math.log(sinpiy) - r9lgmc(y);
532             }
533         }
534         return ans;
535     }
536 
537 
538     //  Series for the interval [0,0.01]
539     private static final double R9LGMC_COEF[] =
540     {
541         .166638948045186324720572965082e0,
542         -.138494817606756384073298605914e-4,
543         .981082564692472942615717154749e-8,
544         -.180912947557249419426330626672e-10,
545         .622109804189260522712601554342e-13,
546         -.339961500541772194430333059967e-15,
547         .268318199848269874895753884667e-17
548     };
549 
550     /*
551      *  Returns the log gamma correction term for argument
552      *  values greater than or equal to 10.0.
553      */
554     static double r9lgmc(double x)
555     {
556         double  ans;
557 
558         if (x < 10.0) {
559             ans = Double.NaN;
560         } else if (x < 9.490626562e+07) {
561             // 9.490626562e+07 = 1/Math.sqrt(EPSILON_SMALL)
562             double y = 10.0/x;
563             ans = csevl(2.0*y*y-1.0, R9LGMC_COEF) /  x;
564         } else if (x < 1.39118e+11) {
565             // 1.39118e+11 = exp(min(log(amach(2) / 12.0), -log(12.0 * amach(1))));
566             // See A&S 6.1.41
567             ans = 1.0/(12.0*x);
568         } else {
569             ans = 0.0; // underflows
570         }
571         return ans;
572     }
573 
574     /***
575      *  Returns the logarithm of the Beta function.
576      *  @param  a   A double value.
577      *  @param  b   A double value.
578      *  @return  The natural logarithm of the Beta function.
579      */
580     static public double logBeta(double a, double b)
581     {
582         double  corr, ans;
583         double  p = Math.min(a, b);
584         double  q = Math.max(a, b);
585 
586         if (p <= 0.0) {
587             ans = Double.NaN;
588         } else if (p >= 10.0) {
589             // P and Q are large;
590             corr = r9lgmc(p) + r9lgmc(q) - r9lgmc(p+q);
591             double temp = dlnrel(-p/(p+q));
592             ans = -0.5*Math.log(q) + 0.918938533204672741780329736406 + corr + (p-0.5)*Math.log(p/(p+q)) + q*temp;
593         } else if (q >= 10.0) {
594             // P is small, but Q is large
595             corr = Sfun.r9lgmc(q) - r9lgmc(p+q);
596             //  Check from underflow from r9lgmc
597             ans = logGamma(p) + corr + p - p*Math.log(p+q) + (q-0.5)*dlnrel(-p/(p+q));
598         } else {
599             // P and Q are small;
600             ans = Math.log(gamma(p)*(gamma(q)/gamma(p+q)));
601         }
602         return ans;
603     }
604 
605     // Series on [-0.375,0.375]
606     final private static double ALNRCS_COEF[] = {
607        .103786935627437698006862677191e1,
608       -.133643015049089180987660415531,
609       .194082491355205633579261993748e-1,
610       -.301075511275357776903765377766e-2,
611       .486946147971548500904563665091e-3,
612       -.810548818931753560668099430086e-4,
613       .137788477995595247829382514961e-4,
614       -.238022108943589702513699929149e-5,
615       .41640416213865183476391859902e-6,
616       -.73595828378075994984266837032e-7,
617       .13117611876241674949152294345e-7,
618       -.235467093177424251366960923302e-8,
619       .425227732760349977756380529626e-9,
620       -.771908941348407968261081074933e-10,
621       .140757464813590699092153564722e-10,
622       -.257690720580246806275370786276e-11,
623       .473424066662944218491543950059e-12,
624       -.872490126747426417453012632927e-13,
625       .161246149027405514657398331191e-13,
626       -.298756520156657730067107924168e-14,
627       .554807012090828879830413216973e-15,
628       -.103246191582715695951413339619e-15,
629       .192502392030498511778785032449e-16,
630       -.359550734652651500111897078443e-17,
631       .672645425378768578921945742268e-18,
632       -.126026241687352192520824256376e-18
633     };
634     
635     /*
636      *  Correction term used by logBeta.
637      */
638     private static double dlnrel(double x)
639     {
640         double ans;
641         
642         if (x <= -1.0) {
643             ans = Double.NaN;
644         } else if (Math.abs(x) <= 0.375) {
645             ans = x*(1.0 - x*Sfun.csevl(x/.375, ALNRCS_COEF));
646         } else {
647             ans = Math.log(1.0 + x);
648         }
649         return ans;
650     }
651 
652 
653     // Series on [0,1]
654     private static final double ERFC_COEF[] = {
655      -.490461212346918080399845440334e-1,
656      -.142261205103713642378247418996e0,
657      .100355821875997955757546767129e-1,
658      -.576876469976748476508270255092e-3,
659      .274199312521960610344221607915e-4,
660      -.110431755073445076041353812959e-5,
661      .384887554203450369499613114982e-7,
662      -.118085825338754669696317518016e-8,
663      .323342158260509096464029309534e-10,
664      -.799101594700454875816073747086e-12,
665      .179907251139614556119672454866e-13,
666      -.371863548781869263823168282095e-15,
667      .710359900371425297116899083947e-17,
668      -.126124551191552258324954248533e-18
669     };
670 
671     // Series on [0.25,1.00]
672     private static final double ERFC2_COEF[] = {
673      -.69601346602309501127391508262e-1,
674      -.411013393626208934898221208467e-1,
675      .391449586668962688156114370524e-2,
676      -.490639565054897916128093545077e-3,
677      .715747900137703638076089414183e-4,
678      -.115307163413123283380823284791e-4,
679      .199467059020199763505231486771e-5,
680      -.364266647159922287393611843071e-6,
681      .694437261000501258993127721463e-7,
682      -.137122090210436601953460514121e-7,
683      .278838966100713713196386034809e-8,
684      -.581416472433116155186479105032e-9,
685      .123892049175275318118016881795e-9,
686      -.269063914530674343239042493789e-10,
687      .594261435084791098244470968384e-11,
688      -.133238673575811957928775442057e-11,
689      .30280468061771320171736972433e-12,
690      -.696664881494103258879586758895e-13,
691      .162085454105392296981289322763e-13,
692      -.380993446525049199987691305773e-14,
693      .904048781597883114936897101298e-15,
694      -.2164006195089607347809812047e-15,
695      .522210223399585498460798024417e-16,
696      -.126972960236455533637241552778e-16,
697      .310914550427619758383622741295e-17,
698      -.766376292032038552400956671481e-18,
699      .190081925136274520253692973329e-18
700     };
701 
702     // Series on [0,0.25]
703     private static final double ERFCC_COEF[] = {
704      .715179310202924774503697709496e-1,
705      -.265324343376067157558893386681e-1,
706      .171115397792085588332699194606e-2,
707      -.163751663458517884163746404749e-3,
708      .198712935005520364995974806758e-4,
709      -.284371241276655508750175183152e-5,
710      .460616130896313036969379968464e-6,
711      -.822775302587920842057766536366e-7,
712      .159214187277090112989358340826e-7,
713      -.329507136225284321486631665072e-8,
714      .72234397604005554658126115389e-9,
715      -.166485581339872959344695966886e-9,
716      .401039258823766482077671768814e-10,
717      -.100481621442573113272170176283e-10,
718      .260827591330033380859341009439e-11,
719      -.699111056040402486557697812476e-12,
720      .192949233326170708624205749803e-12,
721      -.547013118875433106490125085271e-13,
722      .158966330976269744839084032762e-13,
723      -.47268939801975548392036958429e-14,
724      .14358733767849847867287399784e-14,
725      -.444951056181735839417250062829e-15,
726      .140481088476823343737305537466e-15,
727      -.451381838776421089625963281623e-16,
728      .147452154104513307787018713262e-16,
729      -.489262140694577615436841552532e-17,
730      .164761214141064673895301522827e-17,
731      -.562681717632940809299928521323e-18,
732      .194744338223207851429197867821e-18
733     };
734 
735 
736 
737     /***
738      *  Returns the error function of a double.
739      *  @param  x   A double value.
740      *  @return  The error function of x.
741      */
742     static public double erf(double x)
743     {
744         double  ans;
745         double  y = Math.abs(x);
746 
747         if (y <= 1.49012e-08) {
748             // 1.49012e-08 = Math.sqrt(2*EPSILON_SMALL)
749             ans = 2 * x / 1.77245385090551602729816748334;
750         } else if (y <= 1) {
751             ans = x * (1 + csevl(2 * x * x - 1, ERFC_COEF));
752         } else if (y < 6.013687357) {
753             // 6.013687357 = Math.sqrt(-Math.log(1.77245385090551602729816748334 * EPSILON_SMALL))
754             ans = sign(1 - erfc(y), x);
755         } else {
756             ans = sign(1, x);
757         }
758         return ans;
759     }
760 
761 
762     /***
763      *  Returns the complementary error function of a double.
764      *  @param  x   A double value.
765      *  @return  The complementary error function of x.
766      */
767     static public double erfc(double x)
768     {
769         double  ans;
770         double  y = Math.abs(x);
771 
772         if (x <= -6.013687357) {
773             // -6.013687357 = -Math.sqrt(-Math.log(1.77245385090551602729816748334 * EPSILON_SMALL))
774             ans = 2;
775         } else if (y < 1.49012e-08) {
776             // 1.49012e-08 = Math.sqrt(2*EPSILON_SMALL)
777             ans = 1 - 2*x/1.77245385090551602729816748334;
778         } else {
779             double ysq = y*y;
780             if (y < 1) {
781                 ans = 1 - x*(1+csevl(2*ysq-1,ERFC_COEF));
782             } else if (y <= 4.0) {
783                 ans = Math.exp(-ysq)/y*(0.5+csevl((8.0/ysq-5.0)/3.0,ERFC2_COEF));
784                 if (x < 0)  ans = 2.0 - ans;if (x < 0)  ans = 2.0 - ans;
785                 if (x < 0)  ans = 2.0 - ans;
786             } else {
787                 ans = Math.exp(-ysq)/y*(0.5+csevl(8.0/ysq-1,ERFCC_COEF));
788                 if (x < 0)  ans = 2.0 - ans;
789             }
790         }
791         return ans;
792     }
793 
794 }