Repository: josdejong/mathjs (js)

Issue:
Implement erf function
I have tried my best to find erf function, but failed. So， is there one ?


Pull request description:
Implementation of erf()
Includes an implementation of the `erf()` error function that many other languages use, including MATLAB, R, and Python. It's basically a translation of [W. J. Cody's Fortran77 implementation](http://netlib.org/specfun/erf) of the approximation that he outlines in his paper, _[Rational Chebyshev Approximations for the Error Function](http://www.ams.org/journals/mcom/1969-23-107/S0025-5718-1969-0247736-4/S0025-5718-1969-0247736-4.pdf)_.

I put `erf.js` in a new submodule, `function/special/`, because it didn't really seem to fit into any of the other submodules.

The spec tests the accuracy of the function w.r.t. values obtained from python 3's `math.erf()`.

Closes #678


Patch:
diff --git a/lib/function/index.js b/lib/function/index.js
index 236888b5b..fb8173326 100644
--- a/lib/function/index.js
+++ b/lib/function/index.js
@@ -9,6 +9,7 @@ module.exports = [
   require('./matrix'),
   require('./probability'),
   require('./relational'),
+  require('./special'),
   require('./statistics'),
   require('./string'),
   require('./trigonometry'),
diff --git a/lib/function/special/erf.js b/lib/function/special/erf.js
new file mode 100644
index 000000000..37033ab5a
--- /dev/null
+++ b/lib/function/special/erf.js
@@ -0,0 +1,200 @@
+'use strict';
+
+var deepMap = require('../../utils/collection/deepMap');
+var sign = require('../../utils/number').sign;
+
+/**
+ * Upper bound for the first approximation interval, 0 <= x <= THRESH
+ * @constant
+ */
+var THRESH = 0.46875;
+
+/**
+ * Constant used by W. J. Cody's Fortran77 implementation to denote sqrt(pi)
+ * @constant
+ */
+var SQRPI = 5.6418958354775628695e-1;
+
+/**
+ * Coefficients for each term of the numerator sum (p_j) for each approximation
+ * interval (see W. J. Cody's paper for more details)
+ * @constant
+ */
+var P = [[
+  3.16112374387056560e00, 1.13864154151050156e02,
+  3.77485237685302021e02, 3.20937758913846947e03,
+  1.85777706184603153e-1
+], [
+  5.64188496988670089e-1, 8.88314979438837594e00,
+  6.61191906371416295e01, 2.98635138197400131e02,
+  8.81952221241769090e02, 1.71204761263407058e03,
+  2.05107837782607147e03, 1.23033935479799725e03,
+  2.15311535474403846e-8
+], [
+  3.05326634961232344e-1, 3.60344899949804439e-1,
+  1.25781726111229246e-1, 1.60837851487422766e-2,
+  6.58749161529837803e-4, 1.63153871373020978e-2
+]];
+
+/**
+ * Coefficients for each term of the denominator sum (q_j) for each approximation
+ * interval (see W. J. Cody's paper for more details)
+ * @constant
+ */
+var Q = [[
+  2.36012909523441209e01, 2.44024637934444173e02,
+  1.28261652607737228e03, 2.84423683343917062e03
+], [
+  1.57449261107098347e01, 1.17693950891312499e02,
+  5.37181101862009858e02, 1.62138957456669019e03,
+  3.29079923573345963e03, 4.36261909014324716e03,
+  3.43936767414372164e03, 1.23033935480374942e03
+], [
+  2.56852019228982242e00, 1.87295284992346047e00,
+  5.27905102951428412e-1, 6.05183413124413191e-2,
+  2.33520497626869185e-3
+]];
+
+/**
+ * Maximum/minimum safe numbers to input to erf() (in ES6+, this number is
+ * Number.[MAX|MIN]_SAFE_INTEGER). erf() for all numbers beyond this limit will
+ * return 1
+ */
+var MAX_NUM = Math.pow(2, 53);
+
+
+function factory (type, config, load, typed) {
+  var multiply  = load(require('../arithmetic/multiply'));
+  var pow       = load(require('../arithmetic/pow'));
+  var sum       = load(require('../statistics/sum'));
+
+  /**
+   * Compute the erf function of a value using a rational Chebyshev
+   * approximations for different intervals of x.
+   *
+   * This is a translation of W. J. Cody's Fortran implementation from 1987
+   * ( http://www.netlib.org/specfun/erf ). See the AMS publication
+   * "Rational Chebyshev Approximations for the Error Function" by W. J. Cody
+   * for an explanation of this process.
+   *
+   * For matrices, the function is evaluated element wise.
+   *
+   * Syntax:
+   *
+   *    math.erf(x)
+   *
+   * Examples:
+   *
+   *    math.erf(0.2);    // returns 0.22270258921047847
+   *    math.erf(-0.5);   // returns -0.5204998778130465
+   *    math.erf(4);      // returns 0.9999999845827421
+   *
+   * @param {number | Array | Matrix} x   A real number
+   * @return {number | Array | Matrix}    The erf of `x`
+   */
+  var erf = typed('erf', {
+    'number': function (x) {
+      var y = Math.abs(x);
+
+      if (y >= MAX_NUM) {
+        return sign(x);
+      }
+      if (y <= THRESH) {
+        return sign(x) * erf1(y);
+      }
+      if (y <= 4.0) {
+        return sign(x) * (1 - erfc2(y));
+      }
+      return sign(x) * (1 - erfc3(y));
+    },
+
+    // TODO: Not sure if there's a way to guarantee some degree of accuracy here.
+    //  Perhaps it would be best to set the precision of the number to that which
+    //  is guaranteed by erf()
+    'BigNumber': function (n) {
+      return new type.BigNumber(erf(n.toNumber()));
+    },
+
+    'Array | Matrix': function (n) {
+      return deepMap(n, erf);
+    }
+
+    // TODO: For complex numbers, use the approximation for the Faddeeva function
+    //  from "More Efficient Computation of the Complex Error Function" (AMS)
+
+  });
+
+  /**
+   * Approximates the error function erf() for x <= 0.46875 using this function:
+   *               n
+   * erf(x) = x * sum (p_j * x^(2j)) / (q_j * x^(2j))
+   *              j=0
+   */
+  function erf1(y) {
+    var ysq = y * y;
+    var xnum = P[0][4]*ysq;
+    var xden = ysq;
+    var i;
+
+    for (i = 0; i < 3; i += 1) {
+      xnum = (xnum + P[0][i]) * ysq;
+      xden = (xden + Q[0][i]) * ysq;
+    }
+    return y * (xnum + P[0][3]) / (xden + Q[0][3]);
+  }
+
+  /**
+   * Approximates the complement of the error function erfc() for
+   * 0.46875 <= x <= 4.0 using this function:
+   *                       n
+   * erfc(x) = e^(-x^2) * sum (p_j * x^j) / (q_j * x^j)
+   *                      j=0
+   */
+  function erfc2(y) {
+    var xnum = P[1][8] * y;
+    var xden = y;
+    var i;
+
+    for (i = 0; i < 7; i += 1) {
+      xnum = (xnum + P[1][i]) * y;
+      xden = (xden + Q[1][i]) * y;
+    }
+    var result = (xnum + P[1][7]) / (xden + Q[1][7]);
+    var ysq = parseInt(y * 16) / 16;
+    var del = (y - ysq) * (y + ysq);
+    return Math.exp(-ysq*ysq) * Math.exp(-del) * result;
+  }
+
+  /**
+   * Approximates the complement of the error function erfc() for x > 4.0 using
+   * this function:
+   *
+   * erfc(x) = (e^(-x^2) / x) * [ 1/sqrt(pi) +
+   *               n
+   *    1/(x^2) * sum (p_j * x^(-2j)) / (q_j * x^(-2j)) ]
+   *              j=0
+   */
+  function erfc3(y) {
+    var ysq = 1 / (y * y);
+    var xnum = P[2][5] * ysq;
+    var xden = ysq;
+    var i;
+
+    for (i = 0; i < 4; i += 1) {
+      xnum = (xnum + P[2][i]) * ysq;
+      xden = (xden + Q[2][i]) * ysq;
+    }
+    var result = ysq * (xnum + P[2][4]) / (xden + Q[2][4]);
+    result = (SQRPI - result) / y;
+    ysq = parseInt(y * 16) / 16;
+    var del = (y - ysq) * (y + ysq);
+    return Math.exp(-ysq*ysq) * Math.exp(-del) * result;
+  }
+
+  erf.toTex = {1: 'erf\\left(${args[0]}\\right)'};
+
+  return erf;
+}
+
+exports.name = 'erf';
+exports.factory = factory;
diff --git a/lib/function/special/index.js b/lib/function/special/index.js
new file mode 100644
index 000000000..1fd1b35e2
--- /dev/null
+++ b/lib/function/special/index.js
@@ -0,0 +1,3 @@
+module.exports = [
+  require('./erf')
+];


Test patch:
diff --git a/test/function/special/erf.test.js b/test/function/special/erf.test.js
new file mode 100644
index 000000000..74b038526
--- /dev/null
+++ b/test/function/special/erf.test.js
@@ -0,0 +1,158 @@
+var assert          = require('assert')
+  , approx          = require('../../../tools/approx')
+  , error           = require('../../../lib/error/index')
+  , math            = require('../../../index')
+  , bigUtil         = require('../../../lib/utils/index').bignumber
+  , bignumber       = math.bignumber
+  , subtract        = math.subtract
+  , abs             = math.abs
+  , smaller         = math.smaller
+  , erf             = math.erf
+  , actualErfValues = require('./erf.values.json');
+
+var DIFF_THRESH = 5e-16;
+
+describe('erf', function () {
+
+  it('should calculate erf(x), |x| < 0.46875', function () {
+    var i, diff, actual, expected;
+    for (i = -4; i <= 4; i += 1) {
+      actual = erf(i / 10);
+      expected = actualErfValues[(i / 10).toFixed(1)];
+      diff = Math.abs(actual - expected);
+      assert.ok(diff < DIFF_THRESH, [i, actual, expected, diff]);
+    }
+  });
+
+  it('should calculate erf(x), 0.46875 <= |x| <= 4', function () {
+    var i, diff, actual, expected;
+    for (i = -40; i < -4; i += 1) {
+      actual = erf(i / 10);
+      expected = actualErfValues[(i / 10).toFixed(1)];
+      diff = Math.abs(actual - expected);
+      assert.ok(diff < DIFF_THRESH, [i, actual, expected, diff]);
+    }
+    for (i = 5; i <= 40; i += 1) {
+      actual = erf(i / 10);
+      expected = actualErfValues[(i / 10).toFixed(1)];
+      diff = Math.abs(actual - expected);
+      assert.ok(diff < DIFF_THRESH, [i, actual, expected, diff]);
+    }
+  });
+
+  it('should calculate erf(x), |x| > 4', function () {
+    var i, diff, actual, expected;
+    for (i = -70; i < -40; i += 1) {
+      actual = erf(i / 10);
+      expected = actualErfValues[(i / 10).toFixed(1)];
+      diff = Math.abs(actual - expected);
+      assert.ok(diff < DIFF_THRESH, [i, actual, expected, diff]);
+    }
+    for (i = 41; i < 70; i += 1) {
+      actual = erf(i / 10);
+      expected = actualErfValues[(i / 10).toFixed(1)];
+      diff = Math.abs(actual - expected);
+      assert.ok(diff < DIFF_THRESH, [i, actual, expected, diff]);
+    }
+  });
+
+  it('should calculate the erf of a nonpositive integer', function () {
+    assert.ok(Math.abs(erf(-1) - actualErfValues['-1.0']) < DIFF_THRESH);
+    assert.ok(Math.abs(erf(-2) - actualErfValues['-2.0']) < DIFF_THRESH);
+    assert.ok(Math.abs(erf(-3) - actualErfValues['-3.0']) < DIFF_THRESH);
+  });
+
+  it('should calculate the erf of a rational bignumber', function () {
+    assert.ok(smaller(
+      abs(subtract(erf(bignumber(0.1)), bignumber(actualErfValues['0.1']))),
+      DIFF_THRESH
+    ));
+    assert.ok(smaller(
+      abs(subtract(erf(bignumber(0.2)), bignumber(actualErfValues['0.2']))),
+      DIFF_THRESH
+    ));
+    assert.ok(smaller(
+      abs(subtract(erf(bignumber(0.3)), bignumber(actualErfValues['0.3']))),
+      DIFF_THRESH
+    ));
+    assert.ok(smaller(
+      abs(subtract(erf(bignumber(1.5)), bignumber(actualErfValues['1.5']))),
+      DIFF_THRESH
+    ));
+    assert.ok(smaller(
+      abs(subtract(erf(bignumber(2.5)), bignumber(actualErfValues['2.5']))),
+      DIFF_THRESH
+    ));
+
+    var bigmath = math.create({ precision: 15 });
+    assert.ok(smaller(
+      abs(subtract(bigmath.erf(bignumber(-1.5)), bigmath.bignumber(actualErfValues['-1.5']))),
+      DIFF_THRESH
+    ));
+    bigmath.config({ precision: 13 });
+    assert.ok(smaller(
+      abs(subtract(bigmath.erf(bignumber(-1.5)), bigmath.bignumber(actualErfValues['-1.5']))),
+      DIFF_THRESH
+    ));
+  });
+
+  it('should calculate the erf of a boolean (true = 1, false = 0)', function () {
+    assert.ok(Math.abs(erf(true) - actualErfValues['1.0']) < DIFF_THRESH);
+    assert.ok(Math.abs(erf(false) - actualErfValues['0.0']) < DIFF_THRESH);
+  });
+
+  it('should calculate the erf of null (0)', function () {
+    assert.ok(Math.abs(erf(null) - actualErfValues['0.0']) < DIFF_THRESH);
+  });
+
+  it('should calculate the erf of each element in a matrix', function () {
+    math.subtract(
+      erf(math.matrix([0,1,2,3,4,5])),
+      math.matrix([0,1,2,3,4,5].map(function (x) {
+        return actualErfValues[x.toFixed(1)];
+      }))
+    ).forEach(function (diff) {
+      assert.ok(diff < DIFF_THRESH);
+    });
+  });
+
+  it('should calculate the erf of each element in an array', function () {
+    math.subtract(
+      erf(math.matrix([0,1,2,3,4,5])),
+      math.matrix([0,1,2,3,4,5].map(function (x) {
+        return actualErfValues[x.toFixed(1)];
+      }))
+    ).forEach(function (diff) {
+      assert.ok(diff < DIFF_THRESH);
+    });
+  });
+
+  it('should throw en error if called with invalid number of arguments', function() {
+    assert.throws(function () { erf(); });
+    assert.throws(function () { erf(1,3); });
+  });
+
+  it('should throw en error if called with invalid type of argument', function() {
+    assert.throws(function () { erf(new Date()); });
+    assert.throws(function () { erf('a string'); });
+  });
+
+  it('should LaTeX erf', function () {
+    var expression = math.parse('erf(2.5)');
+    assert.equal(expression.toTex(), 'erf\\left(2.5\\right)');
+  });
+
+  it('should return 1 for numbers greater than 2**53 (including Infinity)', function () {
+    assert.equal(erf(Math.pow(2, 53)), 1);
+    assert.equal(erf(Infinity), 1);
+  });
+
+  it('should return -1 for numbers less than -2**53 (including -Infinity)', function () {
+    assert.equal(erf(-Math.pow(2, 53)), -1);
+    assert.equal(erf(-Infinity), -1);
+  });
+
+  // TODO: Test with nums really close to 0
+  // TODO: Once this works for complex numbers and imaginary units, test
+
+});
diff --git a/test/function/special/erf.values.json b/test/function/special/erf.values.json
new file mode 100644
index 000000000..3eff57eac
--- /dev/null
+++ b/test/function/special/erf.values.json
@@ -0,0 +1,46 @@
+{
+  "-7.0": -1.0, "-6.9": -1.0, "-6.8": -1.0, "-6.7": -1.0, "-6.6": -1.0, "-6.5": -1.0,
+  "-6.4": -1.0, "-6.3": -1.0, "-6.2": -1.0, "-6.1": -1.0, "-6.0": -1.0,
+  "-5.9": -0.9999999999999999,  "-5.8": -0.9999999999999998,  "-5.7": -0.9999999999999992,
+  "-5.6": -0.9999999999999977,  "-5.5": -0.9999999999999927,  "-5.4": -0.9999999999999777,
+  "-5.3": -0.9999999999999338,  "-5.2": -0.9999999999998075,  "-5.1": -0.9999999999994507,
+  "-5.0": -0.9999999999984626,  "-4.9": -0.999999999995781,   "-4.8": -0.9999999999886479,
+  "-4.7": -0.9999999999700474,  "-4.6": -0.999999999922504,   "-4.5": -0.9999999998033839,
+  "-4.4": -0.999999999510829,   "-4.3": -0.9999999988065282,  "-4.2": -0.9999999971445058,
+  "-4.1": -0.9999999932999724,  "-4.0": -0.9999999845827421,  "-3.9": -0.9999999652077514,
+  "-3.8": -0.9999999229960725,  "-3.7": -0.9999998328489421,  "-3.6": -0.999999644137007,
+  "-3.5": -0.9999992569016276,  "-3.4": -0.9999984780066371,  "-3.3": -0.9999969422902035,
+  "-3.2": -0.9999939742388483,  "-3.1": -0.9999883513426328,  "-3.0": -0.9999779095030014,
+  "-2.9": -0.9999589021219005,  "-2.8": -0.9999249868053346,  "-2.7": -0.9998656672600594,
+  "-2.6": -0.9997639655834707,  "-2.5": -0.999593047982555,   "-2.4": -0.999311486103355,
+  "-2.3": -0.9988568234026434,  "-2.2": -0.9981371537020182,  "-2.1": -0.997020533343667,
+  "-2.0": -0.9953222650189527,  "-1.9": -0.9927904292352574,  "-1.8": -0.9890905016357308,
+  "-1.7": -0.9837904585907746,  "-1.6": -0.976348383344644,   "-1.5": -0.9661051464753108,
+  "-1.4": -0.9522851197626488,  "-1.3": -0.9340079449406524,  "-1.2": -0.9103139782296352,
+  "-1.1": -0.8802050695740816,  "-1.0": -0.842700792949715,   "-0.9": -0.796908212422832,
+  "-0.8": -0.7421009647076606,  "-0.7": -0.6778011938374184,  "-0.6": -0.6038560908479259,
+  "-0.5": -0.5204998778130465,  "-0.4": -0.42839235504666845, "-0.3": -0.3286267594591274,
+  "-0.2": -0.22270258921047842, "-0.1": -0.11246291601828488,  "0.0": 0.0,
+   "0.1": 0.11246291601828488,   "0.2": 0.22270258921047842,   "0.3": 0.3286267594591274,
+   "0.4": 0.42839235504666845,   "0.5": 0.5204998778130465,    "0.6": 0.6038560908479259,
+   "0.7": 0.6778011938374184,    "0.8": 0.7421009647076606,    "0.9": 0.796908212422832,
+   "1.0": 0.842700792949715,     "1.1": 0.8802050695740816,    "1.2": 0.9103139782296352,
+   "1.3": 0.9340079449406524,    "1.4": 0.9522851197626488,    "1.5": 0.9661051464753108,
+   "1.6": 0.976348383344644,     "1.7": 0.9837904585907746,    "1.8": 0.9890905016357308,
+   "1.9": 0.9927904292352574,    "2.0": 0.9953222650189527,    "2.1": 0.997020533343667,
+   "2.2": 0.9981371537020182,    "2.3": 0.9988568234026434,    "2.4": 0.999311486103355,
+   "2.5": 0.999593047982555,     "2.6": 0.9997639655834707,    "2.7": 0.9998656672600594,
+   "2.8": 0.9999249868053346,    "2.9": 0.9999589021219005,    "3.0": 0.9999779095030014,
+   "3.1": 0.9999883513426328,    "3.2": 0.9999939742388483,    "3.3": 0.9999969422902035,
+   "3.4": 0.9999984780066371,    "3.5": 0.9999992569016276,    "3.6": 0.999999644137007,
+   "3.7": 0.9999998328489421,    "3.8": 0.9999999229960725,    "3.9": 0.9999999652077514,
+   "4.0": 0.9999999845827421,    "4.1": 0.9999999932999724,    "4.2": 0.9999999971445058,
+   "4.3": 0.9999999988065282,    "4.4": 0.999999999510829,     "4.5": 0.9999999998033839,
+   "4.6": 0.999999999922504,     "4.7": 0.9999999999700474,    "4.8": 0.9999999999886479,
+   "4.9": 0.999999999995781,     "5.0": 0.9999999999984626,    "5.1": 0.9999999999994507,
+   "5.2": 0.9999999999998075,    "5.3": 0.9999999999999338,    "5.4": 0.9999999999999777,
+   "5.5": 0.9999999999999927,    "5.6": 0.9999999999999977,    "5.7": 0.9999999999999992,
+   "5.8": 0.9999999999999998,    "5.9": 0.9999999999999999,    "6.0": 1.0, "6.1": 1.0,
+   "6.2": 1.0, "6.3": 1.0, "6.4": 1.0, "6.5": 1.0, "6.6": 1.0, "6.7": 1.0, "6.8": 1.0,
+   "6.9": 1.0
+}