Difference between revisions of "Calculate digits of pi"

Implementations

```------------------------------------------------------------------
--| Pi to 2400 places
------------------------------------------------------------------

procedure Pi is

a, b, c, d, e, g: Integer;
f: array (0 .. 8400) of Integer;

procedure Dec(n : in out Integer) is
begin
n := n - 1;
end Dec;

procedure Print_Digits(i : integer) is
s : String(1 .. 4);
begin
put(s, i);
for i in 1 .. s'Length loop
if s(i) = ' ' then s(i) := '0';
end if;
end loop;
put(s);
end Print_Digits;

begin
a := 10000; b := 0; c := 4200; e := 0;
for b in 0 .. c loop f(b) := a / 5; end loop;
while c > 0 loop
d := 0; g := c*2; b := c;
while b > 0 loop
d := d + f(b) * a; Dec(g); f(b) := d mod g;
d := d / g; Dec(g); Dec(b);
if b > 0 then d := d * b; end if;
end loop;
c := c - 14;
Print_Digits(e + d / a);
e := d mod a;
end loop;
end;
```

BASIC

```'PI.BAS 'Prints PI to number of decimal places desired. 19830623tac.
'        Ref: American Mathematical Monthly vol#45,1938 p657-667.
'        Ref: Best of Micro vol.#1 p85-86.
'
'
'             *       ----                      ----
'        *****         \     16(-1)^(k+1)        \         4(-1)^(k+1)
'       * * *   ---     >   ------------    --    >    ----------------
'         * *   ---    /    (2k-1)5^(2k-1)       /      (2k-1)239^(2k-1)
'         * *         ----                      ----
'        *   *        k=1                       k=1
'
'

main:
TEN = 100
INPUT "Enter number of digits wanted";SIZE:
SIZE = (SIZE + 1) OR 1
PLACES = SIZE - 1

IF SIZE > 200 THEN TEN = 10 ELSE SIZE = (SIZE + 1) / 2

DIM POWER(SIZE),TERM(SIZE),RESULT(SIZE)
K(1) = 25 :K(2) = 239   'Constants
C(1) = 5 :C(2) = 239

'Main loop.
FOR PASS = 1 TO 2

Init:
FOR L = 0 TO SIZE
POWER(L) = 0 :TERM(L) = 0 :IF PASS = 1 THEN RESULT = 0
NEXT L
POWER(0) = 16 / PASS ^2
MODE = 0
DIVISOR = C(PASS)
gosub divide
XPONENT = 1
SIGN = 3 - 2 * PASS
'

copy:
FOR L = 0 TO SIZE :TERM(L) = POWER(L) :NEXT L
'
MODE = 1
DIVISOR = XPONENT
gosub divide

MODE = (SIGN < 0) + ABS(SIGN > 0)

XPONENT = XPONENT + 2
SIGN = -SIGN

MODE = 0
DIVISOR = K(PASS)
gosub divide

IF PASS = 2 THEN gosub divide

IF ZERO <> 0 THEN goto copy

NEXT PASS
'
'Print result.
PRINT
PRINT
PRINT "The value of PI to";PLACES;" decimal places"
PRINT
PRINT RESULT(0);".";
FOR L = 1 TO SIZE
IF (TEN = 100)*(RESULT(L) < 10) THEN
PRINT "0";
PRINT USING "#";RESULT(L);
GOTO nxt
end if

PRINT USING "##";RESULT(L);

nxt:
NEXT L
PRINT
end

divide:
'Division subroutine.   mode0 =  Power(x)/Div else mode1 =  Term(x)/Div.
DIGIT = 0 :ZERO = 0
FOR L1 = 0 TO SIZE
DIGIT = DIGIT + TERM(L1) * MODE + POWER(L1) - POWER(L1) * MODE
QUOTIENT = INT(DIGIT / DIVISOR)
RESIDUE = DIGIT MOD DIVISOR
ZERO = ZERO OR (QUOTIENT + RESIDUE)
IF MODE THEN TERM(L1) = QUOTIENT ELSE POWER(L1) = QUOTIENT
DIGIT = RESIDUE * TEN
NEXT L1
MODE = 0
return

'Add/Subtract subroutine.       mode - 1 =  Subtract  else mode1 = Add.
CARRY = 0
FOR L1 = SIZE TO 0 STEP  - 1
SUM = RESULT(L1) + TERM(L1) * MODE + CARRY * MODE :CARRY = 0
IF (MODE = 1)*(SUM < TEN) + (MODE =  - 1) * (SUM >= 0) THEN goto t1
SUM = SUM + MODE * -TEN :CARRY = 1

t1:
RESULT(L1) = SUM
NEXT L1
MODE = 0
return
```

