Calculate digits of pi

From CodeCodex

Related content:

The code examples below show how to calculate digits of pi in different programming languages.

Implementations[edit]

Ada[edit]

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;

BASIC[edit]

'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 

C[edit]

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[edit]

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[edit]

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));
}

Erlang[edit]

-module(codex).
-export([pi_digits/1]).
-define(SCALE, 10000).
-define(ARRINIT, 2000).

pi_digits(Digits) ->
    lists:foreach(fun(I) -> put(I, ?ARRINIT) end, lists:seq(0, Digits + 1)),
    pi_digits(Digits, 0, 0).

pi_digits(I, _Carry, _Sum) when  I =< 0 ->
    io:nl();
pi_digits(I, Carry, Sum) ->
    Sum2 = pi_digits(I, Sum),
    io:format("~4..0B", [Carry + Sum2 div ?SCALE]),     % zero pad a number
    pi_digits(I - 14, Sum2 rem ?SCALE, Sum2).

pi_digits(J, Sum) when J =< 0 -> Sum;
pi_digits(J, Sum) ->
    Sum2 = Sum * J + ?SCALE * get(J),
    put(J, Sum2 rem (J * 2 - 1)),
    pi_digits(J - 1, Sum2 div (J * 2 - 1)).

Java[edit]

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[edit]

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[edit]

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[edit]

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

Ruby[edit]

SCALE = 10000
ARRINIT = 2000
def pi_digits(digits)
  pi = ""
  carry = 0
  arr = Array.new(digits + 1, ARRINIT)
  digits.step(1, -14) do |i|
    sum = 0
    i.downto(1) do |j|
      sum = sum * j + SCALE * arr[j]
      arr[j] = sum % (j * 2 - 1)
      sum /= j * 2 - 1
    end
    pi << sprintf("%04d", carry + sum / SCALE)
    carry = sum % SCALE;
  end
  pi
end

puts pi_digits(1000)

Seed7[edit]

$ 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[edit]

# 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 {}
}