Rudimentary numerical analysis tools, written for node.js

numerical.js

var exports = module.exports = {};

/**********************************************************
 * @author Alex Nelson
 * @email [email protected]
 * @date 29 September 2014
 **********************************************************/
 
/**
 * My todo list...
 * @todo add unit tests!
 * @todo make sure Newton's method won't cause problems with the dy/y code
 */

/*********************
 * Utility functions *
 *********************/

function computeEpsilon() {
    if (Number.EPSILON) return Number.EPSILON;
    var k = 1.0;
    while (1.0 !== 1.0+k) {
        k = k*0.5;
    }
    return k;
}

/** Machine epsilon, either given as Number.epsilon or computed from
    scratch
    @constant
    @type {number}
    @default
*/
const machineEpsilon = computeEpsilon();

/** The squareroot of machine epsilon
    @constant
    @type {number}
    @default
*/
const sqrtEpsilon = Math.sqrt(machineEpsilon);

// taken from @link {http://stackoverflow.com/a/17156580/1296973}
function getNumberParts(x) {
    var float = new Float64Array(1),
        bytes = new Uint8Array(float.buffer);

    float[0] = x;

    var sign = bytes[7] >> 7,
        exponent = ((bytes[7] & 0x7f) << 4 | bytes[6] >> 4) - 0x3ff;

    bytes[7] = 0x3f;
    bytes[6] |= 0xf0;

    return {
        sign: sign,
        exponent: exponent,
        mantissa: float[0],
    }
}

var abs = function(x) {
    if (x < 0) return -x;
    return x;
}

/**
 * Uses Knuth 2.2's method of comparing floating point numbers.
 *
 * If (x < y),    return -1
 *    (x > y),    return +1
 *    (x "eq" y), return  0
 */
var floatCompare = function(x, y, eps) {
    if (!eps) eps = sqrtEpsilon;
    var exponent = (getNumberParts(abs(x) > abs(y) ? x : y)).exponent;
                   // (Math/get-exponent (if (> (abs a) (abs b)) a b))
    var delta = parseFloat(eps+'e'+exponent); // (Math/scalb eps exponent)
    var diff = x - y;
    if (diff > delta) return 1;
    if (diff < -delta) return -1;
    return 0;
}

/**
 * @param x, y
 * @returns true if they are "equal", false otherwise
 */
var floatEquals = function(x, y, eps) {
    return (floatCompare(x, y, eps) === 0);
}

/**
 * Horner's method will take an array of coefficients, and return a
 * function which evaluates the polynomial at the given point.
 */
var horner = function horner(coefficients) {
    var n = coefficients.length - 1;

    var fn = function(x) {
        var ret = coefficients[n];
        for(var k = n-1; k > -1; k--) {
            ret = coefficients[k] + x*ret;
        }
        return ret;
    }

    return fn;
}

var sgn = function sgn(x) {
    if (x < 0) return -1;
    if (x > 0) return +1;
    return 0;
}

/**********************************
 * Root finding algorithms        *
 **********************************/
var rootFinding = {};

var bisection = function bisection(fn, a, b) {
    if (sgn(fn(a)) === sgn(fn(b)))
        return undefined;

    var midpoint;
    var n;
    for(n = 0; n < 53; n++) {
        midpoint = (a + b)*0.5;
        if(fn(midpoint)===0.0 || (b - a)*0.5 < tolerance)
            return midpoint;
        else if (sgn(fn(a))===sgn(fn(midpoint)))
            a = midpoint;
        else
            b = midpoint;
    }
    return midpoint;
}

var secantMethod = function secantMethod(fn, a, b) {
    var tmp, x, xPrime, y, yPrime;
    x = a;
    if (b)
        xPrime = b;
    else
        xPrime = a + sqrtEpsilon;
    y = fn(x);
    yPrime = fn(xPrime);
    for(var n = 0; n < 13; n++) {
        if (floatEquals(y, 0.0))
            return x;
        // avoid roundoff errors with this approach
        tmp = (x*yPrime - xPrime*y)/(yPrime - y);
        y = yPrime;
        x = xPrime;
        xPrime = tmp;
        yPrime = fn(xPrime);
    }
    return xPrime;
}

var newtonRootMethod = function(fn, df, initialGuess) {
    var tmp, y, dy, x;
    x = initialGuess;
    for(var n = 0; n < 9; n++) {
        dy = df(x);
        y = fn(x);
        // @todo should check that dy/y isn't going to cause problems
        x = x - (dy/y);
    }
    return x;
}

rootFinding.bisection = bisection;
rootFinding.secant = secantMethod;
rootFinding.newton = newtonRootMethod;

exports.machineEpsilon = machineEpsilon;
exports.sqrtEpsilon = sqrtEpsilon;
exports.floatCompare = floatCompare;
exports.floatEquals = floatEquals;
exports.sgn = sgn;
exports.bisection = bisection;
exports.horner = horner;
exports.root = rootFinding;