偶然看到《How to compute digits of ?》中有这么一段:
The following 160 character C program, written by Dik T. Winter at CWI, computes to 800 decimal digits. int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5; for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a, f[b]=d%--g,d/=g--,--b;d*=b);} /*出处: https://cs.uwaterloo.ca/~alopez-o/math-faq/node38.html Alex Lopez-Ortiz Mon Feb 23 16:26:48 EST 1998 注:查了几个国内的网站,看到最早一个转发日期在2001年*/
转C++代码
整理一下格式修改输出语句转成C++代码,在Dev-C++上通过编译,得到一个长度800的字符串:
1.#include <iostream> #include <iomanip> using namespace std; int main(void) { int a=10000,b,c=2800,d,e,f[2801],g; for(;b-c;)f[b++]=a/5; for(;d=0,g=c*2;c-=14,cout<<setw(4)<<setfill('0')<<e+d/a,e=d%a) for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b); return 0; } /* 31415926535897932384626433832795028841971693993751058209749445923078164062862089\ 98628034825342117067982148086513282306647093844609550582231725359408128481117450\ 28410270193852110555964462294895493038196442881097566593344612847564823378678316\ 52712019091456485669234603486104543266482133936072602491412737245870066063155881\ 74881520920962829254091715364367892590360011330530548820466521384146951941511609\ 43305727036575959195309218611738193261179310511854807446237996274956735188575272\ 48912279381830119491298336733624406566430860213949463952247371907021798609437027\ 70539217176293176752384674818467669405132000568127145263560827785771342757789609\ 17363717872146844090122495343014654958537105079227968925892354201995611212902196\ 08640344181598136297747713099605187072113499999983729780499510597317328160963185\ -------------------------------- Process exited after 0.5965 seconds with return value 0 请按任意键继续. . . */
验证准确性
到网上搜了2个pi字串,分别是小数点后2600位和1000位的,用它们验证一下这801位(加了小数点)的准确性:
>>> pi2602=\ '''3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339047802759009946576407895126946839835259570982582262052248940772671947826848260147699090264013639443745530506820349625245174939965143142980919065925093722169646151570985838741059788595977297549893016175392846813826868386894277415599185592524595395943104997252468084598727364469584865383673622262609912460805124388439045124413654976278079771569143599770012961608944169486855584840635342207222582848864815845602850601684273945226746767889525213852254995466672782398645659611635488623057745649803559363456817432411251507606947945109659609402522887971089314566913686722874894056010150330861792868092087476091782493858\ ''' >>> pi1002=\ '''3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989\ ''' >>> pi801=\ '''3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185\ ''' >>> pi1002==pi2602[:1002] True >>> pi801==pi1002[:801] True >>>
改写成Python代码
a,c=10000,2800 b=e=0 f=[0]*2801 while b!=c: f[b]=a//5 b+=1 while c!=0: d=0 g=c*2 b=c while 1: d+=f[b]*a g-=1 f[b]=d%g d//=g g-=1 b-=1 if b!=0: d*=b else: break c-=14 print('%.4d'%(e+d//a),end='') e=d%a
运行结果:
压缩代码行数、字符数
a,c=1e4,2800;b=d=e=0;f=[0]*2801 while b!=c:f[b]=a//5;b+=1 while c!=0: b,g=c,c*2 while 1: d+=f[b]*a;g-=1;f[b]=d%g;d//=g;g-=1;b-=1 if b==0:break d*=b c-=14;print('%.4d'%(e+d//a),end='');e=d%a
Python不需要C的语句结束字符“分号;”,但可以用分号把多行代码连成一行;只是没法减掉必要的缩进,python语法规定使然。另外没压缩前的: if b!=0:d*=b \ else:break 改写成: if b==0:break \ d*=b 也省下5个字符。
统计行数和字符数
压缩后的代码去掉空格不算约177个字符,而原先的C代码没算进“#include <stdio.h>”,如上图的统计压缩后的代码只比C代码输(多)3个字符;输就输在Python没有自增++自减--操作符。另外规范的C代码,main()函数应该有整型声明及整数返回值: int main(){...; return 0;}
改进代码成函数Pi799()
def Pi799(): a,c=1e4,2800;b=d=e=0;f=[0]*2801;s='' while b!=c:f[b]=a//5;b+=1 while c!=0: b,g=c,c*2 while 1: d+=f[b]*a;g-=1;f[b]=d%g;d//=g;g-=1;b-=1 if b==0:break d*=b c-=14;s+=f'{e+d//a:04}';e=d%a return '3.'+s[1:] pi=Pi799() print(pi)
代码多引入一个字符串变量s到最后输出,改变了原代码像慢动作一样的输出;还有未修改前的 '%.4d'%(e+d/a) 可以改写成: str(e+d//a).zfill(4) 如此更好理解但多用了几个字符;最后决定改写成: f'{e+d//a:04}' (版本要求 python 3.6+),字符数刚好一样多,但可读性更好 ^_^
数学原理
数列:f(1)=1,f(n+1)=f(n)*n/(2*n+1)
Pi/2=f(1)+f(2)+f(3)+.....+f(n) n->无穷大
或:Pi = 2 + 2/3 + 2*(1*2)/(3*5) + 2*(1*2*3)/(3*5*7) + 2*(1*2*3*4)/(3*5*7*9)+...
还需要用到使用数组精确计算除法的算法,网上查询到比较可信的推理:
http://www.voidcn.com/article/p-sulbahsi-wq.html
如果就仅在float精度范围内,用上述数列50次循环可计算到pi小数点后14位的精度。
pi,a=0,lambda n:1 if n==1 else a(n-1)*(n-1)/(2*n-1) for i in range(1,51):pi+=2*a(i);print(pi) ''' =============================== RESTART: D:/pi.py ============================== 2 2.6666666666666665 2.933333333333333 3.0476190476190474 3.098412698412698 3.121500721500721 3.132156732156732 3.1371295371295367 3.1394696806461506 3.140578169680336 3.141106021601377 3.1413584725201353 3.1414796489611394 3.1415379931734746 3.1415661593449467 3.1415797881375944 3.1415863960370602 3.141589605588229 3.1415911669915006 3.1415919276751456 3.1415922987403384 3.141592479958223 3.1415925685536337 3.1415926119088344 3.1415926331440347 3.1415926435534467 3.141592648659951 3.14159265116678 3.141592652398205 3.1415926530034817 3.1415926533011587 3.141592653447635 3.141592653519746 3.1415926535552634 3.141592653572765 3.1415926535813923 3.141592653585647 3.141592653587746 3.141592653588782 3.141592653589293 3.1415926535895458 3.1415926535896705 3.1415926535897323 3.141592653589763 3.141592653589778 3.1415926535897856 3.141592653589789 3.141592653589791 3.141592653589792 3.1415926535897922 >>> '''