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

[Werner]

For both accuracy and speed, however, you'd be better off dividing by LN(10) only once, after the loop?

Cheers, Werner

[walter b]

(02-18-2014 10:32 PM)Christoph Giesselink Wrote:  So back to the entire question, what is 254! = "1.31405909214E502".

Hope you enjoy,

Personally, I enjoy using the WP 34S for such problems: what is e.g. 5432! ? Returns: 3.022 553 598 4 E17 931. See the description of LNΓ on p. 102, IIRC.

d:-)

[Gerson W. Barbosa]

(02-19-2014 11:45 AM)walter b Wrote:  Personally, I enjoy using the WP 34S for such problems: what is e.g. 5432! ? Returns: 3.022 553 598 4 E17 931. See the description of LNΓ on p. 102, IIRC.

d:-)

What is the last non-zero digit of 5432! ? My good old HP-33C says it is 2, in about 30 seconds! :-)

(Program in message #80 here)

[Christoph Giesselink]

(02-19-2014 10:45 AM)Werner Wrote:  For both accuracy and speed, however, you'd be better off dividing by LN(10) only once, after the loop?

I agree with the speed reason, but when I wrote the program years ago I got better results adding the LOG instead of the LN numbers.

But your suggested modification is simple for own tests:

Code:

...
  #1+_ONE_DO (DO)
      INDEX@ UNCOERCE%%                 ( convert loop counter to long real )
      %%LN                                                 ( %%LN of number )
      %%+                                                       ( to result )
    LOOP
    OVER %%/                               ( convert result to sum of %%LOG )
    DUP %%FLOOR                                            ( fetch exponent )
...

Christoph

[Dieter]

(02-19-2014 11:45 AM)walter b Wrote:  Personally, I enjoy using the WP 34S for such problems: what is e.g. 5432! ? Returns: 3.022 553 598 4 E17 931. See the description of LNΓ on p. 102, IIRC.

Users without access to a ln Γ function may try one of the Stirling formulas instead. ;-)

BTW, a general caveat regarding methods based on lg n! or ln Γ: Since the exponent 17931 has five digits, at least five mantissa digits of the result are uncertain. On a WP34s in 16-digit SP mode the valid digits are exactly those eleven you posted. Yes, there is a DP option. ;-)

Dieter

[Dieter]

(02-18-2014 10:32 PM)Christoph Giesselink Wrote:  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.

So essentially you're calculating ln n! first and then n! is obtained from this.

There is an easier and much faster method for this: simply use a Stirling approximation. With only two terms, the largest error in log₁₀ n! is less than 1 unit in the 12th significant digit if n ≥ 70. Since n! can be calculated directly for n ≤ 253, we only have to consider n ≥ 254. If evaluated exactly (!), the absolute error of log₁₀ n! is less than 3,3 E-16, i.e. better than the available working precicion. So two terms are all we need:

\(ln n!  =  n · ln n + \frac{1}{2}ln(2 n \pi) - n + \frac{1}{12 n} - \frac{1}{360 n^3} + ...\)

The resulting accuracy is essentially limited by that of the internal 15-digit functions and the order of addition: add the smallest terms first.

Dieter