Pi

From Rosetta Code
Jump to: navigation, search
Task
Pi
You are encouraged to solve this task according to the task description, using any language you may know.
Create a program to continually calculate and output the next digit of π (pi). The program should continue forever (until it is aborted by the user) calculating and outputting each digit in succession. The output should be a decimal sequence beginning 3.14159265 ...


Note: this task is about calculating pi. For information on built-in pi constants see Real constants and functions.

Related Task Arithmetic-geometric mean/Calculate Pi

Contents

[edit] Ada

Works with: Ada 2005
Library: GMP

uses same algorithm as Go solution, from http://web.comlab.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf

pi_digits.adb
with Ada.Command_Line;
with Ada.Text_IO;
with GNU_Multiple_Precision.Big_Integers;
with GNU_Multiple_Precision.Big_Rationals;
use GNU_Multiple_Precision;
 
procedure Pi_Digits is
type Int is mod 2 ** 64;
package Int_To_Big is new Big_Integers.Modular_Conversions (Int);
 
-- constants
Zero : constant Big_Integer := Int_To_Big.To_Big_Integer (0);
One : constant Big_Integer := Int_To_Big.To_Big_Integer (1);
Two : constant Big_Integer := Int_To_Big.To_Big_Integer (2);
Three : constant Big_Integer := Int_To_Big.To_Big_Integer (3);
Four : constant Big_Integer := Int_To_Big.To_Big_Integer (4);
Ten : constant Big_Integer := Int_To_Big.To_Big_Integer (10);
 
-- type LFT = (Integer, Integer, Integer, Integer
type LFT is record
Q, R, S, T : Big_Integer;
end record;
 
-- extr :: LFT -> Integer -> Rational
function Extr (T : LFT; X : Big_Integer) return Big_Rational is
use Big_Integers;
Result : Big_Rational;
begin
-- extr (q,r,s,t) x = ((fromInteger q) * x + (fromInteger r)) /
-- ((fromInteger s) * x + (fromInteger t))
Big_Rationals.Set_Numerator (Item => Result,
New_Value => T.Q * X + T.R,
Canonicalize => False);
Big_Rationals.Set_Denominator (Item => Result,
New_Value => T.S * X + T.T);
return Result;
end Extr;
 
-- unit :: LFT
function Unit return LFT is
begin
-- unit = (1,0,0,1)
return LFT'(Q => One, R => Zero, S => Zero, T => One);
end Unit;
 
-- comp :: LFT -> LFT -> LFT
function Comp (T1, T2 : LFT) return LFT is
use Big_Integers;
begin
-- comp (q,r,s,t) (u,v,w,x) = (q*u+r*w,q*v+r*x,s*u+t*w,s*v+t*x)
return LFT'(Q => T1.Q * T2.Q + T1.R * T2.S,
R => T1.Q * T2.R + T1.R * T2.T,
S => T1.S * T2.Q + T1.T * T2.S,
T => T1.S * T2.R + T1.T * T2.T);
end Comp;
 
-- lfts = [(k, 4*k+2, 0, 2*k+1) | k<-[1..]
K : Big_Integer := Zero;
function LFTS return LFT is
use Big_Integers;
begin
K := K + One;
return LFT'(Q => K,
R => Four * K + Two,
S => Zero,
T => Two * K + One);
end LFTS;
 
-- next z = floor (extr z 3)
function Next (T : LFT) return Big_Integer is
begin
return Big_Rationals.To_Big_Integer (Extr (T, Three));
end Next;
 
-- safe z n = (n == floor (extr z 4)
function Safe (T : LFT; N : Big_Integer) return Boolean is
begin
return N = Big_Rationals.To_Big_Integer (Extr (T, Four));
end Safe;
 
-- prod z n = comp (10, -10*n, 0, 1)
function Prod (T : LFT; N : Big_Integer) return LFT is
use Big_Integers;
begin
return Comp (LFT'(Q => Ten, R => -Ten * N, S => Zero, T => One), T);
end Prod;
 
procedure Print_Pi (Digit_Count : Positive) is
Z : LFT := Unit;
Y : Big_Integer;
Count : Natural := 0;
begin
loop
Y := Next (Z);
if Safe (Z, Y) then
Count := Count + 1;
Ada.Text_IO.Put (Big_Integers.Image (Y));
exit when Count >= Digit_Count;
Z := Prod (Z, Y);
else
Z := Comp (Z, LFTS);
end if;
end loop;
end Print_Pi;
 
N : Positive := 250;
begin
if Ada.Command_Line.Argument_Count = 1 then
N := Positive'Value (Ada.Command_Line.Argument (1));
end if;
Print_Pi (N);
end Pi_Digits;

output:

 3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7 9 5 0 2 8 8 4 1 9 7 1 6 9 3 9 9 3 7 5 1 0 5 8 2 0 9 7 4 9 4 4 5 9 2 3 0 7 8 1 6 4 0 6 2 8 6 2 0 8 9 9 8 6 2 8 0 3 4 8 2 5 3 4 2 1 1 7 0 6 7

[edit] AutoHotkey

Library: MPL

Could be optimized with Ipp functions, but runs fast enough for me as-is. Does not work in AHKLx64.

#NoEnv
#SingleInstance, Force
SetBatchLines, -1
#Include mpl.ahk
dot:=".", i:=0
, MP_SET(q, "1")
, MP_SET(r, "0")
, MP_SET(t, "1")
, MP_SET(k, "1")
, MP_SET(n, "3")
, MP_SET(l, "3")
, MP_SET(ONE, "1")
, MP_SET(TWO, "2")
, MP_SET(THREE, "3")
, MP_SET(FOUR, "4")
, MP_SET(SEVEN, "7")
, MP_SET(TEN, "10")
 
Loop
{
MP_MUL(q4, q, FOUR)
, MP_ADD(q4r, q4, r)
, MP_SUB(q4rt, q4r, t)
, MP_MUL(tn, t, n)
If (MP_CMP(q4rt,tn) = -1)
{
s := MP_DEC(n) . dot
OutputDebug %s%
dot := ""
, i++
, MP_MUL(tn, t, n)
, MP_SUB(rtn, r, tn)
, MP_MUL(nr, rtn, TEN)
, MP_MUL(q3, q, THREE)
, MP_ADD(q3r, q3, r)
, MP_DIV(q3rt, remainder, q3r, t)
, MP_SUB(q3rtn, q3rt, n)
, MP_MUL(n, q3rtn, TEN)
, MP_MUL(tmp, q, TEN)
, MP_CPY(q, tmp)
, MP_CPY(r, nr)
}
Else
{
MP_MUL(q2, q, TWO)
, MP_ADD(q2r, q2, r)
, MP_MUL(nr, q2r, l)
, MP_MUL(k7, k, SEVEN)
, MP_ADD(k72, k7, TWO)
, MP_MUL(qk, q, k72)
, MP_MUL(rl, r, l)
, MP_ADD(qkrl, qk, rl)
, MP_MUL(tl, t, l)
, MP_DIV(nn, remainder, qkrl, tl)
, MP_MUL(tmp, q, k)
, MP_CPY(q, tmp)
, MP_MUL(tmp, t, l)
, MP_CPY(t, tmp)
, MP_ADD(tmp, l, TWO)
, MP_CPY(l, tmp)
, MP_ADD(tmp, k, ONE)
, MP_CPY(k, tmp)
, MP_CPY(n, nn)
, MP_CPY(r, nr)
}
}

[edit] ALGOL 68

Translation of: Pascal
Note: This specimen retains the original Pascal coding style of code.
Works with: ALGOL 68 version Revision 1 - no extensions to language used.
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny.

This codes uses 33 decimals places as a test case. Performance is O(2) based on the number of decimal places required.

#!/usr/local/bin/a68g --script #
 
INT base := 10;
 
MODE YIELDINT = PROC(INT)VOID;
PROC gen pi digits = (INT decimal places, YIELDINT yield)VOID:
BEGIN
INT nine = base - 1;
INT nines := 0, predigit := 0; # First predigit is a 0 #
[decimal places*10 OVER 3]#LONG# INT digits; # We need 3 times the digits to calculate #
FOR place FROM LWB digits TO UPB digits DO digits[place] := 2 OD; # Start with 2s #
FOR place TO decimal places + 1 DO
INT digit := 0;
FOR i FROM UPB digits BY -1 TO LWB digits DO # Work backwards #
INT x := #SHORTEN#(base*digits[i] + #LENG# digit*i);
digits[i] := x MOD (2*i-1);
digit := x OVER (2*i-1)
OD;
digits[LWB digits] := digit MOD base; digit OVERAB base;
nines :=
IF digit = nine THEN
nines + 1
ELSE
IF digit = base THEN
yield(predigit+1); predigit := 0 ;
FOR repeats TO nines DO yield(0) OD # zeros #
ELSE
IF place NE 1 THEN yield(predigit) FI; predigit := digit;
FOR repeats TO nines DO yield(nine) OD
FI;
0
FI
OD;
yield(predigit)
END;
 
main:(
INT feynman point = 762; # feynman point + 4 is a good test case #
# the 33rd decimal place is a shorter tricky test case #
INT test decimal places = UPB "3.1415926.......................502"-2;
 
INT width = ENTIER log(base*(1+small real*10));
 
# iterate throught the digits as they are being found #
# FOR INT digit IN # gen pi digits(test decimal places#) DO ( #,
## (INT digit)VOID: (
printf(($n(width)d$,digit))
)
# OD #);
print(new line)
)

Output:

3141592653589793238462643383279502

[edit] BASIC256

Translation of: Pascal
below, and originally published by Stanley Rabinowitz in [1].
cls
 
n =1000
len = 10*n \ 4
needdecimal = true
dim a(len)
nines = 0
predigit = 0 # {First predigit is a 0}
 
for j = 1 to len
a[j-1] = 2 # {Start with 2s}
next j
 
for j = 1 to n
q = 0
for i = len to 1 step -1
# {Work backwards}
x = 10*a[i-1] + q*i
a[i-1] = x % (2*i - 1)
q = x \ (2*i - 1)
next i
a[0] = q % 10
q = q \ 10
if q = 9 then
nines = nines + 1
else
if q = 10 then
d = predigit+1: gosub outputd
if nines > 0 then
for k = 1 to nines
d = 0: gosub outputd
next k
end if
predigit = 0
nines = 0
else
d = predigit: gosub outputd
predigit = q
if nines <> 0 then
for k = 1 to nines
d = 9: gosub outputd
next k
nines = 0
end if
end if
end if
next j
print predigit
end
 
outputd:
if needdecimal then
if d = 0 then return
print d + ".";
needdecimal = false
else
print d;
end if
return

Output:

3.14159265358979323846264338327950288419716939937510582097494459230781...

[edit] BBC BASIC

[edit] BASIC version

      WIDTH 80
M% = (HIMEM-END-1000) / 4
DIM B%(M%)
FOR I% = 0 TO M% : B%(I%) = 20 : NEXT
E% = 0
L% = 2
FOR C% = M% TO 14 STEP -7
D% = 0
A% = C%*2-1
FOR P% = C% TO 1 STEP -1
D% = D%*P% + B%(P%)*&64
B%(P%) = D% MOD A%
D% DIV= A%
A% -= 2
NEXT
CASE TRUE OF
WHEN D% = 99: E% = E% * 100 + D% : L% += 2
WHEN C% = M%: PRINT ;(D% DIV 100) / 10; : E% = D% MOD 100
OTHERWISE:
PRINT RIGHT$(STRING$(L%,"0") + STR$(E% + D% DIV 100),L%);
E% = D% MOD 100 : L% = 2
ENDCASE
NEXT

[edit] Assembler version

The first 250,000 digits output have been verified.

      DIM P% 32
[OPT 2 :.pidig mov ebp,eax :.pi1 imul edx,ecx : mov eax,[ebx+ecx*4]
imul eax,100 : add eax,edx : cdq : div ebp : mov [ebx+ecx*4],edx
mov edx,eax : sub ebp,2 : loop pi1 : mov eax,edx : ret :]
 
WIDTH 80
M% = (HIMEM-END-1000) / 4
DIM B%(M%) : B% = ^B%(0)
FOR I% = 0 TO M% : B%(I%) = 20 : NEXT
E% = 0
L% = 2
FOR C% = M% TO 14 STEP -7
D% = 0
A% = C%*2-1
D% = USR(pidig)
CASE TRUE OF
WHEN D% = 99: E% = E% * 100 + D% : L% += 2
WHEN C% = M%: PRINT ;(D% DIV 100) / 10; : E% = D% MOD 100
OTHERWISE:
PRINT RIGHT$(STRING$(L%,"0") + STR$(E% + D% DIV 100),L%);
E% = D% MOD 100 : L% = 2
ENDCASE
NEXT

Output:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
65271201909145648566923460348610454326648213393607260249141273724587006606315588
17488152092096282925409171536436789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185480744623799627495673518857527
24891227938183011949129833673362440656643086021394946395224737190702179860943702
77053921717629317675238467481846766940513200056812714526356082778577134275778960
91736371787214684409012249534301465495853710507922796892589235420199561121290219
60864034418159813629774771309960518707211349999998372978049951059731732816096318
....

[edit] bc

The digits of Pi are printed 20 per line, by successively recomputing pi with higher precision. The computation is not accurate to the entire scale (for example, scale = 4; 4*a(1) prints 3.1412 instead of the expected 3.1415), so the program includes two excess digits in the scale. Fixed number of guarding digits will eventually fail because Pi can contain arbitrarily long sequence of consecutive 9s (or consecutive 0s), though for this task it might not matter in practice. The program proceeds more and more slowly but exploits bc's unlimited precision arithmetic.

The program uses three features of GNU bc: long variable names, # comments (for the #! line), and the print command (for zero padding).

Library: bc -l
Works with: GNU bc
Works with: OpenBSD bc
#!/usr/bin/bc -l
 
scaleinc= 20
 
define zeropad ( n ) {
auto m
for ( m= scaleinc - 1; m > 0; --m ) {
if ( n < 10^m ) {
print "0"
}
}
return ( n )
}
 
wantscale= scaleinc - 2
scale= wantscale + 2
oldpi= 4*a(1)
scale= wantscale
oldpi= oldpi / 1
oldpi
while( 1 ) {
wantscale= wantscale + scaleinc
scale= wantscale + 2
pi= 4*a(1)
scale= 0
digits= ((pi - oldpi) * 10^wantscale) / 1
zeropad( digits )
scale= wantscale
oldpi= pi / 1
}

Output:

3.141592653589793238
46264338327950288419
71693993751058209749
44592307816406286208
99862803482534211706
79821480865132823066
47093844609550582231
72535940812848111745
02841027019385211055
59644622948954930381
96442881097566593344
61284756482337867831
65271201909145648566
92346034861045432664
82133936072602491412
73724587006606315588
17488152092096282925
40917153643678925903
60011330530548820466
52138414695194151160
94330572703657595919
....

[edit] Bracmat

Translation of: Icon_and_Unicon
  ( pi
= f,q r t k n l,first
.  !arg:((=?f),?q,?r,?t,?k,?n,?l)
& yes:?first
& whl
' ( 4*!q+!r+-1*!t+-1*!n*!t:<0
& f$!n
& (  !first:yes
& f$"."
& no:?first
|
)
& "compute and update variables for next cycle"
& 10*(!r+-1*!n*!t):?nr
& div$(10*(3*!q+!r).!t)+-10*!n:?n
& !q*10:?q
& !nr:?r
| "compute and update variables for next cycle"
& (2*!q+!r)*!l:?nr
& div$(!q*(7*!k+2)+!r*!l.!t*!l):?nn
& !q*!k:?q
& !t*!l:?t
& !l+2:?l
& !k+1:?k
& !nn:?n
& !nr:?r
)
)
& pi$((=.put$!arg),1,0,1,1,3,3)

