home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / C / JPLC2.ZIP / ERF.C < prev    next >
Encoding:
C/C++ Source or Header  |  1988-07-22  |  3.6 KB  |  117 lines
  1. /* 1.0  04-27-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************/
  7.  
  8. #include "defs.h"
  9. #include "stdtyp.h"
  10. #include "errno.h"
  11. #include "mathcons.h"
  12.  
  13. /*   Coefficients from J. L. Schonfelder, Math. Comp. Vol. 32, No. 144,
  14.  *   Oct. 1978, pp. 1232-1240.  Converted from Chebyshev coefficients
  15.  *   to polynomial coefficients.
  16.  */
  17.  
  18. #define bDEGREE 19
  19. static double b[bDEGREE+1] =
  20. {    /* b[0] */     6.74933236039655050263e-01,
  21.     /* b[1] */    -2.61111860931245377287e-01,
  22.     /* b[2] */     1.19479138609851943329e-01,
  23.     /* b[3] */    -4.86627774491551332176e-02,
  24.     /* b[4] */     1.71283445718166884035e-02,
  25.     /* b[5] */    -5.23487583615739131346e-03,
  26.     /* b[6] */     1.40509142365203239024e-03,
  27.     /* b[7] */    -3.35143533535934984684e-04,
  28.     /* b[8] */     7.18010084384964546189e-05,
  29.     /* b[9] */    -1.39462736935665598139e-05,
  30.     /* b[10] */     2.47580204969781916589e-06,
  31.     /* b[11] */    -4.04509593030131189153e-07,
  32.     /* b[12] */     6.11956771414692047983e-08,
  33.     /* b[13] */    -8.61749835843079723418e-09,
  34.     /* b[14] */     1.13482488917751004919e-09,
  35.     /* b[15] */    -1.40320352244567126036e-10,
  36.     /* b[16] */     1.63314151643890654668e-11,
  37.     /* b[17] */    -1.79944661052623065189e-12,
  38.     /* b[18] */     1.97235204103604890406e-13,
  39.     /* b[19] */    -1.95488364156685071066e-14
  40. };
  41.  
  42. #define cDEGREE 21
  43. static double c[cDEGREE+1] = 
  44. {    /* c[0] */     5.40464821348814759403e-01,
  45.     /* c[1] */    -2.61515522487415932119e-02,
  46.     /* c[2] */    -2.88573438386340392753e-03,
  47.     /* c[3] */    -5.29353396945768734440e-04,
  48.     /* c[4] */    -8.54344544297635788098e-05,
  49.     /* c[5] */    -1.74905253140942193568e-05,
  50.     /* c[6] */    -3.06142196558938128874e-06,
  51.     /* c[7] */    -6.87177932914162287489e-07,
  52.     /* c[8] */    -1.18400244002269115299e-07,
  53.     /* c[9] */    -3.04457664850508328527e-08,
  54.     /* c[10] */    -4.49026736274657305330e-09,
  55.     /* c[11] */    -1.55646576715015340596e-09,
  56.     /* c[12] */    -1.22012057640254963189e-10,
  57.     /* c[13] */    -1.01539218332404736429e-10,
  58.     /* c[14] */     4.97278086303190793842e-12,
  59.     /* c[15] */    -8.76177272209771675989e-12,
  60.     /* c[16] */     3.33713910876554204151e-12,
  61.     /* c[17] */    -1.94148792669164948165e-12,
  62.     /* c[18] */    -2.99117077089093858376e-13,
  63.     /* c[19] */     2.10106516619397560135e-13,
  64.     /* c[20] */     3.40434284957405645400e-13,
  65.     /* c[21] */    -1.73433792163219535723e-13,
  66. };
  67.  
  68. /************************************************************************/
  69.     double
  70. erf(x)        /* return the error function of x, as defined in
  71.            Handbook of Mathematical Functions, National Bureau of
  72.            Standards AMS-55, pp. 297.                */
  73. /*----------------------------------------------------------------------*/
  74. double x;
  75. {
  76.     int minus, k;
  77.     double val, t, erfc();
  78.  
  79.     minus = (x < 0 ? 1 : 0);
  80.     t = (minus ? -x : x);
  81.     if (t > 2.0)
  82.     {    t = 1.0 - erfc(t);
  83.         return (minus ? -t : t);
  84.     }
  85.     else
  86.     {    t = x * x / 2.0 - 1.0;
  87.         val = b[bDEGREE];
  88.         for (k = bDEGREE-1; k >= 0; k--)
  89.             val = val * t + b[k];
  90.         return (x * val);
  91.     }
  92. }
  93.  
  94. /*\p*********************************************************************/
  95.     double
  96. erfc(x)            /* complementary error function, 1-erf(x)    */
  97.  
  98. /*----------------------------------------------------------------------*/
  99. double x;
  100. {
  101.     int k, minus, chrstc();
  102.     double t, val, y, z, erf(), exp();
  103.  
  104.     minus = (x < 0.0 ? 1 : 0);
  105.     z = (minus ? -x : x);
  106.     if (z < 2.0)
  107.         return (1.0 - erf(x));
  108.  
  109.     y = x * x;
  110.     t = (10.5 - y)/(2.5 + y);
  111.     val = c[cDEGREE];
  112.     for (k = cDEGREE-1; k >= 0; k--)
  113.         val = val * t + c[k];
  114.     val = exp(-y) * val / z;
  115.     return (minus ? 2.0 - val : val);
  116. }
  117.