Almkvist Berndt 1988 begins with an investigation of why the agm is such an efficient algorithm, and proves that it converges quadratically. This is an efficient method to calculate 



π


{\displaystyle \pi }

.

With the same notations used in Arithmetic-geometric mean, we can summarize the paper by writing:





π
=



4


a
g
m

(
1
,
1

/



2



)

2




1
−

∑

n
=
1


∞



2

n
+
1


(

a

n


2


−

g

n


2


)





{\displaystyle \pi ={\frac {4\;\mathrm {agm} (1,1/{\sqrt {2}})^{2}}{1-\sum \limits _{n=1}^{\infty }2^{n+1}(a_{n}^{2}-g_{n}^{2})}}}



This allows you to make the approximation, for any large   N:





π
≈



4


a

N


2




1
−

∑

k
=
1


N



2

k
+
1


(

a

k


2


−

g

k


2


)





{\displaystyle \pi \approx {\frac {4\;a_{N}^{2}}{1-\sum \limits _{k=1}^{N}2^{k+1}(a_{k}^{2}-g_{k}^{2})}}}



The purpose of this task is to demonstrate how to use this approximation in order to compute a large number of decimals of 



π


{\displaystyle \pi }

.

from math import sqrt
import times
 
import bignum
 
#---------------------------------------------------------------------------------------------------
 
func sqrt(value, guess: Int): Int =
  result = guess
  while true:
    let term = value div result
    if abs(term - result) <= 1:
      break
    result = (result + term) shr 1
 
#---------------------------------------------------------------------------------------------------
 
func isr(term, guess: Int): Int =
  var term = term
  result = guess
  let value = term * result
  while true:
    if abs(term - result) <= 1:
      break
    result = (result + term) shr 1
    term = value div result
 
#---------------------------------------------------------------------------------------------------
 
func calcAGM(lam, gm: Int; z: var Int; ep: Int): Int =
  var am: Int
  var lam = lam
  var gm = gm
  var n = 1
  while true:
    am = (lam + gm) shr 1
    gm = isr(lam, gm)
    let v = am - lam
    let zi = v * v * n
    if zi < ep:
      break
    dec z, zi
    inc n, n
    lam = am
  result = am
 
#---------------------------------------------------------------------------------------------------
 
func bip(exp: int; man = 1): Int {.inline.} = man * pow(10, culong(exp))
 
#---------------------------------------------------------------------------------------------------
 
func compress(str: string; size: int): string =
  if str.len <= 2 * size: str
  else: str[0..<size] & "..." & str[^size..^1]
 
#———————————————————————————————————————————————————————————————————————————————————————————————————
 
import os
import parseutils
import strutils
 
const DefaultDigits = 25_000
 
var d = DefaultDigits
if paramCount() > 0:
  if paramStr(1).parseInt(d) > 0:
    if d notin 1..999_999:
      d = DefaultDigits
 
let t0 = getTime()
 
let am = bip(d)
let gm = sqrt(bip(d + d - 1, 5), bip(d - 15, int(sqrt(0.5) * 1e15)))
var z = bip(d + d - 2, 25)
let agm = calcAGM(am, gm, z, bip(d + 1))
let pi = agm * agm * bip(d - 2) div z
 
var piStr = $pi
piStr.insert(".", 1)
 
let dt = (getTime() - t0).inMicroseconds.float
let timestr = if dt > 1_000_000:
                (dt / 1e6).formatFloat(ffDecimal, 2) & " s"
              elif dt > 1000:
                $(dt / 1e3).toInt & " ms"
              else:
                $dt.toInt & " µs"
echo "Computed ", d, " digits in ", timeStr
 
echo "π = ", compress(piStr, 20), "..."