Permutations & Combinations for large inputs

[Dave Britten]

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

[Dieter]

(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

[Dave Britten]

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

[Werner]

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

[Dave Britten]

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

[Dieter]

(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

[Werner]

(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

[Werner]

(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

[Dieter]

(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

[Massimo Gnerucci]

(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? :)

[Ángel Martin]

(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

[Dieter]

(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

[Ángel Martin]

(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!