Pi

Pi
You are encouraged to solve this task according to the task description, using any language you may know.

The task is to 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 ...

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: <lang Ada>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;</lang>

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

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. <lang algol68>#!/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 #

1. 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)); 

1. iterate throught the digits as they are being found #
2. FOR INT digit IN # gen pi digits(test decimal places#) DO ( #,

 ## (INT digit)VOID: (
   printf(($n(width)d$,digit))
 )

1. OD #);

 print(new line)

)</lang> Output:

3141592653589793238462643383279502

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

<lang 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

}</lang>

Output:

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

C

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. <lang C>#include <stdio.h>

1. include <stdlib.h>
2. 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; }</lang>

C#

<lang csharp>using System;

namespace RC__Pi {

 internal sealed class Pi
 {
   private const int Scale = 10000;
   private const int ArrInit = 2000;

   private static void Main()
   {
     Console.Write("Enter digits to calculate: (10000) ");
     string inp = Console.ReadLine();
     if (inp==null) return;
     long digits = (inp.Length>0) ? Convert.ToInt64(inp) : 10000;
     long carry = 0;
     var arr = new long[digits+1];

     for (long i = 0; i<=digits; ++i)
       arr[i] = ArrInit;

     for (long i = digits; i>0; i -= 4)
     {
       long sum = 0;
       for (long j = i; j>0; --j)
       {
         sum = sum*j+Scale*arr[j];
         arr[j] = sum%(j*2-1);
         sum /= j*2-1;
       }
       Console.Write((carry+sum/Scale).ToString("D4"));
       carry = sum%Scale;
     }
   }
 }

}</lang>

D

This modified Spigot algorithm does not continue infinitely, because its required memory grow as the number of digits need to print. <lang d>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);

}</lang> Output:

314159265
358979323
846264338
327950288
419716939
937510582
097494459
230781640
628620899
862803482
534211706

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.) <lang go>package main

import (

   "big"
   "fmt"

)

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

}</lang>

Haskell

The code from [1]: <lang haskell>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)</lang>

Icon and Unicon

Translation of: PicoLisp

based on Jeremy Gibbons' Haskell solution.

<lang icon>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</lang>

J

<lang j>pi=:3 :0

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

)</lang>

Example use:

<lang j> pi 3 . 1 4 1 5 9 2 6 5 3 ...</lang>

Java

Translation of: Icon

This is a copy of the icon/unicon implementation translated into java... <lang java>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() ;
 }

}</lang>

Output :

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480 ...

OCaml

The Constructive Real library Creal contains an infinite-precision Pi, so we can just print out its digits. <lang OCaml>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</lang>

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. <lang OCaml>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</lang>

PARI/GP

Uses the built-in Brent-Salamin arithmetic-geometric mean iteration. <lang parigp>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)
 )

};</lang>

Perl

Using Machin's formula. Not the fastest program in the world. <lang Perl>use Math::BigInt; use bigint;

1. Pi/4 = 4 arctan 1/5 - arctan 1/239
2. 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; } }</lang>

PicoLisp

The following script uses the spigot algorithm published by Jeremy Gibbons. Hit Ctrl-C to stop it. <lang PicoLisp>#!/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) )</lang>

Output:

3.14159265358979323846264338327950288419716939937510582097494459 ...

PureBasic

Calculate Pi, limited to ~24 M-digits for memory and speed reasons. <lang PureBasic>#SCALE = 10000

1. 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</lang>

Python

<lang 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</lang>output<lang>3141592653589793238462643383279502884197

1693993751058209749445923078164062862089 9862803482534211706798214808651328230664 7093844609550582231725359408128481117450 2841027019385211055596446229489549303819 6442881097566593344612847564823378678316 5271201909145648566923460348610454326648 2133936072602491412737245870066063155881 7488152092096282925409171536436789259036 0011330530548820466521384146951941511609 4330572703657595919530921861173819326117 ...</lang>

Ruby

Translation of: Icon

<lang ruby>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}</lang>

Tcl

Based on the reference in the D code.

Works with: Tcl version 8.6

<lang tcl>package require Tcl 8.6

1. http://www.cut-the-knot.org/Curriculum/Algorithms/SpigotForPi.shtml
2. 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] }

   }

}</lang> 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. <lang tcl>coroutine piDigit piDigitsBySpigot 10000 fconfigure stdout -buffering none while 1 {

   puts -nonewline [piDigit]

}</lang>