Code formatting; Added credits for erf/erfc implementation

1 files changed, 41 insertions, 37 deletions

diff --git a/src/math.c b/src/math.c
index 56cde6bf5..c6f412055 100644
--- a/src/math.c
+++ b/src/math.c
@@ -17,55 +17,59 @@
 #define atanh(x) (log(1+x) - log(1-x))/2.0
 #define cbrt(x)  pow(x,1.0/3.0)
 
+/* 
+** Implementations of error functions
+** credits to http://www.digitalmars.com/archives/cplusplus/3634.html
+*/
 #define REL_ERROR 1E-12
 /* Implementation of Error function */
 double 
 erf(double x)
 {
   static const double two_sqrtpi=  1.128379167095512574;
-    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);
-        ++j;
-        term*= xsqr/j;
-        sum+= term/(2*j+1);
-        ++j;
-    } while (fabs(term)/sum > REL_ERROR);
-    return two_sqrtpi*sum;
+  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);
+    ++j;
+    term*= xsqr/j;
+    sum+= term/(2*j+1);
+    ++j;
+  } while (fabs(term)/sum > REL_ERROR);
+  return two_sqrtpi*sum;
 }
 
 /* Implementation of complementary Error function */
 double 
 erfc(double x)
 {
-    static const double one_sqrtpi=  0.564189583547756287;        
-    if (fabs(x) < 2.2) {
-        return 1.0 - erf(x);       
-    }                              
-    if (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;
-      } while (fabs(q1-q2)/q2 > REL_ERROR);
-    return one_sqrtpi*exp(-x*x)*q2;
+  static const double one_sqrtpi=  0.564189583547756287;        
+  if (fabs(x) < 2.2) {
+    return 1.0 - erf(x);       
+  }                              
+  if (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;
+  } while (fabs(q1-q2)/q2 > REL_ERROR);
+  return one_sqrtpi*exp(-x*x)*q2;
 }
 
 #endif