Output:

3.1415926535897932384626433832795028841971693993751058209749445923078164062
862089986280348253421170679821480865132823066470938446095505822317253594081
284811174502841027019385211055596446229489549303819644288109756659334461284
756482337867831652712019091456485669234603486104543266482133936072602491412
73724587006606315588174881520...

[edit] C

There are many ways to do this, with quite different performance profiles. A simple measurement of 6 programs:

Digits Spigot 1 Spigot 2 Machin 1 Machin 2 AGM Chudnovsky
1,000 0.008 0.009 0.001 0.001 0.000 0.000
10,000 0.402 0.589 0.020 0.016 0.003 0.002
100,000 39.400 85.600 1.740 1.480 0.084 0.002
1,000,000 177.900 156.800 1.474 0.333
10,000,000 25.420 5.715


Using Machin's formula. The "continuous printing" part is silly: the algorithm really calls for a preset number of digits, so the program repeatedly calculates Pi digits with increasing length and chop off leading digits already displayed. But it's still faster than the unbounded Spigot method by an order of magnitude, at least for the first 100k digits.

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
 
mpz_t tmp1, tmp2, t5, t239, pows;
void actan(mpz_t res, unsigned long base, mpz_t pows)
{
int i, neg = 1;
mpz_tdiv_q_ui(res, pows, base);
mpz_set(tmp1, res);
for (i = 3; ; i += 2) {
mpz_tdiv_q_ui(tmp1, tmp1, base * base);
mpz_tdiv_q_ui(tmp2, tmp1, i);
if (mpz_cmp_ui(tmp2, 0) == 0) break;
if (neg) mpz_sub(res, res, tmp2);
else mpz_add(res, res, tmp2);
neg = !neg;
}
}
 
char * get_digits(int n, size_t* len)
{
mpz_ui_pow_ui(pows, 10, n + 20);
 
actan(t5, 5, pows);
mpz_mul_ui(t5, t5, 16);
 
actan(t239, 239, pows);
mpz_mul_ui(t239, t239, 4);
 
mpz_sub(t5, t5, t239);
mpz_ui_pow_ui(pows, 10, 20);
mpz_tdiv_q(t5, t5, pows);
 
*len = mpz_sizeinbase(t5, 10);
return mpz_get_str(0, 0, t5);
}
 
int main(int c, char **v)
{
unsigned long accu = 16384, done = 0;
size_t got;
char *s;
 
mpz_init(tmp1);
mpz_init(tmp2);
mpz_init(t5);
mpz_init(t239);
mpz_init(pows);
 
while (1) {
s = get_digits(accu, &got);
 
/* write out digits up to the last one not preceding a 0 or 9*/
got -= 2; /* -2: length estimate may be longer than actual */
while (s[got] == '0' || s[got] == '9') got--;
 
printf("%.*s", (int)(got - done), s + done);
free(s);
 
done = got;
 
/* double the desired digits; slows down at least cubically */
accu *= 2;
}
 
return 0;
}

[edit] C#

Translation of: Java

using System;
using System.Numerics;
 
namespace PiCalc {
internal class Program {
private readonly BigInteger FOUR = new BigInteger(4);
private readonly BigInteger SEVEN = new BigInteger(7);
private readonly BigInteger TEN = new BigInteger(10);
private readonly BigInteger THREE = new BigInteger(3);
private readonly BigInteger TWO = new BigInteger(2);
 
private BigInteger k = BigInteger.One;
private BigInteger l = new BigInteger(3);
private BigInteger n = new BigInteger(3);
private BigInteger q = BigInteger.One;
private BigInteger r = BigInteger.Zero;
private BigInteger t = BigInteger.One;
 
public void CalcPiDigits() {
BigInteger nn, nr;
bool first = true;
while (true) {
if ((FOUR*q + r - t).CompareTo(n*t) == -1) {
Console.Write(n);
if (first) {
Console.Write(".");
first = false;
}
nr = TEN*(r - (n*t));
n = TEN*(THREE*q + r)/t - (TEN*n);
q *= TEN;
r = nr;
} else {
nr = (TWO*q + r)*l;
nn = (q*(SEVEN*k) + TWO + r*l)/(t*l);
q *= k;
t *= l;
l += TWO;
k += BigInteger.One;
n = nn;
r = nr;
}
}
}
 
private static void Main(string[] args) {
new Program().CalcPiDigits();
}
}
}

Adopted Version:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
 
