Almkvist-Giullera formula for pi

The Almkvist-Giullera formula for calculating 1/π² 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⁵/3) ∑₀^∞ ((6n)! / (n!⁶))(532n² + 126n + 9) / 1000²ⁿ⁺¹

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:

(2⁵) (6n)! (532n² + 126n + 9) / (3(n!)⁶) (***)

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

Gert Almkvist and Jesús Guillera, Ramanujan-like series for 1/π² and string theory, Experimental Mathematics, 21 (2012), page 2, formula 1.

11l

Translation of: C#

<lang 11l>F isqrt(BigInt =x)

  BigInt q = 1
  BigInt r = 0
  BigInt t
  L q <= x
     q *= 4
  L q > 1
     q I/= 4
     t = x - r - q
     r I/= 2
     I t >= 0
        x = t
        r += q
  R r

F dump(=digs, show)

  V gb = 1
  digs++
  V dg = digs + gb
  BigInt t1 = 1
  BigInt t2 = 9
  BigInt t3 = 1
  BigInt te
  BigInt su = 0
  V t = BigInt(10) ^ (I dg <= 60 {0} E dg - 60)
  BigInt d = -1
  BigInt _fn_ = 1

  V n = 0
  L n < dg
     I n > 0
        t3 *= BigInt(n) ^ 6
     te = t1 * t2 I/ t3
     V z = dg - 1 - n * 6
     I z > 0
        te *= BigInt(10) ^ z
     E
        te I/= BigInt(10) ^ -z
     I show & n < 10
        print(‘#2 #62’.format(n, te * 32 I/ 3 I/ t))
     su += te

     I te < 10
        I show
           digs--
           print("\n#. iterations required for #. digits after the decimal point.\n".format(n, digs))
        L.break

     L(j) n * 6 + 1 .. n * 6 + 6
        t1 *= j
     d += 2
     t2 += 126 + 532 * d

n++

  V s = String(isqrt(BigInt(10) ^ (dg * 2 + 3) I/ su I/ 32 * 3 * BigInt(10) ^ (dg + 5)))
  R s[0]‘.’s[1 .+ digs]

print(dump(70, 1B))</lang>

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

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). <lang csharp>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)); }

}</lang>

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

Erlang

This version uses integer math only (does not resort to a rational number package) Since the denominator is always a power of 10, it's possible to just keep track of the log of the denominator and scale the numerator accordingly; to keep track of the accuracy we get the order of magnitude of the difference between terms by subtracting the log of the numerator from the log of the denominator, so again, no rational arithmetic is needed.

However, Erlang does not have much in the way of calculating with large integers beyond basic arithmetic, so this version implements integer powers, logs, square roots, as well as the factorial function. <lang Erlang> -mode(compile).

% Integer math routines: factorial, power, square root, integer logarithm. % fac(N) -> fac(N, 1). fac(N, A) when N < 2 -> A; fac(N, A) -> fac(N - 1, N*A).

pow(_, N) when N < 0 -> pow_domain_error; pow(2, N) -> 1 bsl N; pow(A, N) -> ipow(A, N).

ipow(_, 0) -> 1; ipow(A, 1) -> A; ipow(A, 2) -> A*A; ipow(A, N) ->

   case N band 1 of
       0 -> X = ipow(A, N bsr 1), X*X;
       1 -> A * ipow(A, N - 1)
   end.

% integer logarithm, based on Zeckendorf representations of integers. % https://www.keithschwarz.com/interesting/code/?dir=zeckendorf-logarithm % we need this, since the exponents get larger than IEEE-754 double can handle.

lognext({A, B, S, T}) -> {B, A+B, T, S*T}. logprev({A, B, S, T}) -> {B-A, A, T div S, S}.

ilog(A, B) when (A =< 0) or (B < 2) -> ilog_domain_error; ilog(A, B) ->

   UBound = bracket(A, {0, 1, 1, B}),
   backlog(A, UBound, 0).

bracket(A, State = {_, _, _, T}) when T > A -> State; bracket(A, State) -> bracket(A, lognext(State)).

backlog(_, {0, _, 1, _}, Log) -> Log; backlog(N, State = {A, _, S, _}, Log) when S =< N ->

   backlog(N div S, logprev(State), Log + A);