C

Complete C99 code:

```#include "stdio.h"
#include "stdlib.h"

#define SCALE 10000
#define ARRINIT 2000

void pi_digits(int digits) {
int carry = 0;
int arr[digits + 1];
for (int i = 0; i <= digits; ++i)
arr[i] = ARRINIT;
for (int i = digits; i > 0; i-= 14) {
int sum = 0;
for (int j = i; j > 0; --j) {
sum = sum * j + SCALE * arr[j];
arr[j] = sum % (j * 2 - 1);
sum /= j * 2 - 1;
}
printf("%04d", carry + sum / SCALE);
carry = sum % SCALE;
}
}

int main(int argc, char** argv) {
int n = argc == 2 ? atoi(argv[1]) : 100;
pi_digits(n);

return 0;
}
```

Common Lisp

This was taken from the debian shootout [1]

```(defconstant +digits-per-line+ 70)

(defun make-digit-generator()
(declare (optimize speed))
(let ((zq 1)
(zr 0)
(zt 1)
(k 0)
(4k+2 2)
(2k+1 1))
(declare (type integer zq zr zt)
(type fixnum k 4k+2 2k+1))
(labels ((extract (j)
(declare (type fixnum j))
(the (integer 0 9) (floor (+ (* zq j) zr) zt)))
(compose (aq ar at bq br bt)
(declare (type integer aq ar at bq br bt))
(setq zq (* aq bq)
zr (+ (* aq br) (* ar bt))
zt (* at bt))))
#'(lambda ()
(let ((y (extract 3)))
(declare (type (integer 0 9) y))
(loop while (not (= y (extract 4)))
do (compose zq zr zt (incf k) (incf 4k+2 4) (incf 2k+1 2))
(setf y (extract 3)))
(compose 10 (* -10 y) 1 zq zr zt)
y)))))

;;; It should be trivial to change this and make it store the result in a string-stream
;;; instead of just printing it to screen.
(defun spigot (digits)
(declare (type fixnum digits))
(let ((digits-printed 0)
(next-digit (make-digit-generator)))
(declare (type fixnum digits-printed))
(loop while (plusp digits)
do (if (>= digits +digits-per-line+)
(progn (loop repeat +digits-per-line+
do (format t "~d" (funcall next-digit)))
(incf digits-printed +digits-per-line+))
(progn (loop repeat digits
do (format t "~d" (funcall next-digit)))
(loop repeat (- +digits-per-line+ digits)
do (format t " "))
(incf digits-printed digits)))
(format t "~a:~d~%" #\Tab digits-printed)
(decf digits +digits-per-line+))))

;;; Usage:
(spigot 1000)

```

D

```import std.stdio: writefln;
import std.conv: toInt;
import std.string: format;

string pi_digits(int digits) {
const int SCALE = 10000;
const int ARRINIT = 2000;
string result;

int carry;
auto arr = new int[digits + 1];
arr[] = ARRINIT;
for (int i = digits; i > 0; i-= 14) {
int sum = 0;
for (int j = i; j > 0; --j) {
sum = sum * j + SCALE * arr[j];
arr[j] = sum % (j * 2 - 1);
sum /= j * 2 - 1;
}
result ~= format("%04d", carry + sum / SCALE);
carry = sum % SCALE;
}
return result;
}

void main(string[] args) {
int n = args.length == 2 ? toInt(args[1]) : 100;
writefln(pi_digits(n));
}
```

Java

This code is very similar to the D one:

```public class PiDigits {
private static final int SCALE = 10000;
private static final int ARRINIT = 2000;

public static String pi_digits(int digits){
StringBuffer pi = new StringBuffer();
int[] arr = new int[digits + 1];
int carry = 0;

for (int i = 0; i <= digits; ++i)
arr[i] = ARRINIT;

for (int i = digits; i > 0; i-= 14) {
int sum = 0;
for (int j = i; j > 0; --j) {
sum = sum * j + SCALE * arr[j];
arr[j] = sum % (j * 2 - 1);
sum /= j * 2 - 1;
}

pi.append(String.format("%04d", carry + sum / SCALE));
carry = sum % SCALE;
}
return pi.toString();
}

public static void main(String[] args) {
int n = args.length > 0 ? Integer.parseInt(args[0]) : 100;
System.out.println(PiDigits.pi_digits(n));
}
}
```

OCaml

```let scale = 10000 and n = 2800
let a = Array.make (n+1) 2000

let rec g j sum =
if j < 1 then sum else
let sum = sum * j + scale * a.(j) in
a.(j) <- sum mod (j * 2 - 1);
g (j - 1) (sum / (j * 2 - 1))

let rec f i carry =
if i = 0 then () else
let sum = g i 0 in
Printf.printf "%04d" (carry + sum / scale);
f (i - 14) (sum mod scale)

let () =
f n 0;
print_newline()
```

PHP

This PHP script is equivalent to the above C program. (This appears to give an incorrect result, at least in some versions of PHP, because of floating point/integer conversions.)

```<?php
\$scale = 10000;
\$maxarr = 2800;
\$arrinit = 2000;
\$carry = 0;
for(\$i=0;\$i<=\$maxarr;\$i++)
\$arr[\$i] = \$arrinit;
for(\$i=\$maxarr;\$i>0;\$i=\$i-14)
{
\$sum = 0;
for(\$j=\$i;\$j>0;\$j--)
{
\$sum = (\$sum * \$j) + (\$scale * \$arr[\$j]);
\$arr[\$j] = \$sum % ((\$j*2)-1);
\$sum = \$sum / ((\$j*2)-1);
}
printf("%04d",(\$carry + (\$sum / \$scale)));
\$carry = \$sum % \$scale;
}
?>
```

Python

This python program calculates the digits of pi. See below for a more straightforward solution.

