HP41C: Factorial (kind of) in MCODE « Next Oldest | Next Newest »

 ▼ Frido Bohn Member Posts: 117 Threads: 16 Joined: Jul 2009 05-24-2012, 08:34 AM 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. ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 05-25-2012, 05:21 AM Way to go Frido, you should post it as an article as well. I noticed you use the "hard-coded" Ln(2) in the final division step. One alternative could've been calling [LN10] with "2" in C, shorter code but of course a tad slower. ```04E C=0 ALL 'Clear C 35C PT= 12 'Pointer at start of mantissa 090 LD@PT- 2 115 ?NC XQ 'calculate natural log 06C ->1B45 [LN10] ``` Since the result of [LN10] is also available as 13-digit form, this method would allow a 13-digit division next - done as multiplication by its inverse: ```239 ?NC XQ 060 ->188E [ON/X13] 238 READ 8(P) 'read the sum of logs 13D ?NC XQ 060 ->184F [MP1_10] 3E0 RTN 'could use a GO in previous step instead) ``` This is an interesting way to overcome inherent limitations of the numeric range in the 41. There are a few other great examples in JM-Baillard's pages, you and others may be interested to check them out: Cheers, ÁM Edited: 25 May 2012, 5:29 a.m. ▼ Frido Bohn Member Posts: 117 Threads: 16 Joined: Jul 2009 05-25-2012, 09:27 AM Dear Ángel, First, I was somewhat disappointed that nobody had a reply on my contribution. Now, the MCODE Guru himself made a comment! :-) Sure, the hard code of ln(10) - this was not ln(2)! - is not that elegant. The way you proposed is more straightforward as it uses fully the abilities of the HP41. The code would be changed as such: ```377 JC -12 'End of loop, see my posting 04E C=0 ALL 'Clear C 35C R=12 'Set pointer to nybble 12 050 LD@R 1 'fill it with 1 39C R=0 'Set pointer to nybble 0 050 LD@R 1 'fill it again with 1 so we have 10 in C 115 ?NC XQ 06C ->1B45 'Perform ln(10) 070 N=C 'Put ln(10) in CPU register N 238 READ 8(P) 'Read the sum of logs into C 0AE A<>C ALL 1BD ?NC XQ 'Divide C by N 040 ->106F 3E0 RTN ``` The solution you offered by calculating the inverse (1/x) of ln(10) and multiplying would of course work too. However, an interesting, but also annoying point about the entry points of HP41 MCODE routines is that one cannot simply figure out where the operands have to be put in. For example, the 1/x function (->188E, [ON/X13]) wants its argument in CPU register B, for the multiplication (->184F, [MP1_10]) these are A and C, and for divide (->106F, [DIVIDE] or [/]) the dividend is in A and the divisor is in N. Then you have to consider where the result will be delivered. Sometimes directly into the stack (this is the case with [DIVIDE]), but also into C only from where you have to do an additional step to write it into the stack. This is not very straightforward, and I wish there would be documentation about this. Also, the preparation of the registers before feeding a function costs a lot of code, so that one has to consider about using an HP41 routine or if it would result in less code designing your own routine. I will put this contribution with a few extensions into the articles section, as you suggested. ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 05-25-2012, 09:59 AM Believe me, I'm quite far from being a MCODE guru - not an appropriate term, I can hardly read 50% of the OS code without stumbling! W.r.t. the input/output conventions for the MATH routines, I think reading the 13-digit writeup I did some time ago will help - available at TOS. The first distinction is whether you're using 10-digit or 13-digit forms, then it all falls logically. Sorry I got it wrong, of course it's not Ln(2). In fact there's one entry point in the OS that returns Ln(10) in 13-digit form to register B - so you could just call it: ```2B5 ?NC XQ 068 ->1AAD [LNC10*] 239 ?NC XQ 060 ->188E [ON/X13] 238 READ 8(P) 13D ?NC GO 061 ->184F [MP1-10] ``` Here using [ON/X13] is just a work-around to keep the 13-digit form of Ln(10). There are a couple of examples about that in the mini-paper. Keep up the good work! ÁM Edited: 25 May 2012, 10:02 a.m. ▼ Frido Bohn Member Posts: 117 Threads: 16 Joined: Jul 2009 05-25-2012, 03:58 PM Ángel, Interestingly, the routine ```2B5 ?NC XQ 068 ->1AAD [LNC10*] ``` gives ln(10)*10^-6 as a result in B. If ```239 ?NC XQ 060 ->188E [ON/X13] ``` is executed directly thereafter the result is 1/ln(10). 1/ln(10) multiplied with a natural logarithm results in the corresponding decadic logarithm (log). It appears that HP designed the [LNC10*] function for this purpose. Thus, log would be calculated primarily by the code for the natural logarithm. This confirms my previous observation that a sum of log is less precise than a sum of ln. Obviously, because the error of the conversion adds up in a chain. Thank you for the link. I will dive into your document and see which new insights it will give to me. Saludos Frido ▼ Frido Bohn Member Posts: 117 Threads: 16 Joined: Jul 2009 05-25-2012, 04:37 PM Quote: Interestingly, the routine ```2B5 ?NC XQ 068 ->1AAD [LNC10*] ``` gives ln(10)*10^-6 as a result in B. What a fool I am, this was a 13-digit precision ln(10). Now, I have to complete reading your paper first, Ángel, before I produce more stupidities. ▼ Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 05-26-2012, 07:19 AM Quote: ... this was a 13-digit precision ln(10). exactly! it´s a little confusing at the beginning but once you get the knack of it is actually easier than the 10-digit routines. In fact I´d go as far as saying that the math routines were designed to be 13-digit from the scratch, and the 10-digit truncation came later as the way to interact with the input/output values. Have fun! ▼ Frido Bohn Member Posts: 117 Threads: 16 Joined: Jul 2009 05-26-2012, 09:18 AM Ángel, The code has been transferred to the articles section. It has been improved using 13-digit precision. Thank you for you very valuable contribution! Frido

 Possibly Related Threads... Thread Author Replies Views Last Post HP-41 MCODE: The Last Function - at last! Ángel Martin 0 676 11-08-2013, 05:11 AM Last Post: Ángel Martin Factorial misbehavior Han 10 2,054 10-09-2013, 02:36 AM Last Post: Ed Hodapp Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand 3 990 06-29-2013, 03:39 PM Last Post: Gjermund Skailand 41-MCODE: Auto XEQ+ALPHA possible? Ángel Martin 5 1,323 05-29-2013, 06:15 AM Last Post: Ángel Martin HP 41 Mcode related Questions Michael Fehlhammer 4 1,261 05-10-2013, 07:09 PM Last Post: Michael Fehlhammer How to check / repair a HP 82106A module for HP41C calculator ? Eduardo Mingo 13 2,374 01-07-2013, 11:05 AM Last Post: Diego Diaz A simple question about the HP41C or HP41CV Antoine M. Couëtte 6 1,583 12-16-2012, 04:06 AM Last Post: Antoine M. Couëtte HP41C repair ? Michael Paul 4 1,174 09-30-2012, 11:26 AM Last Post: db (martinez, ca.) 41-MCODE: Breaking the FAT barrier. Ángel Martin 0 629 09-03-2012, 06:31 AM Last Post: Ángel Martin 41-MCODE: Dr. Jekyll & Mr. Hyde Ángel Martin 9 2,134 07-09-2012, 09:41 AM Last Post: Monte Dalrymple

Forum Jump: