Merge pull request #724 from xxuejie/add_erf

Add erf and erfc implementations

3 files changed, 77 insertions, 0 deletions

diff --git a/AUTHORS b/AUTHORS
index 554d7bfd..a1e995c5 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -40,4 +40,5 @@ a license to everyone to use it as detailed in LICENSE.)
 * Lars Schneider <lars.schneider@autodesk.com> (copyright owned by Autodesk, Inc.)
 * Joel Martin <github@martintribe.org>
 * Manuel Wellmann <manuel.wellmann@autodesk.com> (copyright owned by Autodesk, Inc.)
+* Xuejie Xiao <xxuejie@gmail.com>
 
diff --git a/src/library.js b/src/library.js
index 083ec0f1..837ed71d 100644
--- a/src/library.js
+++ b/src/library.js
@@ -5131,6 +5131,64 @@ LibraryManager.library = {
   atan2f: 'Math.atan2',
   exp: 'Math.exp',
   expf: 'Math.exp',
+
+  // The erf and erfc functions are inspired from
+  // http://www.digitalmars.com/archives/cplusplus/3634.html
+  // and mruby source code at
+  // https://github.com/mruby/mruby/blob/master/src/math.c
+  erfc: function (x) {
+    var MATH_TOLERANCE = 1E-12;
+    var ONE_SQRTPI = 0.564189583547756287;
+    var a = 1;
+    var b = x;
+    var c = x;
+    var d = x * x + 0.5;
+    var n = 1.0;
+    var q2 = b / d;
+    var q1, t;
+
+    if (Math.abs(x) < 2.2) {
+      return 1.0 - _erf(x);
+    }
+    if (x < 0) {
+      return 2.0 - _erfc(-x);
+    }
+    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 (Math.abs(q1 - q2) / q2 > MATH_TOLERANCE);
+    return (ONE_SQRTPI * Math.exp(- x * x) * q2);
+  },
+  erfcf: 'erfcf',
+  erf: function (x) {
+    var MATH_TOLERANCE = 1E-12;
+    var TWO_SQRTPI = 1.128379167095512574;
+    var sum = x;
+    var term = x;
+    var xsqr = x*x;
+    var j = 1;
+
+    if (Math.abs(x) > 2.2) {
+      return 1.0 - _erfc(x);
+    }
+    do {
+      term *= xsqr / j;
+      sum -= term / (2 * j + 1);
+      ++j;
+      term *= xsqr / j;
+      sum += term / (2 * j + 1);
+      ++j;
+    } while (Math.abs(term / sum) > MATH_TOLERANCE);
+    return (TWO_SQRTPI * sum);
+  },
+  erff: 'erf',
   log: 'Math.log',
   logf: 'Math.log',
   sqrt: 'Math.sqrt',
diff --git a/tests/runner.py b/tests/runner.py
index ba5ec148..6404a211 100755
--- a/tests/runner.py
+++ b/tests/runner.py
@@ -1450,6 +1450,24 @@ c5,de,15,8a
         '''
         self.do_run(src, '*3.14,-3.14,1,0,0,0,1,0,1,1,0,2,3,0.0,1.0,0.0,1.0*')
 
+    def test_erf(self):
+        src = '''
+          #include <math.h>
+          #include <stdio.h>
+          int main()
+          {
+            printf("%1.6f, %1.6f, %1.6f, %1.6f, %1.6f, %1.6f\\n",
+                   erf(1.0),
+                   erf(3.0),
+                   erf(-1.0),
+                   erfc(1.0),
+                   erfc(3.0),
+                   erfc(-1.5));
+            return 0;
+          }
+        '''
+        self.do_run(src, '0.842701, 0.999978, -0.842701, 0.157299, 0.000022, 1.966105')
+
     def test_math_hyperbolic(self):
         src = open(path_from_root('tests', 'hyperbolic', 'src.c'), 'r').read()
         expected = open(path_from_root('tests', 'hyperbolic', 'output.txt'), 'r').read()