namespace EnumeratePi {
class Program {
private const int N = 60;
private const string ZS = " +-";
static void Main() {
Console.WriteLine("Digits of PI");
Console.WriteLine(new string('=', N + 13));
 
Console.WriteLine("Decimal  : {0}", string.Concat(PiDigits(10).Take(N).Select(_ => _.ToString("d"))));
Console.WriteLine("Binary  : {0}", string.Concat(PiDigits(2).Take(N).Select(_ => _.ToString("d"))));
Console.WriteLine("Quaternary : {0}", string.Concat(PiDigits(4).Take(N).Select(_ => _.ToString("d"))));
Console.WriteLine("Octal  : {0}", string.Concat(PiDigits(8).Take(N).Select(_ => _.ToString("d"))));
Console.WriteLine("Hexadecimal: {0}", string.Concat(PiDigits(16).Take(N).Select(_ => _.ToString("x"))));
Console.WriteLine("Alphabetic : {0}", string.Concat(PiDigits(26).Take(N).Select(_ => (char) ('A' + _))));
Console.WriteLine("Fun  : {0}", string.Concat(PiDigits(ZS.Length).Take(N).Select(_ => ZS[(int)_])));
 
Console.WriteLine("Nibbles  : {0}", string.Concat(PiDigits(0x10).Take(N/2).Select(_ => string.Format("{0:x1} ", _))));
Console.WriteLine("Bytes  : {0}", string.Concat(PiDigits(0x100).Take(N/3).Select(_ => string.Format("{0:x2} ", _))));
Console.WriteLine("Words  : {0}", string.Concat(PiDigits(0x10000).Take(N/5).Select(_ => string.Format("{0:x4} ", _))));
Console.WriteLine("Dwords  : {0}", string.Concat(PiDigits(0x100000000).Take(N/9).Select(_ => string.Format("{0:x8} ", _))));
 
Console.WriteLine(new string('=', N + 13));
Console.WriteLine("* press any key to exit *");
Console.ReadKey();
}
 
/// <summary>Enumerates the digits of PI.</summary>
/// <param name="b">Base of the Numeral System to use for the resulting digits (default = Base.Decimal (10)).</param>
/// <returns>The digits of PI.</returns>
static IEnumerable<long> PiDigits(long b = 10) {
BigInteger
k = 1,
l = 3,
n = 3,
q = 1,
r = 0,
t = 1
;
 
// skip integer part
var nr = b * (r - t * n);
n = b * (3 * q + r) / t - b * n;
q *= b;
r = nr;
 
for (; ; ) {
var tn = t * n;
if (4 * q + r - t < tn) {
yield return (long)n;
nr = b * (r - tn);
n = b * (3 * q + r) / t - b * n;
q *= b;
} else {
t *= l;
nr = (2 * q + r) * l;
var nn = (q * (7 * k) + 2 + r * l) / t;
q *= k;
l += 2;
++k;
n = nn;
}
r = nr;
}
}
}
}

Output:

Digits of PI
=========================================================================
Decimal  : 141592653589793238462643383279502884197169399375105820974944
Binary  : 001001000011111101101010100010001000010110100011000010001101
Quaternary : 021003331222202020112203002031030103012120220232000313001303
Octal  : 110375524210264302151423063050560067016321122011160210514763
Hexadecimal: 243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e
Alphabetic : DRSQLOLYRTRODNLHNQTGKUDQGTUIRXNEQBCKBSZIVQQVGDMELMUEXROIQIYA
Fun  : + -++ +---- + -++ -+++++ --+----- +++- +-+-+-+- +-++ +
Nibbles  : 2 4 3 f 6 a 8 8 8 5 a 3 0 8 d 3 1 3 1 9 8 a 2 e 0 3 7 0 7 3
Bytes  : 24 3f 6a 88 85 a3 08 d3 13 19 8a 2e 03 70 73 44 a4 09 38 22
Words  : 243f 6a88 85a3 08d3 1319 8a2e 0370 7344 a409 3822 299f 31d0
Dwords  : 243f6a88 85a308d3 13198a2e 03707344 a4093822 299f31d0
=========================================================================
* press any key to exit *

[edit] D

This modified Spigot algorithm does not continue infinitely, because its required memory grow as the number of digits need to print.

import std.stdio, std.conv, std.string;
 
struct PiDigits {
immutable uint nDigits;
 
int opApply(int delegate(ref string /*chunk of pi digit*/) dg){
// Maximum width for correct output, for type ulong.
enum size_t width = 9;
 
enum ulong scale = 10UL ^^ width;
enum ulong initDigit = 2UL * 10UL ^^ (width - 1);
enum string formatString = "%0" ~ text(width) ~ "d";
 
immutable size_t len = 10 * nDigits / 3;
auto arr = new ulong[len];
arr[] = initDigit;
ulong carry;
 
foreach (i; 0 .. nDigits / width) {
ulong sum;
foreach_reverse (j; 0 .. len) {
auto quo = sum * (j + 1) + scale * arr[j];
arr[j] = quo % (j*2 + 1);
sum = quo / (j*2 + 1);
}
auto yield = format(formatString, carry + sum/scale);
if (dg(yield))
break;
carry = sum % scale;
}
return 0;
}
}
 
void main() {
foreach (d; PiDigits(100))
writeln(d);
}

Output:

314159265
358979323
846264338
327950288
419716939
937510582
097494459
230781640
628620899
862803482
534211706

[edit] Alternative version

import std.stdio, std.bigint;
 
void main() {
int ndigits = 0;
auto q = BigInt(1);
auto r = BigInt(0);
auto t = q;
auto k = q;
auto n = BigInt(3);
auto l = n;
 
bool first = true;
while (ndigits < 1_000) {
if (4 * q + r - t < n * t) {
write(n); ndigits++;
if (ndigits % 70 == 0) writeln();
if (first) { first = false; write('.'); }
auto nr = 10 * (r - n * t);
n = ((10 * (3 * q + r)) / t) - 10 * n;
q *= 10;
r = nr;
} else {
auto nr = (2 * q + r) * l;
auto nn = (q * (7 * k + 2) + r * l) / (t * l);
q *= k;
t *= l;
l += 2;
k++;
n = nn;
r = nr;
}
}
}

Output:

3.141592653589793238462643383279502884197169399375105820974944592307816
4062862089986280348253421170679821480865132823066470938446095505822317
2535940812848111745028410270193852110555964462294895493038196442881097
5665933446128475648233786783165271201909145648566923460348610454326648
2133936072602491412737245870066063155881748815209209628292540917153643
6789259036001133053054882046652138414695194151160943305727036575959195
3092186117381932611793105118548074462379962749567351885752724891227938
1830119491298336733624406566430860213949463952247371907021798609437027
7053921717629317675238467481846766940513200056812714526356082778577134
2757789609173637178721468440901224953430146549585371050792279689258923
5420199561121290219608640344181598136297747713099605187072113499999983
7297804995105973173281609631859502445945534690830264252230825334468503
5261931188171010003137838752886587533208381420617177669147303598253490
4287554687311595628638823537875937519577818577805321712268066130019278
76611195909216420198

[edit] Erlang

% Implemented by Arjun Sunel
-module(pi_calculation).
-export([main/0]).
 
main() ->
pi(1,0,1,1,3,3,0).
 
pi(Q,R,T,K,N,L,C) ->
 
if C=:=50 ->
io:format("\n"),
pi(Q,R,T,K,N,L,0) ;
 
true ->
 
if
(4*Q + R-T) < (N*T) ->
io:format("~p",[N]),
P = 10*(R-N*T),
pi(Q*10 , P, T , K , ((10*(3*Q+R)) div T)-10*N , L,C+1);
 
true ->
P = (2*Q+R)*L,
M = (Q*(7*K)+2+(R*L)) div (T*L),
H = L+2,
J =K+ 1,
pi(Q*K, P , T*L ,J,M,H,C)
end
end.
 
Output:
31415926535897932384626433832795028841971693993751
05820974944592307816406286208998628034825342117067
98214808651328230664709384460955058223172535940812
84811174502841027019385211055596446229489549303819
64428810975665933446128475648233786783165271201909
14564856692346034861045432664821339360726024914127
37245870066063155881748815209209628292540917153643
67892590360011330530548820466521384146951941511609
43305727036575959195309218611738193261179310511854
80744623799627495673518857527248912279381830119491
29833673362440656643086021394946395224737190702179
86094370277053921717629317675238467481846766940513
20005681271452635608277857713427577896091736371787
21468440901224953430146549585371050792279689258923
54201995611212902196086403441815981362977477130996
05187072113499999983729780499510597317328160963185
95024459455346908302642522308253344685035261931188
17101000313783875288658753320838142061717766914730
35982534904287554687311595628638823537875937519577
81857780532171226806613001927876611195909216420198
93809525720106548586327886593615338182796823030195
20353018529689957736225994138912497217752834791315
15574857242454150695950829533116861727855889075098
38175463746493931925506040092770167113900984882401
28583616035637076601047101819429555961989467678374
4944825537977472684710404753464620


[edit] F#

Translation of: Haskell
let rec g q r t k n l = seq {
if 4I*q+r-t < n*t
then
yield n
yield! (g (10I*q) (10I*(r-n*t)) t k ((10I*(3I*q+r))/t - 10I*n) l)
else
yield! (g (q*k) ((2I*q+r)*l) (t*l) (k+1I) ((q*(7I*k+2I)+r*l)/(t*l)) (l+2I))
}
 
let π = (g 1I 0I 1I 1I 3I 3I)
 
Seq.take 1 π |> Seq.iter (printf "%A.")
// 6 digits beginning at position 762 of π are '9'
Seq.take 767 (Seq.skip 1 π) |> Seq.iter (printf "%A")
Output:
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066
470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831
652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903
600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527
248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051
320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219
6086403441815981362977477130996051870721134999999

[edit] FunL

The code for compute_pi() is from [2]. The number of digits may be given on the command line as an argument. If there's no argument, the program will run until interrupted.

def compute_pi =
def g( q, r, t, k, n, l ) =
if 4*q + r - t < n*t
n # g( 10*q, 10*(r - n*t), t, k, (10*(3*q + r))\t - 10*n, l )
else
g( q*k, (2*q + r)*l, t*l, k + 1, (q*(7*k + 2) + r*l)\(t*l), l + 2 )
 
g( 1, 0, 1, 1, 3, 3 )
 
if _name_ == '-main-'
print( compute_pi().head() + '.' )
 
if args.isEmpty()
for d <- compute_pi().tail()
print( d )
else
for d <- compute_pi().tail().take( int(args(0)) )
print( d )
 
println()

[edit] Go

Code below is a simplistic translation of Haskell code in Unbounded Spigot Algorithms for the Digits of Pi. This is the algorithm specified for the pidigits benchmark of the Computer Language Benchmarks Game. (The standard Go distribution includes source submitted to the benchmark site, and that code runs stunning faster than the code below.)

package main
 
import (
"fmt"
"math/big"
)
 
type lft struct {
q,r,s,t big.Int
}
 
func (t *lft) extr(x *big.Int) *big.Rat {
var n, d big.Int
var r big.Rat
return r.SetFrac(
n.Add(n.Mul(&t.q, x), &t.r),
d.Add(d.Mul(&t.s, x), &t.t))
}
 
var three = big.NewInt(3)
var four = big.NewInt(4)
 
func (t *lft) next() *big.Int {
r := t.extr(three)
var f big.Int
return f.Div(r.Num(), r.Denom())
}
 
func (t *lft) safe(n *big.Int) bool {
r := t.extr(four)
var f big.Int
if n.Cmp(f.Div(r.Num(), r.Denom())) == 0 {
return true
}
return false
}
 
func (t *lft) comp(u *lft) *lft {
var r lft
var a, b big.Int
r.q.Add(a.Mul(&t.q, &u.q), b.Mul(&t.r, &u.s))
r.r.Add(a.Mul(&t.q, &u.r), b.Mul(&t.r, &u.t))
r.s.Add(a.Mul(&t.s, &u.q), b.Mul(&t.t, &u.s))
r.t.Add(a.Mul(&t.s, &u.r), b.Mul(&t.t, &u.t))
return &r
}
 
