(HP-67/97) Combinatorics - Extended factorial, gamma, permutations, combinations

[Dave Britten]

This program provides implementations of factorial (with support for numbers greater than 69), gamma, combinations (nCr), and permutations (nPr). They run in constant time without any looping (slightly faster for integers no greater than 69 when the built-in factorial can be used). It is 92 steps, and fits comfortably on one side of a card.

The program relies on using Nemes' formula to calculate an approximation of the natural logarithm of gamma without overflowing on large inputs. Large combinations/permutations are generally accurate to around 5 decimal places, but accuracy worsens as the difference between x and y increases (precision is lost when subtracting ln(x!) and ln(y!) when they differ greatly in magnitude).

Keys

A - ln(x!)
B - ln(Γ)
b - Γ(x)
C - ln(yCx)
c - yCx
D - ln(yPx)
d - yPx
E - e^x→x*10^y
e - Round Int

Usage

To use any of these functions, enter arguments in the x and y registers, and press the key for the desired function. The functions that return a natural logarithm will not overflow for large inputs/outputs. The shifted versions of Γ, nCr, and nPr return a decimal result to x, with the natural logarithm in y. Thus if the calculation overflows (9.999999999 99), it's not necessary to repeat the calculation, simply press x><y to get the natural logarithm of the result, and optionally press E to convert the power of e to a decimal mantissa in x, and exponent of 10 in y.

Examples

Calculate the natural log of 144!:

144 A
*** 575.0575 (e^575.0575)


Show a decimal approximation of 144! with the natural log of 144! already in x:
E
*** 5.5503
x><y
*** 249 (5.5503*10^249)


Calculate P(71, 3):

71 ENTER 3 f d
*** 342930.0000 (yCx and yPx are automatically rounded to the nearest integer)


Calculate C(278000, 47000):

278000 ENTER 47000 f c
*** 9.999999999 99 (overflow)
x><y
*** 126317.0360 (e^126317.0360)
E
*** 6.1903
x><y
*** 54858 (6.1903*10^54858)

Code:

LBL b
GSB B
GTO 8
LBL A
7
0
x<=y?
GTO 3
RDown
N!
LN
RTN
LBL 3
RDown
1
+
LBL B
STO 7
1
2
*
.
1
RCL 7
/
-
1/x
RCL 7
+
LN
1
-
RCL 7
*
2
PI
*
RCL 7
/
LN
2
/
+
RTN
LBL c
GSB C
GTO 8
LBL C
x><y
STO 0
GSB A
x><y
STO- 0
GSB A
-
RCL 0
GSB A
-
RTN
LBL d
GSB D
GTO 8
LBL D
x><y
STO 0
x><y
STO- 0
RDown
GSB A
RCL 0
GSB A
-
RTN
LBL E
1
0
LN
/
INT
LASTx
FRAC
10^x
RTN
LBL 8
ENTER
e^x
LBL e
ENTER
FRAC
+
INT
RTN

[Albert Chan]

(06-15-2020 01:42 PM)Dave Britten Wrote:  Large combinations/permutations are generally accurate to around 5 decimal places, but accuracy worsens as the difference between x and y increases (precision is lost when subtracting ln(x!) and ln(y!) when they differ greatly in magnitude).

I think you meant precision is lost when subtracting 2 values similar in magnitude.

For nPr where n ≫ r, ln(nPr) = ln(n!) - ln((n-r)!), both nearly equal.
To avoid catastrophic cancellation, we can start with nPr ≈ n^r, and add a correction.

Reusing code from my log of probability of no repetitions (ln_nr)

Code:

function ln_nr(n,s) return 0.5*s* log1p((s-1)*(s-2-6*n)/(6*n*n)) end
function ln_nPr(n,r) return ln_nr(n,r) + r*log(n) end

Coded in HP12C (26 steps), with assumption n ≫ r

Code:

1  STO 0   ; r
2  X<>Y
3  STO 1   ; n
4  LN
5  *       ; r * ln(n)
6  RCL 0
7  2 
8  -
9  RCL 1 
10 6 
11 * 
12 STO* 1
13 -
14 RCL 1
15 /
16 RCL 0
17 1
18 -
19 *
20 1
21 +
22 LN
23 RCL 0
24 2
25 /
26 *       ; ln_nr(n,r)

Note: register X = ln_nr(n,r), Y = r*log(n). To get ln_nPr(n,r), press "+"

Example, for ln_nPr(1e6, 100)

1e6 ENTER 100
[R/S]           → X = -.004950165034, Y = 1381.551056
+                 → ln_nPr = 1381.546106