# HP Forums

Full Version: Permutations & Combinations for large inputs
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
This program will calculate permutations or combinations (nPr and nCr) for arguments that would normally overflow the FACT function when using the standard forms n!/(n-r)! and n!/(r!(n-r)!). This is a good example of using flag 25 to handle errors. The method employed here is to attempt the standard forms, setting flag 25 before each FACT, then bailing out to a routine that performs repeated multiplication and division if FACT overflowed.

The loop routine will also perform the appropriate division for COMB as the multiplication is being performed, to allow evaluating COMB for arguments that would still overflow PERM even with this method.

For small inputs, execution time should be relatively quick, but larger values will take longer due to the loop operation. The below example takes 33 seconds on my unmodified 41CV.

Inputs:
y: n
x: r

XEQ PERM, or XEQ COMB

Example:
234 ENTER
78 XEQ ALPHA COMB ALPHA
Result: 2.6795468 63

Uses registers 20 and 21, and flags 05 and 25 (error trap). No modules required.

Code:
```01 LBL "PERM" 02 CF 05 03 GTO 00 04 LBL "COMB" 05 SF 05 06 LBL 00 07 STO 21 08 X<>Y 09 STO 20 10 SF 25   //Ignore next error 11 FACT 12 FC?C 25 //Was there an error? 13 GTO 01  //Go to loop routine 14 RCL Y 15 ST- L 16 X<> L 17 SF 25 18 FACT 19 FC?C 25 20 GTO 1 21 / 22 FC? 05 23 RTN 24 RCL Y 25 SF 25 26 FACT 27 FC?C 25 28 GTO 01 29 / 30 RTN 31 LBL 01 32 RCL 20 33 1 34 LBL 02 35 FC? 05 //Permutations? 36 GTO 03 //Skip divide by R21 37 RCL 21 38 / 39 LBL 03 40 RCL Y 41 RCL 21 42 - 43 1 44 + 45 * 46 DSE 21 47 GTO 02 48 RTN 49 END```
(05-24-2015 05:39 PM)Dave Britten Wrote: [ -> ]This program will calculate permutations or combinations (nPr and nCr) for arguments that would normally overflow the FACT function when using the standard forms n!/(n-r)! and n!/(r!(n-r)!). This is a good example of using flag 25 to handle errors. The method employed here is to attempt the standard forms, setting flag 25 before each FACT, then bailing out to a routine that performs repeated multiplication and division if FACT overflowed.

Dave, your method requires just one single overflow test: If n! does not overflow, the smaller values (n–r)! and r! won't either. So lines 17,19 and 20 resp. 25, 27 and 28 can be removed. They should be removed because a set flag 25 may cover other errors, e.g. non-integer or negative arguments.

Just for comparison, here is another version that does not require any data registers and runs faster due to two separate loops and some further improvements within these. For instance two counters are used that are set up before the loops start so that no additional arithmetics within the loops are required. Both loops can be combined to one, but in this case especially nPr would run slower. The program also uses the property Comb(n,r) = Comb(n, n–r), which substantially speeds up the program for r > n/2. Flag 25 is not used – checking whether n>69 will do.

