[Bug-apl] 1000 Digits of Pi

[Mike Duvos]

Here's a nice little floating point intensive benchmark that should run efficiently under GNU APL.

The function ME uses modular exponentiation to calculate N|16*M for positive integer M and N even if 16*M is too large to be exactly representable as a double float.

∇Z←M ME N;⎕IO;I;B;W
 ⎕IO←I←0
 B←∼⌽((1+⌊2⍟M)⍴2)⊤M
 Z←1
 W←N|16
 →B[0]/L1
 Z←N|Z×W
 L1:→((I←I+1)≥⍴B)/0
 W←N|W×W
 →B[I]/L1
 Z←N|Z×W
 →L1
∇

The function PI1000 calculates and returns the first 1000 digits of PI using the Bailey-Borwein-Plouffe formula.

∇Z←PI1000;⎕IO;T32;X;Y;I;C;J;K;T;S;K8
 ⎕IO←0◊T32←⌊2⋆32
 Y←125⍴I←¯1
 L1:→((I←I+1)≥125)/L4
 K←S←0
 J←8ׯ1+I
 L2:→(J≤0)/L3
 S←S+((4×J ME 1+K8)÷1+K8)-((2×J ME 4+K8)÷4+K8)+((J ME 5+K8)÷5+K8)+((J ME 6+K8)÷6+K8←8×K)
 K←K+1◊J←J-1
 →L2
 L3:S←S+T←(16⋆J)×(4÷1+K8)-(2÷4+K8)+(÷5+K8)+÷6+K8←8×K
 K←K+1◊J←J-1
 →(1E¯17<|T)/L3
 Y[I]←⌊T32×1|S
 →L1
 L4:X←1000⍴I←¯1
 L5:→((I←I+1)≥⍴X)/L7
 X[I]←Y[0]◊Y[0]←0◊Y←10×Y
 L6:C←⌊Y÷T32◊Y←(1⌽C)+T32|Y
 →(L5,L6)[∨/0≠C]
 L7:Z←20 10⍴(5×⍳200)+¨⊂⍳5
 (∈Z)←'0123456789'[X[∈Z]]
∇

And of course, our familiar timing function.

∇TIME X;TS
 TS←⎕TS
 ⍎X
 (⍕(24 60 60 1000⊥¯4↑⎕TS-TS)÷1000),' Seconds.'
∇

[IBM APL2]

      TIME 'PI1000'
 31415 92653 58979 32384 62643 38327 95028 84197 16939 93751
 05820 97494 45923 07816 40628 62089 98628 03482 53421 17067
 98214 80865 13282 30664 70938 44609 55058 22317 25359 40812
 84811 17450 28410 27019 38521 10555 96446 22948 95493 03819
 64428 81097 56659 33446 12847 56482 33786 78316 52712 01909
 14564 85669 23460 34861 04543 26648 21339 36072 60249 14127
 37245 87006 60631 55881 74881 52092 09628 29254 09171 53643
 67892 59036 00113 30530 54882 04665 21384 14695 19415 11609
 43305 72703 65759 59195 30921 86117 38193 26117 93105 11854
 80744 62379 96274 95673 51885 75272 48912 27938 18301 19491
 29833 67336 24406 56643 08602 13949 46395 22473 71907 02179
 86094 37027 70539 21717 62931 76752 38467 48184 67669 40513
 20005 68127 14526 35608 27785 77134 27577 89609 17363 71787
 21468 44090 12249 53430 14654 95853 71050 79227 96892 58923
 54201 99561 12129 02196 08640 34418 15981 36297 74771 30996
 05187 07211 34999 99983 72978 04995 10597 31732 81609 63185
 95024 45945 53469 08302 64252 23082 53344 68503 52619 31188
 17101 00031 37838 75288 65875 33208 38142 06171 77669 14730
 35982 53490 42875 54687 31159 56286 38823 53787 59375 19577
 81857 78053 21712 26806 61300 19278 76611 19590 92164 20198
26.39 Seconds.

Not a bad time for my desktop Optiplex, and it even got the right answer.  Almost all of the CPU time should be in the ME function, which is called four times for every term as it adds up the series.

I've attached PI.atf as a file, in case pasting the functions dropped any characters.

Currently, this won't run correctly under GNU APL because of issues with Residue and Floor.

PI.atf
Description: Binary data