Skip to content

Instantly share code, notes, and snippets.

@getify
Created November 16, 2015 22:10
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save getify/f770998e5094bda75ae0 to your computer and use it in GitHub Desktop.
Save getify/f770998e5094bda75ae0 to your computer and use it in GitHub Desktop.
// This is a port of the C-language reference implementation, which can be
// found here: http://www.experimentalmath.info/bbp-codes/piqpr8.c
/*
This program implements the BBP algorithm to generate a few hexadecimal
digits beginning immediately after a given position id, or in other words
beginning at position id + 1. On most systems using IEEE 64-bit floating-
point arithmetic, this code works correctly so long as d is less than
approximately 1.18 x 10^7. If 80-bit arithmetic can be employed, this limit
is significantly higher. Whatever arithmetic is used, results for a given
position id can be checked by repeating with id-1 or id+1, and verifying
that the hex digits perfectly overlap with an offset of one, except possibly
for a few trailing digits. The resulting fractions are typically accurate
to at least 11 decimal digits, and to at least 9 hex digits.
*/
/* David H. Bailey 2006-09-08 */
function getDigit(id) {
var pid, s1, s2, s3, s4;
var NHX = 16;
var chx;
/* id is the digit position. Digits generated follow immediately after id. */
s1 = series(1, id);
s2 = series(4, id);
s3 = series(5, id);
s4 = series(6, id);
pid = 4 * s1 - 2 * s2 - s3 - s4;
pid = pid - Math.trunc(pid) + 1;
chx = ihex(pid, NHX);
// printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n",
// id, pid, chx);
// console.log(chx);
return chx;
}
function ihex(x, nhx) {
/* This returns the first nhx hex digits of the fraction of x. */
var i;
var y;
var hx = "0123456789ABCDEF";
var chx = [];
y = Math.abs(x);
for (i = 0; i < nhx; i++) {
y = 16 * (y - Math.floor(y));
chx[i] = hx[Math.trunc(y)];
}
return chx;
}
function series(m, id) {
/* This routine evaluates the series sum_k 16^(id-k)/(8*k+m)
using the modular exponentiation technique. */
var k;
var ak, eps, p, s, t;
var eps = 1e-17;
s = 0;
/* Sum the series up to id. */
for (k = 0; k < id; k++) {
ak = 8 * k + m;
p = id - k;
t = expm(p, ak);
s = s + t / ak;
s = s - Math.trunc(s);
}
/* Compute a few terms where k >= id. */
for (k = id; k <= id + 100; k++) {
ak = 8 * k + m;
t = Math.pow(16, (id - k)) / ak;
if (t < eps) break;
s = s + t;
s = s - Math.trunc(s);
}
return s;
}
function expm(p, ak) {
/* expm = 16^p mod ak. This routine uses the left-to-right binary
exponentiation scheme. */
var i, j;
var p1, pt, r;
var ntp = 25;
expm.tp = [];
expm.tp1 = 0;
/* If this is the first call to expm, fill the power of two table tp. */
if (expm.tp1 == 0) {
expm.tp1 = 1;
expm.tp[0] = 1;
for (i = 1; i < ntp; i++) expm.tp[i] = 2 * expm.tp[i - 1];
}
if (ak == 1) return 0;
/* Find the greatest power of two less than or equal to p. */
for (i = 0; i < ntp; i++)
if (expm.tp[i] > p) break;
pt = expm.tp[i - 1];
p1 = p;
r = 1;
/* Perform binary exponentiation algorithm modulo ak. */
for (j = 1; j <= i; j++) {
if (p1 >= pt) {
r = 16 * r;
r = r - Math.trunc(r / ak) * ak;
p1 = p1 - pt;
}
pt = 0.5 * pt;
if (pt >= 1) {
r = r * r;
r = r - Math.trunc(r / ak) * ak;
}
}
return r;
}
@username1565
Copy link

username1565 commented May 2, 2019

<script>
//Put this code into txt file, rename to html, open in browser and see console.log() - F12 button in your browser.

//fast modPow
function powmod( base, exp, mod ){
  if (exp == 0) return 1;
  if (exp % 2 == 0){
    return Math.pow( powmod( base, (exp / 2), mod), 2) % mod;
  }
  else {
    return (base * powmod( base, (exp - 1), mod)) % mod;
  }
}

/**
   Bailey-Borwein-Plouffe digit-extraction algorithm for pi 
    <https://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula#BBP_digit-extraction_algorithm_for_.CF.80>
*/

function run(n) {
    var partial = function(d, c) {
        var sum = 0;

        // Left sum
        var k;
        for (k = 0; k <= d - 1; k++) {
            //sum += (Math.pow(16, d - 1 - k) % (8 * k + c)) / (8 * k + c);
            sum += (powmod(16, (d - 1 - k), (8 * k + c))) / (8 * k + c);
        }

        // Right sum. This converges fast...
        var prev = undefined;
        for(k = d; sum !== prev; k++) {
            prev = sum;
            sum += Math.pow(16, d - 1 - k) / (8 * k + c);
        }

        return sum;
    };

    /**
        JavaScript's modulus operator gives the wrong
        result for negative numbers. E.g. `-2.9 % 1`
        returns -0.9, the correct result is 0.1.
    */
    var mod1 = function(x) {
        return x < 0 ? 1 - (-x % 1) : x % 1;
    };

    var s = 0;
    s +=  4 * partial(n, 1);
    s += -2 * partial(n, 4);
    s += -1 * partial(n, 5);
    s += -1 * partial(n, 6);

    s = mod1(s);
    return Math.floor(s * 16).toString(16);
}


var pi_hex = Math.PI.toString(16);
	//TEST
console.log(
				'decimal Pi: ', Math.PI
	,	'\n'+	'hexadecimal Pi: ', pi_hex
);

// Pi in hex is 3.243f6a8885a308d313198a2e037073...
pi_hex = pi_hex.split('.').join('');	//remove "."
for( digit = 0; digit < pi_hex.length ; digit++ ){
	console.log( run( digit ), ( run( digit ) === pi_hex[ digit ] ) ); // 0th hexadecimal digit of pi is the leading 3
	//check 3243f6a8885a3
}
console.log( '...' );

//test 55 hexadecimal digits
var res = '3.';
for( dig = 1; dig <= 55; dig++ ){
	res += run( dig );
}
console.log( res );

console.log( '100501-th digit: '+run( 100501 ) );

</script>

@username1565
Copy link

@getify, In my example, N-th hexadecimal digit of Pi available as run( N ).
And pi-hex-digits.js working too. N-th digit available there as getDigit(N)[0]...
But my example working good, and pi-hex-digits.js using Math.trunc.
To make pi-hex-digits.js compatible with old browsers, need to add there polyfill for Math.trunc() :

//Polyfill Math.trunc:
if (!Math.trunc) {
	Math.trunc = function(v) {
		v = +v;
		return (v - v % 1)   ||   ( ( !isFinite(v) || (v === 0) ) ? v : ( ( v < 0 ) ? -0 : 0 ) );

		// returns:
		//  0        ->  0
		// -0        -> -0
		//  0.2      ->  0
		// -0.2      -> -0
		//  0.7      ->  0
		// -0.7      -> -0
		//  Infinity ->  Infinity
		// -Infinity -> -Infinity
		//  NaN      ->  NaN
		//  null     ->  0
	};
}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment