アルキメデスの方法による円周率の近似値

直径 1 の円の外接正 n 角形の周長を Pn，内接正 n 角形の周長を pn とすれば，

n = 3 * 2 ^ k, k = 0, 1, 2, ...

P2n = 2 * Pn * pn

p2n = sqrt(P2n * pn)

ということで，以下のプログラム

import scipy as sp

def Pi_Archimedes():
    print("{0:>2s} {1:>9s} {2:>17s} {3:>17s}".format("k", "n", "Pn", "pn"))
    k = 0 # 正三角形から始める
    Pn = 3 * sp.sqrt(3)
    pn = Pn / 2
    while True:
        n = 3 * 2 ** k
        print("{0:2d} {1:9d} {2:.15f} {3:.15f}".format(k, n, Pn, pn))
        if Pn == pn:
            return Pn
        k += 1
        P2n = 2 * Pn * pn / (Pn + pn)
        p2n = sp.sqrt(P2n * pn)
        Pn, pn = P2n, p2n

Pi_Archimedes()

 k         n                Pn                pn
 0         3 5.196152422706632 2.598076211353316
 1         6 3.464101615137754 3.000000000000000
 2        12 3.215390309173473 3.105828541230249
 3        24 3.159659942097500 3.132628613281238
 4        48 3.146086215131435 3.139350203046867
 5        96 3.142714599645368 3.141031950890509
 6       192 3.141873049979824 3.141452472285462
 7       384 3.141662747056849 3.141557607911857
 8       768 3.141610176604690 3.141583892148318
 9      1536 3.141597034321526 3.141590463228050
10      3072 3.141593748771351 3.141592105999271
11      6144 3.141592927385096 3.141592516692157
12     12288 3.141592722038614 3.141592619365384
13     24576 3.141592670701998 3.141592645033691
14     49152 3.141592657867844 3.141592651450767
15     98304 3.141592654659306 3.141592653055036
16    196608 3.141592653857171 3.141592653456104
17    393216 3.141592653656637 3.141592653556370
18    786432 3.141592653606503 3.141592653581437
19   1572864 3.141592653593970 3.141592653587703
20   3145728 3.141592653590837 3.141592653589270
21   6291456 3.141592653590053 3.141592653589661
22  12582912 3.141592653589857 3.141592653589759
23  25165824 3.141592653589808 3.141592653589784
24  50331648 3.141592653589796 3.141592653589790
25 100663296 3.141592653589793 3.141592653589791
26 201326592 3.141592653589792 3.141592653589792
27 402653184 3.141592653589792 3.141592653589792

3.141592653589792

正 4 億多角形ですと！

 

R によるプログラムもほとんど同じで，以下のように

Pi.Archimedes = function() {
    Pn = 3 * sqrt(3)
    pn = Pn / 2
    while (TRUE) {
        if (Pn == pn) {
            return(Pn)
        }
        tmp = 2 * Pn * pn /  (Pn + pn)
        pn = sqrt(tmp * pn)
        Pn = tmp
    }
}

> options(digits=16)
> Pi.Archimedes()
[1] 3.141592653589792

追記

Julia では

function piarchimedes()
    Pn = 3sqrt(3)
    pn = Pn / 2
    while true
        Pn == pn && return Pn
        tmp = 2Pn * pn /  (Pn + pn)
        pn = sqrt(tmp * pn)
        Pn = tmp
    end
end

piarchimedes() # 3.141592653589792