21. Computing PI to 10000 decimal places

 

20. Handling big, multiline data using the cross-cellset cursor

According to the PI formula devised by J.Marchin:

imagepng

The nth term is:

imagepng

To compute pi to L decimal places, we need to compute the formula to at least the nth term. The value of n is:

imagepng

The brackets [] in the formula are all Gaussian functions, which are also called rounding functions and which are evaluated as the largest integers that are equal to or less than their actual values.

Use the above formula to compute pi to 10000 decimal places and estimate it takes how long (unit: millisecond) to finish the computation.

Tip:

It is hard to realize the 10000 decimal digits of precision using any type of data. In this case we can use a big integer to handle the computation. To do this we need to prepare a big integer computation script called, say, 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

In this script, parameter calc inputs operators +, -, * or /; parameters a and b are sequences of big integers or ordinary integers:

To help understand the script, A3, A7, A11 and A15 record comments, which are explanations and headed by slash / and won’t be taken care of. A1 calls different subprograms to execute the computation. A sequence parameter is used for the invocation by filling its members in the starting cell of the subprogram. During the computation every 10 digits are treated as one “number” and then considered as a new digit value to generate a 10-digit, radix of 1010 number. The computing process is similar to that for ordinary four arithmetic operations. In order to distinguish jobs of the subprograms, we add comment before each subprogram.

Take the multiplication as an example. The computation starts from the lowest digit “number”. To do this B12 and C12 reverse the order of “number” digits. Variable c is used to record the carry of handling each digit of “number” – compute the product of the current “number” digit and then use integer division to compute the carry while making remainder the radix result. After the loop is finished on a digit, check whether there is a carry and, if there is, add a higher place value to the result.

Specifically, when dealing with big integer computation with the division operation, the divisor must not exceed the radix, and the result should retain only the integer part by discarding the remainder.

Expected results:

SPL code:

Use call() function to invoke the script BigIntegerCalc.splx provided in the Tip part to assist the computation. Between the nth term and the (n+1) th term in the formula, the first part of the latter is 25 times greater than that of the former while there is 239*239 times difference between their second halves. So the divisors in computing the first term should be also adjusted to 25 and 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 defines enough number of “number” digit to store the result – where each “number” is a 10-digit integer – and replaces the real number computation with the computation of a big integer with specified number of digits. B2 and B3 store parameters used by the formula in each round of loop. A4 executes the specified number of loops and after that we can view the result sequence stored by “number” digits. A11 concatenates the final result. Note that if a “number” consists of less than 10 digits, corresponding number of 0s must be added before the existing number.

In handling this task, as repeated calls of *BigIntegerCalc.splx *are needed, we can register the script as a regular function. Change A3’s code to >register(“bigintcalc”,“BigIntegerCalc.splx”), for example, we should also use the registered function when calling the function later. In that case B5’s code will be >B2=bigintcalc(B2,25), and the result is the same after execution. A registered function is valid across all scripts.


22. The simple use of graph plotting
Contents and Exercise Data