func (t *lft) prod(n *big.Int) *lft {
var r lft
r.q.SetInt64(10)
r.r.Mul(r.r.SetInt64(-10), n)
r.t.SetInt64(1)
return r.comp(t)
}
 
func main() {
// init z to unit
z := new(lft)
z.q.SetInt64(1)
z.t.SetInt64(1)
 
// lfts generator
var k int64
lfts := func() *lft {
k++
r := new(lft)
r.q.SetInt64(k)
r.r.SetInt64(4*k+2)
r.t.SetInt64(2*k+1)
return r
}
 
// stream
for {
y := z.next()
if z.safe(y) {
fmt.Print(y)
z = z.prod(y)
} else {
z = z.comp(lfts())
}
}
}

[edit] Groovy

Translation of: Java

Solution:

BigInteger q = 1, r = 0, t = 1, k = 1, n = 3, l = 3
String nn
boolean first = true
 
while (true) {
(nn, first, q, r, t, k, n, l) = (4*q + r - t < n*t) \
? ["${n}${first?'.':''}", false, 10*q, 10*(r - n*t), t , k , 10*(3*q + r)/t - 10*n , l ] \
 : ['' , first, q*k , (2*q + r)*l , t*l, k + 1, (q*(7*k + 2) + r*l)/(t*l), l + 2]
print nn
}

Output (thru first 1000 iterations):

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337

[edit] Haskell

The code from [3]:

pi_ = g(1,0,1,1,3,3) where
g (q,r,t,k,n,l) =
if 4*q+r-t < n*t
then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

[edit] Complete command-line program

Works with: GHC version 7.4.1
#!/usr/bin/runhaskell
 
import Control.Monad
import System.IO
 
pi_ = g(1,0,1,1,3,3) where
g (q,r,t,k,n,l) =
if 4*q+r-t < n*t
then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)
 
digs = insertPoint digs'
where insertPoint (x:xs) = x:'
.':xs
digs'
= map (head . show) pi_
 
main = do
hSetBuffering stdout $ BlockBuffering $ Just 80
forM_ digs putChar
Output:
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198

[edit] Icon and Unicon

Translation of: PicoLisp
based on Jeremy Gibbons' Haskell solution.
procedure pi (q, r, t, k, n, l)
first := "yes"
repeat { # infinite loop
if (4*q+r-t < n*t) then {
suspend n
if (\first) := &null then suspend "."
# compute and update variables for next cycle
nr := 10*(r-n*t)
n := ((10*(3*q+r)) / t) - 10*n
q *:= 10
r := nr
} else {
# compute and update variables for next cycle
nr := (2*q+r)*l
nn := (q*(7*k+2)+r*l) / (t*l)
q *:= k
t *:= l
l +:= 2
k +:= 1
n := nn
r := nr
}
}
end
 
procedure main ()
every (writes (pi (1,0,1,1,3,3)))
end

[edit] J

pi=:3 :0
smoutput"0'3.1'
n=.0 while.n=.n+1 do.
smoutput-/1 10*<.@o.10x^1 0+n
end.
)

Example use:

   pi''
3
.
1
4
1
5
9
2
6
5
3
...

[edit] Java

Translation of: Icon
import java.math.BigInteger ;
 
public class Pi {
final BigInteger TWO = BigInteger.valueOf(2) ;
final BigInteger THREE = BigInteger.valueOf(3) ;
final BigInteger FOUR = BigInteger.valueOf(4) ;
final BigInteger SEVEN = BigInteger.valueOf(7) ;
 
BigInteger q = BigInteger.ONE ;
BigInteger r = BigInteger.ZERO ;
BigInteger t = BigInteger.ONE ;
BigInteger k = BigInteger.ONE ;
BigInteger n = BigInteger.valueOf(3) ;
BigInteger l = BigInteger.valueOf(3) ;
 
public void calcPiDigits(){
BigInteger nn, nr ;
boolean first = true ;
while(true){
if(FOUR.multiply(q).add(r).subtract(t).compareTo(n.multiply(t)) == -1){
System.out.print(n) ;
if(first){System.out.print(".") ; first = false ;}
nr = BigInteger.TEN.multiply(r.subtract(n.multiply(t))) ;
n = BigInteger.TEN.multiply(THREE.multiply(q).add(r)).divide(t).subtract(BigInteger.TEN.multiply(n)) ;
q = q.multiply(BigInteger.TEN) ;
r = nr ;
System.out.flush() ;
}else{
nr = TWO.multiply(q).add(r).multiply(l) ;
nn = q.multiply((SEVEN.multiply(k))).add(TWO).add(r.multiply(l)).divide(t.multiply(l)) ;
q = q.multiply(k) ;
t = t.multiply(l) ;
l = l.add(TWO) ;
k = k.add(BigInteger.ONE) ;
n = nn ;
r = nr ;
}
}
}
 
public static void main(String[] args) {
Pi p = new Pi() ;
p.calcPiDigits() ;
}
}

Output :

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480 ...

[edit] jq

Works with: jq version 1.4

The focus in this section is on the Gibbons spigot algorithm as it is relatively simple and therefore provides a gentle introduction to how such algorithms can be implemented in jq.

Since the Gibbons algorithm quickly fails in the absence of support for large integers, we shall assume BigInt support, such as provided by BigInt.jq.

The jq program presented here closely follows the Groovy and Python examples on this page. The spigot generator is named "next", and is driven by an annotation function, "decorate"; thus the main program is just "S0 | decorate(next)" where S0 is the initial state. One advantage of this approach is that the generator's state is exposed, thus making it easy to restart the stream at any point.

The annotation defined here results in a triple for each digit of pi: [index, digit, space], where "space" is the the sum of the lengths of the strings in the six-dimensional state vector, [q, r, t, k, n, l]. The output shows that the space requirements of the Gibbons spigot grow very slightly more than linearly.

# The Gibbons spigot, in the mold of the [[#Groovy]] and ython]] programs shown on this page.  
# The "bigint" functions
needed are: long_minus long_add long_multiply long_div
 
def pi_spigot:
 
# S is the sixtuple:
# q r t k n l
# 0 1 2 3 4 5
 
def long_lt(x;y): if x == y then false else lessOrEqual(x;y) end;
 
def check:
long_lt(long_minus(long_add(long_multiply("4"; .[0]); .[1]) ; .[2]);
long_multiply(.[4]; .[2]));
 
# state: [d, S] where digit is null or a digit ready to be printed
def next:
.[1] as $S
| $S[0] as $q | $S[1] as $r | $S[2] as $t | $S[3] as $k | $S[4] as $n | $S[5] as $l
| if $S|check
then [$n,
[long_multiply("10"; $q),
long_multiply("10"; long_minus($r; long_multiply($n;$t))),
$t,
$k,
long_minus( long_div(long_multiply("10";long_add(long_multiply("3"; $q); $r)); $t );
long_multiply("10";$n)),
$l ]]
else [null,
[long_multiply($q;$k),
long_multiply( long_add(long_multiply("2";$q); $r); $l),
long_multiply($t;$l),
long_add($k; "1"),
long_div( long_add(long_multiply($q; long_add(long_multiply("7";$k); "2")) ; long_multiply($r;$l));
long_multiply($t;$l) ),
long_add($l; "2") ]]
end;
 
# Input: input to the filter "nextstate"
# Output: [count, space, digit] for successive digits produced by "nextstate"
def decorate( nextstate ):
 
# For efficiency it is important that the recursive
# function have arity 0 and be tail-recursive:
def count:
.[0] as $count
| .[1] as $state
| $state[0] as $value
| ($state[1] | map(length) | add) as $space
| (if $value then [$count, $space, $value] else empty end),
( [if $value then $count+1 else $count end, ($state | nextstate)] | count);
[0, .] | count;
 
# q=1, r=0, t=1, k=1, n=3, l=3
[null, ["1", "0", "1", "1", "3", "3"]] | decorate(next)
;
 
