(19C/29C) (k-)Permutations, Combinations, Linear Regression and Pseudo-random number

[Archilog]

Here is a set of routines in 96 steps to add to your HP-19C, HP-29C or HP-29E.

The three first routines (steps 01 to 37) 1, 2 and 3 compute:

- the factorial of an integer between 0 and 69 - or number of permutations of n elements without repetition(*);
- the number of k-permutations from a set of n elements;
- the number of combinations of k elements from a set of n elements.

In this program factorials are computed into the k-permutations routine therefore the two routines (steps 01 to 25) are needed; for the same reason computing combinations needs the whole set (steps 01 to 37); this program does not suit for large N and small K - see below for other references.
Please note that input values are not checked for validity.

(*): Number of permutations of n elements with repetition is given by n^k.

Code:


01 LBL 1 | Permutations or Factorials
02 1
03 X>Y?
04 RTN   | If x=0, then 0!=1
05 X<->Y
06 ENTER
07 LBL 2 | k-permutation: P(n,k)=n!/(n-k)!
08 FIX 0 | will display a rounded number - for some results, the mantissa may differ
09 1
10 STO 1
11 Rd    | Roll down
12 X<->Y
13 -
14 Lst X
15 +
16 STO-0
17 Lst X
18 LBL 0
19 STO*1
20 1
21 -
22 ISZ
23 GTO 0
24 RCL 1
25 RTN
26 LBL 3 | Combination
27 STO 2
28 GSB 2 | uses P(n,k)=n!/(n-k)!
29 RCL 2
30 X<->Y
31 STO 2
32 X<->Y
33 GSB 1 | computes k!
34 RCL 2
35 X<->Y
36 /     | [n!/(n-k)!]/k!
37 RTN

Registers 0, 1 and 2 are used. It is recommended to clear those registers before operation.

Examples:
[0] [GSB] [1]: 0!=1
[1] [2] [GSB] [1]: 12!=479001600
[8] [ENTER] [3] [GSB] [2]: P(8,3)=336
[7] [8] [ENTER] [6] [GSB] [3]: C(78,6)=256851595


Here is a longer version to compute Factorials, Permutations and Combinations with bigger n and smaller k.
Another interesting way to proceed should be to use common logarithms like in this program for HP-25.

The fourth routine (steps 38 to 88) computes a linear regression (y=a1.x + a0) with a0, a1 and coefficient of determination r² respectively in register 3, 4 and 5 from a set of statistical datas (registers 10 to 15); then the program waits for a new variable x to predict y.

Code:


38 LBL 4 | Linear regression
39 ENG 3
40 RCL.5 | Sum{x*y}
41 RCL.1 | Sum{x}
42 RCL.3 | Sum{y}
43 *
44 RCL.0 | n
45 /
46 -     | Sum{x*y}-Sum{x}*Sum{y}/n, also used to calculate r²
47 STO 3
48 RCL.2 | Sum{x²}
49 RCL.1 | Sum{x}
50 x²
51 RCL.0 | n
52 /
53 -     | Sum{x²}-Sum²{x}/n, also used to calculate r²
54 STO 5
55 /     | a1=[Sum{x*y}-Sum{x}*Sum{y}/n]/[Sum{x²}-Sum²{x}/n]
56 STO 4 |
57 RCL 5 | ending r² calculus
58 RCL.4 | Sum{y²}
59 RCL.3 
60 x²    | Sum²{y}
61 RCL.0 | n
62 /
63 -     | 2nd part of denominator
64 *     | denominator completed
65 RCL 3 | numerator
66 x²
67 X<->Y
68 /
69 STO 5 | r²
70 RCL.3 | a0 calculus 
71 RCL.0 | n
72 /
73 RCL.1
74 RCL.0
75 /
76 RCL 4 | a1
77 *
78 -     | a0
79 STO 3 | R3=a0; R4=a1; R5=r²
80 RCL 5
81 PAUSE | displaying r²
82 RCL 3
83 PAUSE | displaying a0 (y for x=0)
84 RCL 4 | displays a1 and waits for new x
85 R/S
86 *
87 +     | executes y=a1.x+a0
88 RTN

Registers 3, 4 and 5 are used.

Example (from HP25 Applications Programs):
An eccentric professor of numerical analysis wakes up one morning and feels feverish. A search through his medicine cabinet reveals one oral thermometer which, unfortunately, is in degrees centigrade, a scale he is not familiar with. As he stares disconsolately out his window, he spies the outdoor thermometer affixed to the windowframe. This thermometer, however, will not fit comfortably into his mouth. Still, with some ingenuity...

The professor suspects that the relationship is F = a1 C + a0. If he can get a few data pairs for F and C, he can run a linear regression program to find a1 and a0, then convert any reading in °C to °F through the equation. So tossing both thermometers into a sink of lukewarm water, he reads the following pairs of temperatures as the water cools:

Code:

C | 40.5  | 38.6 | 37.9 | 36.2 | 35.1 | 34.6 |
----------------------------------------------
F | 104.5 | 102  | 100  | 97.5 | 95.5 | 94   |

If the relationship is indeed F = a1 C + a0, what are the values for a1 and a0? What is the coefficient of determination?

Solution:
In RUN mode, press [f] [Σ],

104.5 [ENTER] 40.5 [Σ+]
102 [ENTER] 38.6 [Σ+]
100 [ENTER] 37.9 [Σ+]
97.5 [ENTER] 36.2 [Σ+]
95.5 [ENTER] 35.1 [Σ+]
94 [ENTER] 34.6 [Σ+]
[GSB] [4]

The calculator displays in engineering notation:
- r² (961.3 -03);
- a0 (125.8 00);
- a1 (171.5 00), waiting for a new x.

Suppose the professor puts the centigrade thermometer in his mouth and finds he has a temperature of 37°C. Should he be worried?

Press 37 [R/S] ---------------------> 98.65°F (98.65 00)

It looks he is safe.

You can compute another value by pressing [GSB] [4], entering a new variable x and pressing [R/S].

Pseudo-random number is useful in statistics, here is a quick code to fill the program zone and generate some; you only need to enter a seed [0 ; +1] in X register and to press [GSB] [5] (steps 89 to 96):

Code:


89 LBL 5
90 FIX 2
91 PI
92 +
93 5
94 y^x
95 FRAC
96 RTN

X register displays a pseudo-random number which can be used as seed to generate others.

Examples:
0.5 in X register will generate sequence 0.41, 0.58, 0.45, 0.51, 0.11, 0.26, etc.


Edition: two useless steps were removed from first routine, two links were added.

[Gene]

Great set of routines!

One suggestion on LBL 5 would be to consider dropping the FIX 2 step.

Another would be to consider getting away from the PI / Y^X generator as they are fairly poor per some tests done in the PPC Journal for getting a good distribution of outputs. Perhaps this?

Saves a step too in case a change elsewhere needs to be made.

But seriously, I love this combination of routines. Makes a great edition to the HP-19C/29C function set.

Once I get my LP HP-29C setup, this will be in memory.

Code:

91 LBL 5
92 9
93 9
94 7
95 *
96 FRAC
97 RTN

[Archilog]

(01-15-2020 04:10 PM)Gene Wrote:  Great set of routines!

One suggestion on LBL 5 would be to consider dropping the FIX 2 step.

Another would be to consider getting away from the PI / Y^X generator as they are fairly poor per some tests done in the PPC Journal for getting a good distribution of outputs. Perhaps this?

Saves a step too in case a change elsewhere needs to be made.

But seriously, I love this combination of routines. Makes a great edition to the HP-19C/29C function set.

Once I get my LP HP-29C setup, this will be in memory.

Code:

91 LBL 5
92 9
93 9
94 7
95 *
96 FRAC
97 RTN

Thank you for your opinion, Gene. As you imagine, this last routine was almost only intended to complete the program zone - an economic behavior for every HP-29E (or LP 29C) user to benefit from each step of memory. I just copied it from my HP25 Applications Programs, yours is better, I totally agree.

I had decided to strongly distinguish every outputs by using a different display, but feel free to modify that!

Regards

[David Hayden]

In the factorial function, you could remove steps 7 and 8 (GSB 2, RTN) and just let control fall through to LBL 2.

Also, the combinations function looks like it will overflow for large N and small K. You can avoid this by alternately multiplying and dividing. For example, compute 10C3 by doing 8 * 9 / 2 * 10 / 3. This has the remarkable property that the divisions always result in an integer value!

[Archilog]

Thanks for these tips.

(01-15-2020 08:16 PM)David Hayden Wrote:  In the factorial function, you could remove steps 7 and 8 (GSB 2, RTN) and just let control fall through to LBL 2.

