The Museum of HP Calculators

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.

Calculating Factorials in HP41 MCODE

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.4329x10^18
30	7.93	0.82	2.6525E32	2.6525x10^32
50	12.13	0.97	3.0414E64	3.0414x10^64
100	22.47	1.28	9.3326E157	9.3326x10^157
200	44.94	1.68	7.8866E374	7.8865x10^374
500	109.81	3.01	1.2201E1134	1.2201x10^1134
1000	218.21	5.31	4.0239E2567	4.0238x10^2567
5000	1113.30	24.25	4.2285E16325	4.2285x10^16325
10000	n.d.	47.55	2.8461E35659	2.8462x10^35659
100000	n.d.	473.14	2.8029E456573	2.8242x10^456573
A 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).
The precision deteriorates at great numbers. However, the difference is about 1% at 100000!, thus I would rather prefer the absolute difference in dollars if someone asks me.

The program accepts any non-negative integer >0. Theoretically, the sum of logs must be only <10100. Any key can be pressed to abort the program if the user feels it lasts too long (LOGSUM would then render the iterations left to be done).

Listing of MCODE program "LOGSUM"
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.

This little FOCAL program renders a legible result in decimal representation, which however cannot be used for further calculations:

Listing of FOCAL program "XFACT"
02 LOGSUM	 'Calculate the logarithmic sum from 1 to n
03 INT	         'Take the integer. It is the exponential part.
05 FRC	         'Take the fraction, the antilogarithm of which
06 10^X	         'is the mantissa.
07 FIX 4
08 CLA
10 "|-E"	 'This means APPEND "E".
11 FIX 0
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