在技術社區(qū)搜pi跳出來多半是樹莓派、PI控制器、pi agent這類內容真要搜“計算pi小數(shù)點后10000位”反而會掉進一堆年代久遠的代碼片段里有的用C語言全篇宏定義有的只貼出幾千位就說“已算到一萬位”。我自己動手完整做了一遍之后最大的感受是這個題目非常適合當作“任意精度計算”的入門實踐它逼著你把算法收斂速度、中間截斷誤差、浮點數(shù)的精度天花板這些平時被框架掩蓋掉的問題全部面對一遍。文章后面會給出可直接運行的兩套Python實現(xiàn)一套用decimal模塊邏輯直觀一套用純整數(shù)運算速度更好并講清楚驗證、性能、踩坑三個環(huán)節(jié)。無論是把它當面試題、項目引子還是性能基準測試這篇文章都能讓你少走彎路。1. 10000位背后的真實難度一個看似簡單的編程題1.1 目標不只是“算出來”而是“算對這10000位”很多人一上來就寫while True: pi ...跑完把數(shù)字貼出來結果對前50位后面就亂了。這個問題的本質不是循環(huán)次數(shù)而是精度系統(tǒng)的搭建。我們把這個需求拆開看實際上是三個子目標得到至少10000個正確的十進制小數(shù)位而不是一個近似浮點數(shù)計算過程可以被驗證別人能復現(xiàn)你的結果耗時可控不至于讓一次計算變成等待兩小時的煎熬。我把這三個子目標寫進計劃之后才意識到這不是一個“寫個公式就完事”的題目。它橫跨了數(shù)值分析選公式、估計截斷誤差、編程語言的數(shù)值模型float和Decimal的區(qū)別、大整數(shù)運算當數(shù)字變成10^10000級別時普通類型都失效三個層面。從規(guī)模上看10000位小數(shù)是double可表示精度的600倍左右。double在大多數(shù)語言里只有53位二進制有效數(shù)字換算成十進制大約是15到17位。也就是說用原生浮點類型你連第18位都保證不了。這個是后面所有坑的總根源。1.2 浮點數(shù)的“精度天花板”到底在哪IEEE 754規(guī)定C/C的double、Java的double、Python的float都使用64位存儲其中1位符號、11位指數(shù)、52位尾數(shù)加上隱含位可以看作53位。53位二進制對應的十進制精度是log10(2^53)≈15.95所以通常說“double有16位有效數(shù)字”。用這樣的類型去算π就算你鍵盤敲冒煙屏幕上永遠只會顯示3.141592653589793再往后都是噪音。要突破這個天花板只有兩條路一是引入任意精度庫比如GMP、MPFR、Java的BigDecimal、Python的decimal二是在整數(shù)空間里做運算把小數(shù)部分放大到10的N次冪全程用整數(shù)加減乘除最后再把小數(shù)點插回去。后面的整數(shù)版本用的就是第二條路。2. 算法選型哪些公式能撐起一萬位圓周率公式在數(shù)學史上非常多但真正適合編程計算的就那么幾類。我篩選時首先放棄的不是“錯”的公式而是“收斂太慢導致物理意義上不可能”的公式。2.1 蒙特卡洛與萊布尼茨級數(shù)入門可以上萬位不行蒙特卡洛法往正方形里隨機撒點靠面積比估計π。隨機采樣的誤差收斂速度是O(1/√N)也就是說要把誤差壓到10^(-10)需要10^20次采樣。放到10000位需要的采樣次數(shù)是10^20000級別宇宙毀滅都算不完。萊布尼茨級數(shù)π4(1-1/31/5-1/7...)就更夸張了。它是一個交替調和級數(shù)誤差衰減速度是1/(2k1)。每算一項小數(shù)點后的有效位數(shù)增加約0.3位。想靠它算到10000位需要約3×10^10000項同樣不可能。這兩個例子說明一個關鍵判斷標準當你要沖擊極高精度時級數(shù)的項與精度的關系必須是對數(shù)級的或者至少是冪級數(shù)中收斂極快的否則就是死路。2.2 馬青公式中等精度繞不開的經(jīng)典馬青公式Machin formula是1706年發(fā)現(xiàn)的π 16 arctan(1/5) - 4 arctan(1/239)把arctan展開成泰勒級數(shù)arctan(x) x - x^3/3 x^5/5 - ...代入x1/5和x1/239之后每一項的大小分別按1/25和1/57121的比例衰減。1/25的log10是-1.39794也就是說arctan(1/5)的級數(shù)每迭代一項大約能多1.4位小數(shù)而arctan(1/239)的每項衰減是4.756位。要達到10000位精度考慮上截斷余量arctan(1/5)需要大約7200項arctan(1/239)需要大約2200項。這個計算量非常溫和現(xiàn)代CPU毫秒級就能跑完。這也是為什么馬青公式是“萬位級精度”最實用的選擇。2.3 再看一眼Chudnovsky精度更高但復雜度也更高Chudnovsky算法是1989年提出的公式長這樣π 426880 √10005 / Σ_{k0}^∞ ( (6k)! (13591409 545140134k) ) / ( (3k)! (k!)^3 (-262537412640768000)^k )它的優(yōu)點非常嚇人每一項貢獻約14.18位十進制有效數(shù)字。算10000位只需要約710項算1億位也就600多萬項。但代價是每一項都要做超大整數(shù)的階乘、乘方和除法還涉及高精度的平方根計算。要用好它通常需要配合二進制分割binary splitting技術代碼復雜度直接上一個臺階。對于10000位這個精度馬青公式和Chudnovsky差距并不大。我的建議是如果你把這次任務當作算法學習馬青公式足夠如果你打算以后沖擊百萬位、千萬位那直接學Chudnovsky更值。2.4 我最終選型先馬青再用整數(shù)優(yōu)化我的最終方案分成兩步先用馬青公式的Decimal版本把邏輯跑通驗證前幾百位正確再切換成整數(shù)運算版本把速度提上去。這樣的好處是兩個實現(xiàn)互為參照算出來的結果還可以交叉驗證一旦有一個出問題立刻能發(fā)現(xiàn)。3. 從公式到代碼兩種可落地的實現(xiàn)方案我用的語言是Python 3。先聲明一點Python內置的float完全不參與這次計算核心是decimal模塊和大整數(shù)。3.1 Decimal版本最容易讀懂的實現(xiàn)Python的decimal模塊提供了任意精度的十進制浮點數(shù)核心是把精度上下文getcontext().prec設成目標位數(shù)。下面是完整實現(xiàn)from decimal import Decimal, getcontext def arctan_inv_decimal(x, n): 計算 arctan(x) 的泰勒級數(shù)x 必須是 Decimal total Decimal(0) term x xx x * x sign 1 for k in range(1, 2 * n, 2): total sign * term / k term * xx sign -sign return total def calc_pi_decimal(ndigits10000): # 留出20位余量避免中間舍入污染最后一位 getcontext().prec ndigits 20 # 迭代次數(shù)粗略估算arctan(1/5) 需要約 ndigits/1.397 項 n int(ndigits / 1.3) 300 a arctan_inv_decimal(Decimal(1) / Decimal(5), n) b arctan_inv_decimal(Decimal(1) / Decimal(239), n) pi 16 * a - 4 * b return str(pi)[:ndigits 2] if __name__ __main__: print(calc_pi_decimal(10000))注意兩個關鍵點getcontext().prec必須在做除法之前設置。如果在默認精度28下先算Decimal(1) / Decimal(239)得到的是一個只有28位有效數(shù)字的數(shù)后續(xù)無論怎么加精度誤差已經(jīng)埋進去了。迭代次數(shù)n不用算得特別精確取大一點不虧最多多跑幾千次循環(huán)但取小了最后若干位就是錯的。3.2 整數(shù)運算版本更快、更可控Decimal版本容易理解但每次循環(huán)都做Decimal除法本質上是模擬十進制浮點運算開銷不低。更貼近底層、也更快的方式是把整個結果放大10^prec倍用純整數(shù)來算泰勒級數(shù)。思路是這樣的arctan(1/d)的第k項是(-1)^(k) / ((2k1) * d^(2k1))我先把分子固定為10^prec用一個整數(shù)term表示當前項放大后的值def arctan_int(den, prec): 計算 arctan(1/den) * 10^prec 的整數(shù)近似值 total 0 term 10 ** prec // den k 1 sign 1 den2 den * den while term: total sign * (term // k) term // den2 k 2 sign -sign return total def calc_pi_int(ndigits10000): prec ndigits 20 # 余量留大一點更穩(wěn) a arctan_int(5, prec) b arctan_int(239, prec) pi_int 16 * a - 4 * b return pi_int這里term // den2的作用是讓當前項從1/5^(2k-1)過渡到1/5^(2k1)每一步只需要一次大整數(shù)除法。整個過程中所有的數(shù)都是整數(shù)不存在浮點舍入誤差只來自每一次整除的向下取整。由于我留了20位余量向下取整帶來的損失會被控制在最后十幾位以內不會污染前10000位。3.3 輸出格式與運行效果整數(shù)版本算出來的pi_int是一個大約有10020位數(shù)字的大整數(shù)第一位是3后面跟著10019位小數(shù)部分。輸出時只需要把它轉成字符串然后在第一位后面插入小數(shù)點def pi_to_string(pi_int, ndigits): s str(pi_int) # 防止某些極端情況下整數(shù)位數(shù)不夠先補零 if len(s) ndigits 1: s s.zfill(ndigits 1) return s[0] . s[1:ndigits 1] pi_int calc_pi_int(10000) print(pi_to_string(pi_int, 10000))在我的筆記本上跑一遍前幾行輸出是3.14159265358979323846264338327950288419716939937510 58209749445923078164062862089986280348253421170679 ...第一眼看到這個結果我就知道整個流程跑通了。但“看到了π”和“確認這一萬位全對”是兩碼事我單獨把驗證環(huán)節(jié)拎出來說。4. 驗證結果算出來的10000位怎么保證沒錯很多人算出結果就結束了但如果你真要把這個結果用于基準測試、算法對比或者教學演示一定要做驗證。這里分享幾種我用下來覺得靠譜的方式。4.1 前綴對比前100位一眼定勝負π的前100位是公開常數(shù)隨手可查3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679我在代碼里固定存了一段前綴字符串算完后直接startswith檢查。這一步能過濾掉90%的明顯錯誤公式抄錯、泰勒展開符號錯、小數(shù)點位置錯基本都逃不過這雙火眼金睛。4.2 交叉驗證用兩個獨立實現(xiàn)互算我前面特意保留了兩套實現(xiàn)Decimal版和整數(shù)版它們各有各的舍入來源。用同一個馬青公式分別算10000位再把字符串做一次全量對比如果完全一致那基本可以判定正確。交叉驗證里有個容易被忽略的細節(jié)兩套實現(xiàn)要盡量獨立不要復制同一份代碼。我的Decimal版和整數(shù)版從數(shù)據(jù)結構、循環(huán)方式到誤差來源都不一樣交叉驗證才有意義。如果你只是改改變量名驗證就是自欺欺人。4.3 分段切片核對與哈希校驗前綴對比只能證明開頭對交叉驗證能證明兩套代碼一致但還不能證明“兩套代碼一起錯了”這種極端情況。為了徹底打消疑慮可以引入第三方結果。方法很簡單找一個與你的代碼完全無關的高精度計算工具比如gmpy2.const_pi()、mpmath的mp.dps10000; mpmath.pi或者從OEIS、可信的開源倉庫下載標準π文本文件然后在隨機位置分段切片做對比。我習慣的做法是抽查三處第1000位附近取第990到1010位第5000位附近取第4990到5010位第9990位附近取第9970到10000位。如果這三段都能對上那基本可以確認算到了第10000位。更進一步把整個10000位字符串做一次SHA256哈希與官方文本的哈希比對一旦對上連“中間某處錯一位”的可能也被排除。在Python里做這個只是幾行代碼的事。5. 性能實測與優(yōu)化路徑5.1 三個精度檔位的耗時對比我在自己的筆記本Intel i58GB內存Python 3.11上分別跑了1000位、10000位、100000位耗時量級大致如下目標位數(shù)Decimal版耗時整數(shù)版耗時1,000約0.05秒約0.02秒10,000約0.9秒約0.2秒100,000約60秒約7秒這個數(shù)據(jù)不是精確基準不同機器差異很大但量級關系是穩(wěn)定的Decimal版在10萬位時有明顯吃力感整數(shù)版快了近一個數(shù)量級卻也開始逼近秒級。5.2 性能瓶頸到底在哪里馬青公式的計算量由兩部分組成第一是級數(shù)項數(shù)前面算過10000位需要約72002200項100000位就需要約7200022000項項數(shù)和精度成正比。第二是每一項操作的大數(shù)規(guī)模。整數(shù)版里的term有10^prec量級也就是10萬位時需要處理一個十萬位的整數(shù)每做一次整除開銷跟大數(shù)的字節(jié)長度成正比。項數(shù)乘上每次操作的大數(shù)長度總復雜度大致是O(n^2)。這就是為什么1000位時感覺不到時間100000位時明顯卡頓。Python的大整數(shù)雖然有C語言底層優(yōu)化但O(n^2)的曲線擺在那里位數(shù)每翻10倍時間要翻約100倍。5.3 更進一步從O(n^2)往O(n log n)走如果只是算到10000位優(yōu)化空間已經(jīng)不大。但如果想繼續(xù)沖更高精度可以考慮三條路用gmpy2庫替換Python原生int和decimal。gmpy2.mpz底層是GMP大整數(shù)乘除法比Python原生快數(shù)倍到數(shù)十倍同樣的馬青公式代碼改成gmpy2之后10萬位能壓進1秒以內。換Chudnovsky公式加二進制分割。二進制分割能把階乘和級數(shù)求和變成分治形式總復雜度降到接近O(n log n)這是目前百萬位以上的主流做法。如果只是臨時驗證直接用gmpy2.const_pi(prec)它會調用MPFR把π算到任意精度一行代碼速度還快得離譜。不過對于10000位這個目標我的結論是沒必要為了性能引入復雜方案馬青整數(shù)運算已經(jīng)是性價比最高的組合。6. 踩坑記錄幾個容易讓結果悄悄出錯的地方最后這部分是這次實操里最值錢的經(jīng)驗。下面每個坑我都實際踩過或者看著它們讓結果悄悄出錯。6.1 Decimal精度設置在“操作之前”而不是“表達式之前”Python的decimal精度上下文是全局狀態(tài)不是表達式屬性。最容易犯的錯誤是getcontext().prec 28 x Decimal(1) / Decimal(239) # 這里的除法已經(jīng)在28位精度下算完了 getcontext().prec 10050 # 改晚了 pi 16 * x - ...你以為后面把精度調到10050x已經(jīng)是一個28位精度的數(shù)后續(xù)計算結果的有效位數(shù)最多28位。正確做法是先把getcontext().prec設為目標精度再執(zhí)行任何除法、開方等會產(chǎn)生舍入的操作。6.2 整數(shù)版本的整除截斷誤差與余量設計整數(shù)版本每一步term // den2都會丟掉一點余數(shù)項數(shù)越多向下取整的累計誤差越大。我最初用prec ndigits去跑結果第9990位開始就和參考值對不上。后來把余量從10加到20尾端才穩(wěn)定下來。所以整數(shù)版本的余量不能省。prec ndigits 20是我實測夠用的值但如果你的機器環(huán)境不同建議算完后抽查尾部100位。6.3 迭代次數(shù)不足比超跑更可怕Decimal版本里我把迭代次數(shù)設成int(ndigits / 1.3) 300這個經(jīng)驗值夠用。但如果你圖省力寫nndigits也能過寫nndigits//2就會在很靠后的位置出現(xiàn)錯誤。重要的是理解估算邏輯arctan(1/5)每項增加約1.4位arctan(1/239)每項增加約4.8位按精度需求反推項數(shù)再留出10%左右的余量就不會踩坑。6.4 字符串輸出時的長度與補零整數(shù)版本里π乘以10^prec后整數(shù)部分有prec1位字符串長度足夠一般不用補零。但如果你把prec設成ndigits20str長度是ndigits21切片時[1:ndigits1]會正確拿到10000位。如果你在別的公式里遇到首項特別大或特別小的情況建議還是加一句zfill兜底。這種細節(jié)在現(xiàn)場跑數(shù)據(jù)時最磨人寧可多寫一行防御代碼。我個人做這類項目習慣先寫一個簡單的驗證函數(shù)把前綴、中間段、尾部段三處斷言寫進去每次改完代碼立刻全量自檢。這樣即使后面迭代了很多版本也不會在某個深夜把一段錯誤的結果當成“正確的一萬位”發(fā)布出去。這個計算任務表面上是玩數(shù)字實際上把數(shù)值穩(wěn)定性、大整數(shù)運算、算法復雜度分析全練了一遍。如果你也想動手試試建議從馬青公式的整數(shù)版本開始一步步把代碼寫出來再親手踩一遍精度余量的坑。等你能穩(wěn)定輸出并驗證10000位時后面再接觸Chudnovsky、二進制分割這些高階技巧會輕松很多。