Digits of pi calculation
From CodeCodex
| Image:Lightbulb.png Related content: |
| <DPL>
category=Math shownamespace=false </DPL> |
Contents |
[edit] Implementations
[edit] Ada
------------------------------------------------------------------
--| Pi to 2400 places
--| See also http://home.att.net/~srschmitt/pi-spigot.html
------------------------------------------------------------------
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
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;
[edit] 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) gosub addsub 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 addsub: '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
[edit] 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;
}
[edit] 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));
}
[edit] 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));
}
}
[edit] 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()
[edit] 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;
}
?>
[edit] 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
[edit] 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: [1]
[edit] Common Lisp
This was taken from the debian shootout [2]
(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)
