home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / erfl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-26  |  6.8 KB  |  326 lines
  1. /*                            erf.c
  2.  *
  3.  *    Error function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, erf();
  10.  *
  11.  * y = erf( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * The integral is
  18.  *
  19.  *                           x 
  20.  *                            -
  21.  *                 2         | |          2
  22.  *   erf(x)  =  --------     |    exp( - t  ) dt.
  23.  *              sqrt(pi)   | |
  24.  *                          -
  25.  *                           0
  26.  *
  27.  * The magnitude of x is limited to 9.231948545 for DEC
  28.  * arithmetic; 1 or -1 is returned outside this range.
  29.  *
  30.  * For 0 <= |x| < 1, erf(x) = x * P6(x^2)/Q6(x^2); otherwise
  31.  * erf(x) = 1 - erfc(x).
  32.  *
  33.  *
  34.  *
  35.  * ACCURACY:
  36.  *
  37.  *                      Relative error:
  38.  * arithmetic   domain     # trials      peak         rms
  39.  *    IEEE      0,1          3000       1.8e-19     6.5e-20
  40.  *
  41.  */
  42. /*                            erfc.c
  43.  *
  44.  *    Complementary error function
  45.  *
  46.  *
  47.  *
  48.  * SYNOPSIS:
  49.  *
  50.  * double x, y, erfc();
  51.  *
  52.  * y = erfc( x );
  53.  *
  54.  *
  55.  *
  56.  * DESCRIPTION:
  57.  *
  58.  *
  59.  *  1 - erf(x) =
  60.  *
  61.  *                           inf. 
  62.  *                             -
  63.  *                  2         | |          2
  64.  *   erfc(x)  =  --------     |    exp( - t  ) dt
  65.  *               sqrt(pi)   | |
  66.  *                           -
  67.  *                            x
  68.  *
  69.  *
  70.  * For small x, erfc(x) = 1 - erf(x); otherwise rational
  71.  * approximations are computed.
  72.  *
  73.  *
  74.  *
  75.  * ACCURACY:
  76.  *
  77.  *                      Relative error:
  78.  * arithmetic   domain     # trials      peak         rms
  79.  *    IEEE      0,8         1000       1.8e-18     6.1e-19
  80.  * For x > 1, error is dominated by the calculation of exp(-x^2)
  81.  * and is of order 3e-19 x.
  82.  *
  83.  *
  84.  * ERROR MESSAGES:
  85.  *
  86.  *   message         condition              value returned
  87.  * erfc underflow    x > sqrt(MAXLOGL)          0.0
  88.  *
  89.  *
  90.  */
  91.  
  92. /*                            ndtr.c
  93.  *
  94.  *    Normal distribution function
  95.  *
  96.  *
  97.  *
  98.  * SYNOPSIS:
  99.  *
  100.  * double x, y, ndtr();
  101.  *
  102.  * y = ndtr( x );
  103.  *
  104.  *
  105.  *
  106.  * DESCRIPTION:
  107.  *
  108.  * Returns the area under the Gaussian probability density
  109.  * function, integrated from minus infinity to x:
  110.  *
  111.  *                            x
  112.  *                             -
  113.  *                   1        | |          2
  114.  *    ndtr(x)  = ---------    |    exp( - t /2 ) dt
  115.  *               sqrt(2pi)  | |
  116.  *                           -
  117.  *                          -inf.
  118.  *
  119.  *             =  ( 1 + erf(z) ) / 2
  120.  *             =  erfc(z) / 2
  121.  *
  122.  * where z = x/sqrt(2). Computation is via the functions
  123.  * erf and erfc.
  124.  *
  125.  *
  126.  * ACCURACY:
  127.  *
  128.  * See erfl, erfc.
  129.  *
  130.  * ERROR MESSAGES:
  131.  *
  132.  *   message         condition         value returned
  133.  * erfc underflow    x > sqrt(MAXLOGL)      0.0
  134.  *
  135.  */
  136.  
  137. /*
  138. Cephes Math Library Release 2.2:  June, 1992
  139. Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  140. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  141. */
  142.  
  143. #include "mathl.h"
  144. extern long double __polevll ( long double, long double *, int );
  145. extern long double __p1evll ( long double, long double *, int );
  146.  
  147. #define SQRTH 7.07106781186547524401E-1L
  148. #define MAXLOGL 1.1356523406294143949492E4L
  149.  
  150. /*
  151. erfc(x) = exp(-x^2) P(z)/Q(z)
  152. z(x) = 1/x
  153. Relative error
  154. n=9, d=10
  155. Peak error =  7.761600916141239627567679213879910992123E-21
  156. Relative error spread =  1.023351469044268241908780645010281568783E0
  157. */
  158. static long double P[] = {
  159.  1.257754157740948157833312785512722262237E9L,
  160.  2.544170037359209505076053758183547658874E9L,
  161.  2.546959074195068787444430942851029462365E9L,
  162.  1.605230257662887752768505542850382838437E9L,
  163.  6.899137924684927504585242368865475647796E8L,
  164.  2.066152777382002343549858319788089708173E8L,
  165.  4.227432144169649646138430076248102333111E7L,
  166.  5.463192951525719736944609969792412137906E6L,
  167.  3.509752864577976833796176700600960407890E5L,
  168. -1.071020848327048365754654097516174838747E-5L,
  169. };
  170. static long double Q[] = {
  171. /* 1.000000000000000000000000000000000000000E0 */
  172.  1.257754146471137811029618866219716318753E9L,
  173.  3.963393686787501126271033323119385984564E9L,
  174.  5.761415509088700168863599482318704228527E9L,
  175.  5.089047658881201390929483867868625212361E9L,
  176.  3.023468726254039362316680583964735625632E9L,
  177.  1.259993446190802606495833425577455019099E9L,
  178.  3.710577088683661371497272954319350172686E8L,
  179.  7.524032571027728454099705963391636990057E7L,
  180.  9.683257455518460133999913804480458079007E6L,
  181.  6.220874963793103874137444426727337364012E5L
  182. };
  183. /*
  184. erfc(x) = exp(-x^2) z P(z**2)/Q(z**2)
  185. z(x) = 1/x
  186. x >= 8
  187. Relative error
  188. n=4, d=5
  189. Peak error =  1.82e-21
  190. Relative error spread =  3.9e-3
  191. */
  192. static long double R[] = {
  193.  3.621771504143418323727691824727645990561E0L,
  194.  7.176028075353789451273277986715563076155E0L,
  195.  3.446811092743012309113505238337669421091E0L,
  196.  5.541324340261065929875426053951866015250E-1L,
  197.  2.699905677541423267238580904507582955709E-2L
  198. };
  199. static long double S[] = {
  200. /* 1.000000000000000000000000000000000000000E0 */
  201.  1.073074561751358600517229966140297947840E1L,
  202.  1.534256653809200505130664147433856235472E1L,
  203.  6.576473386900615376978244623190582091578E0L,
  204.  1.006101457677419240008726711451145446254E0L,
  205.  4.785458215239962067727484634993830260306E-2L
  206. };
  207. /*
  208. erf(x)  z P(z**2)/Q(z**2)
  209. 0 <= x <= 1
  210. z(x) = x
  211. Relative error
  212. n=6, d=6
  213. Peak error =  7.640688643888757419028899590163681314010E-23
  214. Relative error spread =  9.108314657350825251106158661154911170107E-3
  215. */
  216. static long double T[] = {
  217.  1.097461701666048418572764102524978668857E-1L,
  218.  5.403114568277600506272199159607654730741E0L,
  219.  2.871755797503643801098000314688440010017E2L,
  220.  2.677455035557940040284806618777105723130E3L,
  221.  4.825836895686958085365801683371047641820E4L,
  222.  1.549871261553367266536226453640515782388E5L,
  223.  1.104347959782514112319886330985904796384E6L
  224. };
  225. static long double U[] = {
  226. /* 1.000000000000000000000000000000000000000E0 */
  227.  4.525736349514749653000420739995094197946E1L,
  228.  9.715176200212976593594958594375455394272E2L,
  229.  1.245878656969449881132817817129763145725E4L,
  230.  9.942693068079099139656602031366763195331E4L,
  231.  4.635880633067639826809494193933792409347E5L,
  232.  9.787028970280835391789815008496258590638E5L
  233. };
  234.  
  235. #define UTHRESH 37.519379347
  236.  
  237. /*
  238. double ndtr(a)
  239. double a;
  240. {
  241. double x, y, z;
  242. double fabs(), erf(), erfc();
  243.  
  244. x = a * SQRTH;
  245. z = fabs(x);
  246.  
  247. if( z < SQRTH )
  248.     y = 0.5 + 0.5 * erf(x);
  249.  
  250. else
  251.     {
  252.     y = 0.5 * erfc(z);
  253.  
  254.     if( x > 0 )
  255.         y = 1.0 - y;
  256.     }
  257.  
  258. return(y);
  259. }
  260. */
  261.  
  262. long double erfcl(long double a)
  263. {
  264. long double p,q,x,y,z;
  265.  
  266.  
  267. if( a < 0.0L )
  268.     x = -a;
  269. else
  270.     x = a;
  271.  
  272. if( x < 1.0L )
  273.     return( 1.0L - erfl(a) );
  274.  
  275. z = -a * a;
  276.  
  277. if( z < -MAXLOGL )
  278.     {
  279. under:
  280. /* Some kind of error flagging needed. */
  281. /*    mtherr( "erfcl", UNDERFLOW ); */
  282.     if( a < 0 )
  283.         return( 2.0L );
  284.     else
  285.         return( 0.0L );
  286.     }
  287.  
  288. z = expl(z);
  289.  
  290. if( x < 8.0L )
  291.     {
  292.       y = 1.0/a;
  293.     p = __polevll( y, P, 9 );
  294.     q = __p1evll( y, Q, 10 );
  295.     }
  296. else
  297.     {
  298.     y = 1.0/(a*a);
  299.     p = __polevll( y, R, 4 ) / a;
  300.     q = __p1evll( y, S, 5 );
  301.     }
  302. y = (z * p)/q;
  303.  
  304. if( a < 0 )
  305.     y = 2.0L - y;
  306.  
  307. if( y == 0.0L )
  308.     goto under;
  309.  
  310. return(y);
  311. }
  312.  
  313.  
  314.  
  315. long double erfl(long double x)
  316. {
  317. long double y, z;
  318.  
  319. if( fabsl(x) > 1.0L )
  320.     return( 1.0L - erfcl(x) );
  321.  
  322. z = x * x;
  323. y = x * __polevll( z, T, 6 ) / __p1evll( z, U, 6 );
  324. return( y );
  325. }
  326.