Almkvist-Giullera formula for pi: Difference between revisions

Content added Content deleted

Inline

Revision as of 21:29, 12 October 2020

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.

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

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

REXX

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

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

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