```#!/usr/bin/python

from __future__ import division
import copy
import sys
from types import IntType,FloatType,LongType,ListType,TupleType,InstanceType,ClassType

# set error output to stdout instead
#__errout = sys.stderr
__errout = sys.stdout

#
# Define class contfrac
# Continued Fraction object version 1
#
class contfrac:

def __init__(self,numerator,denominator=0,limit=256,epsilon=2e-9):
self.classname="contfracv1"
# if called by A=contfrac(3.1416)
if denominator==0:
if type(numerator)==IntType or type(numerator)==LongType:
self.num=numerator
self.den=1
self.cf=[numerator]
elif type(numerator)==FloatType:
self.cf = decimal2contfrac(numerator,limit,epsilon)
(self.num,self.den) = contfrac2rational(self.cf)
elif type(numerator)==ListType:
self.cf=copy.copy(numerator)
self.num,self.den=contfrac2rational(numerator)
elif is_contfrac(numerator):
self.num=numerator.num
self.den=numerator.den
self.cf=copy.copy(numerator.cf)
else:
# if all else fails its the devil's number
self.num=666
self.den=1
self.cf=[666]
# if called by A=(2,3)
else:
self.cf = rational2contfrac(long(numerator),long(denominator))
(self.num,self.den) = contfrac2rational(self.cf)

def __len__(self):
return len(self.cf)

def lid(cfobj):
return lazy_identity(cfobj)

#
# Define class lazy_identity
# Lazy connector for contfrac objects
#
class lazy_identity:

def __init__(self,cf):
self.classname="lazy_identity"
self.cfobj=cf
self.n=0
self.a=0
self.b=1
self.c=1
self.d=0

def egest(self,debugflag=False):
# First check if we can egest
# if we can then egest a digit
# otherwise loop until we can
if debugflag : print >> __errout,"a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
while 1 :
if self.c == 0 and self.d == 0 :
break
if self.c <> 0 and self.d <> 0 and (self.a//self.c) == (self.b//self.d) :
if debugflag : print >> __errout,"We can egest"
break
# Now we must ingest one term from cfobj
if debugflag : print >> __errout,"n=%d len(cfobj)=%d" % (self.n,len(self.cfobj))
if self.n < len(self.cfobj) :
p=self.cfobj.cf[self.n]
if debugflag : print >> __errout,"term ingested is ",p
self.n=self.n+1
# Now calculate the new a,b,c,d
self.a,self.b,self.c,self.d=self.b,self.a+self.b*p,self.d,self.c+self.d*p
else :
if debugflag : print >> __errout,"term ingested is INF"
self.n=self.n+1
self.a,self.b,self.c,self.d=self.b,self.b,self.d,self.d
if debugflag : print >> __errout,"after ingestion a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
# out of the while loop
if self.c== 0 and self.d == 0 :
if debugflag : print >> __errout,"egest the number INF"
return 'INF'
else:
q=self.a//self.c
if debugflag : print >> __errout,"egest the number ",q
# Now calculate the new a,b,c,d
self.a,self.b,self.c,self.d=self.c,self.d,self.a-self.c*q,self.b-self.d*q
if debugflag : print >> __errout,"after egestion a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
return q

#
# Define class lazy_rcf_4onpi
# Lazy connector for 4/pi in regular continued fractions format
# [(1,1),(3,4),(5,9),(7,16),(9,25),(11,36),...]
#        (2n+1,(n+1)^2)
#
class lazy_rcf_4onpi:

def __init__(self):
self.classname="lazy_rcf_4onpi"
self.n=0

def egest(self,debugflag=False):
a=  self.n * 2 + 1
b= (self.n + 1) * (self.n + 1)
self.n = self.n + 1
return (a,b)

#
# Define class lazy_rcf2cf
# Lazy connector for turning regular continued fractions
# into standardize continued fractions
#
class lazy_rcf2cf:

def __init__(self,rcf):
self.classname="lazy_rcf2cf"
self.cfobj=rcf
self.a=0
self.b=1
self.c=1
self.d=0

def egest(self,debugflag=False):
# First check if we can egest
# if we can then egest a digit
# otherwise loop until we can
if debugflag : print >> __errout,"a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
while 1 :
if self.c == 0 and self.d == 0 :
break
if self.c <> 0 and self.d <> 0 and (self.a//self.c) == (self.b//self.d) :
if debugflag : print >> __errout,"We can egest"
break
# Now we must ingest one term from cfobj
(p,q)=self.cfobj.egest()
if debugflag : print >> __errout,"term ingested is ",p
# Now calculate the new a,b,c,d
self.a,self.b,self.c,self.d=self.b*q,self.a+self.b*p,self.d*q,self.c+self.d*p
if debugflag : print >> __errout,"after ingestion a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
# out of the while loop
if self.c== 0 and self.d == 0 :
if debugflag : print >> __errout,"egest the number INF"
return 'INF'
else:
q=self.a//self.c
if debugflag : print >> __errout,"egest the number ",q
# Now calculate the new a,b,c,d
self.a,self.b,self.c,self.d=self.c,self.d,self.a-self.c*q,self.b-self.d*q
if debugflag : print >> __errout,"after egestion a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
return q

#
# Define class lazy_decimal
# Lazy connector for contfrac objects
#
class lazy_decimal:

def __init__(self,lazycf):
self.classname="lazy_decimal"
self.cfobj=lazycf
self.a=0
self.b=1
self.c=1
self.d=0

def egest(self,debugflag=False):
# First check if we can egest
# if we can then egest a digit
# otherwise loop until we can
if debugflag : print >> __errout,"a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
while 1 :
if self.c == 0 and self.d == 0 :
break
if self.c <> 0 and self.d <> 0 and (self.a//self.c) == (self.b//self.d) :
if debugflag : print >> __errout,"We can egest"
break
# Now we must ingest one term from cfobj
p=self.cfobj.egest()
if p == 'INF' :
if debugflag : print >> __errout,"term ingested is INF"
self.a,self.b,self.c,self.d=self.b,self.b,self.d,self.d
else :
if debugflag : print >> __errout,"term ingested is ",p
# Now calculate the new a,b,c,d
self.a,self.b,self.c,self.d=self.b,self.a+self.b*p,self.d,self.c+self.d*p
if debugflag : print >> __errout,"after ingestion a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
# out of the while loop
if self.c== 0 and self.d == 0 :
if debugflag : print >> __errout,"egest the number INF"
return 'INF'
else:
r=self.a//self.c
if debugflag : print >> __errout,"egest the number ",r
# Now calculate the new a,b,c,d
self.a,self.b,self.c,self.d=self.c,self.d,self.a-self.c*r,self.b-self.d*r
if debugflag : print >> __errout,"after egestion a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
# Now calculate the egestion of r2=0 and s2=10
# new_a = c * s2
# new_b = d * s2
# new_c = a - c * r2
# new_d = b - d * r2
self.a,self.b,self.c,self.d=self.c*10,self.d*10,self.a,self.b
if debugflag : print >> __errout,"after renormalisation a=%d,b=%d,c=%d,d=%d" % (self.a,self.b,self.c,self.d)
return r

#
# Define class lazy_base
# Lazy base for lazy operators
#
class lazy_base:

def egest(self,debugflag=False):
# First check if we can egest
# if we can then egest a ter
```

