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"

08D "M"

015 "U"

013 "S"

007 "G"

00F "O"

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"

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

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.**