为什么这段代码计算pi以及它基于什么原理？

Question

我在邮件的签名中偶然发现了这个片段:

int f[9814],b,c=9814,g,i;long a=1e4,d,e,h;main(){for(;b=c,c-=14;i=printf(
"%04d",e+d/a),e=d%a)while(g=--b*2)d=h*b+a*(i?f[b]:a/5),h=d/--g,f[b]=d%g;}

这计算pi(π)的输出:

3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354886230577456498035593634568174324112515076069479451096596094025228879710893145669136867228748940560101503308617928680920874760917824938589009714909675985261365549781893129784821682998948722658804857564014270477555132379641451523746234364542858444795265867821051141354735739523113427166102135969536231442952484937187110145765403590279934

我对此有一些问题:

• 它基于哪种数学原理？
• 这是什么样的迭代近似(哪个系列)？
• 为什么输出长2800个字符？

Answer 1

这似乎是基于Dik T. Winter的160字节C程序来计算pi的前800位数,扩展到更高的精度.

他的原始代码:

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);}

因此可以更容易地重写:

#include <stdio.h>

int main() {
    int r[2800 + 1];
    int i, k;
    int b, d;
    int c = 0;

    for (i = 0; i < 2800; i++) {
        r[i] = 2000;
    }

    for (k = 2800; k > 0; k -= 14) {
        d = 0;

        i = k;
        for (;;) {
            d += r[i] * 10000;
            b = 2 * i - 1;

            r[i] = d % b;
            d /= b;
            i--;
            if (i == 0) break;
            d *= i;
        }
        printf("%.4d", c + d / 10000);
        c = d % 10000;
    }

    return 0;
}

在k循环期间,这是计算以下数字:

由Beeler等人提出.(1972)是pi的近似值:

有关详细信息,请参见此处

• https://crypto.stanford.edu/pbc/notes/pi/code.html

Answer 2

我会回答最后一个问题:"为什么输出长2800个字符？"

在代码的初始化行中,我们看到了

c=9814

然后在循环条件中for我们看到循环的表达式

c-=14

因此,循环从c每个循环中减去14 ,从9814开始,并在c到达零时停止.但我们看到的 9814 / 14是701.由于该终止条件在第一个循环开始之前执行,因此循环执行700次数.

for循环的最后一部分,通常会增加循环变量,而不是执行

printf("%04d",e+d/a)

因此4,在每个循环中打印整数的字符.这些4字符,打印700时间,导致2800总共打印字符.