Almkvist-Giullera formula for pi

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/π² 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:

Failed to parse (syntax error): {\displaystyle \frac{1}{pi^2} = \frac{2^5}{3}\sum_{0}^n\frac{6n}{n!^6(532n^2 + 126n + 9)\frac{1}{1000^2n+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:

(2^5) (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.

Julia

<lang julia>using Formatting

setprecision(BigFloat, 72)

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.1415926535897932384628340155181824844277116426383145153522491455078125
Computer π is     3.1415926535897932384628340155181824844277116426383145153522491455078125

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