Code:
```01 LBL"NPR" // Permutations 02 CF 05 03 GTO 00 04 LBL"NCR" // Combinations 05 SF 05 06 LBL 00 07 x<>y 08 69 09 x<y?     // n>69? 10 GTO 01   // then use iterative method 11 RDN      // else compute result directly 12 FACT     // n! 13 LastX 14 RCL Z 15 - 16 FACT     // (n-r)! 17 /        // nPr = n!/(n-r)! 18 FC?C 05  // Permutations? 19 RTN      // then quit 20 R↑       // else recall r from stack 21 FACT     // r! 22 /        // nCr = nPr / r! 23 RTN 24 LBL 01   // iterative methods start here 25 FS?C 05  // Combinations? 26 GTO 03   // then continue at LBL 03 27 SIGN     // x:=1 28 RCL Z 29 x<>y 30 x>y?     // r < 1 (i.e. = 0)? 31 RTN      // then return 1 32 x<>y 33 RDN 34 LBL 02   // start nPr loop 35 RCL Y    // multiply by n, n-1, n-2, ... 36 * 37 DSE Y    // decrement n 38 LBL 00   // dummy label = NOP 39 DSE Z    // decrement counter 40 GTO 02 41 RTN 42 LBL 03   // prepare nCr loop 43 RDN 44 RCL X 45 RCL Z 46 - 47 LastX 48 x<y? 49 x<>y    // leave min(r, n-r) in y 50 SIGN    // x:=1 51 x>y?    // r resp. n-r < 1 (i.e. = 0)? 52 RTN     // then return 1 53 LBL 04 54 RCL Z   // multiply by n, n-1, n-2, ... 55 * 56 RCL Y   // divide by r, r-1, r-2, ..., 1 57 / 58 DSE Z   // decrement n 59 LBL 00  // dummy label = NOP 60 DSE Y   // decrement r (= counter) 61 GTO 04 62 END```

I sincerely hope you do not mind my suggestion – it includes some basic ideas that may be helpful for other programs as well.

Edit: program modified to handle r=0 (and n–r=0 for nCr) correctly.

Dieter
(05-26-2015 09:26 AM)Dieter Wrote: [ -> ]I sincerely hope you do not mind my suggestion – it includes some basic ideas that may be helpful for other programs as well.

Not at all. Half the reason I post stuff here is to find better ways to do stuff.

I find that if you want to ask an engineer the best way to do something, it's better to simply show him a poor way of doing it.
Common pitfalls for a COMB routine:
- C(8,3) calculated as 8/3*7/2*6 returns 56.00000001
solution : calculate it as 8/1*7/2*6/3, guaranteeing integers all along, but then...:
- C(335,167) overflows while the result is < 1e100
solution: calculate 0.001*n/1*(n-1)/2*(n-2)/3... *1000
1e3 is large enough because we are looking for the largest p for which C(n,p)<1e100,
and that is 168 (C(336,168) = 6e99)
The largest intermediate number formed is p*C(n,p).
- C(n,0) = 1
- C(n,n-1) = C(n,1)
solution: calculate C(n,min(p,n-p))

The following routine avoids these, but may return slightly inaccurate
results > 1e10 due to roundoff. Nothing to be done about that.

Code:
```                L       X       Y       Z       T 00 { 45-Byte Prgm } 01>LBL"COMB"            p       n       z 02 ST- Y                p       n-p     z       - 03 X>Y? 04 X<>Y                 m       n-m     z       - 05 +            m       n       z       -       - 06 1 E-3        m       ,001    n       z       - 07 STx L        0,m     ,001    n       z       - 08 X<>Y         0,m     n       ,001    z       - 09 ISG L        1,m     n       ,001    z       - 10 X=0?            Always False Filler 11 GTO 00 12>LBL 02       1,m     n       10^-3   z       - 13 STx Y 14 LASTX 15 INT 16 ST/ Z 17 Rv 18 DSE X 19 ISG L 20 GTO 02 21>LBL 00        -        -        c/10^3        z        - 22 Rv 23 1 E3 24 x 25 END```
(05-26-2015 09:26 AM)Dieter Wrote: [ -> ]Dave, your method requires just one single overflow test: If n! does not overflow, the smaller values (n–r)! and r! won't either.

Yeah, those can probably go. I'm just hard-wired for defensive coding when it comes to error trapping.

(05-26-2015 09:26 AM)Dieter Wrote: [ -> ]So lines 17,19 and 20 resp. 25, 27 and 28 can be removed. They should be removed because a set flag 25 may cover other errors, e.g. non-integer or negative arguments.

In this case, it shouldn't be terribly harmful. The worst that will happen is it will detect some other error condition and run the looping algorithm instead. But on a related note, you should (nearly) always test flag 25 with FC?C so as not to accidentally leave it set and mask an error later.

