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