HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
Calculation of factorials using logarithms (HP41 MCODE)
Posted by Frido Bohn on 26 May 2012, 9:15 a.m.
This is a routine in MCODE I developed to overcome the situation that an HP41C can only calculate factorials of up to 69 as this is <10100. Jean-Marc Baillard created also some great routines dealing with this "problem".
The idea of this solution is to calculate the sum of logarithms from 1 to n. Thus the result of "LOGSUM" is the decimal logarithm of the factorial of n. The program takes advantage of 13-digit precision for the calculation of the natural logarithm (ln) and its sum. Eventually ln is converted to the decimal logarithm (log). The intermediate results are stored in scratch registers using MCODE routines of the operating system designed to handle 13-digit precision ([STCR] and [RCSCR]). The loop counter is also kept in a scratch register.
As the result is not that vivid, a little FOCAL program is presented which converts the logarithm into an approximation of the natural number using the alpha registers. For example LOGUM 10 results in 6.5598 which is log(10!), and XFACT transforms it to 3.6288E6 (10! = 3628800).
This program runs well on every flavor of HP41. However, as the calculations are based on an iterative process the larger the number is the longer it will last. Therefore, HP41CL in TURBO mode and emulators on fast processors have a noticeable advantage.
The execution times of XFACT at normal speed and 50x turbo and results are given below. As a reference, the truncated outputs of Wolfram Alpha are shown.
(XFACT was compiled and transferred to HEPRAM. The HP-IL module was attached for printing.)
XFACT 1x Speed 50x Speed Result Wolfram Alpha (n) (sec) (sec) (output) (truncated) 10 3.82 0.75 3.6288E6 36288 x10^6 20 5.81 0.76 2.4329E18 2.4329…x10^18 30 7.93 0.82 2.6525E32 2.6525…x10^32 50 12.13 0.97 3.0414E64 3.0414…x10^64 100 22.47 1.28 9.3326E157 9.3326…x10^157 200 44.94 1.68 7.8866E374 7.8865…x10^374 500 109.81 3.01 1.2201E1134 1.2201…x10^1134 1000 218.21 5.31 4.0239E2567 4.0238…x10^2567 5000 1113.30 24.25 4.2285E16325 4.2285…x10^16325 10000 n.d. 47.55 2.8461E35659 2.8462…x10^35659 100000 n.d. 473.14 2.8029E456573 2.8242…x10^456573A linear regression of the time elapsed and the number entered shows a very good correlation coefficient RR >0.9999. The ratio of the slopes of the regression line resulted to be pretty exactly 47x times larger at 50x TURBO if compared to normal speed. The small difference from 50 may be attributed to the overhead (key input, FOCAL code, display and HP-IL output).
08D "M" 015 "U" 013 "S" 007 "G" 00F "O" 00C "L" 'Name of function 2A0 SETDEC 'Start (to be linked to FAT) 1A0 A=B=C=0 'Clears CPU registers A, B and C 089 ?NC XQ 'so the scratch registers Q and + get cleared to 064 ->1922 '[STSCR] 0F8 READ 3(X) 'Reads the number to be factorialized 228 WRIT 8(P) 'into register P to act as loop counter (Start of loop) 10E A=C ALL 'and into register A where MS and S&X is kept 0FA C<>B M 'while register B gets the mantissa 121 ?NC XQ 'So the 13-digit function of ln is fed correctly 06C ->1B48 '[LN13] 0D1 ?NC XQ 'Recall scratch registers Q and + 064 ->1934 '[RCSCR] 031 ?NC XQ 'add ln to the scratch using 13-digit precision 060 ->180C '[AD2-13] 089 ?NC XQ 'and store the result back to the scratch 064 ->1922 '[STSCR] 238 READ 8(P) 'Read the loop counter into register C 10E A=C ALL 'transfer it to register A 04E C=0 ALL 'clear C 130 LDI S&X 001 CON: 1 'Load 1 into S&X of C and right shift twice to get a 1 23C RCR 2 'make CHS transformation to -1 2BE C=-C-1 MS 'As C is wrapped around the carry is set and the entry 01D ?C XQ 'point for a 10-digit addition can be accessed properly. 061 ->1807 '[AD2-10] 3CC KEY ? 'If a key is pressed the program is interrupted and shows 05F JC +&0B 'the number of loops missing for completion 2FA ?C#0 M 'Otherwise check if the counter is zero, if not go back to 34F JC -&17 'loop start 2B5 ?NC XQ 'Perform ln(10) 068 ->1AAD '[LNC10*] 239 ?NC XQ 'Make 1/ln(10) of it 060 ->188E '[1/X13] 0D1 ?NC XQ 'Call back the scratch content - the sum of ln 064 ->1934 '[RCSCR] 149 ?NC XQ 'and multiply it with 1/ln(10) to obtain the decimal log 060 ->1852 '[MP2-13] 0E8 WRIT 3(X) 'Put the result into stack register X 3E0 RTN 'and finish.
01*LBL "XFACT" 02 LOGSUM 'Calculate the logarithmic sum from 1 to n 03 INT 'Take the integer. It is the exponential part. 04 LASTX 05 FRC 'Take the fraction, the antilogarithm of which 06 10^X 'is the mantissa. 07 FIX 4 08 CLA 09 ARCL X 10 "|-E" 'This means APPEND "E". 11 FIX 0 12 ARCL Y 13 AVIEW 'Output as a.aaaaEbbbbb in the alpha registers 14 FIX 4 15 END
[ Return to the Message Index ]
Go back to the main exhibit hall