Multiple Precision in JavaScript: Rational Arithmetic III

Elementary Functions: Second Try

Not that I found the code for the rational exponential function I wrote (that’s lost) but I found it unfair to say “too slow” without giving a reason, so i did it again, for exp() and log().

Exponential Function

Let’s take the standard series for the exponential function, to simplify the thing a bit, together with some argument reduction and evaluate it with Horner’s scheme. We use the slightly different function E = \exp(x) - 1 because it has some properties useful for later implementation with Bigfloats.

We can implement the series very simply with the help of Bigint and Bigrational:

function expoz(n,z){
  var r = z,t,t2;
  t = z.divInt(n).addInt(1);
  for(var i = n-1;i>= 2;i--){
     t2 = z.divInt(i)//.fastround(BIGRATIONAL_PRECISION);
     t = t2.mul(t).addInt(1).fastround(BIGRATIONAL_PRECISION);
  }
  return t.mul(z)//.fastround(BIGRATIONAL_PRECISION);
}

The places to do a rounding are different for different input; the rounding function itself has a significant overhead. The placement above is for fractions using the full precision; with a precision of 10 decimal digits for example, the fraction \frac{1}{2} doesn’t use the full precision but the number \frac{14862109042}{5467464369} does.

We take the quite random fraction

{\displaystyle{ \frac{2993558589961767975520115124024319199289207105579416583949678276150555438468529121279331718975468284}{1089360959377383732084311481992855973031093237661551378735236638866020386130681850286603456932617083} }}

It has a hundred decimal digit precision.
The argument reduction gets done by dividing the number by d = 2^{15} and a later expansion by the iteration E_{n+1} = 2E_n + E_n^2 done \log_2 d times, the number of series terms is 20 (twenty);

BIGRATIONAL_PRECISION = 332; // 100 decimal digits
var num = "2993558589961767975520115124024319199289207105579416583949678276150555438468529121279331718975468284";
var den = "1089360959377383732084311481992855973031093237661551378735236638866020386130681850286603456932617083";
num = num.toBigint();
den = den.toBigint();
var N = 15;
den.lShiftInplace(N);
var ret = expoz(20,new Bigrational(num,den))
for(var k = 0;k< N;k++ ){
    var tmp = ret.mulInt(2);
    var tmp2 = ret.sqr()//.fastround(BIGRATIONAL_PRECISION);
    ret=tmp.add(tmp2).fastround(BIGRATIONAL_PRECISION);
}
ret = ret.addInt(1);
//console.log(ret.toString())

It needs about 9 (nine) seconds on my old and slow machine. That is just not acceptable.

Logarithmic Function to the Base e

The logarithm is even worse.
If we try the Padé approximants (continued fraction) from Jörg Arndt’s book together with a moderate argument reduction:

BIGRATIONAL_PRECISION = 332;
var num = "2993558589961767975520115124024319199289207105579416583949678276150555438468529121279331718975468284";
var den = "1089360959377383732084311481992855973031093237661551378735236638866020386130681850286603456932617083";
num = num.toBigint();
den = den.toBigint();
var N = 16; // 2^4
var BR = new Bigrational(num,den);
var Start = performance.now();
BR = BR.nthroot(N).subInt(1).fastround(BIGRATIONAL_PRECISION);
var stop = performance.now();
console.log(stop - Start);
var ret =  log_pade(28,BR);
stop = performance.now();
console.log(stop - Start);
Start = performance.now();
ret = ret.mulInt(16).fastround(BIGRATIONAL_PRECISION);
stop = performance.now();
console.log(stop - Start);
/*
   My output:
13938.746253002435
68769.60310699977
432.42878299951553
"8314754741434890005606506915947915692960187527007493539396437290363369820043517533705476450148602249987/8225332259772227587959495044322232660558536920237514246081714698244856299137514215595359255233659725457"
*/

More than one and a half minute. The Padé approximants produces about one decimal digit per round without reduction and needs for the resulting 101 rounds more than four and a half minutes. Reducing the argument further is of no use because the nthroot() doubles its runtime for every doubling of the exponent (the algorithm is ok, the implementation: not so 😉 ).
Replacing it with consecutive square-roots and using 2^{64} brings it down to about 75 seconds.

AGM

For higher precisions and floating point numbers the AGM method is the fastest (for now). If it is also faster for rational numbers needs to be evaluated but it would need \pi in full precision.
We could ask Google for some but the one I liked most is an obfuscated C-program by Dik Winter and Achim Flammenkamp but Google also gave me the information that it is based on the book “Unbounded Spigot Algorithms for the Digits of Pi”, by Jeremy Gibbons, Math. Monthly, April 2006, pages 318-328 so I had no need to decipher it. But I took the whole idea, so honor to whom honor is due.

function obfuscatedjspi(n) {
    var a;
    var b = 0;
    var ab;
    var c;
    var d = 0,
        e = 0;
    var f = Math.pow(10,4);
    var f5 = Math.floor(f/5);
    var g = 1,
        h = 0;
    var out, chunk;
    var digits;
    var m = 14;
    digits = Math.floor( (n + 40)/40 ) ;
    a = new Array(digits * 10 * m);
    c = digits * 10 * m;
    out = new Bigint();
    b = c;
    while (b) {
        d %= f;
        e = d;
        b--;
        g = b *2;
        while(g){
            ab = (h == 0) ? f5 : a[b];
            d = d * b + f * ab;
            g--;
            a[b] = d % g;
            d = Math.floor(d/g);
            b--;
            g = b * 2;
        }
        chunk = e + Math.floor(d / f);
        h = Math.floor(Math.log(chunk) / Math.log(10)) + 1;
        if(h != 0) {
        out = out.mulInt(f);
        out = out.addInt(chunk);
        }
        c -= m;
        b = c;
    }
    return out;
}

But such things are for moderate precisions only, higher precisions should be computed with the AGM, too.

The hypergeometric series might be a better idea, too. It allows to compute several useful constants with one small algorithm, albeit for low precisions only. See
http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.pdf.
An example of it is implemented in their online-calculator.

The AGM for the rational numbers makes not much sense as you can find out for yourself:

// a minute for 200 digits (actually produced: 348)
function agmpi(d) {
    var a = new Bigrational(1, 1);
    var b = new Bigrational(2, 1);
    var s = new Bigrational(1, 4);
    var p = new Bigrational(1, 1);
    var n = Math.ceil(Math.log(d) / Math.log(2));
    var k, c, P, t;
    b = b.sqrt();
    b.inverse();
    for (k = 0; k < n; k++) {
        c = a.add(b).divInt(2).fastround(BIGRATIONAL_PRECISION);
        b = a.mul(b).fastround(BIGRATIONAL_PRECISION).sqrt();
        t = c.sub(a).sqr().fastround(BIGRATIONAL_PRECISION);
        t = p.mul(t) //.fastround(BIGRATIONAL_PRECISION);
        s = s.sub(t) //.fastround(BIGRATIONAL_PRECISION);
        p = p.mulInt(2) //.fastround(BIGRATIONAL_PRECISION);
        a = c;
    }
    P = a.sqr().fastround(BIGRATIONAL_PRECISION).div(s).round(
        BIGRATIONAL_PRECISION);
    return P;
}

If you thought the algorithms above were slow with rationals… 😉

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s