Talk:Pi with the BBP formula (Python)

From LiteratePrograms

It'd be interesting to measure the performance of a straight C port of the integer version, especially if 64-bit integers were used on 64-bit hardware. I can't try that out right now; maybe soon... Fredrik 14:46, 30 October 2006 (PST)

Interesting idea. I'd like to try a C implementation. I have some 64-bit hardware lying around at work that I could try it on. Deco 21:53, 1 November 2006 (PST)

Gross error in Python code

I noticed that when I run this in Python 2.5.1 on Windows, that in the example pibbp.py for function S, that the Left hand sum, s is always zero. This is because pow(16, n-k, r) is an integer and pow(16, n-k, r)/r in all-integer arithmetic will always be zero. To fix this problem, I added a float() function around the call to modular pow().--75.37.8.1 20:29, 20 September 2008 (EDT)

Excellent. Thanks! -- Derek Ross | Talk 00:50, 21 September 2008 (EDT)

Code in Java

FYI: There are a few implementation nits in Java. In this implementation, I simply did a printf() did not bother with main args or with the sprintf()-like function, which requires java.util.Formatter .--75.36.169.116 02:38, 14 October 2008 (EDT)

import java.lang.Math;
class bbp {
	public static void main(String arg[]){
		bbp mybbp = new bbp();
		mybbp.pi(1);
	}
	void pi(long n){
		double x;
		String ret = new String();
		n -= 1;
		x = (4*S(1, n) - 2*S(4, n) - S(5, n) - S(6, n)) % 1.0;
		// Java floating modulo works differently than Python
		// Convert explicitly to positive number
		if(x < 0.0){
			x = x + 1.0;
		}
		System.out.printf("%014x\n" , (long)(x * Math.pow(16,14) ) );
		return;
	}
	double S(long j, long n){
		double s, t, newt;
		long k, r;
		// Left sum
		s = 0.0;
		k = 0;
		while (k <= n){
			r = 8*k+j;
			s = (s + (double)modpow(16, n-k, r) / r) % 1.0;
			k += 1;
		}
		// Right sum
		t = 0.0;
		k = n + 1;
		while (true){
			// This pow() is double, not long
			newt = t + Math.pow(16, n-k) / (8*k+j);
			// Iterate until t no longer changes
			if (t == newt)
				break;
			else
				t = newt;
			k += 1;
		}
		return s + t;
	}
	public static final long modpow(long base, long exponent, long modulus) {
		long result = 1;
		while (exponent > 0) {
			if ((exponent & 1) == 1) {
				result = (result * base) % modulus;
			}
			exponent = exponent >> 1;
			base = (base * base) % modulus;
		}
		return result;
	}
}