backlog(N, State, Log) -> backlog(N, logprev(State), Log).

isqrt(N) when N < 0 -> isqrt_domain_error; isqrt(N) ->

   X0 = pow(2, ilog(N, 2) div 2),
   iterate(N, newton(X0, N), N).

iterate(A, B, _) when A =< B -> A; iterate(_, B, N) -> iterate(B, newton(B, N), N).

newton(X, N) -> (X + N div X) div 2.

% With this out of the way, we can get down to some serious calculation. % term(N) -> { % returns numerator and log10 of the denominator.

   (fac(6*N)*(N*(532*N + 126) + 9) bsl 5) div (3*pow(fac(N), 6)),
   6*N + 3
   }.

neg_term({N, D}) -> {-N, D}. abs_term({N, D}) -> {abs(N), D}.

add_term(T1 = {_, D1}, T2 = {_, D2}) when D1 > D2 -> add_term(T2, T1); add_term({N1, D1}, {N2, D2}) ->

   Scale = pow(10, D2 - D1),
   {N1*Scale + N2, D2}.

calculate(Prec) -> calculate(Prec, {0, 0}, 0). calculate(Prec, T0, K) ->

   T1 = add_term(T0, term(K)),
   {N, D} = abs_term(add_term(neg_term(T1), T0)),
   Accuracy = D - ilog(N, 10),
   if
       Accuracy < Prec -> calculate(Prec, T1, K + 1);
       true -> T1
   end.

get_pi(Prec) ->

   {N0, D0} = calculate(Prec),
   % from the term, t = n0/10^d0, calculate 1/√t
   % if the denominator is an odd power of 10, add 1 to the denominator and multiply the numerator by 10.
   {N, D} = case D0 band 1 of
       0 -> {N0, D0};
       1 -> {10*N0, D0 + 1}
   end,
   [Three|Rest] = lists:sublist(
           integer_to_list(pow(10, D) div isqrt(N)), Prec),
   [Three, $. | Rest].

show_term({A, Decimals}) ->

   Str = integer_to_list(A),
   [$0, $.] ++ lists:duplicate(Decimals - length(Str), $0) ++ Str.

main(_) ->

   Terms = [term(N) || N <- lists:seq(0, 9)],
   io:format("The first 10 terms as scaled decimals are:~n"),
   [io:format("    ~s~n", [show_term(T)]) || T <- Terms],
   io:format("~nThe sum of these terms (pi^-2) is ~s~n",
               [show_term(lists:foldl(fun add_term/2, {0, 0}, Terms))]),
   Pi70 = get_pi(71),
   io:format("~npi to 70 decimal places:~n"),
   io:format("~s~n", [Pi70]).

</lang>

Output:

The first 10 terms as scaled decimals are:
    0.096
    0.005122560
    0.000190722470400
    0.000007574824857600000
    0.000000312546150372456000000
    0.000000013207874703225491420651520
    0.000000000567273919793089083292259942400
    0.000000000024650600248172987140112763715584000
    0.000000000001080657854354639453670407474439566400000
    0.000000000000047701779391594966287470570490839978880000000

The sum of these terms (pi^-2) is 0.101321183642335555356499725503850584160514406378880000000

pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164

Factor

Works with: Factor version 0.99 2020-08-14

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

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

Translation of: Wren

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

}</lang>

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

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

</lang>

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

Mathematica/Wolfram Language

