Π day

[DavidM]

Back in March of 2004, Olivier De Smet posted this fairly compact RPL version of a pi digits program (see post #17):

Code:

%%HP: T(3)A(R)F(.);
\<< -105. CF \-> N
  \<< 280000 N ALOG * DUP 'R' STO 2 \-> y
    \<<
      DO y * y 1 + IQUOT 50 IQUOT DUP 'R' STO+ 2 'y' STO+
      UNTIL DUP 0 ==
      END DROP
    \>> 30336 N ALOG * DUP 'R' STO+ 2 \-> y
    \<<
      DO 9 * y * 6250 IQUOT y 1 + IQUOT DUP 'R' STO+ 2 'y' STO+
      UNTIL DUP 0 ==
      END DROP
    \>>
  \>> R
\>>

It uses exact integers for intermediate computations, thus providing digit counts in excess of the inherent approximate limitations. He attributed it as an RPL translation of a program originally posted by Katie Wasserman which was written for the HP 32SII. Katie, in turn, cited the algorithm as essentially being the Euler/Machin approach (Pi = 20 * arctan(1/7) + 8 * arctan(3/79)).

As soon as I saw Olivier's version some time last year, I wanted to try a SysRPL implementation based on his code that would avoid the inherent performance penalty of IQUOT. This seemed like a good time to try it.

This is essentially Olivier's program, but coded in SysRPL and providing a result rounded to the specified number of digits (Olivier's provides 6 additional digits). The code will compile on a 49g+/50g with the standard HP extable installed (use 256.06 MENU ASM). I've used some DCs and DEFINEs to simplify the code listing a bit.

Code:

!NO CODE
!ASM
DC Z0_      273AB
DC Z1_      273B6
DC Z2_      273C2
DC Z9_      27416
!RPL
DEFINE      z*                FPTR2 ^QMul
DEFINE      z+                FPTR2 ^QAdd
DEFINE      zInc1             Z1_ z+
DEFINE      zInc2             Z2_ z+
DEFINE      zIntDiv           FPTR2 ^ZDIVext DROP
DEFINE      zDupZero?         FPTR2 ^DupQIsZero?
DEFINE      zALOG             FPTR 001 00E1
::
   CK1NoBlame
   CK&DISPATCH1
   real ::
      ( convert digit count to BINT )
      COERCE
      
      ( "bad argument value" if digit count is 0 )
      DUP#0=case SETSIZEERR
      
      ( setup locals )
      FPTR2 ^#>Z zALOG                       ( mag - antilog of digit count )
      Z2_                                    ( num - initially 2 )
      OVER ZINT 280000 z* DUP4UNROLL         ( result - initially mag*280000 )
      {{ result num mag }}
      
      ( initial result already on stack from above )
      BEGIN
         num z*
         num zInc1 zIntDiv
         ZINT 50 zIntDiv
         DUP result z+ !result
         num zInc2 !num
         zDupZero?   
      UNTIL
      DROP
      
      ZINT 30336 mag z*
      DUP result z+ !result
      Z2_ !num
      
      BEGIN
         Z9_ num z* z*
         ZINT 6250 zIntDiv
         num zInc1 zIntDiv
         DUP result z+ !result
         num zInc2 !num
         zDupZero?   
      UNTIL
      DROP

      ( retrieve result, abandon LAMs )
      1GETABND
      
      ( round the result to the number of digits specified )
      ZINT 500000 z+
      ZINT 1000000 zIntDiv
   ;
;
@

Performance improved, as expected. Generating 1000 digits with Olivier's original RPL program takes about 231 seconds (you actually get 1006 digits, but the last 4 are not accurate). Generating 1006 digits using the SysRPL version takes about 111 seconds, and all resulting digits appear to be correct (at least they match Wolfram Alpha's output for 1006 digits).

The usual advantages of SysRPL apply here, but a significant part of the improvement is in bypassing the overhead of IQUOT and going straight to the internal exact integer division code.

I'm sure NewRPL on the 50g and Python on the Prime would perform much better, but I thought this was reasonable performance for the task at hand on the original 50g platform.