HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
nCr and nPr - a 41-MCODE implementation
Posted by Ángel Martin on 27 Nov 2009, 10:46 a.m.
You can refer to this forum treads for a more detailed discussion of this topic:
http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv019.cgi?read=160031#160031
http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv013.cgi?read=39655
Here's the revised code with the following additions:
- Gets the integer part of the input values, forcing them to be positive
- Checks that neither one is Zero, and that n>r
- Uses the minimum of {r, (n-r)} to expedite the calculation time
- The Out Of Range condition is checked at every multiplication, so it's determined as soon as possible
- The chain of multiplication proceeds right-to-left, with the largest quotients first.
- The algorithm works within the nummeric range of the 41. Example: C(335,167) is calculated without problems.
- It doesn't perform any rounding on the results. Partial divisions are done as opposed to calculating first nPr and dividing it by r!
Incidentally, the same checking routines are used by the GCD and LCM functions in the SandMath as well... so there's some nice code reuse.
Hope you enjoy it - pls. send comments and improvements!
Best wishes, ÁM
-----------------------------------------------
Header A49D 092 "R" nCr= nPr / r! Header A49E 003 "C" nCr = PROD{(n-k)/(r-k)} Header A49F 00E "N" k=0,1,…,(r-2),(r-1) NCR A4A0 108 SETF 8 Ángel Martin NCR A4A1 379 PORT DEP: Naturalize inputs NCR A4A2 03C XQ Integer & Positive NCR A4A3 0F4 ->A4F4 [NATXY] NCR A4A4 070 N=C ALL r NCR A4A5 0B8 READ 2(Y) n NCR A4A6 10E A=C ALL NCR A4A7 085 ?NC XQ Calculates (n-r) NCR A4A8 050 ->1421 [Y_MINUS_X] NCR A4A9 2EE ?C#0 ALL is n=r? NCR A4AA 053 JNC +10d case r=n NCR A4AB 128 WRIT 4(L) see if n>r NCR A4AC 23E C=C+1 MS Carry set if n<r NCR A4AD 30D ?C GO NCR A4AE 063 ->18C3 [ERR0] NCR A4AF 0F8 READ 3(X) NCR A4B0 10E A=C ALL r NCR A4B1 085 ?NC XQ Calculates r-(n-r) NCR A4B2 050 ->1421 [Y_MINUS_X] NCR A4B3 23E C=C+1 MS r-(n-r) = 2r-n NCR A4B4 05F JC +11d Carry set if (2r-n)<0, or r<(n-r) NCR A4B5 138 READ 4(L) NCR A4B6 0E8 WRIT 3(X) NCR A4B7 043 JNC +08 Header A4B8 092 "R" nPr = PROD{n-k} Header A4B9 010 "P" k=0,1,…,(r-2),(r-1) Header A4BA 00E "N" Ángel Martin NPR A4BB 104 CLRF 8 NPR A4BC 379 PORT DEP: Naturalize inputs NPR A4BD 03C XQ Integer & Positive NPR A4BE 0F4 ->A4F4 [NATXY] NPR A4BF 04E C=0 ALL NPR A4C0 35C PT= 12 builds "1" in C NPR A4C1 050 LD@PT- 1 NPR A4C2 268 WRIT 9(Q) initial product NPR A4C3 0F8 READ 3(X) min[r,(n-r)] NPR A4C4 10E A=C ALL r NPR A4C5 04E C=0 ALL NPR A4C6 2DC PT= 13 Build "-1" in C NPR A4C7 250 LD@PT- 9 - NPR A4C8 050 LD@PT- 1 1 NPR A4C9 01D ?NC XQ Adds normalized A and C NPR A4CA 060 ->1807 [AD2_10] NPR A4CB 070 N=C ALL (r-k), index counter NPR A4CC 2BE C=-C-1 MS -(r-k) NPR A4CD 3C4 ST=0 Clears Carry if set NPR A4CE 128 WRIT 4(L) -(r-k) NPR A4CF 10E A=C ALL -(r-k) NPR A4D0 0B8 READ 2(Y) n NPR A4D1 01D ?NC XQ Adds normalized A and C NPR A4D2 060 ->1807 [AD2_10] NPR A4D3 10C ?FSET 8 NPR A4D4 063 JNC +12d skip if nPr case NPR A4D5 228 WRIT 8(P) save [n-(r-k)] it in P NPR A4D6 0F8 READ 3(X) r NPR A4D7 10E A=C ALL NPR A4D8 138 READ 4(L) -(r-k) NPR A4D9 01D ?NC XQ Adds normalized A and C NPR A4DA 060 ->1807 [AD2_10] NPR A4DB 10E A=C ALL r-(r-k)=k NPR A4DC 238 READ 8(P) :-( NPR A4DD 0AE A<>C ALL NPR A4DE 261 ?NC XQ Divides A over C NPR A4DF 060 ->1898 [DV2_10] NPR A4E0 10E A=C ALL n-[(r-k)-1] NPR A4E1 278 READ 9(Q) partial product PP NPR A4E2 135 ?NC XQ it *uses* M !! NPR A4E3 060 ->184D [MP2_10] NPR A4E4 2F6 ?C#0 XS See if we overflow? NPR A4E5 289 ?C GO "Out of Range" NPR A4E6 003 ->00A2 [ERROF] NPR A4E7 268 WRIT 9(Q) PP*{n-[(r-k)-1)]} NPR A4E8 0B0 C=N ALL (r-k)-1 NPR A4E9 2EE ?C#0 ALL NPR A4EA 2D7 JC -38d loop back NPR A4EB 278 READ 9(Q) NPR A4EC 0EE B<>C ALL final solution NPR A4ED 391 ?NC GO (REQUIRES Value in B) NPR A4EE 002 ->00E4 [DROPST] Header A4EF 099 "Y" Header A4F0 018 "X" Header A4F1 014 "T" Naturalize inputs Header A4F2 001 "A" INT, ABS, #0 Header A4F3 00E "N" NATXY A4F4 0F8 READ 3(X) r NATXY A4F5 0B8 READ 2(Y) n NATXY A4F6 351 ?NC XQ (includes SETDEC) NATXY A4F7 050 ->14D4 [CHK_NO_S1] NATXY A4F8 088 SETF 5 Integer NATXY A4F9 0ED ?NC XQ Take Integer part NATXY A4FA 064 ->193B [INTFRC] NATXY A4FB 05E C=0 MS make it positive >0 NATXY A4FC 2EE ?C#0 ALL NATXY A4FD 03B JNC +07 NATXY A4FE 0A8 WRIT 2(Y) n (integer, >0) NATXY A4FF 0F8 READ 3(X) r NATXY A500 0ED ?NC XQ Take Integer part NATXY A501 064 ->193B [INTFRC] NATXY A502 05E C=0 MS make it positive >0 NATXY A503 2EE ?C#0 ALL Carry set if NOT Zero NATXY A504 30D ?NC GO NATXY A505 062 ->18C3 [ERR0] NATXY A506 0E8 WRIT 3(X) r (integer, >0) NATXY A507 3E0 RTN
Edited: 25 Oct 2010, 1:28 p.m.