(05-26-2015 09:26 AM)Dieter Wrote: [ -> ]Just for comparison, here is another version that does not require any data registers and runs faster due to two separate loops and some further improvements within these. For instance two counters are used that are set up before the loops start so that no additional arithmetics within the loops are required. Both loops can be combined to one, but in this case especially nPr would run slower. The program also uses the property Comb(n,r) = Comb(n, n–r), which substantially speeds up the program for r > n/2. Flag 25 is not used – checking whether n>69 will do.

Good idea with doing multiple loops. It briefly crossed my mind, but I wasn't sure which would be faster on the 41, and hadn't tested it yet.

(05-26-2015 12:22 PM)Werner Wrote: [ -> ]Common pitfalls for a COMB routine:
- C(8,3) calculated as 8/3*7/2*6 returns 56.00000001
solution : calculate it as 8/1*7/2*6/3, guaranteeing integers all along, but then...:
- C(335,167) overflows while the result is < 1e100
solution: calculate 0.001*n/1*(n-1)/2*(n-2)/3... *1000
1e3 is large enough because we are looking for the largest p for which C(n,p)<1e100,
and that is 168 (C(336,168) = 6e99)
The largest intermediate number formed is p*C(n,p).
- C(n,0) = 1
- C(n,n-1) = C(n,1)
solution: calculate C(n,min(p,n-p))

Mine deals with the first two situations, since it tries factorials first, and switches to looped division/multiplication if that overflows. As for r=0 or r>n, I haven't done anything to detect those. Sometimes you just have to accept garbage in garbage out.
(05-26-2015 12:22 PM)Werner Wrote: [ -> ]Common pitfalls for a COMB routine:
- C(8,3) calculated as 8/3*7/2*6 returns 56.00000001
solution : calculate it as 8/1*7/2*6/3, guaranteeing integers all along, but then...:

That's the way I do it in higher-level programming languages where count-up loops are much more relaxed than on calculators like the HP41. ;-) However, small input like in your example is evaluated directly via factorials in the program I posted. Since the program calculates n!/(n–r)! first (which always is an integer) the results should be exact here.

(05-26-2015 12:22 PM)Werner Wrote: [ -> ]- C(335,167) overflows while the result is < 1e100

OK. But at least the version I posted correctly retuns C(336, 168) as 6,0887 E+99. #-)

(05-26-2015 12:22 PM)Werner Wrote: [ -> ]solution: calculate 0.001*n/1*(n-1)/2*(n-2)/3... *1000
1e3 is large enough because we are looking for the largest p for which C(n,p)<1e100,
and that is 168 (C(336,168) = 6e99)
The largest intermediate number formed is p*C(n,p).

That's an interesting approach.

(05-26-2015 12:22 PM)Werner Wrote: [ -> ]- C(n,0) = 1
- C(n,n-1) = C(n,1)
solution: calculate C(n,min(p,n-p))

That's what my version implements as well. However, the C(n, 0) case has to be handled separately. I now posted an update.

(05-26-2015 12:22 PM)Werner Wrote: [ -> ]The following routine avoids these, but may return slightly inaccurate
results > 1e10 due to roundoff. Nothing to be done about that.

Yes, you can't expect 10-digit accuracy with merely 10-digit precision.

Dieter
(05-26-2015 05:14 PM)Dieter Wrote: [ -> ]Since the program calculates n!/(n–r)! first (which always is an integer) the results should be exact here.
Only for n<16 I'm afraid.
eg. PERM(16,5) = 16!/(11!) = 524160.0001. There's 110 cases like this for n<70.
All of them can be corrected with 0.1 + INT.
Worst case is PERM(64,5)=

Then there are a few cases with a result of 10 significa
(05-26-2015 05:14 PM)Dieter Wrote: [ -> ]Since the program calculates n!/(n–r)! first (which always is an integer) the results should be exact here.
Only for n<16 I'm afraid.
eg. PERM(16,5) = 16!/(11!) = 524160.0001. There's 110 cases like this for n<70.
These can be corrected with 0.1 + INT.
Worst case is PERM(64,5) = 914941440.4