pi_spigot
Output:
$ jq -M -n -c -f pi.bigint.jq
[0,9,"3"]
[1,14,"1"]
[2,29,"4"]
[3,36,"1"]
[4,51,"5"]
[5,69,"9"]
[6,80,"2"]
[7,95,"6"]
[8,115,"5"]
[9,125,"3"]
[10,142,"5"]
[11,167,"8"]
[12,181,"9"]
[13,197,"7"]
[14,226,"9"]
[15,245,"3"]
[16,263,"2"]
[17,276,"3"]
[18,300,"8"]
[19,320,"4"]
[20,350,"6"]
[21,363,"2"]
[22,383,"6"]
[23,408,"4"]
[24,429,"3"]
[25,442,"3"]
[26,475,"8"]
[27,502,"3"]
[28,510,"2"]
[29,531,"7"]
[30,563,"9"]
[31,611,"5"]
[32,613,"0"]
[33,628,"2"]
[34,649,"8"]
[35,676,"8"]
[36,711,"4"]
[37,720,"1"]
[38,748,"9"]
[39,783,"7"]
[40,792,"1"]
[41,814,"6"]
[42,849,"9"]
[43,870,"3"]
[44,886,"9"]
[45,923,"9"]
[46,939,"3"]
[47,967,"7"]
[48,1004,"5"]
[49,1041,"1"]
[50,1043,"0"]
[51,1059,"5"]
[52,1103,"8"]
[53,1133,"2"]
[54,1135,"0"]
[55,1165,"9"]
[56,1195,"7"]
[57,1212,"4"]
[58,1242,"9"]
[59,1273,"4"]
[60,1297,"4"]
[61,1313,"5"]
[62,1358,"9"]
[63,1375,"2"]
[64,1421,"3"]
[65,1423,"0"]
[66,1447,"7"]
[67,1493,"8"]
[68,1501,"1"]
[69,1533,"6"]
[70,1579,"4"]
[71,1581,"0"]
[72,1613,"6"]
[73,1630,"2"]
[74,1662,"8"]
[75,1701,"6"]
[76,1733,"2"]
[77,1735,"0"]
[78,1781,"8"]
[79,1792,"9"]
[80,1816,"9"]
[81,1849,"8"]
[82,1889,"6"]
[83,1898,"2"]
[84,1961,"8"]
[85,1963,"0"]
[86,1988,"3"]
[87,2013,"4"]
[88,2054,"8"]
[89,2071,"2"]
[90,2104,"5"]
[91,2129,"3"]
[92,2162,"4"]
[93,2195,"2"]
[94,2220,"1"]
[95,2230,"1"]
[96,2287,"7"]
[97,2289,"0"]
[98,2314,"6"]
[99,2340,"7"]
[100,2373,"9"]
[101,2414,"8"]
[102,2448,"2"]
[103,2458,"1"]
[104,2484,"4"]
[105,2534,"8"]
[106,2536,"0"]
[107,2569,"8"]
[108,2602,"6"]
[109,2645,"5"]
[110,2662,"1"]
[111,2696,"3"]
[112,2707,"2"]
[113,2756,"8"]
[114,2775,"2"]
[115,2825,"3"]
[116,2827,"0"]
[117,2853,"6"]
[118,2887,"6"]
[119,2914,"4"]
[120,2964,"7"]
[121,2966,"0"]
[122,3008,"9"]
[123,3027,"3"]
[124,3061,"8"]
[125,3088,"4"]
[126,3114,"4"]
[127,3165,"6"]
[128,3167,"0"]
[129,3202,"9"]
[130,3237,"5"]
[131,3287,"5"]
[132,3289,"0"]
[133,3316,"5"]
[134,3360,"8"]
[135,3387,"2"]
[136,3414,"2"]
[137,3456,"3"]
[138,3466,"1"]
[139,3510,"7"]
[140,3529,"2"]
[141,3564,"5"]
[142,3583,"3"]
[143,3610,"5"]
[144,3653,"9"]
[145,3697,"4"]
[146,3699,"0"]
[147,3752,"8"]
[148,3770,"1"]
[149,3789,"2"]
[150,3825,"8"]
[151,3852,"4"]
[152,3905,"8"]
[153,3933,"1"]
[154,3960,"1"]
[155,3970,"1"]
[156,4006,"7"]
[157,4033,"4"]
[158,4102,"5"]
[159,4104,"0"]
[160,4124,"2"]
[161,4159,"8"]
[162,4203,"4"]
[163,4248,"1"]
[164,4250,"0"]
[165,4269,"2"]
[166,4348,"7"]
[167,4350,"0"]
[168,4361,"1"]
[169,4405,"9"]
[170,4424,"3"]
[171,4460,"8"]
[172,4497,"5"]
[173,4542,"2"]
[174,4569,"1"]
[175,4605,"1"]
[176,4607,"0"]
[177,4644,"5"]
[178,4672,"5"]
[179,4691,"5"]
[180,4727,"9"]
[181,4764,"6"]
[182,4792,"4"]
[183,4820,"4"]
[184,4865,"6"]
[185,4893,"2"]
[186,4913,"2"]
[187,4949,"9"]
[188,4968,"4"]
[189,5005,"8"]
[190,5042,"9"]
[191,5070,"5"]
[192,5098,"4"]
[193,5144,"9"]
[194,5198,"3"]
[195,5200,"0"]
[196,5219,"3"]
[197,5266,"8"]
[198,5276,"1"]
[199,5313,"9"]
[200,5350,"6"]
[201,5387,"4"]
[202,5416,"4"]
[203,5435,"2"]
[204,5471,"8"]
[205,5526,"8"]
[206,5556,"1"]
[207,5558,"0"]
[208,5594,"9"]
[209,5632,"7"]
[210,5660,"5"]
[211,5689,"6"]
[212,5726,"6"]
[213,5746,"5"]
[214,5792,"9"]
[215,5821,"3"]
[216,5849,"3"]
[217,5887,"4"]
[218,5906,"4"]
[219,5961,"6"]
[220,5981,"1"]
[221,6002,"2"]
[222,6038,"8"]
[223,6068,"4"]
[224,6096,"7"]
[225,6134,"5"]
[226,6163,"6"]
[227,6191,"4"]
[228,6238,"8"]
[229,6267,"2"]
[230,6296,"3"]
[231,6316,"3"]
[232,6344,"7"]
[233,6383,"8"]
[234,6411,"6"]
[235,6440,"7"]
[236,6487,"8"]
[237,6525,"3"]
[238,6545,"1"]
[239,6574,"6"]
[240,6621,"5"]
[241,6641,"2"]
[242,6688,"7"]
[243,6717,"1"]
[244,6782,"2"]
[245,6784,"0"]
[246,6795,"1"]
[247,6852,"9"]
[248,6854,"0"]
[249,6910,"9"]
[250,6929,"1"]
[251,6959,"4"]
[252,6988,"5"]
[253,7027,"6"]
[254,7046,"4"]
[255,7085,"8"]
[256,7115,"5"]
[257,7153,"6"]
[258,7181,"6"]
[259,7229,"9"]
[260,7258,"2"]
[261,7288,"3"]
[262,7317,"4"]
[263,7383,"6"]
[264,7385,"0"]
[265,7415,"3"]
[266,7435,"4"]
[267,7474,"8"]
[268,7530,"6"]
[269,7569,"1"]
[270,7571,"0"]
[271,7609,"4"]
[272,7639,"5"]
[273,7678,"4"]
[274,7716,"3"]
[275,7736,"2"]
[276,7766,"6"]
[277,7805,"6"]
[278,7826,"4"]
[279,7873,"8"]
[280,7912,"2"]
[281,7933,"1"]
[282,7971,"3"]
[283,7991,"3"]
[284,8030,"9"]
[285,8060,"3"]
[286,8118,"6"]
[287,8120,"0"]
[288,8168,"7"]
[289,8189,"2"]
[290,8264,"6"]
[291,8266,"0"]
[292,8287,"2"]
[293,8317,"4"]
[294,8374,"9"]
[295,8395,"1"]
[296,8443,"4"]
[297,8464,"1"]
[298,8485,"2"]
[299,8524,"7"]
[300,8544,"3"]
[301,8593,"7"]
[302,8623,"2"]
...
 

[edit] Julia

Julia comes with built-in support for computing π in arbitrary precision (using the GNU MPFR library). This implementation computes π at precisions that are repeatedly doubled as more digits are needed, printing one digit at a time and never terminating (until it runs out of memory) as specified:

prec = get_bigfloat_precision()
spi = ""
digit = 1
while true
if digit > length(spi) - 6
prec *= 2
set_bigfloat_precision(prec)
spi = string(big(π))
end
print(spi[digit])
digit += 1
end

Output:

3.141592653589793238462643383279502884195e69399375105820974944592307816406286198e9862803482534211706798214808651328230664709384460955058223172535940812848115e450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997e0631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526357e8277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201...

[edit] Lasso

Based off Dik T. Winter's C implementation of Beeler et al. 1972, Item 120.

#!/usr/bin/lasso9
 
