jiggzson / nerdamer

a symbolic math expression evaluator for javascript
http://www.nerdamer.com
MIT License
511 stars 81 forks source link

404: precision loss not found #404

Open Happypig375 opened 6 years ago

Happypig375 commented 6 years ago

0. TL;DR

This proposal is my take on supporting arbitrary precise decimals, through a new BigDecimal data type along with library support for its operations. Bonuses included. Also fixes buildFunction().

1. Problem

Nerdamer uses native JavaScript Numbers which restricts Nerdamer's ability to achieve a high decimal precision. This proposal attempts to solve that using a BigDecimal data type which can store arbitrary precise decimals.

2. Objectives

3. Separation of non-mutating and mutating operations

i.e. Use non-mutating operations where possible; only use mutating operations where intermediate results are ignored or in performance-tight loops.

I propose the following in the Parser: a. Rename existing functions add, subtract, multiply, divide, pow and expand to mutatingAdd, mutatingSubtract, mutatingMultiply, mutatingDivide, mutatingPow and mutatingExpand. b. Investigate each line where these functions are used and determine if the mutating operations are really needed. If so, change them to use mutatingAdd etc. c. Insert new non-mutating functions, e.g.

add: function(a, b) {
    return _.mutatingAdd(a.clone(), b.clone());
}

The SAFE flag of Settings and blocks of SAFE code can now be removed.

4. The BigDecimal structure

function BigDecimal(mantissa, exponent) {
    if(!this instanceof BigDecimal) return new BigDecimal(mantissa, exponent);
    this.mantissa = bigInt(mantissa);
    this.exponent = bigInt(exponent);
}

This is the classical structure of every real number data type out there. The mantissa and the exponent store bigInts to achieve arbitrary precision.

The precision of calculations will be specified at Settings.PRECISION, with negative meaning number of significant figures, and number of decimal places otherwise.

Here are some functions to help out inside the arithmetic functions:

//To avoid rounding errors as much as possible
BigDecimal.getSafePrecision = function() { return Settings.PRECISION + (Settings.PRECISION < 0 ? -3 : 3); };
//For reuse in Parser methods
BigDecimal.prototype.clone = function() { return new BigDecimal(this.mantissa, this.exponent); };

BigDecimal.scaleExponent = function(n, scale) {
    if(scale > 0) {
        n.exponent = n.exponent.add(scale);
        var tenToScale = bigInt[10].pow(scale);
        var divmod = n.mantissa.divmod(tenToScale);
        //round 0.5~1 to 1
        n.mantissa = divmod.quotient.add(divmod.remainder.compareAbs(tenToScale.divide(2)) === -1 ? 0 : 1);
    } else {
        n.exponent = n.exponent.add(scale);
        n.mantissa = n.mantissa.multiply(bigInt[10].pow(bigInt(scale).negate()));
    }
};
BigDecimal.scaleToExponent = function(n, exponent) { BigDecimal.scaleExponent(n, n.exponent.subtract(exponent).negate()); };

//aka get mantissa length
BigDecimal.getSigFigs = function(n) {
    //Directly working with the fields of bigInts
    var valueLength = n.mantissa.value.length, length = 0;
    for(var i = 0; i < valueLength - 1; i++) length += 7;
    length += (n.mantissa.value[i] || n.mantissa.value).toString().length;
    return length;
};
//aka set mantissa length
BigDecimal.setSigFigs = function(n, figures) { BigDecimal.scaleExponent(n, BigDecimal.getSigFigs(n) - figures); };

//To be called at the end of all BigDecimal operations (add, subtract etc)
BigDecimal.simplify = function(n) {
    if(n.mantissa.equals(0)) {
        n.exponent = bigInt.zero;
    } else {
        var zeroes = /0+$/.exec(n.mantissa.toString())[0];
        if(!zeroes) return;
        BigDecimal.scaleExponent(n, zeroes.length);
    }
};

//For add, subtract and divide
BigDecimal.equateExponents = function(a, b) {
    if(a.exponent.lt(b.exponent)) BigDecimal.scaleToExponent(b, a.exponent);
    else BigDecimal.scaleToExponent(a, b.exponent);
};

//Equality and comparison
BigDecimal.areCloseEnough = function(n1, n2, p /*precision*/) {
    var a = n1.clone(), b = n2.clone();
    if(p < 0) { //significant figures
        BigDecimal.setSigFigs(a, -p);
        BigDecimal.setSigFigs(b, -p);
    } else { //decimal places
        BigDecimal.scaleToExponent(a, p);
        BigDecimal.scaleToExponent(b, p);
    }
    return a.equals(b);
};

{
    var comparison = ['eq', 'equals', 'neq', 'gt', 'lt', 'geq', 'leq'];
    for(x in comparison)
        BigDecimal.prototype[comparison[x]] = function(n) {
             var a = this.clone(), b = n.clone();
             BigDecimal.equateExponents(a, b);
             return a.mantissa[comparison[x]](b.mantissa);
        };
}

//For the following...
BigDecimal.isNegligible = function(n, ref) {
    n = n.clone();
    var p = BigDecimal.getSafePrecision();
    if(p < 0 && +ref) { //significant figures
        BigDecimal.scaleToExponent(n, ref.exponent.add(p));
    } else { //decimal places
        BigDecimal.scaleToExponent(n, p);
    }
    return n.mantissa.equals(0);
};

5. Reusing Parser operations

Currently the Parser operations support Symbols, Vectors and Matrices. I propose that BigDecimals (, Fracs and bigInts)[OPTIONAL] should also be supported.

For example:

        this.multiply = function(a, b) { 
            return _.mutatingMultiply(a.clone(), b.clone());
        };

5.1. BigDecimals

Add checks for BigDecimals, e.g.

        this.mutatingMultiply = function(a, b) { 
+           if (a instanceof BigDecimal) {
+               if (b instanceof Symbol) b = b.toDecimal(Settings.PRECISION);
+               if (b instanceof BigDecimal) {
+                   a.mantissa = a.mantissa.multiply(b.mantissa);
+                   a.exponent = a.exponent.multiply(b.exponent);
+                   return a;
+               }
+           }
@@ @@ Omitted code
        }

Here are a few helper functions:

BigDecimal.integer = function(int) { return new BigDecimal(int, bigInt.zero); };
    function Symbol(obj) { 
        var isInfinity = obj === 'Infinity';
        //this enables the class to be instantiated without the new operator
        if(!(this instanceof Symbol)) { 
            return new Symbol(obj); 
-       };
+       }
+       if(obj instanceof Symbol) return obj;
    Symbol.integer = function(int) {
        var s = new Symbol();
        s.group = N;
        s.value = CONST_HASH; 
        s.multiplier = Frac.integer(int);
        return s;
    };
    Frac.integer = function(int) {
        var f = new Frac();
        f.num = bigInt(int);
        f.den = bigInt.one;
        return f;
    };

By the way, all calls like _.parse('1') or new Symbol(1) should be replaced with Symbol.integer(1).

5.2. Fracs [OPTIONAL]

Move all instance functions to Parser, and convert them to BigDecimals if the other argument is BigDecimal.

        this.mutatingMultiply = function(a, b) { 
            if (a instanceof BigDecimal) {
-               if (b instanceof Symbol) b = b.toDecimal(Settings.PRECISION);
+               if (b instanceof Symbol || b instanceof Frac) b = b.toDecimal(Settings.PRECISION);
                if (b instanceof BigDecimal) {
                    a.mantissa = a.mantissa.multiply(b.mantissa);
                    a.exponent = a.exponent.multiply(b.exponent);
                    return a;
                }
            }
+           if (a instanceof Frac) {
+               if (b instanceof BigDecimal) return _.mutatingMultiply(b, a);
+               if (b instanceof Frac) {
+                   var n1 = a.den, n2 = b.den;
+                   var a = a.num, b = b.num;
+                   if(n1.equals(n2)) {
+                       a.num = a.add(b);
+                   }
+                   else {
+                       a.num = a.multiply(n2).add(b.multiply(n1));
+                       a.den = n1.multiply(n2);
+                   }
+                   return a.simplify();
+               }
+               if (b instanceof Symbol) return _.mutatingMultiply(a, new Symbol(b));
+           }
@@ @@ Frac.prototype = {
-       neg: function() {
-           this.num = this.num.multiply(-1);
-           return this;
-       },
-       add: function(m) { 
-           var n1 = this.den, n2 = m.den, c = this.clone();
-           var a = c.num, b = m.num;
-           if(n1.equals(n2)) {
-               c.num = a.add(b);
-           }
-           else {
-               c.num = a.multiply(n2).add(b.multiply(n1));
-               c.den = n1.multiply(n2);
-           }
-           return c.simplify();
-       },
@@ @@
    function Symbol(obj) { 
        var isInfinity = obj === 'Infinity';
        //this enables the class to be instantiated without the new operator
        if(!(this instanceof Symbol)) { 
            return new Symbol(obj); 
        };
        //define numeric symbols
-       if(/^(\-?\+?\d+)\.?\d*e?\-?\+?\d*/i.test(obj)) { 
+       if(obj instanceof Frac || /^(\-?\+?\d+)\.?\d*e?\-?\+?\d*/i.test(obj)) { 
            this.group = N;
            this.value = CONST_HASH; 
            this.multiplier = new Frac(obj);
        }

5.3. bigInts [OPTIONAL]

May not need to actually mutate the arguments, see Section 5.4.

5.4. Determinism of the mutated argument of mutating operations

If we want the mutating behaviour be deterministic (preferred): The argument which is mutated should be determined by:

  1. The broadest numeric data type: bigInt [OPTIONAL] < Number < Frac [OPTIONAL] < Symbol (with no variables) < BigDecimal | bigInt [OPTIONAL] < Number < Frac [OPTIONAL] < Symbol (all groups)
  2. The first argument if both arguments are of the same type.

For Numbers and bigInts, (choose one) a. When the arguments are NOT mutated despite using the mutating operations: The argument to be mutated is of type bigInt or Number. In this case, variable = _.mutatingMultiply(variable, ...) cannot be reduced to _.mutatingMultiply(variable, ...). (preferred) b. Force mutation of bigInts by copying the implementation inside BigInteger.js, and force mutation of Numbers by using Number wrappers and converting every === on Numbers to ==:

var a = { value: 1, valueOf: function() { return this.value; } };
a == 1; //true

Pros and Cons: Pros: Further reduce the performance cost of clone()s by changing _.add to only clone when necessary. Enables the Parser operations to be further restructured to e.g. _.add(a, b, mutate). Cons: Need to remember the rules of mutation (can reference here though); One small change in these operations can result in a lot of failing unit tests (is normal in a development process)

If we want the mutating behaviour to be non-deterministic: All uses of mutating operations must be in the form variable = _.mutatingMultiply(variable, ...). Pros: No need to remember the rules of mutation. Cons: Cannot effectively remove the performance cost of clone()ing.

5.5. BigDecimal with variables

When a BigDecimal is operated with a Symbol with variables, (choose one) a. the Symbol is mutated with its group set to D (a new group) and disallowed to operate with other symbols; or b. the Symbol is mutated with its group set to DN (e.g. 2.33e2), DS (e.g. x^2.33e2), DEX (e.g. 2.332^x), etc. and are dealt separately; or c. a new DecimalSymbol is introduced as the return type and disallowed to operate with other symbols, not mutating any of the arguments; or d. the multiplier or power of the Symbol is assigned a BigDecimal and its group changed accordingly (need to go through checking that all code looking for Symbols will work).

5.6. Alternatives

Note: the mutating operations may need to be rethought if one of the following is adopted.

Alternative 1

Define all methods on Symbol as mutating and all methods on Parser as non-mutating. After restructuring, a.add(b) would mutate a and _.add(a, b) would return a new Symbol for Symbols a and b.

Alternative 2

Like Alternative 1 but the other way around. a.add(b) would return a new Symbol while _.add(a, b) would mutate a.

Note that both alternatives require rewriting of all code from Section 5 onwards.

5.7. Numbers [ABANDONED]

~~These are easy: Just use JavaScript operators or use Math functions. For conversion to Fracs, use new Frac(number); for Symbols, use new Symbol(new Frac(number)); for BigDecimals, use~~

var number = ...;
var exponent = 0;
if(typeof number !== 'number' || !isFinite(number)) throw new TypeError('number is not a valid decimal!');
while(true) {
    if(number % 1 === 0) return BigDecimal.simplify(new BigDecimal(number, exponent));
    number *= 10;
    exponent--;
} 

~~. May not need to actually mutate the arguments, see Section 5.4.~~

6. Constructing via infinite series or infinite products

Since all elementary functions and most special functions are analytic, they can be defined as an infinite series. Therefore, BigDecimals should be able to be constructed from an infinite series or an infinite product. Construction from an infinite product is included in case it is needed.

Requires Frac.toDecimal to return a BigDecimal and this:

        //forces the symbol to be returned as a decimal
-       toDecimal: function(prec) {
+       toDecimal: function(prec, forceNumeric) {
+           if(forceNumeric)
+               switch(this.group) {
+                   case N:
+                       return this.multiplier.toDecimal(prec);
+                   case P:
+                       return _.mutatingPow(this.multiplier.toDecimal(prec), this.power.toDecimal(prec));
+                   default:
+                       throw new Error('Symbol is not numeric!');
+               }
+           var prevPrec = Settings.PRECISION;
-           Settings.precision = prec;
+           Settings.PRECISION = prec;
            var dec = text(this.symbol, 'decimals');
-           Settings.precision = undefined;
+           Settings.PRECISION = prevPrec;
            return dec;
        },

Frac.toDecimal should also have the second argument specified in the argument list to avoid misaligned new arguments later on.

Frac.toDecimal should support negative precision, aka significant figures.

/**
 * Creates a BigDecimal from an infinite series.
 * 
 * @param {Function} generalTerm - The function representing the general term with one bigInt parameter and one BigDecimal result
 * @param {bigInt} i - The lower bound of the infinite series, defaults to 0
 * @returns {BigDecimal}
 */
BigDecimal.fromInfiniteSeries = function(generalTerm, i) {
    i = i || bigInt.zero;
    var current, decimal = new BigDecimal(bigInt.zero, bigInt.zero);
    do {
        //add current term
        _.mutatingAdd(decimal, current = generalTerm(i));

        //advance index by one
        i = i.next();
    } while(!BigDecimal.isNegligible(current, decimal));
}

/**
 * Creates a BigDecimal from an infinite product.
 * 
 * @param {Function} generalTerm - The function representing the general term with one bigInt parameter and one BigDecimal result
 * @param {bigInt} i - The lower bound of the infinite product, defaults to 1
 * @returns {BigDecimal}
 */
BigDecimal.fromInfiniteProduct = function(generalTerm, i) {
    i = i || bigInt.one;
    var current, decimal = new BigDecimal(bigInt.one, bigInt.zero);
    do {
        //multiply current term
        _.mutatingMultiply(decimal, current = generalTerm(i));

        //advance index by one
        i = i.next();
    } while(!BigDecimal.isNegligible(current, decimal) || !BigDecimal.isNegligible(_.mutatingSubtract(current, 1), decimal));
}

7. Optimizations for constants

For constants like pi and e, I suggest a BigDecimalGenerator structure:

function BigDecimalGenerator(generalTerm, i, partial, preciseTo) {
    function err() {
        throw new Error('The constructor of BigDecimalGenerator was called directly!');
    }
    if(!this instanceof BigDecimalGenerator) err();

    this.generalTerm = generalTerm;
    this.i = i;
    this.partial = partial;
    this.preciseTo = preciseTo;
    this.generate = err;
}

/**
 * Creates a BigDecimalGenerator from an infinite series with a partial sum.
 * 
 * @param {Function} generalTerm - The function representing the general term with one bigInt parameter and one Frac result, to be called if a more accurate value is requested
 * @param {bigInt} i - The index of the last term of the partial sum
 * @param {Frac} partialSum - The current partial sum of the infinite series
 * @param {Number} preciseTo - The number of decimals that this partial sum is accurate to
 * @returns {BigDecimalGenerator}
 */
BigDecimalGenerator.fromInfiniteSeries = function(generalTerm, i, partialSum, preciseTo) {
    var g = new BigDecimalGenerator(generalTerm, i, partialSum, preciseTo);
    g.generate = function() {
        var prec = Settings.PRECISION;

        if(this.preciseTo < prec) {
            var p = BigDecimal.getSafePrecision(), current;
            do {
                //add current term
                this.partial = this.partial.add(current = this.generalTerm(this.i));
                //Or if mutating Parser operations support Fracs: 
                //_.mutatingAdd(this.partial, current = this.generalTerm(this.i));

                //advance index by one
                this.i = this.i.next();
            } while(!BigDecimal.isNegligible(current.toDecimal(p), this.partial.toDecimal(p)))
            this.preciseTo = prec;
        }
       return this.partial.toDecimal(p);
    };
    return g;
}

/**
 * Creates a BigDecimalGenerator from an infinite product with a partial product.
 * 
 * @param {Function} generalTerm - The function representing the general term with one bigInt parameter and one Frac result, to be called if a more accurate value is requested
 * @param {bigInt} i - The index of the last term of the partial product
 * @param {Frac} partialProduct - The current partial product of the infinite product
 * @param {Number} preciseTo - The number of decimals that this partial product is accurate to
 * @returns {BigDecimalGenerator}
 */
BigDecimalGenerator.fromInfiniteProduct = function(generalTerm, i, partialProduct, preciseTo) {
    var g = new BigDecimalGenerator(generalTerm, i, partialProduct, preciseTo);
    g.generate = function() {
        if(this.preciseTo < Settings.PRECISION) {
            var p = BigDecimal.getSafePrecision(), current;
            do {
                //multiply current term
                this.partial = this.partial.multiply(current = this.generalTerm(this.i));
                //Or if mutating Parser operations support Fracs: 
                //_.mutatingMultiply(this.partial, current = this.generalTerm(this.i));

                //advance index by one
                this.i = this.i.next();
            } while(!BigDecimal.isNegligible(current.toDecimal(p), this.partial.toDecimal(p)) || !BigDecimal.isNegligible(_.mutatingSubtract(current.toDecimal(p), BigDecimal.integer(1)), this.partial.toDecimal(p)))
            this.preciseTo = Settings.PRECISION;
        }
       return this.partial.toDecimal(p);
    };
    return g;
}

/**
 * Creates a BigDecimalGenerator from an infinite series where each term depends on the previous term and a partial sum.
 * 
 * @param {Function} generalTerm - The function representing the general term with one bigInt parameter and one Frac result, to be called if a more accurate value is requested
 * @param {Frac} partialProduct - The current partial sum of the infinite series
 * @param {Number} preciseTo - The number of decimals that this partial sum is accurate to
 * @returns {BigDecimalGenerator}
 */
BigDecimalGenerator.fromBootstrapSeries(generalTerm, partialSum, preciseTo) {
    var g = new BigDecimalGenerator(generalTerm, undefined, partialSum, preciseTo);
    g.generate = function() {
        var prec = Settings.PRECISION;

        if(this.preciseTo < prec) {
            var p = BigDecimal.getSafePrecision(), current;
            do {
                //add current term
                this.partial = this.partial.add(current = this.generalTerm(this.partial));
                //Or if mutating Parser operations support Fracs: 
                //_.mutatingAdd(this.partial, current = this.generalTerm(this.partial));
            } while(!BigDecimal.isNegligible(current.toDecimal(p), this.partial.toDecimal(p)))
            this.preciseTo = prec;
        }
       return this.partial.toDecimal(p);
    };
    return g;
}

The decimal generators of constants like pi and e will be stored in the Parser. Whenever nerdamer evaluates an expression, the generate() method of each generator will be called and each of their values will be stored in the respective fields in MathD.

Like this: (Initialization)

        _.PI = BigDecimalGenerator.fromInfiniteSeries(function(k){return _.ifactorial(k).multiply(2).divide(_.idfactorial(k.multiply(2).next()))}, bigInt[100], Frac.quick('69857222584256563248421626232855718905618239601272387446360031947601767526320570368','22236244569910437892880691136109856086794490658698290303471308035343578465879432175'), 30),
        _.E = BigDecimalGenerator.fromInfiniteSeries(function(k){return Frac.quick(bigInt[3].minus(k.square().multiply(4)),_.ifactorial(k.multiply(2).next()))}, bigInt[30], Frac.quick('41810487026297450678630938607562589344797971920562745147744661141255717033116978111','15381218602340145418207782187170461431091046332970787972709270893363200000000000000'), 83),
        _.SQRT2 = BigDecimalGenerator.fromInfiniteSeries(function(k){return Frac.quick(_.ifactorial(k.multiply(2).next()),bigInt[2].pow(k.multiply(3).next()).multiply(_.ifactorial(k).square())}, bigInt[105], Frac.quick('11799772793954915567626029419139823535513556097082996259637480088230036296369382361445481881759','8343699359066055009355553539724812947666814540455674882605631280555545803830627148527195652096'), 30),
        _.SQRT1_2 = BigDecimalGenerator.fromInfiniteSeries(Frac.quick(_.ifactorial(k.multiply(2).next()),bigInt[2].pow(k.multiply(3).add(2)).multiply(_.ifactorial(k).square())), bigInt[105], Frac.quick('11799772793954915567626029419139823535513556097082996259637480088230036296369382361445481881759','16687398718132110018711107079449625895333629080911349765211262561111091607661254297054391304192'), 31)
        // similarly for LN2, LN10, LOG2E and LOG10E

(Whenever an expression is evaluated)

        MathD.PI = _.PI.generate();
        MathD.E = _.E.generate();
        MathD.SQRT2 = _.SQRT2.generate();
        MathD.SQRT1_2 = _.SQRT1_2.generate();
        // similarly for LN2, LN10, LOG2E and LOG10E

8. An arbitrarily precise Math object

The Parser currently both parses symbols and abstracts math operations. I propose to separate the math functions inside the Parser that can be made to only use abstract operations into Math2.

Therefore, the Parser must have functionality of all JavaScript math operators, and a new MathD object must be made to have functionality equivalent to the built-in Math object (and nothing more).

For functions that cannot return finite decimals, infinite series or infinite products should be used to provide decimal values when PARSE2NUMBER is set or evaluate() is called.

I propose that all functions in Math2 that cannot depend on Parser and MathD entirely be moved out into the Parser. This avoids the current near-spaghetti code that mixes logic when an integer is encountered and fraction is encountered.

The Math2 object should be entirely dependent on the Parser and MathD, not referencing any implementation details like Number, bigInt, Frac, Symbol or BigDecimal.

Question: Should the symbolic part of trig functions etc. be moved entirely into MathD or only the numeric part will be in there?

For example:

+       //More precise than built-in Math object
+       MathD = {
+           exp: function(x) {
+               return BigDecimal.fromInfiniteSeries(function(k){return _.divide(_.pow(x,Symbol.integer(k)), _.ifactorial(k)));
+           }
+       },
        //This object holds additional functions for nerdamer. Think of it as an extension of the Math object.
        //I really don't like touching objects which aren't mine hence the reason for Math2. The names of the 
        //functions within are pretty self-explanatory.
        Math2 = {
@@ -904 @@
            },
-           //factorial
-           bigfactorial: function(x) {
-               var retval = new Frac(1);
-               for (var i = 2; i <= x; i++) 
-                   retval = retval.multiply(new Frac(i));
-               return retval;
-           },
            //https://en.wikipedia.org/wiki/Logarithm#Calculation
@@ -927 @@
                    retval = retval.add(r);

                }
                return retval.multiply(new Frac(2));
            },
-           //the factorial function but using the big library instead
-           factorial: function(x) {
-               if(x < 0)
-                   throw new Error('factorial not defined for negative numbers');
-               var retval=1;
-               for (var i = 2; i <= x; i++) retval = retval * i;
-               return retval;
-           },
-           //double factorial
-           dfactorial: function(x) {
-               var even = x % 2 === 0;
-               // If x = even then n = x/2 else n = (x-1)/2
-               var n = even ? x/2 : (x+1)/2; 
-               //the return value
-               var r = new Frac(1);
-               //start the loop
-               if(even)
-                   for(var i=1; i<=n; i++)
-                       r = r.multiply(new Frac(2).multiply(new Frac(i)));
-               else
-                   for(var i=1; i<=n; i++)
-                       r = r.multiply(new Frac(2).multiply(new Frac(i)).subtract(new Frac(1)));
-               //done
-               return r;
-            },
            GCD: function() {
                var args = arrayUnique([].slice.call(arguments)
@@ @@
+          //integer factorial
+           _.ifactorial = function(x) {
+               x = +x; //because x must be small or this will lag
+               if(x < 0)
+                   throw new Error('factorial not defined for negative numbers');
+               var retval = bigInt.one;
+               for (var i = 2; i <= x; i++) retval = retval.multiply(i);
+               return retval;
+           };
+           //integer double factorial
+           _.idfactorial = function(x) {
+               x = +x; //because x must be small or this will lag
+               var even = x % 2 === 0;
+               // If x = even then n = x/2 else n = (x-1)/2
+               var n = even ? x/2 : (x+1)/2; 
+               //the return value
+               var r = bigInt.one;
+               //start the loop
+               if(even)
+                   for(var i=1; i<=n; i++)
+                       r = r.multiply(bigInt[2]).multiply(i);
+               else
+                   for(var i=1; i<=n; i++)
+                       r = r.multiply(bigInt[2]).multiply(i).prev();
+               //done
+               return r;
+           }
+           _.icombination = function(n, r) { return _.ifactorial(n).divide(_.ifactorial(n.subtract(r)).multiply(_.ifactorial(r))); };

The Math2 object also should not use any mutating methods nor clone() for easier implementation of Section 12.

All functions in the following two sections would work on BigDecimals and have arbitrary precision, and never mutate the arguments.

8.1. JavaScript operations/functions to implement

Each implementation may rely on bigInteger.js, BigDecimal methods, fields of mathematical constants in MathD or functions that came before it. For simplicity's sake, BigDecimal.fromInfiniteSeries(function(k){...},0) is simplified to sum(...).

All the implementations should be used for PARSE2NUMBER scenarios in their symbolic equivalent.

Some functions may involve BigDecimal raised to an integer power. The method goes like this (No BigDecimal.simplify to increase speed): function dpowi(x,k) { var c = x.clone(); c.mantissa = c.mantissa.pow(k); c.exponent = c.exponent.pow(k); return c; } In my pseudo-implementations, I may mix JavaScript notation with math notation. The ^ operator always stands for _.pow and not _.xor.

JavaScript Operation Equivalent Pseudo-implementation with BigDecimals
x+=y _.mutatingAdd(x,y) (ASSUMING type of x is wider than y) BigDecimal.equateExponents(x, y); BigDecimal.simplify(new BigDecimal(x.mantissa.add(y.mantissa), x.exponent))
x-=y _.mutatingSubtract(x,y) (ASSUMING type of x is wider than y) BigDecimal.equateExponents(x, y); BigDecimal.simplify(new BigDecimal(x.mantissa.subtract(y.mantissa), x.exponent))
x*=y _.mutatingMultiply(x,y) (ASSUMING type of x is wider than y) BigDecimal.simplify(new BigDecimal(x.mantissa.multiply(y.mantissa), x.exponent.multiply(y.exponent)))
x/=y _.mutatingDivide(x,y) (ASSUMING type of x is wider than y) BigDecimal.equateExponents(x, y); new Frac(x.mantissa, y.mantissa).toDecimal(BigDecimal.getSafePrecision())
x%=y _.mutatingMod(x,y) (ASSUMING type of x is wider than y) BigDecimal.equateExponents(x, y); BigDecimal.simplify(new BigDecimal(x.mantissa.mod(y.mantissa), x.exponent))
x++ (mutating) .mutatingInc(x) or \.mutatingIncrement(x) _.mutatingAdd(x,1)
x-- (mutating) .mutatingDec(x) or \.mutatingDecrement(x) _.mutatingSubtract(x,1)
x=-x _.mutatingNegate(x) x.mantissa=x.mantissa.negate();x
x<<=y .mutatingLShift(x,y) or \.mutatingLeftShift(x,y) (ASSUMING type of x is wider than y) (See Section 8.1.1)
x>>=y .mutatingRShift(x,y) or \.mutatingRightShift(x,y) (ASSUMING type of x is wider than y) (See Section 8.1.1)
x&=y _.mutatingAnd(x,y) (ASSUMING type of x is wider than y) BigDecimal.scaleToExponent(x,0);BigDecimal.scaleToExponent(y,0);x.mantissa.and(y.mantissa)
x|=y _.mutatingOr(x,y) (ASSUMING type of x is wider than y) BigDecimal.scaleToExponent(x,0);BigDecimal.scaleToExponent(y,0);x.mantissa.or(y.mantissa)
x^=y _.mutatingXor(x,y) (ASSUMING type of x is wider than y) BigDecimal.scaleToExponent(x,0);BigDecimal.scaleToExponent(y,0);x.mantissa.xor(y.mantissa)
x=~x _.mutatingNot(x) BigDecimal.scaleToExponent(x,0);BigDecimal.scaleToExponent(y,0);x.mantissa.not(y.mantissa)
x+y _.add(x,y) _.mutatingAdd(x.clone(), y.clone())
x-y _.subtract(x,y) _.mutatingSubtract(x.clone(), y.clone())
x*y _.multiply(x,y) _.mutatingMultiply(x.clone(), y.clone())
x/y _.divide(x,y) _.mutatingDivide(x.clone(), y.clone())
x%y _.mod(x,y) _.mutatingMod(x.clone(), y.clone())
x++ (NOT mutating) .inc(x) or \.increment(x) _.mutatingInc(x.clone())
x-- (NOT mutating) .dec(x) or \.decrement(x) _.mutatingDec(x.clone())
+x _.parse(x) (already implemented)
-x _.negate(x) _.mutatingNegate(x.clone())
x<<y .lshift(x,y) or \.leftshift(x,y) _.mutatingLShift(x.clone(),y.clone())
x>>y .rshift(x,y) or \.rightshift(x,y) _.mutatingRShift(x.clone(),y.clone())
x&y _.and(x,y) _.mutatingAnd(x.clone(), y.clone())
x|y _.or(x,y) _.mutatingOr(x.clone(), y.clone())
x^y _.xor(x,y) _.mutatingXor(x.clone(), y.clone())
~x _.not(x) _.mutatingNot(x.clone())
x==y .equals(x,y) or \.eq(x,y) x.equals(y)
x!=y .notEquals or \.neq(x,y) x.neq(y)
x<y .lesser(x,y) or \.lt(x,y) x.lt(y)
x>y .greater(x,y) or \.gt(x,y) x.gt(y)
x<=y .lesserOrEquals(x,y) or \.leq(x,y) x.leq(y)
x>=y .greaterOrEquals(x,y) or \.geq(x,y) x.geq(y)
x=y x=y x=y
x===y x===y x===y
x!==y x!==y x!==y
Math.exp(x) MathD.exp(x) sum(x^k/k!)
Math.log(x) MathD.log(x) 2(x-1)/(x+1)sum(1/(2k+1)((x-1)^2/(x+1)^2)^k)
Math.log10(x) MathD.log10(x) MathD.log(x)/MathD.LN10
Math.log2(x) MathD.log2(x) MathD.log(x)/MathD.LN2
x**y or Math.pow(x,y) _.pow(x,y) MathD.exp(y*MathD.log(x))
x**=y _.mutatingPow(x,y) (ASSUMING type of x is wider than y) x=_.pow(x,y)
Math.abs(x) MathD.abs(x) var n = x.clone(); n.mantissa = n.mantissa.abs(); n
Math.sqrt(x) MathD.sqrt(x) MathD.exp(MathD.log(x)/2)
Math.cbrt(x) MathD.cbrt(x) MathD.exp(MathD.log(x)/3)
Math.asin(x) MathD.asin(x) sum(_.icombination(2*k,k)x^(2k+1)/(4^k(2k+1)))
Math.acos(x) MathD.acos(x) pi/2-asin(x)
Math.atan(x) MathD.atan(x) sum(2^(2k)(k!)^2/(2k+1)!x^(2k+1)/(1+x^2)^(k+1))
Math.atan2(y, x) MathD.atan2(y, x) x>0 ? 2*atan(y/(sqrt(x^2+y^2)+x)) : x<=0 && y!=0 ? 2*atan(y/(sqrt(x^2+y^2)+x)) ? x<0 && y=0 ? pi : throw new UndefinedError("atan2 is undefined for y="+y+" and x="+x)
Math.asinh(x) MathD.asinh(x) MathD.log(x+sqrt(x^2+1))
Math.acosh(x) MathD.acosh(x) MathD.log(x+sqrt(x^2-1))
Math.atanh(x) MathD.atanh(x) MathD.log((1+x)/(1-x))/2
Math.floor(x) MathD.floor(x) x - _.mod(x,1)
Math.ceil(x) MathD.ceil(x) x + _.mod(-x,1)
Math.sin(x) MathD.sin(x) sum((-1)^k*x^(2k+1)/(2k+1)!)
Math.cos(x) MathD.cos(x) sum((-1)^k*x^(2k)/(2k)!)
Math.tan(x) MathD.tan(x) sin(x)/cos(x)
Math.sinh(x) MathD.sinh(x) sum(x^(2k+1)/(2k+1)!)
Math.cosh(x) MathD.cosh(x) sum(x^(2k)/(2k)!)
Math.tanh(x) MathD.tanh(x) sinh(x)/cosh(x)
Math.hypot([x[, y[, …]]]) MathD.hypot([x[, y[, …]]]) sqrt([x^2+[y^2[+ …]]])
x>>>y _.uRightShift(x,y) (See Section 8.1.1)
x>>>=y _.mutatingURightShift(x,y) (See Section 8.1.1)
Math.max([x[, y[, …]]]) MathD.max([x[, y[, …]]]) var ret = arguments[0]; for(var i = 1; i < arguments.length; i++) if(ret < arguments[i]) ret = arguments[i]; ret
Math.min([x[, y[, …]]]) MathD.min([x[, y[, …]]]) var ret = arguments[0]; for(var i = 1; i < arguments.length; i++) if(ret > arguments[i]) ret = arguments[i]; ret
Math.random() MathD.random() new BigDecimal(bigInt.randBetween("-1e"+(Math.abs(BigDecimal.getSafePrecision())), "1e"+(Math.abs(BigDecimal.getSafePrecision()))), bigInt(-Math.abs(BigDecimal.getSafePrecision()))
Math.round(x) MathD.round(x) ceil(x - 1/2) + (floor(1/4 (2 x + 1)) == x ? 1 : 0)
Math.sign(x) MathD.sign(x) x.mantissa.isZero() ? 0 : x.mantissa.isPositive() ? 1 : -1
Math.trunc(x) MathD.trunc(x) x <= 0 ? ceil(x) : floor(x)

8.1.1. Problem with shifting operators

Should we allow decimals in shifting operators? This will resolve into one of the following:

a. Disallow decimals This is in accordance to most programming languages. Attempting to use decimals with shifting operators will result in an Error.

b. Silently convert to integers All decimals will be floored to an integer. 0.4<<1 will result in 0 while 2.2>>1 will result in 1.

c. Implement as multiplication or division of 2 0.4<<1 will be 0.8 while 2.2>>1 will be 1.1. This means that 1>>1 will be 0.5 instead of 0. Wolfram Alpha uses this option but it computes (0.1*10)>>1 as 0.5 while 1>>1 results in 0.

d. Treat decimals differently from integers 1>>1 will be 0 while 0.5<<1 will be 1. This also means that (0.5<<1)>>1 will be 0. So, (x<<1)>>1 cannot be simplified.

e. Leave as unevaluated Similarly to how Wolfram Alpha handles BitAnd[2.2, 1].

8.2. Functions to implement for the sake of communicating with external math functions (see Section 12)

JavaScript Operation Equivalent Pseudo-implementation with BigDecimals
Math.imul(x, y) MathD.imul(x, y) See polyfill
Math.clz32(x) MathD.clz32(x) See polyfill
Math.expm1(x) MathD.expm1(x) exp(x)-1
Math.fround(x) MathD.fround(x) See polyfill
Math.log1p(x) MathD.log1p(x) MathD.log(1+x)

8.3. Function not worth implementing (LOL)

JavaScript Operation Equivalent Pseudo-implementation with BigDecimals
Math.toSource() MathD.toSource() Math

8.4. Function for parsing a decimal

Maybe this can be called by _.parse when PARSE2NUMBER is true. This might be useless though.

_.parseDecimal = function(str) {
    var i = str.indexOf('.');
    //str is an integer if it contains no .
    if (i === -1) return BigDecimal.integer(str);
    var l = str.length;
    str = str.replace(/^0+|0+$/g, '');
    return new BigDecimal(str.slice(0,i)+str.slice(i+1),-l+i+1);
};

8.5. Functions to import from BigInteger.js [BONUS]

There are a number of functions in BigInteger.js that are not part of nerdamer. Aside from implementing the functions from the built-in Math object, why not reference some from BigInteger.js too? :wink:

These should probably be able to be used from nerdamer itself, e.g. nerdamer("isEven(2)")

9. Exposing the API to Expression

The Expression class currently exposes add, subtract, multiply, divide and pow. Other math operations should be exposed too, including all of the functions in nerdamer that match with the ones in the JavaScript Math class (Sections 8.1, 8.2), along with the ones in Section 8.5.

The toDecimal method should return an Expression with a decimal symbol which behavior is as chosen in Section 5.5. Since the decimal symbol handles all the decimal-symbol interoperations, the Expression class can be used as is.

10. Text and LaTeX representation

Scientific notation for BigDecimals is defined as mantissa×10ᵉˣᵖᵒᵑᵉᵑᵗ, where a dot is inserted after the leading digit of the mantissa and the number of digits after the dot is added to the exponent.

Representations of decimal symbols should follow the scientific notation if its exponent after addition is out of the range [-5, 5], and display as normal if not.

The LaTeX for × is \times while the text for × is *.

For example (under precision over 6 decimal places or 7 significant figures):

Mantissa Exponent Text LaTeX LaTeX Image
22 0 22 22 image
22 -1 2.2 2.2 image
2 6 2*10^6 2 \times 10^6 image
1234567 0 1.234567*10^6 1.234567 \times 10^6 image
123456 -6 0.123456 0.123456 image
1 -5 0.00001 0.00001 image
-1 -6 -10^-6 -10^{-6} image

When the mantissa is 1 or -1 and scientific notation is used, the 1 and the multiplication sign will be omitted. When the mantissa is 0, the entire BigDecimal will display as a good ol' egg.

Remember to bracket the decimal if it becomes the base of exponentiation and requires scientific notation to display, or the exponent of exponentiation when displayed as text.

Example Text LaTeX LaTeX Image
_.pow(new BigDecimal(2,6), new Symbol('x')) (2*10^6)^x (2 \times 10^6)^x image
_.pow(new Symbol('x'), new BigDecimal(1,6)) x^(10^6) x^{10^6} image

The text and LaTeX described in this Section will be used for toString(), text(), toTeX() and latex() of decimal symbols and Expressions with them.

11. External math functions but with arbitrary precision

Let's create a new function called importFunction() and buff replaceFunction(). importFunction() takes the same arguments as replaceFunction() and exposes the given function to nerdamer("...") and nerdamer["..."] after the following enhancements.

Whenever a function is imported with one of the aforementioned two functions, it is parsed and all the operations will be replaced with calls to the Parser.

For example, function(x,y){ return x+y; } will be converted to function(x,y){ return _.add(x,y); } before being consumed.

This enables the following amazing scenario:

function someFunc(x, y) {
    //complicated logic
    x = x + y;
    y = y + x;
    var z = Math.sin(y) - Math.cos(x) + 0.245;
    y = z++;
    return ++z + x / y;
}
nerdamer.importFunction('someFunc', someFunc, 2);
nerdamer.someFunc(1,2).text(); //works!
nerdamer("someFunc(1,2)").text(); //works too!

The function is transformed to:

function someFunc(x, y) {
    //complicated logic
    x = _.add(x, y);
    y = _.add(y, x);
    var z = _.add(_.subtract(_.trig.sin(y), _.trig.cos(x)), _.parse("0.245"));
    y = z;
    _.mutatingNext(z);
    return _.add(_.mutatingNext(z), _.divide(x, y));
}

This can potentially be simplified more by throwing new Symbol('x') and new Symbol('y') into it and observing the outcome.

Now, we can enjoy symbolic and arbitrary-precise numeric output simply via nerdamer("someFunc(1,2)"). This can surely become a selling point for nerdamer.

12. Interoperability with buildFunction

Now we can fix buildFunction(). Previously, if you did nerdamer("sinc(x)").buildFunction(), you'd get function(x) {var sinc = function (x) { if(x.equals(0)) return 1; return Math.sin(x)/x; }; return sinc(x); }. Notice the equals(0).

If we rewrite all the functions a bit to require Parser methods to be used, then our _.equals can be translated to == under the tables of Section 8.

The tables of Section 8 should provide enough information to replace Parser methods to JavaScript native operations. Neat!

13. Related issues

221

269

357

394

14. Postscript

GG. Two months of work in one proposal (March 7 ~ May 5, 2018). This is totally unexpected by me :sob:

This one proposal has totally drained my motivation for new proposals/PRs. I am going to post this as the #404 Not Found special, as well as pre-celebration for my first year anniversary in here (August 10, 2018) :wink: I am going to post my to-do list as #405 then take a temporary break from nerdamer. I will still respond to new issues and conversations though. Goodbye :wave:

jiggzson commented 6 years ago

@Happypig375, thanks for the insightful review. A few points:

  1. Item 3 is only a temporary fix. I guess it would improve readability but it still doesn't eliminate clone being called unnecessarily. The ideal solution would be just that, to determine when clone is needed. My initial thought was that cloning wasn't need at all but I quickly learned otherwise and this patch has been the running solution since.

  2. Item 6 would be another one that would be ideal but due to limitations in JavaScript numbers and the need to use a BigNumber library, the results haven't been the greatest. The problem is worsened by the fact that nerdamer uses the Frac library which uses a bigInt for both the denominator and the numerator. Just imagine what happens when we get a really tiny number for instance. The denominator get HUGE quickly in an infinite series. Even with a low number of iterations.

  3. Item 7 is on the TODO list for the parser rewrite and has given me much grief. I suspect that it might take more than one iteration to get this just right but let's see.

I'm going to be re-reading this issue since there's much to consider and it's pretty long :smile:

Happypig375 commented 6 years ago
  1. The idea here is to follow https://github.com/jiggzson/nerdamer/pull/337#issuecomment-353820872 and only minimize clone()s when a notable performance difference is present.
  2. This proposal attempts to implement BigNumbers (BigDecimals) in terms of BigIntegers instead of strings, and BigInteger.js is pretty fast for arithmetic operations, and are probably much faster than string operations. For really tiny numbers that can potentially lag the library (say more than a few thousand), it will most probably be well under the precisions for everyday use. The benchmark has shown that BigInteger.js can still cope with them though.
  3. I don't think Item 7 without all previous items can work, hence the sequential order. Without a proper big number data type, it is going to be hard.

I'm going to be re-reading this issue since there's much to consider and it's pretty long 😄

And the amount of consideration have taken me 2 months to write this 😛

ibreaker14 commented 5 years ago

Hi @jiggzson do you have any updates for this issue?

jiggzson commented 5 years ago

@ibreaker14, yes. I'm going to be starting this shortly.