For something a little more compact, here's a direct translation of the PHP script above (which does not suffer from the floating point/integer issues).

```from sys import stdout

scale = 10000
maxarr = 2800
arrinit = 2000
carry = 0
arr = [arrinit] * (maxarr + 1)
for i in xrange(maxarr, 1, -14):
total = 0
for j in xrange(i, 0, -1):
total = (total * j) + (scale * arr[j])
arr[j] = total % ((j * 2) - 1)
total = total / ((j * 2) - 1)
stdout.write("%04d" % (carry + (total / scale)))
carry = total % scale
```

Seed7

```\$ include "seed7_05.s7i";

# Newtons formula for PI is:
# PI / 2 = sum_n_from_0_to_inf(n! / (2 * n + 1)!!)

# This can be written as:
# PI / 2 = 1 + 1/3 * (1 + 2/5 * (1 + 3/7 * (1 + 4/9 * (1 + ... ))))

# This algorithm puts 2 * 1000 on the right side and computes everything from inside out.

const integer: SCALE is 10000;
const integer: MAXARR is 3500;
const integer: ARRINIT is 2000;

const proc: main is func
local
var integer: i is 0;
var integer: j is 0;
var integer: carry is 0;
var array integer: arr is MAXARR times ARRINIT;
var integer: sum is 0;
begin
for i range MAXARR downto 1 step 14 do
sum := 0;
for j range i downto 1 do
sum := sum*j + SCALE*arr[j];
arr[j] := sum rem pred(j*2);
sum := sum div pred(j*2);
end for;
write(replace(carry + sum div SCALE lpad 4, " ", "0"));
carry := sum rem SCALE;
end for;
end func;
```

Original source: [2]

Tcl

```# 2400 digits of Pi
# More details: A spigot algorithm for the digits of pi,
# Stanley Rabinowitz and Stan Wagon, American Mathematical Monthly,
# March 1995, pp195-203

proc pi {} {
set n 2400
set e 0
set f {}
for {set b 0} {\$b <= 8400} {incr b} {lappend f 2000}
for {set c 8400} {\$c > 0} {incr c -14} {
set d 0
set g [expr {\$c * 2}]
set b \$c
while 1 {
incr d [expr { [lindex \$f \$b] * 10000 }]
lset f \$b [expr {\$d % [incr g -1]}]
set d [expr {\$d / \$g}]
incr g -1
if {[incr b -1] == 0} break
set d [expr {\$d * \$b}]
}
puts -nonewline [format %04d [expr {\$e+\$d / 10000}]]
flush stdout
set e [expr {\$d % 10000}]
}
puts {}
}
```