Then there are a few cases with a wrong result of exactly 10 significant digits:
PERM(31,9) = 7.315688014e12 (exact 7.315688016e12)
PERM(35,8) = 8.489642628e11 (exact 8.489642624e11)
PERM(38,7) = 6.360609025e10 (exact 6.360609024e10)
39 7
45 7
52 7
54 7
59 6
61 6
64 6
68 6

All of these are larger than 1e10 though, so you're at least able to deliver correct results less than 1e10
(05-27-2015 12:27 PM)Werner Wrote: [ -> ]Only for n<16 I'm afraid.
eg. PERM(16,5) = 16!/(11!) = 524160.0001.

Yes, deviations like this can occur if one of the operands cannot be represented exactly with 10 digits.
Here 16! is 20 922 789 888 000 which is rounded to 2,092278989 E+13.

(05-27-2015 12:27 PM)Werner Wrote: [ -> ]There's 110 cases like this for n<70.
These can be corrected with 0.1 + INT.

That's the lazy solution. Too bad there is no ROUNDI on the '41. ;-)

(05-27-2015 12:27 PM)Werner Wrote: [ -> ]Worst case is PERM(64,5) = 914941440.4

Then there are a few cases with a wrong result of exactly 10 significant digits:
...
All of these are larger than 1e10 though, so you're at least able to deliver correct results less than 1e10

I suggested an upgrade to the combinations routine of the Binomial program that is discussed in another thread in this forum. Here the denominator loop counts up so that the intermediate results are integers.

On the other hand, there is no chance to reliably get a result with 10 exact digits if the calculator works with the same 10-digit precision. So a slight error in the last digit might be acceptable.

Dieter
(05-30-2015 01:07 PM)Dieter Wrote: [ -> ]That's the lazy solution. Too bad there is no ROUNDI on the '41. ;-)

See this thread about CEIL/FLOOR. Do you recognize anybody among the authors? :)
(05-27-2015 12:27 PM)Werner Wrote: [ -> ]
(05-26-2015 05:14 PM)Dieter Wrote: [ -> ]Since the program calculates n!/(n–r)! first (which always is an integer) the results should be exact here.
Only for n<16 I'm afraid.
eg. PERM(16,5) = 16!/(11!) = 524160.0001. There's 110 cases like this for n<70.
These can be corrected with 0.1 + INT.
Worst case is PERM(64,5) = 914941440.4

The NCR and NPR functions on the SandMath can be a good reference. Since they use internal 13-digit accuracy the results are pretty much right on, always integers (no, there isn't any rounding done by the function in case you wonder).

But I'm positive Dieter will find a counter-example to this rule as soon as he starts using the SandMath on V41 ;-)

Cheers,
'AM
(05-31-2015 07:37 AM)Ángel Martin Wrote: [ -> ]But I'm positive Dieter will find a counter-example to this rule as soon as he starts using the SandMath on V41 ;-)

I wonder if the quotient of two numbers rounded to 13 or even 12 significant digits can have an error that would affect the 10-digit result. #-)

BTW, the Sandmath manual I got is revision 44_E ("this compilation revision 4.44.44 (really!)") from November 2012. Is this the current version that matches the Sandmath .mod file you sent?

I'll have to look a bit deeper into this. At the moment I get some strange results, e.g. floor(pi)=1 or ceil(0,5)=–1.

Dieter
(05-31-2015 04:16 PM)Dieter Wrote: [ -> ]...the Sandmath manual I got is revision 44_E ("this compilation revision 4.44.44 (really!)") from November 2012. Is this the current version that matches the Sandmath .mod file you sent?

No it isn't. The latest one is rev. 5.79.89 from February 2015 - I'll email you a copy.

(05-31-2015 04:16 PM)Dieter Wrote: [ -> ]I'll have to look a bit deeper into this. At the moment I get some strange results, e.g. floor(pi)=1 or ceil(0,5)=–1.

What did I say - you just found another bug right off the shoe... thanks!
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :