I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Almkvist-Giullera formula for pi

From Rosetta Code
Almkvist-Giullera formula for pi is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The Almkvist-Giullera formula for calculating 1/π2 is based on the Calabi-Yau differential equations of order 4 and 5, which were originally used to describe certain manifolds in string theory. The formula is:

1/π2 = (25/3) ∑0 ((6n)! / (n!6))(532n2 + 126n + 9) / 10002n+1

This formula can be used to calculate the constant π-2, and thus to calculate π.

Note that, because the product of all terms but the power of 1000 can be calculated as an integer, the terms in the series can be separated into a large integer term:

(25) (6n)! (532n2 + 126n + 9) / (3(n!)6)     (***)

multiplied by a negative integer power of 10:

10-(6n + 3)
Task
  • Print the integer portions (the starred formula, which is without the power of 1000 divisor) of the first 10 terms of the series.
  • Use the complete formula to calculate and print π to 70 decimal digits of precision.
Reference



C#[edit]

A little challenging due to lack of BigFloat or BigRational. Note the extended precision integers displayed for each term, not extended precision floats. Also features the next term based on the last term, rather than computing each term from scratch. And the multiply by 32, divide by 3 is reserved for final sum, instead of each term (except for the 0..9th displayed terms).

using System;
using BI = System.Numerics.BigInteger;
using static System.Console;
 
class Program {
 
static BI isqrt(BI x) { BI q = 1, r = 0, t; while (q <= x) q <<= 2; while (q > 1) {
q >>= 2; t = x - r - q; r >>= 1; if (t >= 0) { x = t; r += q; } } return r; }
 
  static string dump(int digs, bool show = false) {
int gb = 1, dg = ++digs + gb, z;
BI t1 = 1, t2 = 9, t3 = 1, te, su = 0,
t = BI.Pow(10, dg <= 60 ? 0 : dg - 60), d = -1, fn = 1;
for (BI n = 0; n < dg; n++) {
if (n > 0) t3 *= BI.Pow(n, 6);
te = t1 * t2 / t3;
if ((z = dg - 1 - (int)n * 6) > 0) te *= BI.Pow (10, z);
else te /= BI.Pow (10, -z);
if (show && n < 10)
WriteLine("{0,2} {1,62}", n, te * 32 / 3 / t);
su += te; if (te < 10) {
if (show) WriteLine("\n{0} iterations required for {1} digits " +
"after the decimal point.\n", n, --digs); break; }
for (BI j = n * 6 + 1; j <= n * 6 + 6; j++) t1 *= j;
t2 += 126 + 532 * (d += 2);
}
string s = string.Format("{0}", isqrt(BI.Pow(10, dg * 2 + 3) /
su / 32 * 3 * BI.Pow((BI)10, dg + 5)));
return s[0] + "." + s.Substring(1, digs); }
 
static void Main(string[] args) {
WriteLine(dump(70, true)); }
}
Output:
 0  9600000000000000000000000000000000000000000000000000000000000
 1   512256000000000000000000000000000000000000000000000000000000
 2    19072247040000000000000000000000000000000000000000000000000
 3      757482485760000000000000000000000000000000000000000000000
 4       31254615037245600000000000000000000000000000000000000000
 5        1320787470322549142065152000000000000000000000000000000
 6          56727391979308908329225994240000000000000000000000000
 7           2465060024817298714011276371558400000000000000000000
 8            108065785435463945367040747443956640000000000000000
 9              4770177939159496628747057049083997888000000000000

53 iterations required for 70 digits after the decimal point.

3.1415926535897932384626433832795028841971693993751058209749445923078164

Factor[edit]

Works with: Factor version 0.99 2020-08-14
USING: continuations formatting io kernel locals math
math.factorials math.functions sequences ;
 
:: integer-term ( n -- m )
32 6 n * factorial * 532 n sq * 126 n * + 9 + *
n factorial 6 ^ 3 * / ;
 
: exponent-term ( n -- m ) 6 * 3 + neg ;
 
: nth-term ( n -- x )
[ integer-term ] [ exponent-term 10^ * ] bi ;
 
! Factor doesn't have an arbitrary-precision square root afaik,
! so make one using Heron's method.
 
: sqrt-approx ( r x -- r' x ) [ over / + 2 / ] keep ;
 
:: almkvist-guillera ( precision -- x )
0 0 :> ( summed! next-add! )
[
100,000,000 <iota> [| n |
summed n nth-term + next-add!
next-add summed - abs precision neg 10^ <
[ return ] when
next-add summed!
] each
] with-return
next-add ;
 
CONSTANT: 1/pi 113/355  ! Use as initial guess for square root approximation
 
: pi ( -- )
1/pi 70 almkvist-guillera 5 [ sqrt-approx ] times
drop recip "%.70f\n" printf ;
 
! Task
"N Integer Portion Pow Nth Term (33 dp)" print
89 CHAR: - <repetition> print
10 [
dup [ integer-term ] [ exponent-term ] [ nth-term ] tri
"%d  %44d  %3d  %.33f\n" printf
] each-integer nl
"Pi to 70 decimal places:" print pi
Output:
N                               Integer Portion  Pow  Nth Term (33 dp)
-----------------------------------------------------------------------------------------
0                                            96   -3  0.096000000000000000000000000000000
1                                       5122560   -9  0.005122560000000000000000000000000
2                                  190722470400  -15  0.000190722470400000000000000000000
3                              7574824857600000  -21  0.000007574824857600000000000000000
4                         312546150372456000000  -27  0.000000312546150372456000000000000
5                    13207874703225491420651520  -33  0.000000013207874703225491420651520
6                567273919793089083292259942400  -39  0.000000000567273919793089083292260
7           24650600248172987140112763715584000  -45  0.000000000024650600248172987140113
8      1080657854354639453670407474439566400000  -51  0.000000000001080657854354639453670
9  47701779391594966287470570490839978880000000  -57  0.000000000000047701779391594966287

Pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164

Go[edit]

Translation of: Wren
package main
 
import (
"fmt"
"math/big"
"strings"
)
 
func factorial(n int64) *big.Int {
var z big.Int
return z.MulRange(1, n)
}
 
var one = big.NewInt(1)
var three = big.NewInt(3)
var six = big.NewInt(6)
var ten = big.NewInt(10)
var seventy = big.NewInt(70)
 
func almkvistGiullera(n int64, print bool) *big.Rat {
t1 := big.NewInt(32)
t1.Mul(factorial(6*n), t1)
t2 := big.NewInt(532*n*n + 126*n + 9)
t3 := new(big.Int)
t3.Exp(factorial(n), six, nil)
t3.Mul(t3, three)
ip := new(big.Int)
ip.Mul(t1, t2)
ip.Quo(ip, t3)
pw := 6*n + 3
t1.SetInt64(pw)
tm := new(big.Rat).SetFrac(ip, t1.Exp(ten, t1, nil))
if print {
fmt.Printf("%d  %44d  %3d  %-35s\n", n, ip, -pw, tm.FloatString(33))
}
return tm
}
 
func main() {
fmt.Println("N Integer Portion Pow Nth Term (33 dp)")
fmt.Println(strings.Repeat("-", 89))
for n := int64(0); n < 10; n++ {
almkvistGiullera(n, true)
}
 
sum := new(big.Rat)
prev := new(big.Rat)
pow70 := new(big.Int).Exp(ten, seventy, nil)
prec := new(big.Rat).SetFrac(one, pow70)
n := int64(0)
for {
term := almkvistGiullera(n, false)
sum.Add(sum, term)
z := new(big.Rat).Sub(sum, prev)
z.Abs(z)
if z.Cmp(prec) < 0 {
break
}
prev.Set(sum)
n++
}
sum.Inv(sum)
pi := new(big.Float).SetPrec(256).SetRat(sum)
pi.Sqrt(pi)
fmt.Println("\nPi to 70 decimal places is:")
fmt.Println(pi.Text('f', 70))
}
Output:
N                               Integer Portion  Pow  Nth Term (33 dp)
-----------------------------------------------------------------------------------------
0                                            96   -3  0.096000000000000000000000000000000
1                                       5122560   -9  0.005122560000000000000000000000000
2                                  190722470400  -15  0.000190722470400000000000000000000
3                              7574824857600000  -21  0.000007574824857600000000000000000
4                         312546150372456000000  -27  0.000000312546150372456000000000000
5                    13207874703225491420651520  -33  0.000000013207874703225491420651520
6                567273919793089083292259942400  -39  0.000000000567273919793089083292260
7           24650600248172987140112763715584000  -45  0.000000000024650600248172987140113
8      1080657854354639453670407474439566400000  -51  0.000000000001080657854354639453670
9  47701779391594966287470570490839978880000000  -57  0.000000000000047701779391594966287

Pi to 70 decimal places is:
3.1415926535897932384626433832795028841971693993751058209749445923078164

Julia[edit]

using Formatting
 
setprecision(BigFloat, 300)
 
function integerterm(n)
p = BigInt(532) * n * n + BigInt(126) * n + 9
return (p * BigInt(2)^5 * factorial(BigInt(6) * n)) ÷ (3 * factorial(BigInt(n))^6)
end
 
exponentterm(n) = -(6n + 3)
 
nthterm(n) = integerterm(n) * big"10.0"^exponentterm(n)
 
println(" N Integer Term Power of 10 Nth Term")
println("-"^90)
for n in 0:9
println(lpad(n, 3), lpad(integerterm(n), 48), lpad(exponentterm(n), 4),
lpad(format("{1:22.19e}", nthterm(n)), 35))
end
 
function AlmkvistGuillera(floatprecision)
summed = nthterm(0)
for n in 1:10000000
next = summed + nthterm(n)
if abs(next - summed) < big"10.0"^(-floatprecision)
return next
end
summed = next
end
end
 
println("\nπ to 70 digits is ", format(big"1.0" / sqrt(AlmkvistGuillera(70)), precision=70))
 
println("Computer π is ", format(π + big"0.0", precision=70))
 
Output:
  N                       Integer Term              Power of 10     Nth Term
------------------------------------------------------------------------------------------
  0                                              96  -3          9.6000000000000000000e-02
  1                                         5122560  -9          5.1225600000000000000e-03
  2                                    190722470400 -15          1.9072247040000000000e-04
  3                                7574824857600000 -21          7.5748248576000000000e-06
  4                           312546150372456000000 -27          3.1254615037245600000e-07
  5                      13207874703225491420651520 -33          1.3207874703225491421e-08
  6                  567273919793089083292259942400 -39          5.6727391979308908329e-10
  7             24650600248172987140112763715584000 -45          2.4650600248172987140e-11
  8        1080657854354639453670407474439566400000 -51          1.0806578543546394537e-12
  9    47701779391594966287470570490839978880000000 -57          4.7701779391594966287e-14

π to 70 digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164
Computer π is     3.1415926535897932384626433832795028841971693993751058209749445923078164

Phix[edit]

requires("0.8.2") 
include mpfr.e
mpfr_set_default_prec(-70)
 
function almkvistGiullera(integer n, bool bPrint)
mpz {t1,t2,ip} = mpz_inits(3)
mpz_fac_ui(t1,6*n)
mpz_mul_si(t1,t1,32) -- t1:=2^5*(6n)!
mpz_fac_ui(t2,n)
mpz_pow_ui(t2,t2,6)
mpz_mul_si(t2,t2,3) -- t2:=3*(n!)^6
mpz_mul_si(ip,t1,532*n*n+126*n+9) -- ip:=t1*(532n^2+126n+9)
mpz_fdiv_q(ip,ip,t2) -- ip:=ip/t2
integer pw := 6*n+3
mpz_ui_pow_ui(t1,10,pw) -- t1 := 10^(6n+3)
mpq tm = mpq_init_set_z(ip,t1) -- tm := rat(ip/t1)
if bPrint then
string ips = mpz_get_str(ip),
tms = mpfr_sprintf("%.50Rf",mpfr_init_set_q(tm))
tms = trim_tail(tms,"0")
printf(1,"%d  %44s  %3d  %s\n", {n, ips, -pw, tms})
end if
return tm
end function
 
constant hdr = "N --------------- Integer portion ------------- Pow ----------------- Nth term (50 dp) -----------------"
printf(1,"%s\n%s\n",{hdr,repeat('-',length(hdr))})
for n=0 to 9 do
{} = almkvistGiullera(n, true)
end for
 
mpq {res,prev,z} = mpq_inits(3),
prec = mpq_init_set_str(sprintf("1/1%s",repeat('0',70)))
integer n = 0
while true do
mpq term := almkvistGiullera(n, false)
mpq_add(res,res,term)
mpq_sub(z,res,prev)
mpq_abs(z,z)
if mpq_cmp(z,prec) < 0 then
exit
end if
mpq_set(prev,res)
n += 1
end while
mpq_inv(res,res)
mpfr pi = mpfr_init_set_q(res)
mpfr_sqrt(pi,pi)
printf(1,"\nCalculation of pi took %d iterations using the Almkvist-Giullera formula.\n\n",n)
printf(1,"Pi to 70 d.p.: %s\n",mpfr_sprintf("%.70Rf",pi))
mpfr_const_pi(pi)
printf(1,"Pi (builtin) : %s\n",mpfr_sprintf("%.70Rf",pi))
Output:
N --------------- Integer portion -------------  Pow  ----------------- Nth term (50 dp) -----------------
----------------------------------------------------------------------------------------------------------
0                                            96   -3  0.096
1                                       5122560   -9  0.00512256
2                                  190722470400  -15  0.0001907224704
3                              7574824857600000  -21  0.0000075748248576
4                         312546150372456000000  -27  0.000000312546150372456
5                    13207874703225491420651520  -33  0.00000001320787470322549142065152
6                567273919793089083292259942400  -39  0.0000000005672739197930890832922599424
7           24650600248172987140112763715584000  -45  0.000000000024650600248172987140112763715584
8      1080657854354639453670407474439566400000  -51  0.0000000000010806578543546394536704074744395664
9  47701779391594966287470570490839978880000000  -57  0.00000000000004770177939159496628747057049083997888

Calculation of pi took 52 iterations using the Almkvist-Giullera formula.

Pi to 70 d.p.: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Pi (builtin) : 3.1415926535897932384626433832795028841971693993751058209749445923078164

Python[edit]

import mpmath as mp
 
with mp.workdps(72):
 
def integer_term(n):
p = 532 * n * n + 126 * n + 9
return (p * 2**5 * mp.factorial(6 * n)) / (3 * mp.factorial(n)**6)
 
def exponent_term(n):
return -(mp.mpf("6.0") * n + 3)
 
def nthterm(n):
return integer_term(n) * mp.mpf("10.0")**exponent_term(n)
 
 
for n in range(10):
print("Term ", n, ' ', int(integer_term(n)))
 
 
def almkvist_guillera(floatprecision):
summed, nextadd = mp.mpf('0.0'), mp.mpf('0.0')
for n in range(100000000):
nextadd = summed + nthterm(n)
if abs(nextadd - summed) < 10.0**(-floatprecision):
break
 
summed = nextadd
 
return nextadd
 
 
print('\nπ to 70 digits is ', end='')
mp.nprint(mp.mpf(1.0 / mp.sqrt(almkvist_guillera(70))), 71)
print('mpmath π is ', end='')
mp.nprint(mp.pi, 71)
 
Output:
Term  0    96
Term  1    5122560
Term  2    190722470400
Term  3    7574824857600000
Term  4    312546150372456000000
Term  5    13207874703225491420651520
Term  6    567273919793089083292259942400
Term  7    24650600248172987140112763715584000
Term  8    1080657854354639453670407474439566400000
Term  9    47701779391594966287470570490839978880000000

π to 70 digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164
mpmath π is       3.1415926535897932384626433832795028841971693993751058209749445923078164

Raku[edit]

# 20201013 Raku programming solution
 
use BigRoot;
use Rat::Precise;
use experimental :cached;
 
BigRoot.precision = 75 ;
my $Precision = 70 ;
my $AGcache = 0 ;
 
sub postfix:<!>(Int $n --> Int) is cached { [*] 1 .. $n }
 
sub Integral(Int $n --> Int) is cached {
(2*(6*$n)! * (532*$n² + 126*$n + 9)) div (3*($n!))
}
 
sub A-G(Int $n --> FatRat) is cached { # Almkvist-Giullera
Integral($n).FatRat / (10**(6*$n + 3)).FatRat
}
 
sub Pi(Int $n --> Str) {
(1/(BigRoot.newton's-sqrt: $AGcache += A-G $n)).precise($Precision)
}
 
say "First 10 integer portions : ";
say $_, "\t", Integral $_ for ^10;
 
my $target = Pi my $Nth = 0;
 
loop { $target eq ( my $next = Pi ++$Nth ) ?? ( last ) !! $target = $next }
 
say "π to $Precision decimal places is :\n$target"
Output:
First 10 integer portions :
0       96
1       5122560
2       190722470400
3       7574824857600000
4       312546150372456000000
5       13207874703225491420651520
6       567273919793089083292259942400
7       24650600248172987140112763715584000
8       1080657854354639453670407474439566400000
9       47701779391594966287470570490839978880000000
π to 70 decimal places is :
3.1415926535897932384626433832795028841971693993751058209749445923078164

REXX[edit]

/*REXX program uses the Almkvist─Giullera formula for   1 / (pi**2)     [or  pi ** -2]. */
numeric digits length(pi() ) + 1; w= 58
say $( , 3) $( , 44) $('power', 7) $( , w)
say $('N', 3) $('integer term', 44) $('of 10', 7) $('Nth term', w)
say $( , 3, "─") $( , 44, "─") $( , 7, "─") $( , w, "─")
s= 0 /*initialize S (the sum) to zero. */
do n=0 until old=s; old= s /*use the "older" value of S for OLD.*/
a= 2**5 *  !(6*n) * (532 * n**2 + 126*n + 9) / (3 * !(n)**6)
z= 10 ** -(6*n + 3)
s= s + a * z
if n>9 then iterate
say $(n, 3) right(a, 44) $(powX(z), 7) right( lowE( format(a*z, 1, w-6, 2, 0)), w)
end /*n*/
say
say 'calculation of pi took ' n " iterations using" subword( sourceLine(1), 4, 3).
say
numeric digits length(pi() ) - length(.)
say 'calc pi to' digits() - length(.) "fractional decimal digits is" sqrt(1/s)
say 'true pi to' digits() - length(.) "fractional decimal digits is" pi()
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
$: procedure; parse arg text,width,fill; return center(text, width, left(fill, 1) )
!: procedure; parse arg x; !=1; do j=2 to x;  != !*j; end; return !
lowE: procedure; parse arg x; return translate(x, 'e', "E")
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164
powX: procedure; parse arg p; return right( format( p, 1, 3, 2, 0), 3) + 0
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
m.=9; numeric form; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_ % 2
do j=0 while h>9; m.j= h; h= h % 2 + 1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g= (g + x/g) * .5; end /*k*/
numeric digits d; return g/1
output   when using the internal default input:
                                                  power
 N                  integer term                  of 10                           Nth term
─── ──────────────────────────────────────────── ─────── ──────────────────────────────────────────────────────────
 0                                            96   -3    9.6000000000000000000000000000000000000000000000000000e-02
 1                                       5122560   -9    5.1225600000000000000000000000000000000000000000000000e-03
 2                                  190722470400   -15   1.9072247040000000000000000000000000000000000000000000e-04
 3                              7574824857600000   -21   7.5748248576000000000000000000000000000000000000000000e-06
 4                         312546150372456000000   -27   3.1254615037245600000000000000000000000000000000000000e-07
 5                    13207874703225491420651520   -33   1.3207874703225491420651520000000000000000000000000000e-08
 6                567273919793089083292259942400   -39   5.6727391979308908329225994240000000000000000000000000e-10
 7           24650600248172987140112763715584000   -45   2.4650600248172987140112763715584000000000000000000000e-11
 8      1080657854354639453670407474439566400000   -51   1.0806578543546394536704074744395664000000000000000000e-12
 9  47701779391594966287470570490839978880000000   -57   4.7701779391594966287470570490839978880000000000000000e-14

calculation of pi took  54  iterations using the Almkvist─Giullera formula.

calc pi to 70 fractional decimal digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164
true pi to 70 fractional decimal digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164

Visual Basic .NET[edit]

Translation of: C#
Imports System, BI = System.Numerics.BigInteger, System.Console
 
Module Module1
 
Function isqrt(ByVal x As BI) As BI
Dim t As BI, q As BI = 1, r As BI = 0
While q <= x : q <<= 2 : End While
While q > 1 : q >>= 2 : t = x - r - q : r >>= 1
If t >= 0 Then x = t : r += q
End While : Return r
End Function
 
Function dump(ByVal digs As Integer, ByVal Optional show As Boolean = False) As String
digs += 1
Dim z As Integer, gb As Integer = 1, dg As Integer = digs + gb
Dim te As BI, t1 As BI = 1, t2 As BI = 9, t3 As BI = 1, su As BI = 0, t As BI = BI.Pow(10, If(dg <= 60, 0, dg - 60)), d As BI = -1, fn As BI = 1
For n As BI = 0 To dg - 1
If n > 0 Then t3 = t3 * BI.Pow(n, 6)
te = t1 * t2 / t3 : z = dg - 1 - CInt(n) * 6
If z > 0 Then te = te * BI.Pow(10, z) Else te = te / BI.Pow(10, -z)
If show AndAlso n < 10 Then WriteLine("{0,2} {1,62}", n, te * 32 / 3 / t)
su += te : If te < 10 Then
digs -= 1
If show Then WriteLine(vbLf & "{0} iterations required for {1} digits " & _
"after the decimal point." & vbLf, n, digs)
Exit For
End If
For j As BI = n * 6 + 1 To n * 6 + 6
t1 = t1 * j : Next
d += 2 : t2 += 126 + 532 * d
Next
Dim s As String = String.Format("{0}", isqrt(BI.Pow(10, dg * 2 + 3) _
/ su / 32 * 3 * BI.Pow(CType(10, BI), dg + 5)))
Return s(0) & "." & s.Substring(1, digs)
End Function
 
Sub Main(ByVal args As String())
WriteLine(dump(70, true))
End Sub
 
End Module
Output:
 0  9600000000000000000000000000000000000000000000000000000000000
 1   512256000000000000000000000000000000000000000000000000000000
 2    19072247040000000000000000000000000000000000000000000000000
 3      757482485760000000000000000000000000000000000000000000000
 4       31254615037245600000000000000000000000000000000000000000
 5        1320787470322549142065152000000000000000000000000000000
 6          56727391979308908329225994240000000000000000000000000
 7           2465060024817298714011276371558400000000000000000000
 8            108065785435463945367040747443956640000000000000000
 9              4770177939159496628747057049083997888000000000000

53 iterations required for 70 digits after the decimal point.

3.1415926535897932384626433832795028841971693993751058209749445923078164

Wren[edit]

Library: Wren-big
Library: Wren-fmt
import "/big" for BigInt, BigRat
import "/fmt" for Fmt
 
var factorial = Fn.new { |n|
if (n < 2) return BigInt.one
var fact = BigInt.one
for (i in 2..n) fact = fact * i
return fact
}
 
var almkvistGiullera = Fn.new { |n, print|
var t1 = factorial.call(6*n) * 32
var t2 = 532*n*n + 126*n + 9
var t3 = factorial.call(n).pow(6)*3
var ip = t1 * t2 / t3
var pw = 6*n + 3
var tm = BigRat.new(ip, BigInt.ten.pow(pw))
if (print) {
Fmt.print("$d $44i $3d $-35s", n, ip, -pw, tm.toDecimal(33))
} else {
return tm
}
}
 
System.print("N Integer Portion Pow Nth Term (33 dp)")
System.print("-" * 89)
for (n in 0..9) {
almkvistGiullera.call(n, true)
}
 
var sum = BigRat.zero
var prev = BigRat.zero
var prec = BigRat.new(BigInt.one, BigInt.ten.pow(70))
var n = 0
while(true) {
var term = almkvistGiullera.call(n, false)
sum = sum + term
if ((sum-prev).abs < prec) break
prev = sum
n = n + 1
}
var pi = BigRat.one/sum.sqrt(70)
System.print("\nPi to 70 decimal places is:")
System.print(pi.toDecimal(70))
Output:
N                               Integer Portion  Pow  Nth Term (33 dp)
-----------------------------------------------------------------------------------------
0                                            96   -3  0.096                              
1                                       5122560   -9  0.00512256                         
2                                  190722470400  -15  0.0001907224704                    
3                              7574824857600000  -21  0.0000075748248576                 
4                         312546150372456000000  -27  0.000000312546150372456            
5                    13207874703225491420651520  -33  0.00000001320787470322549142065152 
6                567273919793089083292259942400  -39  0.000000000567273919793089083292260
7           24650600248172987140112763715584000  -45  0.000000000024650600248172987140113
8      1080657854354639453670407474439566400000  -51  0.000000000001080657854354639453670
9  47701779391594966287470570490839978880000000  -57  0.000000000000047701779391594966287

Pi to 70 decimal places is:
3.1415926535897932384626433832795028841971693993751058209749445923078164