erfc.cpp

   1 //-----------------------------------------------------------------------------
   2 //
   3 //      Author : philippe.billet@noos.fr
   4 //
   5 //      erfc(x)
   6 //
   7 //      GW      Added erfc() from Numerical Recipes in C
   8 //
   9 //-----------------------------------------------------------------------------
  10
  11 #include "stdafx.h"
  12 #include "defs.h"
  13 static void yyerfc(void);
  14
  15 void
  16 eval_erfc(void)
  17 {
  18         push(cadr(p1));
  19         eval();
  20         yerfc();
  21 }
  22
  23 void
  24 yerfc(void)
  25 {
  26         save();
  27         yyerfc();
  28         restore();
  29 }
  30
  31 static void
  32 yyerfc(void)
  33 {
  34         double d;
  35
  36         p1 = pop();
  37
  38         if (isdouble(p1)) {
  39                 d = erfc(p1->u.d);
  40                 push_double(d);
  41                 return;
  42         }
  43
  44         push_symbol(ERFC);
  45         push(p1);
  46         list(2);
  47         return;
  48 }
  49
  50 // from Numerical Recipes in C
  51
  52 #ifndef LINUX
  53 double
  54 erfc(double x)
  55 {
  56         double t, z, ans;
  57         z = fabs(x);
  58         t = 1.0 / (1.0 + 0.5 * z);
  59
  60         ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
  61         t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
  62         t*(-0.82215223+t*0.17087277)))))))));
  63
  64         return x >= 0.0 ? ans : 2.0-ans;
  65 }
  66 #endif
  67
  68 #if SELFTEST
  69
  70 static char *s[] = {
  71
  72         "erfc(a)",
  73         "erfc(a)",
  74
  75         "erfc(0.0)",
  76         "1",
  77
  78         "float(erfc(0))",
  79         "1",
  80 #if 0
  81         "float(erfc(1))",
  82         "0.157299",
  83 #endif
  84 };
  85
  86 void
  87 test_erfc(void)
  88 {
  89         test(__FILE__, s, sizeof s / sizeof (char *));
  90 }
  91
  92 #endif