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

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()

(原始碼)稍微修改，也能夠計算其他的交錯級數。