Pi: Difference between revisions

Content added Content deleted

Inline

Revision as of 12:12, 19 April 2011

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

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 ;

enum Width = 9 ; // maximum width for correct output, for type ulong enum Scale = 10UL ^^ Width ; enum InitDigit = 2UL * 10UL ^^ (Width - 1) ; enum Format = "%0" ~ to!string(Width) ~ "d" ;

struct PiDigit {

   immutable uint ndigit ;
   int opApply(int delegate(ref string /* chunk of pi digit*/) dg) {
       uint len = 10 * ndigit / 3 ;
       ulong[] arr = new ulong[](len) ;
       ulong carry ;

       arr[] = InitDigit ;

       foreach(i;0..ndigit /Width) {
           auto sum = 0UL ;
           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(Format, carry + sum/Scale) ;
           if(dg(yield))
               break ;
           carry = sum % Scale ;
       }
       return 0 ;
   }

}

void main() {

   foreach(d ; PiDigit(100))
       write(d) ;

}</lang> Output:

31415926535897932384626433832795028841971693993751
0582097494459230781640628620899862803482534211706

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>

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:

PureBasic

Calculate Pi, limited to ~24 M-digits for memory and speed reasons. <lang PureBasic>#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</lang>

Tcl

Based on the reference in the D code.

Works with: Tcl version 8.6

<lang tcl>package require Tcl 8.6

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>