21. 计算 PI 的小数后 10000 位
根据下面的 J.Marchin 的圆周率计算公式：
其第 n 项为：
如果需要计算到 pi 的小数点后第 L 位，则至少需要计算至第 n 项，其中 n 为：
以上的中括号 [] 均为高斯函数，即取整函数，等于不大于该数值的最大整数。
下面，请利用以上公式，将 pi 推算到小数点后 10000 位，并估算计算所用毫秒数。
提示：
由于小数点后 10000 位的精度是各种类型的数据都难以满足的，因此可以考虑使用大实数计算来处理，为此需要先准备大整数计算的脚本 BigIntegerCalc.splx：
| A | B | C | D | |
|---|---|---|---|---|
| 1 | return case(calc,“+”:func(A4,[a,b]),“-”:func(A8,[a,b]),“*”:func(A12,[a,b]),“/”:func(A16,[a,b])) | |||
| 2 | 10000000000L | return A1 | ||
| 3 | /add | |||
| 4 | func | =A4(1).rvs() | =A4(2).rvs() | >c=0 |
| 5 | =B4.((sum=~+C4(#)+c,if(sum>A2,(c=1,sum-A2),(c=0,sum)))) | |||
| 6 | >if(c>0,B5=B5|long(c)) | return B5.rvs() | ||
| 7 | /minus | |||
| 8 | func | =A8(1).rvs() | =A8(2).rvs() | >c=0 |
| 9 | =B8.((differ=~-C8(#)-c,if(differ<0,(c=1,differ+A2),(c=0,differ)))) | |||
| 10 | return B9.rvs() | |||
| 11 | /multiply | |||
| 12 | func | =A12(1).rvs() | =A12(2) | >c=0 |
| 13 | =B4.((prodct=~+C12(#)+c,c=long(product/A2),long(product%A2))) | |||
| 14 | >if(c>0,B13=B13|long(c)) | return B13.rvs() | ||
| 15 | /divide | |||
| 16 | func | =A16(1) | =A16(2) | >r=0 |
| 17 | =B16.((dividend=~+r,(r=(dividend%C16)*A2, dividend\C16))) | |||
| 18 | return B17 | |||
其中，calc 为参数输入运算符号 +-*/，a 和 b 为参与运算的大实数序列或普通整数：
为了易于理解，A3,A7 等格子中记录了注释信息，以 / 开头的字符不做任何处理，仅用于记录说明文字。A1 中根据不同的符号，调用不同的子程序执行计算，调用时用序列参数，这种使用方法会使得调用时，参数都填入子程序的起始格中。计算时，每 10 位为一个“数字”，而后将每 10 位视为一个新的位数值，形成一个 10 位的 1010进制的数字，而计算时的处理类似于普通四则运算。为了区分各个子程序的功能，在每个子程序之前添加了注释信息来说明。
以乘法为例，由最低位“数字”开始计算，为此在计算前在 B12 和 C12 中将各个“数字”位倒转。用变量 c 记录每位“数字”计算时的进位，计算出当前“数字”位的乘积 product 后，用整数除法计算出进位，余数作为本位结果。每一位循环计算完毕后，需判断是否仍然有进位，若有需要在结果中添加一个高位。
特别的，由于这里是在处理大整数计算，因此在计算除法时，要求除数不超过进制单位，而结果仅保存整数部分，除不尽的余数舍弃。
参考答案：
解答：
使用提示中的 BigIntegerCalc.splx 用 call 函数辅助计算。由于公式中的第 n 项和第 n+1 项中，前半部分相差 25 倍，后半部分相差 239*239 倍，因此计算首项时同样把除数调整为 25 和 239*239。
| A | B | C | |
|---|---|---|---|
| 1 | =now() | 10000 | =int(B1/10)+3 |
| 2 | =C1*[0] | =C1*[0] | >B2(1)=16*5 |
| 3 | BigIntegerCalc.splx | =C1*[0] | >B3(1)=4*239 |
| 4 | for int(C1/(2*lg(5))+1) | ||
| 5 | >B2=call(A3,“/”,B2,25) | ||
| 6 | >B3=call(A3,“/”,B3,239*239) | ||
| 7 | =call(A3,“-”,B2,B3) | ||
| 8 | =call(A3,“/”,B7,2*#A4-1) | ||
| 9 | if #A4%2==1 | >A2=call(A3,“+”,A2,B8) | |
| 10 | else | >A2=call(A3,“-”,A2,B8) | |
| 11 | =string(A2(1))+“.”+A2.to(2,).(string(~,“0000000000”)).concat() | =interval@ms(A1,now()) |
C1 中准备足够的“数字”位来存储结果，其中每位“数字”为 10 位整数，用指定位数的大整数计算来替代实数计算。B2 和 B3 中存储每次循环时，公式中所需的计算参数，在 A4 中执行所需的循环次数。执行后在 A2 中能够查看按“数字”位存储的结果序列，A11 中拼接出最终结果，需要注意的是，如果某个“数字”不满 10 位需要在前方补 0。
在本题中，需要频繁调用脚本 BigIntegerCalc.splx，为此也可以将它登记为常规函数，如将 A3 中代码修改为 >register(“bigintcalc”,“BigIntegerCalc.splx”)，此后调用脚本时也需要改为直接使用登记的函数，以 B5 中代码为例，需要修改为 >B2=bigintcalc(B2,25)，其执行效果是一致的。登记过后，这个函数将对所有脚本都有效。
英文版