Python 算圓週率 這次算到無窮位

Posted by tjwei on 星期日, 5月 04, 2008 with No comments
python 的 generator 很好用,理論上來說,可以表示無窮小數。而無窮小數的四則運算,可以用前面有限位來計算到一定的位數,所以,如果用 generator 生成的十進位展開,配合 operator +-*/ 等等,我們在 python 上可以弄出一個無窮精確的 real number class。
這樣,其實有點用有限的電腦來存放無窮位的資料的感覺。
這樣,我們可以很優雅的算出圓週率近似值到無窮位。
當然,速度一定很慢。反正太匆忙本來就不太優雅。
之前有看過利用類似概念來算圓周率的 python 程式,實在不錯。不過那個程式太長了,而且年代久遠,沒有用 operator 和 generator。
另一個好像叫做 accurate-math 的東西,也是可以做無窮位計算,不過呢,他的方法比較類似先固定算前一千位,然後等到需要一千零一位時,從頭開始算前兩千位。
像這種沒有實用價值的東西,用這麼醜的方法是不行的。
於是,我嘗試了一下寫一般性的 real number class ,才發現很麻煩(主要是 0.999...這類的問題)。而且, python 有 recursion 的 limit ,雖然可以改大一點,但再怎麼大也不是無窮位,除非用 stackless python。當然,可以避開不用 recursion,但程式結構就沒有那麼漂亮了。
所以就放棄,改成以算圓週率為目標。這樣,只需要簡單的除法、減法以及交錯級數就行。
結果需要三十幾行,跟之前只算固定位的三行比起來,差太遠了,不夠漂亮。
所以算是失敗,不過還是可以看看就是了。雖然超級慢,但是只要記憶體夠,可以一直跑到無窮位的圓週率。(這個程式碼的風格很糟,原本希望能夠用短短幾句話,很清楚的描述四則運算和交錯級數的計算,可惜失敗)

from itertools import *
import sys
# a/b 的無窮位小數展開
def
ldiv(a,b,p=0):
f=lambda x:x.append((x[0]%b)*10) or x.pop(0)/b
return iter(f.__get__([a*10**p]), -1)
# 比較上界ub 下界lb,然後將相同的位數傳回(p 是用來補前面的 0)
def digits(lb,ub,p):
u,l="%d"%ub, "%d"%lb
n,m=len(l),len(u)
common=takewhile(lambda i:l[i]==u[i], range(n))
return [] if m-n else chain(repeat(0, p-n), (int(l[i]) for i in common))
# 交錯級數計算(只算一個位數,要代入後面的模版)
def altseries(ub,lb, num, prec, a):
if num:
s, d=sum(a(i).next() for i in range(num-1)), a(num-1).next()
ub, lb= ub+s+max(d,0), lb+s+min(d,0)
for i in count(num):
d=a(i, prec).next()
if not d:
return ub, lb, i
ub,lb = (lb+d, lb) if d>0 else (ub, ub+d)
# 減法(只算一個位數,要代入後面的模版)
def subop(ub, lb, num, prec, a,b):
ub=lb+a.next()
return ub, ub-b.next(),1
# 減法以及交錯級數計算的模版
# 真正的上下界是 ub+num 和 lb-num
def opgen(ub,lb, num, precs, op, *a):
for prec in count():
ub,lb,num=op(ub*10,lb*10,num, prec, *a)
for x in digits(lb-num/2-1, ub+num/2+1, prec-precs):
yield x
precs+=1
ub,lb=[x%(10**(prec-precs)) for x in (ub,lb)]
# Pi 的交錯級數(其實是 pi(n) 是數列的第 n 項)
s=lambda k,p:opgen(0,0,0,-1,subop, ldiv(16, k*5**k, p), ldiv(4, k*239**k, p))
def pi(n, p=0, d={}):
return d.setdefault(n,imap(lambda x:-x, s(2*n+1,p)) if n%2 else s(2*n+1,p))
# 開始利用 generator 展開小數
for n in opgen(0,0,0,-1,altseries, pi):
sys.stdout.write("%d"%n)
sys.stdout.flush()
(原始碼)稍微修改,也能夠計算其他的交錯級數。
Categories: