home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / comp / unix / bsd / 8497 < prev    next >
Encoding:
Text File  |  1992-11-06  |  2.7 KB  |  130 lines
  1. Path: sparky!uunet!ukma!cs.widener.edu!dsinc!ub!acsu.buffalo.edu!usenet
  2. From: jones@acsu.buffalo.edu (Terry A. Jones)
  3. Newsgroups: comp.unix.bsd
  4. Subject: Re: Where are the functions erf(), and erfc()?
  5. Message-ID: <Bx9vtF.2Bo@acsu.buffalo.edu>
  6. Date: 6 Nov 92 02:20:02 GMT
  7. References: <Bwz4qG.9p6@acsu.buffalo.edu>
  8. Sender: nntp@acsu.buffalo.edu
  9. Organization: State University of New York at Buffalo/Comp Sci
  10. Lines: 117
  11. Nntp-Posting-Host: hyades.cs.buffalo.edu
  12.  
  13.  
  14.     Well I did not get any feedback on this one so I implemented
  15. the functions myself.  I couldn't wait to get Spice compiled.  Cut into
  16. erf.c, compile and add to your libm.a if you like.
  17.  
  18.     Terry Jones        jones@cs.buffalo.edu
  19.  
  20. ---------------------------cut here-----------------------------------------
  21. #include "mathimpl.h"
  22.  
  23. /*
  24.  *  Predefined terms in expansion series for the error function.
  25.  */
  26. #define T2    0.6666667
  27. #define T3    0.2666667
  28. #define T4    0.07619048
  29. #define T5    0.01693122
  30. #define T6    3.078403e-3
  31. #define T7    4.736005e-4
  32. #define T8    6.314673e-5
  33. #define T9    7.429027e-6
  34. #define T10    7.820028e-7
  35. #define T11    7.447646e-8
  36. #define T12    6.476214e-9
  37.  
  38. /*
  39.  *  Internal version of the error function.
  40.  */
  41. static double errorf( x )
  42.      double x;
  43. {
  44.   double x2;
  45.  
  46.   x2 = x * x;
  47.  
  48.   return( ( 2.0 * exp( -x2 ) / sqrt( M_PI ) ) * ( x * ( 1 + x2 * ( T2 + x2 *
  49.         ( T3 + x2 *( T4 + x2 *( T5 + x2 *( T6 + x2 *( T7 + x2 *( T8 + x2 *
  50.         ( T9 + x2 *( T10 + x2 *( T11 + x2 * T12)))))))))))));
  51. }
  52.  
  53.  
  54. /*
  55.  *  Compute the complement to the error function.
  56.  */
  57. static double errorf_complement( x )
  58.      double x;
  59. {
  60.   double x2;
  61.   double v;
  62.  
  63.   x2 = x * x;
  64.   v = 1.0 / ( 2.0 * x2 );
  65.  
  66.   return( 1.0 / ( exp( x2 ) * x * sqrt( M_PI ) * ( 1 + v / ( 1 + 2 * v /
  67.           ( 1 + 3 * v / ( 1 + 4 * v / ( 1 + 5 * v / ( 1 + 6 * v / 
  68.           ( 1 + 7 * v / ( 1 + 8 * v / ( 1 + 9 * v / ( 1 + 10 * v / 
  69.           ( 1 + 11 * v / ( 1 + 12 * v / ( 1 + 13 * v / ( 1 + 14 * v /
  70.           ( 1 + 15 * v / ( 1 + 16 * v / ( 1 + 17 * v )))))))))))))))))));
  71. }
  72.  
  73. /*
  74.  *  Global entry point for the error function.
  75.  */
  76. double xerf( x )
  77.      double x;
  78. {
  79.   double sign;
  80.  
  81.   if ( 0.0 == x )
  82.     return( 0.0 );
  83.  
  84.   /*
  85.    *  Handle case of negative argument.
  86.    */
  87.   if ( x < 0.0 )
  88.     {
  89.       sign = -1.0;
  90.       x = -x;
  91.     }
  92.  
  93.   else
  94.     sign = 1.0;
  95.  
  96.   if ( x < 1.5 )
  97.     return( errorf( x ) * sign );
  98.   else
  99.     return( ( 1.0 - errorf_complement( x )) * sign );
  100.  
  101. }
  102.  
  103. /*
  104.  * Global entry point for the error function complement.
  105.  */
  106. double xerfc( x )
  107.      double x;
  108. {
  109.  
  110.   if ( 0.0 == x )
  111.     return( 1.0 );
  112.  
  113.   if ( x < 0.0 )
  114.     {
  115.       x = -x;
  116.       if ( x < 1.5 )
  117.     return( 1.0 + errorf( x ) );
  118.       else
  119.     return( 2.0 - errorf_complement( x ) );
  120.     }
  121.  
  122.   else
  123.     {
  124.       if ( x < 1.5 )
  125.     return( 1.0 - errorf( x ) );
  126.       else
  127.     return( errorf_complement( x ) );
  128.     }
  129. }
  130.