True, I don't see any issue.

Also, the combinations function looks like it will overflow for large N and small K. You can avoid this by alternately multiplying and dividing. For example, compute 10C3 by doing 8 * 9 / 2 * 10 / 3. This has the remarkable property that the divisions always result in an integer value!

I have to try this; will the program fit? That is the question.

[Gene]

(01-15-2020 06:57 PM)Archilog Wrote:  Thank you for your opinion, Gene. As you imagine, this last routine was almost only intended to complete the program zone - an economic behavior for every HP-29E (or LP 29C) user to benefit from each step of memory. I just copied it from my HP25 Applications Programs, yours is better, I totally agree.

I had decided to strongly distinguish every outputs by using a different display, but feel free to modify that!

Regards

Gene: Ah, now I understand the FIX 2. Makes perfect sense. And the 997x RN generator is simply different. It has its own quirks, and is certainly not mine :-) but I suspect any RNs generated will be better distributed with the 997 than the PI generator.

Great set of routines. What OTHER 98 step routine filled programs can we put together for the 29C to fill in some gaps? I'm looking forward to my LP 29C with room for 100 of these programs.

[David Hayden]

(01-15-2020 08:36 PM)Archilog Wrote:  I have to try this; will the program fit? That is the question.

I wrote a version a few weeks ago using the ideas that I mentioned. It fits in a 29C but it's much longer than yours. I just edited my post to include a link to your program.
https://www.hpmuseum.org/forum/thread-14141.html

Dave

[Archilog]

(01-16-2020 03:26 AM)David Hayden Wrote:  

(01-15-2020 08:36 PM)Archilog Wrote:  I have to try this; will the program fit? That is the question.

I wrote a version a few weeks ago using the ideas that I mentioned. It fits in a 29C but it's much longer than yours. I just edited my post to include a link to your program.
https://www.hpmuseum.org/forum/thread-14141.html

Dave

Hello Dave, you are right: actually your better algorithm wouldn't fit in the set.
I modified the first routine as you suggested and added a link to you program. Look in the first post, there is also another way to do...

[Archilog]

Hello Gene,

To briefly answer your question, I already had the idea to write an 'as-complete-as-possible' set of routines in order to study a function (extremum, derivative, integral area, etc.). I think it would be quite challenging...
I also have other ideas more specific, in astronomy for instance.

Now about the distribution: I wrote probably an error saying that the "A(n+1) =frac(997∗A(n))" algorithm is better than the "A(n+1) =frac((π+A(n))^5" algorithm.

I checked this morning the first one in the program, and when I tried to lazily modify the example ("0.5 in X register will generate sequence 0.41, 0.58, 0.45, 0.51, 0.11, 0.26, etc."), to my surprise the sequence became... 0.5, 0.5, 0.5 (ad lib), which is obvious since every odd integer multiplied by 0.5 (or divided by 2) gives a fractional part equaling 0.5.

So I remembered an 'old' article I highly recommend about the subject. Despite some good statistical appearances (see below), a graphical representation of their distribution may show an anormal - or too oriented - result (see picture).

Code:


Seed: 0.123456                         A(n+1) = frac(997∗A(n)) || A(n+1) = frac((π+A(n))^5
MEAN                                      0.50124              ||             0.48313 
VARIANCE                                  0.08247              ||             0.08221
AUTOCORRELATION COEF                      0.24780              ||             0.23630

   

Here is the article I refer to (see page 7). Sorry, it is in French, but I am more French than statistician ;)). However one can find therein good and quick codes in BASIC or for the HP-39GII easy to understand to compare some PRN generators.

I don't know to which PPC articles you refer, be sure I would like to read them, please. I own a DVD of Jake Schwartz if needed.

In conclusion, I didn't modify the last routine. But I am open-minded, I could do so if necessary. :)

Regards!

[Gene]

Ok, let's do this...

This is the SOFTWARE forum and I'm guilty here as well.

Let's move any future discussion of the 19c/29c routines to the general forum.

[Archilog]

(01-16-2020 10:46 PM)Gene Wrote:  Ok, let's do this...

This is the SOFTWARE forum and I'm guilty here as well.

Let's move any future discussion of the 19c/29c routines to the general forum.

Well said!
And I hope that our oversight board will raise the warning level of all of us in order to avoid any further violation of the rules and teach us all good manners!