Merged MSVC branch

1 files changed, 30 insertions, 21 deletions

diff --git a/src/math.c b/src/math.c
index c6f412055..cbe773c58 100644
--- a/src/math.c
+++ b/src/math.c
@@ -17,6 +17,10 @@
 #define atanh(x) (log(1+x) - log(1-x))/2.0
 #define cbrt(x)  pow(x,1.0/3.0)
 
+/* Declaration of complementary Error function */
+double
+erfc(double x);
+
 /* 
 ** Implementations of error functions
 ** credits to http://www.digitalmars.com/archives/cplusplus/3634.html
@@ -26,18 +30,20 @@
 double 
 erf(double x)
 {
-  static const double two_sqrtpi=  1.128379167095512574;
+  static const double two_sqrtpi =  1.128379167095512574;
+  double sum  = x;
+  double term = x;
+  double xsqr = x*x;
+  int j= 1;
   if (fabs(x) > 2.2) {
     return 1.0 - erfc(x);
   }
-  double sum= x, term= x, xsqr= x*x;
-  int j= 1;
   do {
-    term*= xsqr/j;
-    sum-= term/(2*j+1);
+    term *= xsqr/j;
+    sum  -= term/(2*j+1);
     ++j;
-    term*= xsqr/j;
-    sum+= term/(2*j+1);
+    term *= xsqr/j;
+    sum  += term/(2*j+1);
     ++j;
   } while (fabs(term)/sum > REL_ERROR);
   return two_sqrtpi*sum;
@@ -48,26 +54,29 @@ double
 erfc(double x)
 {
   static const double one_sqrtpi=  0.564189583547756287;        
+  double a = 1; 
+  double b = x;               
+  double c = x; 
+  double d = x*x+0.5;         
+  double q1, q2;                  
+  double n = 1.0;
+  double t;
   if (fabs(x) < 2.2) {
     return 1.0 - erf(x);       
   }                              
-  if (signbit(x)) {              
+  if (x < 0.0) { /*signbit(x)*/              
     return 2.0 - erfc(-x);     
   }                              
-  double a=1, b=x;               
-  double c=x, d=x*x+0.5;         
-  double q1,q2;                  
-  double n= 1.0, t;
   do {
-    t= a*n+b*x;
-    a= b;
-    b= t;
-    t= c*n+d*x;
-    c= d;
-    d= t;
-    n+= 0.5;
-    q1= q2;
-    q2= b/d;
+    t  = a*n+b*x;
+    a  = b;
+    b  = t;
+    t  = c*n+d*x;
+    c  = d;
+    d  = t;
+    n += 0.5;
+    q1 = q2;
+    q2 = b/d;
   } while (fabs(q1-q2)/q2 > REL_ERROR);
   return one_sqrtpi*exp(-x*x)*q2;
 }