|HP41C: Factorial (kind of) in MCODE|
Message #1 Posted by Frido Bohn on 24 May 2012, 8:34 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 <10^100.
It may serve as a neat example how to handle MCODE in conjunction with the HP41CL in TURBO mode or with emulators on fast processors. The calculation of big numbers in an original HP41 may last a little too long. As far as I know, this "problem" has not been considered yet on a HP41 system. If so, I would like to get informed as I would be curious to learn about it.
With a little trick the numeric ability of the factorial problem can be expanded enormously by the use of logarithms. We can overcome the limitations of natural numeric expression as the antilogarithm of a sum of logarithms from 1 to n corresponds to the factorial of n.
The MCODE routine named LOGSUM takes a number (integer) from stack register X and calculates its natural logarithm starting with the number itself and looping down to 1. The result is then converted to a base-10 logarithm. I experimented with a sum of the base-10 logarithm but found it less reliable than the natural logarithm. Interestingly, the results of the latter were more precise in the HP41.
The main challenge of the MCODE program was to store intermediate results. Therefore, I used the status registers P for the storage of the intermediate log sums, and Q for the loop counter. The use of the CPU registers would increase the performance, however, they are used during the execution of the internal functions for the calculation of the natural logarithm (LN10, entry point 1B45) and the addition (AD2-10, entry point 1807) so intermediate results would be disturbed. In addition, it is more practical than the use of CPU 56-bit registers B, N and M as they must be accessed over C and this would mean much more code.
Eventually, the logarithm can be expressed as a natural number (better said: its approximation) by the use of the alphanumeric capabilities of the calculator. The result is shown as a.aaaaEbbbb in the alpha registers. For this conversion, a FOCAL program is used.
Listing of MCODE program "LOGSUM"
00C "L" 'Function name: LOGSUM
2A0 SETDEC 'set HP-41 to decimal mode
04E C=0 ALL
00E A=0 ALL 'Emptying of CPU arithmetic registers
228 WRIT 8(P) 'and Register P
0F8 READ 3(X) 'Store on stack in C and save it into
268 WRIT 9(Q) 'Register Q to serve as loop counter
115 ?NC XQ 'Start of loop
06C ->1B45 'Calculate natural log
10E A=C ALL 'Copy log to A
238 READ 8(P) 'Register P will contain the sums of the logs
01D ?NC XQ
060 ->1807 'Add Register P and most recent log
228 WRIT 8(P) 'Store the result back into P
278 READ 9(Q) 'Read the loop counter back into C
10E A=C ALL 'and transfer it to A
04E C=0 ALL
130 LDI S&X 'The mantissa of C will become 1 by placing it
001 CON: 1 'first into S&X (exponent part) and then rotating
23C RCR 2 'it right twice
2BE C=-C-1 MS 'Works like a CHS operation on the mantissa
01D ?C XQ 'and the negative value is added resulting
061 ->1807 'in a subtraction of the loop counter by "1"
268 WRIT 9(Q) 'Store the loop counter back into P
2FA ?C#0 M 'and check whether the mantissa of C=0
377 JC -12 'If not 0, carrier set and jump back to loop start
238 READ 8(P) 'If 0, read the sum of logs
10E A=C ALL 'Store the sum into A
04E C=0 ALL 'Clear C
35C PT= 12 'Pointer at start of mantissa
090 LD@PT- 2 'Load mantissa with 2.302585093 which will be
0D0 LD@PT- 3 'the divisor of the sum of natural logs to
010 LD@PT- 0 'get the 10-based logarithm (ln(10)=2.3025…)
090 LD@PT- 2
150 LD@PT- 5
210 LD@PT- 8
150 LD@PT- 5
010 LD@PT- 0
250 LD@PT- 9
0D0 LD@PT- 3
070 N=C ALL
1BD ?NC XQ
040 ->106F 'Perform "DIVIDE" function and put the result
3E0 RTN 'of logsum(n) into Register X
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
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
I entered a set of numbers and calculated the factorial using XFACT at normal speed and 50x turbo on a 41CL. XFACT was compiled and transferred to HEPRAM.
XFACT 1x Speed 50x Speed Result
(n) (sec) (sec) (output)
10 4.01 1.31 3.6288E6
20 6.23 1.31 2.4329E18
30 8.19 1.33 2.6525E32
50 12.35 1.42 3.0414E64
100 22.55 1.59 9.3326E157
200 44.85 1.98 7.8866E374
500 109.35 3.35 1.2201E1134
1000 216.85 5.58 4.0238E2567
5000 1104.16 23.36 4.2289E16325
10000 n.d. 45.32 2.8462E35659
100000 n.d. 446.79 2.7855E456573
A linear regression resulted in RR of >0.9999 and the 1x-speed slope was about 48-fold larger than the 50x-speed slope. By elimination of the overhead (y-axis intersection), supposedly the key read out and output onto display of about 1 second, I got a 49.5-fold difference of the slope which matches pretty exactly 50x TURBO.
Disclaimer: This code may be used at your own responsibility. I will not be responsible for any kind of inconveniences caused by its use.