Calculate large factorials on a HP48 (SysRPL)

[Christoph Giesselink]

It's nice to see that some of you want learn or want to have a look inside SysRPL. One idea was using SysRPL for speed up calculation speed, here's another approach using SysRPL. This is an example from my fund I wrote years ago.

What is 253! = 5.17346099264E499

And what is 254! = ! Error: Overflow

The easy way to solve this, instead of multiplying the numbers, simply add the logarithm to the base 10 of each number. This is the UserRPL code doing this:

Code:


%%HP: T(3)A(R)F(.);
\<< 0 1 ROT
  FOR y y LOG +
  NEXT DUP FP ALOG
\->STR "E" + SWAP IP
\->STR +
\>>

But calculate 10! with this program and you see the disadvantage of this method: inaccuracy!

A solution for this is using the internal data type "long real". This SysRPL program using the same algorithm like before. Because the logarithm to the base 10 don't exist for long real I use the natural logarithm function %%LN. Here's the source code for compiling with HPTOOLS:

Code:


TITLE FAC, calculate factorial of a integer number
* integer argument : Real Number

ASSEMBLE
  NIBASC /HPHP48-E/

RPL
*
* start of secondary
*
::
  CK1NoBlame
  CK&DISPATCH1                                        ( test Stack1 on real )
  real
  ::
    DUP %0< IT :: #305 DO#EXIT ;                     ( < 0, Infinite Result )

    COERCE                                                 ( integer number )
    %%10 %%LN                                      ( pre calculate %%LN[10] )
    %%0                                                        ( for result )
    ROT                                                           ( counter )
    #1+_ONE_DO (DO)
      INDEX@ UNCOERCE%%                 ( convert loop counter to long real )
      %%LN 3PICK %%/                         ( %%LOG, not %%LN for accuracy )
      %%+                                                       ( to result )
    LOOP
    DUP %%FLOOR                                            ( fetch exponent )
    DUP4UNROLL %%-                                        ( fractional part )
    %%* %%EXP                                                      ( %%ALOG )
    2%%>% a%>$                                          ( convert to string )
    CHR_E >T$                                                  ( append "E" )
    SWAP a%>$ &$                                          ( append exponent )
  ;
;

So back to the entire question, what is 254! = "1.31405909214E502".

For single step execution of SysRPL programs on the HP48 I prefer Jazz.

Hope you enjoy,

Christoph