<lang Mathematica>ClearAll[numerator, denominator] numerator[n_] := (2^5) ((6 n)!) (532 n^2 + 126 n + 9)/(3 (n!)^6) denominator[n_] := 10^(6 n + 3) numerator /@ Range[0, 9] val = 1/Sqrt[Total[numerator[#]/denominator[#] & /@ Range[0, 100]]]; N[val, 70]</lang>

Output:

{96,5122560,190722470400,7574824857600000,312546150372456000000,13207874703225491420651520,567273919793089083292259942400,24650600248172987140112763715584000,1080657854354639453670407474439566400000,47701779391594966287470570490839978880000000}
3.141592653589793238462643383279502884197169399375105820974944592307816

Nim

Library: nim-decimal

Derived from Wren version with some simplifications. <lang Nim>import strformat, strutils import decimal

proc fact(n: int): DecimalType =

 result = newDecimal(1)
 if n < 2: return
 for i in 2..n:
   result *= i

proc almkvistGiullera(n: int): DecimalType =

 ## Return the integer portion of the nth term of Almkvist-Giullera sequence.
 let t1 = fact(6 * n) * 32
 let t2 = 532 * n * n + 126 * n + 9
 let t3 = fact(n) ^ 6 * 3
 result = t1 * t2 / t3

let One = newDecimal(1)

setPrec(78) echo "N Integer portion" echo repeat('-', 47) for n in 0..9:

 echo &"{n}  {almkvistGiullera(n):>44}"

echo()

echo "Pi to 70 decimal places:" var

 sum = newDecimal(0)
 prev = newDecimal(0)
 prec = One.scaleb(newDecimal(-70))
 n = 0

while true:

 sum += almkvistGiullera(n) / One.scaleb(newDecimal(6 * n + 3))
 if abs(sum - prev) < prec: break
 prev = sum.clone
 inc n

let pi = 1 / sqrt(sum) echo ($pi)[0..71]</lang>

Output:

N                               Integer portion
-----------------------------------------------
0                                            96
1                                       5122560
2                                  190722470400
3                              7574824857600000
4                         312546150372456000000
5                    13207874703225491420651520
6                567273919793089083292259942400
7           24650600248172987140112763715584000
8      1080657854354639453670407474439566400000
9  47701779391594966287470570490839978880000000

Pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164

Perl

Translation of: Raku

<lang perl>use strict; use warnings; use feature qw(say); use Math::AnyNum qw(:overload factorial);

sub almkvist_giullera_integral {

   my($n) = @_;
   (32 * (14*$n * (38*$n + 9) + 9) * factorial(6*$n)) / (3*factorial($n)**6);

}

sub almkvist_giullera {

   my($n) = @_;
   almkvist_giullera_integral($n) / (10**(6*$n + 3));

}

sub almkvist_giullera_pi {

   my ($prec) = @_;

   local $Math::AnyNum::PREC = 4*($prec+1);

   my $sum = 0;
   my $target = ;

   for (my $n = 0; ; ++$n) {
       $sum += almkvist_giullera($n);
       my $curr = ($sum**-.5)->as_dec;
       return $target if ($curr eq $target);
       $target = $curr;
   }

}

say 'First 10 integer portions: '; say "$_ " . almkvist_giullera_integral($_) for 0..9;

my $precision = 70;

printf("π to %s decimal places is:\n%s\n",

   $precision, almkvist_giullera_pi($precision));</lang>

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

Phix

with javascript_semantics
requires("1.0.0") 
include mpfr.e
mpfr_set_default_precision(-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_get_fixed(mpfr_init_set_q(tm),50)
        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_get_fixed(pi,70))
mpfr_const_pi(pi)
printf(1,"Pi (builtin) : %s\n",mpfr_get_fixed(pi,70))

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

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

</lang>

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

Quackery

<lang Quackery> [ $ "bigrat.qky" loadfile ] now!

 [ 1 swap times [ i^ 1+ * ] ] is !       ( n --> n   )

 [ dup dup 2 ** 532 *
   over 126 * + 9 +
   swap 6 * ! * 32 *
   swap ! 6 ** 3 * / ]        is intterm ( n --> n   )

 [ dup intterm 
   10 rot 6 * 3 + ** 
   reduce ]                   is vterm   ( n --> n/d )

 10 times [ i^ intterm echo cr ] cr
 
 0 n->v 
 53 times [ i^ vterm v+ ]
 1/v 70 vsqrt drop 
 70 point$ echo$ cr</lang>

Output:

96
5122560
190722470400
7574824857600000
312546150372456000000
13207874703225491420651520
567273919793089083292259942400
24650600248172987140112763715584000
1080657854354639453670407474439566400000
47701779391594966287470570490839978880000000

3.1415926535897932384626433832795028841971693993751058209749445923078164

Raku

<lang perl6># 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"</lang>

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

<lang rexx>/*REXX program uses the Almkvist─Giullera formula for 1 / (pi**2) [or pi ** -2]. */ numeric digits length( pi() ) + length(.); w= 102 say $( , 3) $( , w%2) $('power', 5) $( , w) say $('N', 3) $('integer term', w%2) $('of 10', 5) $('Nth term', w) say $( , 3, "─") $( , w%2, "─") $( , 5, "─") $( , 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>10  then do;  do 3*(n==11);  say ' .';  end;  iterate;  end
    say $(n, 3) right(a, w%2)  $(powX(z), 5)  right( lowE( format(a*z, 1, w-6, 2, 0)), w)
    end   /*n*/

say say 'The calculation of pi took ' n " iterations with " digits() ,

   " decimal digits precision using"   subword( sourceLine(1), 4, 3).

say numeric digits length( pi() ) - length(.); d= digits() - length(.); @= ' ↓↓↓ ' say center(@ 'calculated pi to ' d " fractional decimal digits (below) is "@, d+4, '─') say ' 'sqrt(1/s); say say ' 'pi(); @= ' ↑↑↑ ' say center(@ 'the true pi to ' d " fractional decimal digits (above) is" @, d+4, '─') 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") powX: procedure; parse arg p; return right( format( p, 1, 3, 2, 0), 3) + 0 /*──────────────────────────────────────────────────────────────────────────────────────*/ pi: pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286208,

     ||9986280348253421170679821480865132823066470938446095505822317253594081284811174503
     return pi

/*──────────────────────────────────────────────────────────────────────────────────────*/ 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</lang>

output when using the internal default input:

(Shown at two─thirds size.)

                                                        power
 N                     integer term                     of 10                                                Nth term
─── ─────────────────────────────────────────────────── ───── ──────────────────────────────────────────────────────────────────────────────────────────────────────
 0                                                   96  -3   9.600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-02
 1                                              5122560  -9   5.122560000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-03
 2                                         190722470400  -15  1.907224704000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-04
 3                                     7574824857600000  -21  7.574824857600000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-06
 4                                312546150372456000000  -27  3.125461503724560000000000000000000000000000000000000000000000000000000000000000000000000000000000e-07
 5                           13207874703225491420651520  -33  1.320787470322549142065152000000000000000000000000000000000000000000000000000000000000000000000000e-08
 6                       567273919793089083292259942400  -39  5.672739197930890832922599424000000000000000000000000000000000000000000000000000000000000000000000e-10
 7                  24650600248172987140112763715584000  -45  2.465060024817298714011276371558400000000000000000000000000000000000000000000000000000000000000000e-11
 8             1080657854354639453670407474439566400000  -51  1.080657854354639453670407474439566400000000000000000000000000000000000000000000000000000000000000e-12
 9         47701779391594966287470570490839978880000000  -57  4.770177939159496628747057049083997888000000000000000000000000000000000000000000000000000000000000e-14
10    2117262852373157207626265529989139651218848358400  -63  2.117262852373157207626265529989139651218848358400000000000000000000000000000000000000000000000000e-15
 .
 .
 .

The calculation of pi took  122  iterations with  163  decimal digits precision using the Almkvist─Giullera formula.

────────────────────────────────────────────── ↓↓↓  calculated pi to  160  fractional decimal digits (below) is  ↓↓↓ ───────────────────────────────────────────────
 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174503

 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174503
────────────────────────────────────────────── ↑↑↑  the  true  pi to  160  fractional decimal digits (above) is  ↑↑↑ ───────────────────────────────────────────────

Sidef

<lang ruby>func almkvist_giullera(n) {

   (32 * (14*n * (38*n + 9) + 9) * (6*n)!) / (3 * n!**6)

}

func almkvist_giullera_pi(prec = 70) {

   local Num!PREC = (4*(prec+1)).numify

   var sum = 0
   var target = -1

   for n in (0..Inf) {
       sum += (almkvist_giullera(n) / (10**(6*n + 3)))
       var curr = (sum**-.5).as_dec
       return target if (target == curr)
       target = curr
   }

}

say 'First 10 integer portions: '

10.of {|n|

   say "#{n} #{almkvist_giullera(n)}"

}

with(70) {|n|

   say "π to #{n} decimal places is:"
   say almkvist_giullera_pi(n)

}</lang>

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

Visual Basic .NET

Library: System.Numerics

Translation of: C#

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

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

Library: Wren-big

Library: Wren-fmt

<lang ecmascript>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))</lang>

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