define generatePi => {
yield currentCapture
 
local(r = array(), i, k, b, d, c = 0, x)
with i in generateSeries(1, 2800)
do #r->insert(2000)
with k in generateSeries(2800, 1, -14)
do {
#d = 0
#i = #k
while(true) => {
#d += #r->get(#i) * 10000
#b = 2 * #i - 1
#r->get(#i) = #d % #b
#d /= #b
#i--
 !#i ? loop_abort
#d *= #i
}
#x = (#c + #d / 10000)
yield (#k == 2800 ? ((#x * 0.001)->asstring(-precision = 3)) | #x->asstring(-padding=4, -padChar='0'))
#c = #d % 10000
}
}
 
local(pi_digits) = generatePi
loop(200) => {
stdout(#pi_digits())
}

Output (first 100 places):

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067

[edit] Liberty BASIC

Pretty slow if you run for over 100 digits...

    ndigits = 0
 
q = 1
r = 0
t = q
k = q
n = 3
L = n
 
first = 666 ' ANY non-zero =='true' in LB.
 
while ndigits <100
if ( 4 *q +r -t) <( n *t) then
print n;
ndigits =ndigits +1
if not( ndigits mod 40) then print: print " ";
if first =666 then first = 0: print ".";
nr =10 *( r -n *t)
n =int( ( (10 *( 3 *q +r)) /t) -10 *n)
q =q *10
r =nr
else
nr =( 2 *q +r) *L
nn =(q *( 7 *k +2) +r *L) /( t *L)
q =q *k
t =t *L
L =L +2
k =k +1
n =int( nn)
r =nr
end if
scan
wend
 
end
3.141592653589793238462643383279502884197
1693993751058209749445923078164062862089
98628034825342117067

[edit] Lua

Translation of: Pascal
a = {}
n = 1000
len = math.modf( 10 * n / 3 )
 
for j = 1, len do
a[j] = 2
end
nines = 0
predigit = 0
for j = 1, n do
q = 0
for i = len, 1, -1 do
x = 10 * a[i] + q * i
a[i] = math.fmod( x, 2 * i - 1 )
q = math.modf( x / ( 2 * i - 1 ) )
end
a[1] = math.fmod( q, 10 )
q = math.modf( q / 10 )
if q == 9 then
nines = nines + 1
else
if q == 10 then
io.write( predigit + 1 )
for k = 1, nines do
io.write(0)
end
predigit = 0
nines = 0
else
io.write( predigit )
predigit = q
if nines ~= 0 then
for k = 1, nines do
io.write( 9 )
end
nines = 0
end
end
end
end
print( predigit )
03141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086 ...

[edit] Mathematica / Wolfram Language

User can interrupt computation using "Alt+." or "Cmd+." on a Mac.

N[Pi, 1000000!]

[edit] MATLAB / Octave

Matlab / Octave use double precesion numbers per default, and pi is a builtin constant value. Arbitrary precision is only implemented in some additional toolboxes (e.g. symbolic toolbox).

pi
>> pi
ans =  3.1416
> printf('%.60f\n',pi)
3.141592653589793115997963468544185161590576171875000000000000>> format long
Unfortunately this is not the correct value!
3.14159265358979323846264338327950288419716939937510582
=================??????????????????????????????????????

[edit] NetRexx

Translation of: Java
/* NetRexx */
options replace format comments java crossref symbols binary
import java.math.BigInteger
 
runSample(arg)
return
 
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
method runSample(arg) private static
parse arg places .
if places = '' then places = -1
 
TWO = BigInteger.valueOf(2)
THREE = BigInteger.valueOf(3)
FOUR = BigInteger.valueOf(4)
SEVEN = BigInteger.valueOf(7)
 
q_ = BigInteger.ONE
r_ = BigInteger.ZERO
t_ = BigInteger.ONE
k_ = BigInteger.ONE
n_ = BigInteger.valueOf(3)
l_ = BigInteger.valueOf(3)
 
nn = BigInteger
nr = BigInteger
 
first = isTrue()
digitCt = 0
loop forever
if FOUR.multiply(q_).add(r_).subtract(t_).compareTo(n_.multiply(t_)) == -1 then do
digitCt = digitCt + 1
if places > 0 & digitCt - 1 > places then leave
say n_'\-'
if first then do
say '.\-'
first = isFalse()
end
nr = BigInteger.TEN.multiply(r_.subtract(n_.multiply(t_)))
n_ = BigInteger.TEN.multiply(THREE.multiply(q_).add(r_)).divide(t_).subtract(BigInteger.TEN.multiply(n_))
q_ = q_.multiply(BigInteger.TEN)
r_ = nr
end
else do
nr = TWO.multiply(q_).add(r_).multiply(l_)
nn = q_.multiply((SEVEN.multiply(k_))).add(TWO).add(r_.multiply(l_)).divide(t_.multiply(l_))
q_ = q_.multiply(k_)
t_ = t_.multiply(l_)
l_ = l_.add(TWO)
k_ = k_.add(BigInteger.ONE)
n_ = nn
r_ = nr
end
end
say
 
return
 
method isTrue() private static returns boolean
return (1 == 1)
method isFalse() private static returns boolean
return \isTrue()
 
Output:
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679...

[edit] Nimrod

Library: bigints
import strutils, unsigned, bigints
 
var
tmp1, tmp2, tmp3, acc, k, dd = initBigInt(0)
den, num, k2 = initBigInt(1)
 
proc extractDigit(): int32 =
if num > acc:
return -1
 
tmp3 = num shl 1
tmp3 += num
tmp3 += acc
tmp2 = tmp3 mod den
tmp1 = tmp3 div den
tmp2 += num
 
if tmp2 >= den:
return -1
 
result = int32(tmp1.limbs[0])
 
proc eliminateDigit(d: int32) =
acc -= den * d
acc *= 10
num *= 10
 
proc nextTerm() =
k += 1
k2 += 2
tmp1 = num shl 1
acc += tmp1
acc *= k2
den *= k2
num *= k
 
var i = 0
 
while true:
var d: int32 = -1
while d < 0:
nextTerm()
d = extractDigit()
 
stdout.write chr(ord('0') + d)
inc i
if i == 40:
echo ""
i = 0
eliminateDigit d

Output:

3141592653589793238462643383279502884197
1693993751058209749445923078164062862089
9862803482534211706798214808651328230664
7093844609550582231725359408128481117450
...

[edit] OCaml

The Constructive Real library Creal contains an infinite-precision Pi, so we can just print out its digits.

open Creal;;
 
let block = 100 in
let segment n =
let s = to_string pi (n*block) in
String.sub s ((n-1)*block) block in
let counter = ref 1 in
while true do
print_string (segment !counter);
flush stdout;
incr counter
done

However that is cheating if you want to see an algorithm to generate Pi. Since the Spigot algorithm is already used in the pidigits program, this implements Machin's formula.

open Num
 
(* series for: c*atan(1/k) *)
class atan_sum c k = object
val kk = k*/k
val mutable n = 0
val mutable kpow = k
val mutable pterm = c*/k
val mutable psum = Int 0
val mutable sum = c*/k
method next =
n <- n+1; kpow <- kpow*/kk;
let t = c*/kpow//(Int (2*n+1)) in
pterm <- if n mod 2 = 0 then t else minus_num t;
psum <- sum;
sum <- sum +/ pterm
method error = abs_num pterm
method bounds = if pterm </ Int 0 then (sum, psum) else (psum, sum)
end;;
 
let inv i = (Int 1)//(Int i) in
let t1 = new atan_sum (Int 16) (inv 5) in
let t2 = new atan_sum (Int (-4)) (inv 239) in
let base = Int 10 in
let npr = ref 0 in
let shift = ref (Int 1) in
let d_acc = inv 10000 in
let acc = ref d_acc in
let shown = ref (Int 0) in
while true do
while t1#error >/ !acc do t1#next done;
while t2#error >/ !acc do t2#next done;
let (lo1, hi1), (lo2, hi2) = t1#bounds, t2#bounds in
let digit x = int_of_num (floor_num ((x -/ !shown) */ !shift)) in
let d, d' = digit (lo1+/lo2), digit (hi1+/hi2) in
if d = d' then (
print_int d;
if !npr = 0 then print_char '.';
flush stdout;
shown := !shown +/ ((Int d) // !shift);
incr npr; shift := !shift */ base;
) else (acc := !acc */ d_acc);
done

[edit] PARI/GP

Uses the built-in Brent-Salamin arithmetic-geometric mean iteration.

pi()={
my(x=Pi,n=0,t);
print1("3.");
while(1,
if(n>=default(realprecision),
default(realprecision,default(realprecision)*2);
x=Pi
);
print1(floor(x*10^n++)%10)
)
};

[edit] Pascal

Works with: Free_Pascal

With minor editing changes as published by Stanley Rabinowitz in [4]. Minor improvement of <user>Mischi</user> { speedup ~2 ( n=100000 , rumtime 4s-> 1,44s fpc 2.6.4 -O3 }, by calculating only necessary digits up to n.

Program Pi_Spigot;
const
n = 1000;
len = 10*n div 3;
 
var
j, k, q, nines, predigit: integer;
a: array[0..len] of longint;
 
function OneLoop(i:integer):integer;
var
x: integer;
begin
{Only calculate as far as needed }
{+16 for security digits ~5 decimals}
i := i*10 div 3+16;
IF i > len then
i := len;
result := 0;
repeat {Work backwards}
x := 10*a[i] + result*i;
result := x div (2*i - 1);
a[i] := x - result*(2*i - 1);//x mod (2*i - 1)
dec(i);
until i<= 0 ;
end;
 
begin
 
for j := 1 to len do
a[j] := 2; {Start with 2s}
nines := 0;
predigit := 0; {First predigit is a 0}
 
for j := 1 to n do
begin
q := OneLoop(n-j);
a[1] := q mod 10;
q := q div 10;
if q = 9 then
nines := nines + 1
else
if q = 10 then
begin
write(predigit+1);
for k := 1 to nines do
write(0); {zeros}
predigit := 0;
nines := 0
end
else
begin
write(predigit);
predigit := q;
if nines <> 0 then
begin
for k := 1 to nines do
write(9);
nines := 0
end
end
end;
writeln(predigit);
end.

Output:

% ./Pi_Spigot
0314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692345909920711885680651674519425285434096427182587759486723722781805263647624000125396676958803786501460828202281029373153209312571960818250880554547473856719984185863746093965326651520597323020047032736774801454941806997686761324327197543884752285691446935454123565860783675210974816763774404244435759134410295498726734853215534522607487460726127229475668144011717114616426232946428024519416553103668226495999762474794964971942580962672394353639091253676024666416611752768663496859511953379074492301907334533332168653316848429072409790153107727405217592705578811026822246918375946146833893295768948223454910618813165693638675020776147309271039650546608291414446931319129743860775232869818000001373465107828774083538328579350606183962012466

[edit] Perl

Perl being what it is, there are many ways to do this with many variations. With a fixed number of digits and the Math::BigInt::GMP library installed, the [[Arithmetic-geometric mean/Calculate Pi code will be much faster than any of these methods other than some of the modules. If Math::GMP is installed, then replacing "use bigint" with "use Math::GMP qw/:constant/" in either the Perl6 spigot or Machin methods below will be pretty fast. They are not too bad if the Math::BigInt::GMP library is installed. With the default Math::BigInt backend, the AGM code isn't very fast and the Perl6 spigot and Machin methods are very slow.

[edit] Simple Spigot

This takes a numer-of-digits argument, but we can make it large (albeit using memory and some startup time). Unlike the other two, this uses no modules and does not require bigints so is worth showing.

sub pistream {
my $digits = shift;
my(@out, @a);
my($b, $c, $d, $e, $f, $g, $i, $d4, $d3, $d2, $d1);
my $outi = 0;
 
$digits++;
$b = $d = $e = $g = $i = 0;
$f = 10000;
$c = 14 * (int($digits/4)+2);
@a = (20000000) x $c;
print "3.";
while (($b = $c -= 14) > 0 && $i < $digits) {
$d = $e = $d % $f;
while (--$b > 0) {
$d = $d * $b + $a[$b];
$g = ($b << 1) - 1;
$a[$b] = ($d % $g) * $f;
$d = int($d / $g);
}
$d4 = $e + int($d/$f);
if ($d4 > 9999) {
$d4 -= 10000;
$out[$i-1]++;
for ($b = $i-1; $out[$b] == 1; $b--) {
$out[$b] = 0;
$out[$b-1]++;
}
}
$d3 = int($d4/10);
$d2 = int($d3/10);
$d1 = int($d2/10);
$out[$i++] = $d1;
$out[$i++] = $d2-$d1*10;
$out[$i++] = $d3-$d2*10;
$out[$i++] = $d4-$d3*10;
print join "", @out[$i-15 .. $i-15+3] if $i >= 16;
}
# We've closed the spigot. Print the remainder without rounding.
print join "", @out[$i-15+4 .. $digits-2], "\n";
}

[edit] Perl6 spigot

As mentioned earlier, replacing "use bigint" with "use Math::GMP qw/:constant/" will result in many orders of magnitude faster performance.

Translation of: Perl 6
use bigint try=>"GMP";
sub stream {
my ($next, $safe, $prod, $cons, $z, $x) = @_;
$x = $x->();
sub {
while (1) {
my $y = $next->($z);
if ($safe->($z, $y)) {
$z = $prod->($z, $y);
return $y;
} else {
$z = $cons->($z, $x->());
}
}
}
}
 
sub extr {
use integer;
my ($q, $r, $s, $t) = @{shift()};
my $x = shift;
($q * $x + $r) / ($s * $x + $t);
}
 
sub comp {
my ($q, $r, $s, $t) = @{shift()};
my ($u, $v, $w, $x) = @{shift()};
[$q * $u + $r * $w,
$q * $v + $r * $x,
$s * $u + $t * $w,
$s * $v + $t * $x];
}
 
my $pi_stream = stream
sub { extr shift, 3 },
sub { my ($z, $n) = @_; $n == extr $z, 4 },
sub { my ($z, $n) = @_; comp([10, -10*$n, 0, 1], $z) },
\&comp,
[1, 0, 0, 1],
sub { my $n = 0; sub { $n++; [$n, 4 * $n + 2, 0, 2 * $n + 1] } },
;
$|++;
print $pi_stream->(), '.';
print $pi_stream->() while 1;

[edit] Machin's Formula

Here is an original Perl 5 code, using Machin's formula. Not the fastest program in the world. As with the previous code, using either Math::GMP or Math::BigInt::GMP instead of the default bigint Calc backend will make it run thousands of times faster.

use bigint try=>"GMP"
 
# Pi/4 = 4 arctan 1/5 - arctan 1/239
# expanding it with Taylor series with what's probably the dumbest method
 
my ($ds, $ns) = (1, 0);
my ($n5, $d5) = (16 * (25 * 3 - 1), 3 * 5**3);
my ($n2, $d2) = (4 * (239 * 239 * 3 - 1), 3 * 239**3);
 
sub next_term {
my ($coef, $p) = @_[1, 2];
$_[0] /= ($p - 4) * ($p - 2);
$_[0] *= $p * ($p + 2) * $coef**4;
}
 
my $p2 = 5;
my $pow = 1;
 
$| = 1;
for (my $x = 5; ; $x += 4) {
($ns, $ds) = ($ns * $d5 + $n5 * $pow * $ds, $ds * $d5);
 
next_term($d5, 5, $x);
$n5 = 16 * (5 * 5 * ($x + 2) - $x);
 
while ($d5 > $d2) {
($ns, $ds) = ($ns * $d2 - $n2 * $pow * $ds, $ds * $d2);
$n2 = 4 * (239 * 239 * ($p2 + 2) - $p2);
next_term($d2, 239, $p2);
$p2 += 4;
}
 
my $ppow = 1;
while ($pow * $n5 * 5**4 < $d5 && $pow * $n2 * $n2 * 239**4 < $d2) {
$pow *= 10;
$ppow *= 10;
}
 
if ($ppow > 1) {
$ns *= $ppow;
#FIX? my $out = $ns->bdiv($ds); # bugged?
my $out = $ns / $ds;
$ns %= $ds;
 
$out = ("0" x (length($ppow) - length($out) - 1)) . $out;
print $out;
}
 
if ( $p2 % 20 == 1) {
my $g = Math::BigInt::bgcd($ds, $ns);
$ds /= $g;
$ns /= $g;
}
}

[edit] Modules

While no current CPAN module does continuous printing, there are (usually fast) ways to get digits of Pi. Examples include:

 
use ntheory qw/Pi/;
say Pi(10000);
 
use Math::Pari qw/setprecision Pi/;
setprecision(10000);
say Pi;
 
use Math::MPFR;
my $pi = Math::MPFR->new();
Math::MPFR::Rmpfr_set_prec($pi, int(10000 * 3.322)+40);
Math::MPFR::Rmpfr_const_pi($pi, 0);
say Math::MPFR::Rmpfr_get_str($pi, 10, 10000, 0);
 
use Math::BigFloat try=>"GMP"; # Slow without Math::BigInt::GMP installed
say Math::BigFloat::bpi(10000); # For over ~2k digits, slower than AGM
 
use Math::Big qw/pi/; # Very slow
say pi(10000);
 

[edit] Perl 6

# based on http://www.mathpropress.com/stan/bibliography/spigot.pdf
 
sub stream(&next, &safe, &prod, &cons, $z is copy, @x) {
gather loop {
$z = safe($z, my $y = next($z)) ??
prod($z, take $y) !!
cons($z, @x[(state $)++])
}
}
 
sub extr([$q, $r, $s, $t], $x) {
($q * $x + $r) div ($s * $x + $t)
}
 
sub comp([$q,$r,$s,$t], [$u,$v,$w,$x]) {
[$q * $u + $r * $w,
$q * $v + $r * $x,
$s * $u + $t * $w,
$s * $v + $t * $x]
}
 
my @pi :=
stream -> $z { extr($z, 3) },
-> $z, $n { $n == extr($z, 4) },
-> $z, $n { comp([10, -10*$n, 0, 1], $z) },
&comp,
<1 0 0 1>,
(1..*).map: { [$_, 4 * $_ + 2, 0, 2 * $_ + 1] }
 
loop {
print @pi.shift;
once print '.'
}

[edit] PicoLisp

The following script uses the spigot algorithm published by Jeremy Gibbons. Hit Ctrl-C to stop it.

#!/usr/bin/picolisp /usr/lib/picolisp/lib.l
 
(de piDigit ()
(job '((Q . 1) (R . 0) (S . 1) (K . 1) (N . 3) (L . 3))
(while (>= (- (+ R (* 4 Q)) S) (* N S))
(mapc set '(Q R S K N L)
(list
(* Q K)
(* L (+ R (* 2 Q)))
(* S L)
(inc K)
(/ (+ (* Q (+ 2 (* 7 K))) (* R L)) (* S L))
(+ 2 L) ) ) )
(prog1 N
(let M (- (/ (* 10 (+ R (* 3 Q))) S) (* 10 N))
(setq Q (* 10 Q) R (* 10 (- R (* N S))) N M) ) ) ) )
 
(prin (piDigit) ".")
(loop
(prin (piDigit))
(flush) )

Output:

3.14159265358979323846264338327950288419716939937510582097494459 ...

[edit] PL/I

/* Uses the algorithm of S. Rabinowicz and S. Wagon, "A Spigot Algorithm */
/* for the Digits of Pi". */
(subrg, fofl, size):
Pi_Spigot: procedure options (main); /* 21 January 2012. */
declare (n, len) fixed binary;
 
n = 1000;
len = 10*n / 3;
begin;
declare ( i, j, k, q, nines, predigit ) fixed binary;
declare x fixed binary (31);
declare a(len) fixed binary (31);
 
a = 2; /* Start with 2s */
nines, predigit = 0; /* First predigit is a 0 */
do j = 1 to n;
q = 0;
do i = len to 1 by -1; /* Work backwards */
x = 10*a(i) + q*i;
a(i) = mod (x, (2*i-1));
q = x / (2*i-1);
end;
a(1) = mod(q, 10); q = q / 10;
if q = 9 then nines = nines + 1;
else if q = 10 then
do;
put edit(predigit+1) (f(1));
do k = 1 to nines;
put edit ('0')(a(1)); /* zeros */
end;
predigit, nines = 0;
end;
else
do;
put edit(predigit) (f(1)); predigit = q;
do k = 1 to nines; put edit ('9')(a(1)); end;
nines = 0;
end;
end;
put edit(predigit) (f(1));
end; /* of begin block */
end Pi_Spigot;

output:

03141592653589793238462643383279502884197169399375105820974944592307816406286208
99862803482534211706798214808651328230664709384460955058223172535940812848111745
02841027019385211055596446229489549303819644288109756659334461284756482337867831
65271201909145648566923460348610454326648213393607260249141273724587006606315588
17488152092096282925409171536436789259036001133053054882046652138414695194151160
94330572703657595919530921861173819326117931051185480744623799627495673518857527
24891227938183011949129833673362440656643086021394946395224737190702179860943702
77053921717629317675238467481846766940513200056812714526356082778577134275778960
91736371787214684409012249534301465495853710507922796892589235420199561121290219
60864034418159813629774771309960518707211349999998372978049951059731732816096318
59502445945534690830264252230825334468503526193118817101000313783875288658753320
83814206171776691473035982534904287554687311595628638823537875937519577818577805
32171226806613001927876611195909216420198

[edit] PureBasic

Calculate Pi, limited to ~24 M-digits for memory and speed reasons.

#SCALE = 10000
#ARRINT= 2000
 
Procedure Pi(Digits)
Protected First=#True, Text$
Protected Carry, i, j, sum
Dim Arr(Digits)
For i=0 To Digits
Arr(i)=#ARRINT
Next
i=Digits
While i>0
sum=0
j=i
While j>0
sum*j+#SCALE*arr(j)
Arr(j)=sum%(j*2-1)
sum/(j*2-1)
j-1
Wend
Text$ = RSet(Str(Carry+sum/#SCALE),4,"0")
If First
Text$ = ReplaceString(Text$,"3","3.")
First = #False
EndIf
Print(Text$)
Carry=sum%#SCALE
i-14
Wend
EndProcedure
 
If OpenConsole()
SetConsoleCtrlHandler_(?Ctrl,#True)
Pi(24*1024*1024)
EndIf
End
 
Ctrl:
PrintN(#CRLF$+"Ctrl-C was pressed")
End

[edit] Python

def calcPi():
q, r, t, k, n, l = 1, 0, 1, 1, 3, 3
while True:
if 4*q+r-t < n*t:
yield n
nr = 10*(r-n*t)
n = ((10*(3*q+r))//t)-10*n
q *= 10
r = nr
else:
nr = (2*q+r)*l
nn = (q*(7*k)+2+(r*l))//(t*l)
q *= k
t *= l
l += 2
k += 1
n = nn
r = nr
 
import sys
pi_digits = calcPi()
i = 0
for d in pi_digits:
sys.stdout.write(str(d))
i += 1
if i == 40: print(""); i = 0
output
3141592653589793238462643383279502884197
1693993751058209749445923078164062862089
9862803482534211706798214808651328230664
7093844609550582231725359408128481117450
2841027019385211055596446229489549303819
6442881097566593344612847564823378678316
5271201909145648566923460348610454326648
2133936072602491412737245870066063155881
7488152092096282925409171536436789259036
0011330530548820466521384146951941511609
4330572703657595919530921861173819326117
...

[edit] Racket

Utilizing Jeremy Gibbons spigot algorithm and racket generator:

 
#lang racket
(require racket/generator)
 
(define pidig
(generator ()
(let loop ([q 1] [r 0] [t 1] [k 1] [n 3] [l 3])
(if (< (- (+ r (* 4 q)) t) (* n t))
(begin (yield n)
(loop (* q 10) (* 10 (- r (* n t))) t k
(- (quotient (* 10 (+ (* 3 q) r)) t) (* 10 n))
l))
(loop (* q k) (* (+ (* 2 q) r) l) (* t l) (+ 1 k)
(quotient (+ (* (+ 2 (* 7 k)) q) (* r l)) (* t l))
(+ l 2))))))
 
(for ([i (in-naturals)])
(display (pidig))
(when (zero? i) (display "." ))
(when (zero? (modulo i 80)) (newline)))
 

Output:

 
3.14159265358979323846264338327950288419716939937510...
 

[edit] REXX

Calculate digits of   π   using John Machin's formula.

It should be noted that this mechanism spits out the next (new) digits of   π,   not just a single digit.
It will spit out as many (new) digits of   π   that it finds.

The following REXX program uses the formula:

                    ┌─   ─┐                ┌─     ─┐
  π                 │  1  │                │   1   │           John
 ───  =   4 ∙ arctan│ ─── │     -    arctan│ ───── │             Machin's
  4                 │  5  │                │  239  │               formula
                    └─   ─┘                └─     ─┘

 which expands into:

     ┌─                                                                      ─┐
     │    1         1          1          1          1           1            │
4 ∙  │   ───  -  ──────  +  ──────  -  ──────  +  ──────  -  ────────  +  ... │
     │     1         3          5          7          9           11          │
     │  1∙5       3∙5        5∙5        7∙5        9∙5        11∙5            │
     └─                                                                      ─┘


     ┌─                                                                      ─┐
     │    1         1          1          1          1           1            │
 -   │   ───  -  ──────  +  ──────  -  ──────  +  ──────  -  ────────  +  ... │
     │      1         3          5          7          9           11         │
     │ 1∙239     3∙239      5∙239      7∙239      9∙239      11∙239           │
     └─                                                                      ─┘
/*REXX program spits out digits of pi (one at a time) until  Ctrl-Break.*/
arg digs .; if digs=='' then digs=1e6 /*allow the specification of digs*/
fn = 'PI_DIGITS.OUT' /*file used for output: PI digits*/
numeric digits digs /*big digs, the slower the spits.*/
pi=0; s=16; r=4; v=5; vs=v*v; g=239; gg=g*g; j=1; spit=0; old=
call time 'Reset' /*reset the REXX wall-clock timer*/
/*───calculate PI with increasing*/
do n=1 by 2 /*───accuracy (up to DIGS digits)*/
pi=pi + s/(n*v) - r/(n*g) /*───using John Machin's formula.*/
if pi==old then leave /*have exceeded DIGITS accuracy. */
s=-s; r=-r; v=v*vs; g=g*gg /*set some variable for shortcuts*/
if n\==1 then do j=spit+1 to compare(pi,old) /*spit out some π digs.*/
spit=substr(pi,j,1) /*obtain a digit of π to spit out*/
call charout ,spit /*spit out one (new) digit of pi.*/
call charout fn,spit /* ···and also echo it to a file.*/
end /*j*/
spit=j-1 /*adjust for DO index increment.*/
old=pi /*use the "OLD" value next time. */
end /*n*/
 
say; say n%2+1 'iterations took' format(time("Elapsed"),,2) 'seconds.'
/*stick a fork in it, we're done.*/

output (until the   Ctrl-Break   key was pressed):

3.141592653589793238462643383279402884197179499374105820974944592307816406286108998628034825341117067982148086513282306647093844609550582231725359408128481117450284002701938521105559644622948954920381
96442981097566593344612847564823388678316527110190914564856692346034861045432664821339360726024914127372459700660631568817488152092096282925409171536436799269035001132053054982046652138414695194151160
94330572703657595919530921861173819326118931051185480744623899627595673518857527249912279381830119491398336733624406566420860213949463952247371907021898609437027705392171762931767523846748184676794051
31000568127145263560827785771342757799609173637178721468430801225953430146559585370050892279789258923542029956012128021960863034418159813629774771319960518707211349999998372978049941059731732816096218
59402445945534690830264251230825334468503526193118817000000313783875288658753310838142061717766914730359825349042875546873116956286388235388769375295778185878053217122680661200192787661119590921642019
89380952572000654858632788659361533818289682303019510353018529699958736225994138912597217752834791315155748572424541506969508395330168617278558890750983817546374649393192550604009277016711390098488240
12858361603563707650104700181942955596199946768837459448255389774726847104047534646208046684259079491293313677029989152104752162056966024057038150193511253382420035587640247496473263914299272604269922
79678235478163500934172164121992458631502028618297455570675983850549458858692799569092720079750920295532116534498720275595023648066549912988183479775356636980742654252786255181841757467289098777289380
00816470600161452491921732172147723401414419735685481613611573525521334757418494684385233239073941433345477624168625189835695855620992192221842725502542568976717905946016534668059886272327917860857843
83827967976681454000953883786360950680064225125205117393984896084128488626945604241965285022200661186306744278622049194944047123713886950956364371917287467764657573962413990865832646995813390478027590
19946576407895126946839835259560982582262052248940772671947826848250147699090264013639443745530506820349625245174939965143142980919065925093722169646151571985838741069788595977297549892016175392846813
82686838789427741569918559252469539594310599725246808459872736446958486538367362226260991246080512438843904412441365597627807977156914359977001296160894416958685558484063534210722258285886481584560284
06016842749452267467678995252138522559954666727823986456595116354986220577456498035593634568174324112515076069479441096596094025228889710893145669136867228749940560001403208617928680920874760917824938
58900971490967598526136554978199312978481168299994872265880485756301427047755513237964145152374623436454285844489526586781105114135473574952311342716610213596953622144295248493718711014576540369027993
44037420073105785390622983974478084785896833114457138697519435064302184521900484810053706146806749192781912979399510614196634287544406437451237181921899983900159195618146751426912497489409071864942319
61567945208095146550225231603881930141093762137856956639938787083039079791077346722182562599661401421403068038447734549202605414665925101497442850732518665002132434088190700486331734649651453905796268
56000550810665879799816357473638405257145910299706413010097120628043904976951567715770041033786993600723055976317636942187212514712053292819182618612586722158929841484882916446060957527069572209175671
16722900981690915280173506712748583222871835209354965725121083579151379881091444200067500334671003141267011379908658516498315019701651511685171437657618351556508849099998599823873455283216355076489185
3589322618548963213293308985706410467525907091548141654985946163718027098199431992448895757128289069232332609729971208443357326549938239129325974636672058350414281388303103824903758985243745

[edit] Ruby

Translation of: Icon
def pi
q, r, t, k, n, l = 1, 0, 1, 1, 3, 3
dot = nil
loop do
if 4*q+r-t < n*t
yield n
if dot.nil?
yield '.'
dot = '.'
end
nr = 10*(r-n*t)
n = ((10*(3*q+r)) / t) - 10*n
q *= 10
r = nr
else
nr = (2*q+r) * l
nn = (q*(7*k+2)+r*l) / (t*l)
q *= k
t *= l
l += 2
k += 1
n = nn
r = nr
end
end
end
 
pi {|digit| print digit; $stdout.flush}

[edit] Scala

object Pi {
class PiIterator extends Iterable[BigInt]{
var r:BigInt=0
var q, t, k:BigInt=1
var n, l:BigInt=3
var nr, nn:BigInt=0
 
def iterator: Iterator[BigInt]=new Iterator[BigInt]{
def hasNext=true
def next():BigInt={
while((4*q+r-t) >= (n*t)) {
nr = (2*q+r)*l
nn = (q*(7*k)+2+(r*l))/(t*l)
q = q * k
t = t * l
l = l + 2
k = k + 1
n = nn
r = nr
}
val ret=n
nr = 10*(r-n*t)
n = ((10*(3*q+r))/t)-(10*n)
q = q * 10
r = nr
ret
}
}
}
 
def main(args: Array[String]): Unit = {
val it=new PiIterator
println((it head) + "." + (it take 300 mkString))
}
}

Output:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998
62803482534211706798214808651328230664709384460955058223172535940812848111745028410
27019385211055596446229489549303819644288109756659334461284756482337867831652712019
09145648566923460348610454326648213393607260249141273

[edit] Seed7

$ include "seed7_05.s7i";
include "bigint.s7i";
 
const proc: main is func
local
var bigInteger: q is 1_;
var bigInteger: r is 0_;
var bigInteger: t is 1_;
var bigInteger: k is 1_;
var bigInteger: n is 3_;
var bigInteger: l is 3_;
var bigInteger: nn is 0_;
var bigInteger: nr is 0_;
var boolean: first is TRUE;
begin
while TRUE do
if 4_ * q + r - t < n * t then
write(n);
if first then
write(".");
first := FALSE;
end if;
nr := 10_ * (r - n * t);
n := 10_ * (3_ * q + r) div t - 10_ * n;
q *:= 10_;
r := nr;
flush(OUT);
else
nr := (2_ * q + r) * l;
nn := (q * (7_ * k + 2_) + r * l) div (t * l);
q *:= k;
t *:= l;
l +:= 2_;
incr(k);
n := nn;
r := nr;
end if;
end while;
end func;

Original source: [5]

[edit] Sidef

__USE_INTNUM__
 
func pi(callback) {
var (q, r, t, k, n, l) = (1, 0, 1, 1, 3, 3);
while (true) {
if ((4*q + r - t) < n*t) {
callback(n);
static dot = (
callback('.'); dot;
);
var nr = 10*(r-n*t);
n = (((10*(3*q + r)) / t) - 10*n);
q *= 10;
r = nr;
}
else {
var nr = ((2*q + r) * l);
var nn = ((q*(7*k + 2) + r*l) / (t*l));
q *= k;
t *= l;
l += 2;
k += 1;
n = nn;
r = nr;
}
}
}
 
$| = 1; # autoflush
pi(func(digit){ print digit });

[edit] Tcl

Based on the reference in the D code.

Works with: Tcl version 8.6
package require Tcl 8.6
 
# http://www.cut-the-knot.org/Curriculum/Algorithms/SpigotForPi.shtml
# http://www.mathpropress.com/stan/bibliography/spigot.pdf
proc piDigitsBySpigot n {
yield [info coroutine]
set A [lrepeat [expr {int(floor(10*$n/3.)+1)}] 2]
set Alen [llength $A]
set predigits {}
while 1 {
set carry 0
for {set i $Alen} {[incr i -1] > 0} {} {
lset A $i [expr {
[set val [expr {[lindex $A $i] * 10 + $carry}]]
% [set modulo [expr {2*$i + 1}]]
}]
set carry [expr {$val / $modulo * $i}]
}
lset A 0 [expr {[set val [expr {[lindex $A 0]*10 + $carry}]] % 10}]
set predigit [expr {$val / 10}]
if {$predigit < 9} {
foreach p $predigits {yield $p}
set predigits [list $predigit]
} elseif {$predigit == 9} {
lappend predigits $predigit
} else {
foreach p $predigits {yield [incr p]}
set predigits [list 0]
}
}
}

The pi digit generation requires picking a limit to the number of digits; the bigger the limit, the more digits can be safely computed. A value of 10k yields values relatively rapidly.

coroutine piDigit piDigitsBySpigot 10000
fconfigure stdout -buffering none
while 1 {
puts -nonewline [piDigit]
}

[edit] zkl

Uses the GMP big int library. Same algorithm as many of the others on this page. Uses in place ops to cut down on big int generation (eg add vs +). Unless GC is given some hints, it will use up 16 gig quickly as it outruns the garbage collector.

var [const] BN=Import("zklBigNum"),
one=BN(1), two=BN(2), three=BN(3), four=BN(4), seven=BN(7), ten=BN(10);
 
fcn calcPiDigits{
reg q=BN(1), r=BN(0), t=BN(1), k=BN(1), n=BN(3), l=BN(3);
first:=True; N:=0;
while(True){ if((N+=1)==1000){ GarbageMan.collect(); N=0; } // take a deep breath ...
if(four*q + r - t < n*t){
n.print(); if(first){ print("."); first=False; }
nr:=(r - n*t).mul(ten); // 10 * (r - n * t);
n=(three*q).add(r).mul(ten) // ((10*(3*q + r))/t) - 10*n;
.div(t).sub(ten*n);
q.mul(ten); // q *= 10;
r=nr;
}else{
nr:=(two*q).add(r).mul(l); // (2*q + r)*l;
nn:=(q*seven).mul(k).add(two) // (q*(7*k + 2) + r*l)/(t*l);
.add(r*l).div(t*l);
q.mul(k); t.mul(l); // q*=k; t*=l;
l.add(two); k.add(one); // l+=2; k++;
n=nn; r=nr;
}
}
}();

Runs until ^C hit, the first 1000 digits match the D output.

Output:
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox