[VA] SRC #012a - Then and Now: Probability

[Fernando del Rey]

Here’s the raw code to get quick results using Emu71/Win, with no attemps to optimize execution speed or to minimize memory usage. I have ideas for how to do both, but I wanted just to test the algorithm and get quick results.

The idea of the algorith is simple. Wherever the man starts, in our case at position (1,1) in the grid, he has a certain probability to move to the adjacent cells.

If the man is located at any if the 3 corner cells of the grid, he has a 1/2 probability to move to any of their two adjacent cells.

If he is located in a border cell, he has a 1/4 probability to move to each of the 4 adjacent cells.

If he is located at any of the remaining inner cells, he has a 1/6 probability to move to each of its 6 adjacent cells.

The program calculates a probability matrix after each step (matrix B) from the results of the probability matrix resulting from the previous step (matrix A).

After some trivial initilization (lines 10 to 50), the program iterates M times (M number of steps) in the big loop ranging from lines 60 to 240.

Line 70 copies matrix B to A and clears matrix B to start the calculation of a new probability matrix at each step.

Lines of code 80, 90 and 100 calculate the probability propagation of the corner cells.

Lines 120 and 130 calculate the probability propagation of the left border cells. Lines 140 and 150 calculate the propagation of the right border cells. Lines 160 and 170 calculate the propagation of the lower border cells.

Lines 190 to 230 calculate the probability propagation of the inner cells.

Finally, to obtain the probabilty of the man finishing at the lower line of the grid, we just need to add the probabilities of all cells in the lower line of the grid (line 250).

If we would like to know the probability of the man finishing at any other location of the grid after M steps, we just need to look at the final probability map in matrix B.

Also, if we want to know what happens if the man starts at an (X,Y) location different from the upper corner, we just need to change line 50 to B(X,Y)=1.

Please note the the program listing is a manual transcription from the code in Emu71/Win. I have tried to be careful to avoid any errors, but I cannot be 100% sure. I don’t know how a to do a copy/paste of the program code in Emu71/Win. If anyone could guide me how to do it, I’d be most grateful.

In a physical 71B, this program will take about 108 minutes to resolve the 30/60 case. As mentioned earlier, it can surely be optimized for speed, but I was looking for simplicity and code clarity just to test the algorithm, not speed or memory usage optimization.

You may have noticed that I am a complete newbe with the 71B. I had to do a quick read of the manual to be able to produce this code as I have practically no experience with the 71B. But I must admit it has been real fun!

Here's the code:

10 DESTROY ALL @ OPTION BASE 1 @ STD
20 INPUT “Grid Size? “;N
30 INPUT “# Steps? “;M
40 REAL A(N,N), B(N,N), F
50 B(1,1)=1
60 FOR K=1 TO M
70 FOR I=1 TO N @ FOR J=1 TO I @ A(I,J)=B(I,J) @ B(I,J)=0 @ NEXT J @ NEXT I
80 F=A(1,1)/2 @ B(2,1)=F @ B(2,2)=F
90 F=A(N,1)/2 @ B(N-1,1)=F @ B(N,2)=F
100 F=A(N,N)/2 @ B(N-1,N-1)=F @ B(N,N-1)=F
110 FOR I=2 TO N-1
120 F=A(I,1)/4 @ B(I-1,1)=B(I-1,1)+F
130 B(I,2)=B(I,2)+F @ B(I+1,1)=B(I+1,1)+F @ B(I+1,2)=B(I+1,2)+F
140 F=A(I,I)/4 @ B(I-1,I-1)=B(I-1,I-1)+F
150 B(I,I-1)=B(I,I-1)+F @ B(I+1,I)=B(I+1,I)+F @ B(I+1,I+1)=B(I+1,I+1)+F
160 F=A(N,I)/4 @ B(N,I-1)=B(N,I-1)+F
170 B(N-1,I-1)=B(N-1,I-1)+F @ B(N-1,I)=B(N-1,I)+F @ B(N,I+1)=B(N,I+1)+F
180 NEXT I
190 FOR I=3 TO N-1 @ FOR J=2 TO I-1
200 F=A(I,J)/6 @ B(I-1,J-1)=B(I-1,J-1)+F @ B(I-1,J)=B(I-1,J)+F
210 B(I,J-1)=B(I,J-1)+F @ B(I,J+1)=B(I,J+1)+F
220 B(I+1,J)=B(I+1,J)+F @ B(I+1,J+1)=B(I+1,J+1)+F
230 NEXT J @ NEXT I
240 NEXT K
250 F=0 @ FOR I=1 TO N @ F=F+B(N,I) @ NEXT I @ DISP “Pr.=”;F

[J-F Garnier]

I just finished to debug my code.
Valentin's hint to check that the sum of the probability of all cells should be 1 was very helpful to discover several coding errors :-)

The principle of my code is, I believe, similar to Fernando's solution, but the triangle cells are computed in a different order, row by row.
Also I took advantage of the symmetry of the starting condition.
Note that I used the fact that the HP BASIC doesn't enter a FOR loop if the final index is smaller than the starting index. This is not true in all BASIC, if needed enable back lines 160 and 220.
Note also that I used some HP-71 Math ROM matrix statements (and for the fun, one statement of the Math 2B version :-)

My code is longer, partially because I avoided multi-statement lines for clarity.
The result is very slightly different: 9.51234350205E-6 [updated, was 9.51234350213E-6 in the previous version]

Here is my code [updated: lines 90, 140 and 240 replaced by lines 91, 141 and 241]:

 10 OPTION BASE 1
 20 N=30 ! rows
 30 S=60 ! steps
 40 DIM A(N,N),B(N,N)
 50 ! set up A
 60 MAT A=ZER @ MAT B=ZER
 70 A(1,1)=1
 80 FOR X=1 TO S
 90 ! B(1,1)=(A(2,1)+A(2,2))/4 ! row 1
 91   B(1,1)=A(2,1)/4+A(2,2)/4 ! row 1
100   B(2,1)=A(1,1)/2+A(2,2)/4+A(3,1)/4+A(3,2)/6
120   B(2,2)=B(2,1) ! row 2
130   B(3,1)=A(2,1)/4+A(3,2)/6+A(4,1)/4+A(4,2)/6
140 ! B(3,2)=(A(2,1)+A(2,2)+A(3,1)+A(3,3))/4+(A(4,2)+A(4,3))/6
141   B(3,2)=A(2,1)/4+A(2,2)/4+A(3,1)/4+A(3,3)/4+A(4,2)/6+A(4,3)/6
150   B(3,3)=B(3,1) ! row 3
160   ! IF N<6 THEN 310
170   ! rows 4..N-2
180   FOR I=4 TO N-2
190     B(I,1)=A(I-1,1)/4+A(I,2)/6+A(I+1,1)/4+A(I+1,2)/6
210     B(I,2)=A(I-1,1)/4+A(I-1,2)/6+A(I,1)/4+A(I,3)/6+A(I+1,2)/6+A(I+1,3)/6
220     ! IF I<5 THEN 270
230     FOR J=3 TO INT((I+1)/2)
240 !     B(I,J)=(A(I-1,J-1)+A(I-1,J)+A(I,J-1)+A(I,J+1)+A(I+1,J)+A(I+1,J+1))/6
241       B(I,J)=A(I-1,J-1)/6+A(I-1,J)/6+A(I,J-1)/6+A(I,J+1)/6+A(I+1,J)/6+A(I+1,J+1)/6
250       B(I,I+1-J)=B(I,J)
260     NEXT J
270     B(I,I-1)=B(I,2)
280     B(I,I)=B(I,1)
290   NEXT I
310   ! row N-1
320   B(N-1,1)=A(N-2,1)/4+A(N-1,2)/6+A(N,1)/2+A(N,2)/4
330   B(N-1,2)=A(N-2,1)/4+A(N-2,2)/6+A(N-1,1)/4+A(N-1,3)/6+A(N,2)/4+A(N,3)/4
340   FOR J=3 TO INT(N/2)
350     B(N-1,J)=A(N-2,J-1)/6+A(N-2,J)/6+A(N-1,J-1)/6+A(N-1,J+1)/6+A(N,J)/4+A(N,J+1)/4
360     B(N-1,N-J)=B(N-1,J)
370   NEXT J
380   B(N-1,N-1)=B(N-1,1)
390   B(N-1,N-2)=B(N-1,2)
400   ! row N
410   B(N,1)=A(N-1,1)/4+A(N,2)/4
420   B(N,2)=A(N,1)/2+A(N,3)/4+A(N-1,1)/4+A(N-1,2)/6
430   FOR J=3 TO INT((N+1)/2)
440     B(N,J)=A(N,J-1)/4+A(N,J+1)/4+A(N-1,J-1)/6+A(N-1,J)/6
450     B(N,N+1-J)=B(N,J)
460   NEXT J
470   B(N,N-1)=B(N,2)
480   B(N,N)=B(N,1)
490   !
500   MAT A=B
550 NEXT X
555 ! output result
560 DIM C(N)
570 MAT C=RSUM(A)
580 DISP "PROB.=";C(N)

>RUN
PROB.= 9.512343502105-6


Now, I would be curious to see RPN/RPL solutions :-)

J-F

[rawi]

Very nice and complicated problem. Thank you very much, Valentin.

Since I was not able to find the exact solution (like Fernando was, chapeau!) I did what statisticians do if brain is not enough: They replace brain by computer power and make simulations of the proplem. So I did with the DM42, which is not really vintage but vintage in programming. So I hope it is within the rules.
I did 250.000 simulations (no typo, it took about 9 hours with power supply attached) and got in total 2 cases where the man ended at the bottom line of the triangle with a triangle size of 30 and 60 movements. This is 8E-6 which is not bad compared to the exact solution of 9.5 E-6. Here is the code:

00 { 203-Byte Prgm }
01 LBL „VA“
02 „ROWS?“
03 PROMPT
04 STO 00
05 „MOVES?“
06 PROMPT
07 STO 01
08 „SIMUL?“
09 PROMPT
10 STO 04
11 „SEED?“
12 PROMPT
13 SEED
14 0
15 STO 05
16 STO 06
17 LBL 00
18 RCL 01
19 1000
20 /
21 STO11
22 1
23 STO 02
24 STO 03
25 STO+06
26 LBL 07
27 0
28 X<>F
29 RCL 03
30 1
31 X=Y?
32 SF 02
33 X=Y?
34 SF 03
35 X<>Y
36 RCL 02
37 X=Y?
38 SF 06
39 X=Y?
40 SF01
41 RCL 00
42 X=Y?
43 SF 04
44 X=Y?
45 SF 05
46 1
47 ENTER
48 ENTER
49 ENTER
50 FS? 01
51 +
52 FS? 02
53 +
54 FS? 03
55 +
56 FS? 04
57 +
58 FS? 05
59 +
60 FS? 06
61 +
62 +/-
63 7
64 +
65 1/X
66 RAN
67 X<>Y
68 /
69 IP
70 1
71 +
72 STO 08
73 0
74 STO 09
75 STO 10
76 LBL 08
77 1
78 STO+ 09
79 FS? IND 09
80 GTO 08
81 STO+ 10
82 RCL 10
83 RCL 08
84 X>Y?
85 GTO 08
86 1
87 XEQ IND 09
88 ISG 11
89 GTO 07
90 RCL 00
91 RCL 02
92 X=Y?
93 XEQ 09
94 RCL 06
95 RCL 04
96 X>Y?
97 GTO 00
98 RCL 05
99 RCL 04
100 /
101 RTN
102 LBL 09
103 1
104 STO+ 05
105 RTN
106 LBL 01
107 STO- 02
108 RTN
109 LBL 02
110 STO- 02
111 STO- 03
112 RTN
113 LBL 03
114 STO- 03
115 RTN
116 LBL 04
117 STO+ 02
118 RTN
119 LBL 05
120 STO+ 02
121 STO+ 03
122 RTN
123 LBL 06
124 STO+ 03
125 RTN
126 END

Some explanations:
Number: Register, L+Number: Line in program
ROWS = Size of the triangle, e.g. 30
MOVES = Steps through the triangle, e.g. 60
SIMUL = Number of Simulations, e.g. 10000
SEED = Starting value of random number generator
LBL 00: Start of simulation
01 Number of moves
02 Actual row position in the triangle
03 Actual column position in the triangle
06 Count of simulation
LBL 07: Start of random walk through triangle, starting at the top (row 1, column 1)
L28: Deleting flags; flags are used to denote
those directions that are blocked starting with FLAG 1 for the upper left direction and then FlAGS 2 to 6 in clockwise direction.
L46 – L64: Getting the number of directions the man can go (depending on position)
L66: Random number for simulated walk
08: Random walk; 1: first possible, starting with upper left direction, 2-6 clockwise
LBL 08: Getting the direction of the random walk
L77-80: Skipping walks that are not allowed because of actual position
09: Direction 1-6 with starting 1 at upper left and 2-6 in a clockwise direction
L90-93: Add 1 to Reg 05 if man ends at bottom
LBL 09: Counting number of cases when random walk ends at the button of the triangle
LBL 1 to LBL 6: Changing position in triangle.
Row from top to bottom, columns from right to left

[Vincent Weber]

Hi Valentin,

I renounced Monte-Carlo simulations and managed to get a direct solving which I find quite elegant. I managed to find correct results.

But I have a dilmena: I did so by using recursion, and the vanilla 17B/27S calculators don't have user defined functions, let alone recursive ones. I had to use Plus42, which extends the HP solver with such functions.

Can I still post the code, or is it too much cheating, which I would perfectly understand?

Best regards

[Valentin Albillo]

.
Hi, Vincent,

(10-11-2022 10:19 PM)Vincent Weber Wrote:  Hi Valentin,

I renounced Monte-Carlo simulations and managed to get a direct solving which I find quite elegant. I managed to find correct results.

But I have a dilmena: I did so by using recursion, and the vanilla 17B/27S calculators don't have user defined functions, let alone recursive ones. I had to use Plus42, which extends the HP solver with such functions.

Can I still post the code, or is it too much cheating, which I would perfectly understand?

Sorry, Vincent, but using Plus42 with extensions not available in vintage HP calcs completely defeats the purpose of this SRC #012, which is to demonstrate that vintage HP calcs (or at least "very old") can still solve nowadays recently proposed non-trivial problems intended to be solved in modern PCs.

Using Plus42 extensions or Python or Matlab or C++ may be pretty satisfying and elegant but demonstrates nothing of the sort and has no place here.

I'd suggest you post your surely-very-interesting solution to a parallel thread, as suggested by pier4r (post #13 above) and Super Moderator rprosperi (post #14 above) for posts like yours that would not fit my stated rules.

Would you, please ?

Thanks a lot for asking and best regards.
V.

[Vincent Weber]

Hi Valentin,

I fully understand and respect that, in fact I was expecting this answer

Before posting my code in another thread, let me try to remember how to convert a recursive algorithm into an iterative one

Best regards

[Werner]

(10-11-2022 10:48 PM)Vincent Weber Wrote:  Hi Valentin,

I fully understand and respect that, in fact I was expecting this answer

Before posting my code in another thread, let me try to remember how to convert a recursive algorithm into an iterative one

Best regards

I'd like to see the recursive one ;-)
Werner
I have a solution much along the lines of the '71 programs above, but for the 42S - but the memory requirements are too high for the 42S ;-) and Free42 is Not Allowed ;-)
(on the other hand, 32K 42S's exist..)

[PeterP]

Hi Valentin,

I was not able to find a solution that can work on a standard HP41 so I worked on a very limited edge case of steps equal to rows minus one, like the example you provided (btw - it was the provision of that example that allowed little ol’ me to engage and learn, perfectly chosen, thank you!)

My code does deliver the correct result for R = 5, but I dont have a good way (especially right now on a plane and my work computer has no simulators installed…) to check if it is correct for R = 30, S=29. (It comes out to 1.311095094 e-13). And given this is such a specific edge case of your wonderful problem I assume posting code and explanation here would not be in your spirit anyway.

However I wanted to thank you for a couple hours of respite from powerpoints with wonderful recreational math and programming and thinking, its been a loon time since I had this pleasure.

And, as always, I learned a lot in the process, which is the most enjoyable part for me.

[Werner]

We're going about this the wrong way, I'm sure, and as usual, Valentin has given us a few clues.
The first being that he has a simple program that can calculate the desired answer for any number of rows and steps, quickly. So our (exponential) approach of summing up all point probabilities in each step is correct, but not feasible on our vintage calculators, for lack of speed or memory, or both.
The second being the test case for R=5 and S=4, with the number of steps just enough to reach the last row. Now, this probability is a lot easier to calculate as each successive row probability only depends on the previous row, and there's no need to keep the whole triangle.
And, the third that the sum of probabilities over the triangle necessarily has to be 1.
So, we have to ask ourselves: how does the last-row probability for R rows and R-1 steps change if we take an extra step?
The new last-row probability is the sum of
- the probability of staying in the last row, which is easy, as it is half of the total probability of being in the row at step R-1
and
- the probability of coming down from row R-1
and here I'm stuck, as that is not easy.. or I'm missing the obvious.
Coming down from the edge has a probability of 1/2 and coming down from any middle point has a probability of 1/3, but you'd still need to know all the point probabilities as well.

Hope it gave someone else an idea..
Cheers, Werner

[J-F Garnier]

More comments on Fernando's solution and mine.

The basic principle is the same iterative method, starting from the probability matrix at step N, it calculates the probability at step N+1.
Fernando scanned the previous matrix cells (step N-1) and distributes the probabilities to the next step matrix (N), with several statements like:
F=A(I,I)/6 @ B(I+1,I)=B(I+1,I)+F ...
On my side, I scanned the next matrix cells (step N) and computes the cell probability (in one go) from the previous cells (step N-1) with something like:
B(I,J)=(A(I-1,J-1)+A(I-1,J)+A(I,J-1)+A(I,J+1)+A(I+1,J)+A(I+1,J+1))/6

Both methods are equivalent, but Fernando's code is shorter.
However, I would have thought that my method was more accurate since most of the cell probabilities (the inner cells) are computed in one go, with only one division by 6.

But it's the contrary: Fernando's result is exact to 12 places, whereas my result had an "error" of 6 ULP.
To figure out why my answer was less accurate, I tried to mimic Fernando's calculation by not factoring the divisor term but distribute it to each cell term, that is replacing:
B(I,J)=(A(I-1,J-1)+A(I-1,J)+A(I,J-1)+A(I,J+1)+A(I+1,J)+A(I+1,J+1))/6
by
B(I,J)=A(I-1,J-1)/6+A(I-1,J)/6+A(I,J-1)/6+A(I,J+1)/6+A(I+1,J)/6+A(I+1,J+1)/6
In this way our methods sum exactly the same terms, but in a different order.

My result is now 9.51234350205E-6 with a difference of only 2 ULP (instead of 6) from the "exact" value.

I updated my code above accordingly.

Now, the question is: why is it better to distribute the divisions to each term, instead of factoring it?
Is it better just by chance?

J-F

[Werner]

Hi J-F,
Can't say why it seems better in your case to do the divisions straight away, but Fernando's code calculates the edges first, and those are the smallest numbers.
It's like summing a row of numbers: if you sum them from small to large, the result will be more accurate.
Cheers, Werner

[Fernando del Rey]

I must admit it’s pure chance that my method seems to be more accurate than J-F’s. I started by summing the corners first and borders next as it seemed easier to program, but I was not considering that those cells would have smaller values and thus give a more accurate result when adding values as commented by Werner.

In my first code I had introduced a silly error in one of the corner cells (I forgot to divide it by 2 in the propagation). I was checking that the addition of all values in the probability matrix B was close to 1, and I noticed that it was a bit higher than 1, but initially I thought the difference was coming from the accumulation of rounding errors.

Later, I discovered the error in the code, and after correcting it I got a value of 1 accurate to 2 ULP, which gave me confidence that the code was now correct.
I have some ideas to improve the program for speed and memory usage. For speed, you can consider that if you are in iteration M all cells below row M+1 will be zero, and you can save the divisions and additions of zeros.

For memory usage, as the matrices A and B are not used to the right of their diagonal, you could work with a single matrix of size Nx(N+1) using the left side of the diagonal for holding results of iteration M and the right side of the diagonal to calculate results of iteration M+1. But I’m afraid it would make the code much less readable.

And you could also consider the symmetry of the solution if the man is starting at cell (1,1), calculating only half of the grid. But then the algorithm would not be valid for a starting position which is not located in the central column of the grid, which is therefore not symmetrical.

In the end, getting the result of the 30/60 case in less than two hours on a plain 71B (with no Math ROM), might still have been considered acceptable at the times of the 71B.

But I am convinced that Valentin will show us a much cleverer and faster method to arrive at the correct solution

[C.Ret]

(10-13-2022 08:14 AM)J-F Garnier Wrote:  I updated my code above accordingly.

I get exactly the same value. But it's not surprising, even if our codes don't look alike, I actually do the calculation in the same direction as Jean-François.

It takes 5708.4 seconds (nearly 2 hours) to display the result using this compact code that is more readable by using line breaks and indents to show all nested loops:
10  DESTROY ALL @ OPTION BASE 1 @ DELAY 0
  @ INPUT "[VA]SRC012a R,S=";R,S @ T0=TIME
20  REAL T,W(R,R),P(R,R),Q(R,R) @ INTEGER A,B,I,J,K,M 
30  DISP "Init"
  @ FOR I=1 TO R
  @   FOR J=1 TO I
  @     W(I,J)=2*(3-(J=1)-(I=J)-(I=R))
  @   NEXT J
  @ NEXT I
40 DISP "Comp"
  @ Q(1,1)=1 @ M=2
  @ FOR K=1 TO S
  @   FOR I=1 TO M
50      FOR J=1 TO CEIL(I/2)
  @       P(I,J)=0 
60        FOR A=MAX(1,I-1) TO MIN(I+1,M)
  @         FOR B=MAX(1,J-(A<=I)) TO MIN(J+(I<=A),A)
70            IF A<>I OR J<>B THEN P(I,J)=P(I,J)+Q(A,B)/W(A,B)
80          NEXT B
  @       NEXT A
  @       P(I,1+I-J)=P(I,J)
  @     NEXT J
  @   NEXT I
  @   IF M<R THEN M=M+1
90    DISP K;TIME-T0
  @   INVERSE 0,131*K*M/S/R
  @   MAT Q=P
  @ NEXT K
100 T=0
  @ FOR J=1 TO R @ T=T+P(R,J) @ NEXT J @ DISP TIME-T0;R;STR$(S);T



To save time, I use the symmetry of the problem by calculating only the right half of each row of the equilateral triangle. In the I-th row, columns J and 1+I-J are symmetrical positions with equal probability P(I,J) = P(I,1+I-J).

Also, after M steps, the man goes further than the M+1 row. This leaves all the probabilities of the following rows at zero, which limits calculations and saves a little time.

My code uses two instructions from specific modules.
* The first is the INVERSE instruction of the JPC ROM module which show the progress of the calculation by inverting the display. Suffice to say that it does not speed things up and that we can easily do without it. It's just a gadget to make the user wait.
* The second is a MAT Q=P instruction which allows you to 'quickly' copy the newly calculated probabilities from table P into the table Q. I was hoping for a smarter use of my new module. but, for the moment, I have not found a way to perform the calculation more directly. By tha way,

Note the power of the HP-71B, while other systems do not have enough registers to memorize all the points of the triangle, the HP-71B allows, for simplicity, to use half of two matrices.
Moreover, I use a third W array to store the weights of each point. So much memory!



But I am not satisfied with this version. I'm like Werner, I think we're on the wrong track and there's a smarter way that avoids calculating all the iterations and probabilities of every point in the triangle.

So maybe my MATH module can be better exploited. Besides, I discover with horror that my 1A version does not have an RSUM command or any other means of obtaining the sum of the probabilities of the last row in one instruction.

Perhaps calculating iteratively, the probabilities of each point of the triangle will help us for the following questions?

[Valentin Albillo]

.
Hi, all,

A few comments before I post my original solution in a few days ...

Vincent Weber Wrote:Before posting my code in another thread, let me try to remember how to convert a recursive algorithm into an iterative one

It will be quite a sight for sure, please do.

PeterP Wrote:[...] like the example you provided (btw - it was the provision of that example that allowed little ol’ me to engage and learn, perfectly chosen, thank you!)

I'm glad it was so helpful to you. I believe that including a 'toy' version of a problem helps immensely to detect and iron out bugs in one's program, as well as eficiency issues, e.g. if the toy case takes too long then the real McCoy will be hopeless so one has better improve the algorithm instead.

PeterP Wrote:I dont have a good way (especially right now on a plane and my work computer has no simulators installed…) to check if it is correct for R = 30, S=29. (It comes out to 1.311095094 e-13).

Looks pretty good for the 10-digit 41C, I get 1.31109509664e-13 in the 12-digit 71B.

PeterP Wrote:And given this is such a specific edge case of your wonderful problem I assume posting code and explanation here would not be in your spirit anyway.

On the contrary, go ahead with the 41C code, the more the merrier.

PeterP Wrote:However I wanted to thank you for a couple hours of respite from powerpoints with wonderful recreational math and programming and thinking, its been a loon time since I had this pleasure. And, as always, I learned a lot in the process, which is the most enjoyable part for me.

Wow, what can I say ... many thanks and I'm extremely glad you enjoyed it and even learned while dealing with it, that's the idea !

Fernando del Rey Wrote:But I am convinced that Valentin will show us a much cleverer and faster method to arrive at the correct solution.

You're such a good friend, Fernando, but let's not overhype my abilities lest disillusionment ensues, ok ?


As a general remark, I'm somewhat mystified that no RPL solutions have been posted or even discussed so far. I know that writing RPL code adds an enormous layer of sheer incomprehensibleness to the task but still ...

Best regards.
V.

[Albert Chan]

(10-13-2022 08:14 AM)J-F Garnier Wrote:  B(I,J)=A(I-1,J-1)/6+A(I-1,J)/6+A(I,J-1)/6+A(I,J+1)/6+A(I+1,J)/6+A(I+1,J+1)/6

When we divide by 6, we almost always get an error of 1/3 ULP
It may be better to go for scaled probability (by 6^S), then unscale it.

10 OPTION BASE 1 ! JF scaled P version
20 N=30 ! rows
30 S=60 ! steps
40 T0=TIME @ DIM A(N,N),B(N,N)
50 ! set up A
60 MAT A=ZER @ MAT B=ZER
70 A(1,1)=1
80 FOR X=1 TO S
90 B(1,1)=(A(2,1)+A(2,2))*1.5 ! row 1
100 B(2,1)=A(3,2) + (A(2,2)+A(3,1))*1.5 + A(1,1)*3
120 B(2,2)=B(2,1) ! row 2
130 B(3,1)=(A(3,2)+A(4,2)) + (A(2,1)+A(4,1))*1.5
140 B(3,2)=(A(4,2)+A(4,3)) + (A(2,1)+A(2,2)+A(3,1)+A(3,3))*1.5
150 B(3,3)=B(3,1) ! row 3
160 ! IF N<6 THEN 310
170 ! rows 4..N-2
180 FOR I=4 TO N-2
190 B(I,1)=(A(I,2)+A(I+1,2)) + (A(I-1,1)+A(I+1,1))*1.5
210 B(I,2)=(A(I-1,2)+A(I,3)+A(I+1,2)+A(I+1,3)) + (A(I-1,1)+A(I,1))*1.5
220 ! IF I<5 THEN 270
230 FOR J=3 TO INT((I+1)/2)
240 B(I,J)=(A(I-1,J-1)+A(I-1,J)+A(I,J-1)+A(I,J+1)+A(I+1,J)+A(I+1,J+1))
250 B(I,I+1-J)=B(I,J)
260 NEXT J
270 B(I,I-1)=B(I,2)
280 B(I,I)=B(I,1)
290 NEXT I
310 ! row N-1
320 B(N-1,1)=A(N-1,2) + (A(N-2,1)+A(N,2))*1.5 + A(N,1)*3
330 B(N-1,2)=(A(N-1,3)+A(N-2,2)) + (A(N-2,1)+A(N-1,1)+A(N,2)+A(N,3))*1.5
340 FOR J=3 TO INT(N/2)
350 B(N-1,J)=(A(N-2,J-1)+A(N-2,J)+A(N-1,J-1)+A(N-1,J+1)) + (A(N,J)+A(N,J+1))*1.5
360 B(N-1,N-J)=B(N-1,J)
370 NEXT J
380 B(N-1,N-1)=B(N-1,1)
390 B(N-1,N-2)=B(N-1,2)
400 ! row N
410 B(N,1)=(A(N-1,1)+A(N,2))*1.5
420 B(N,2)=A(N-1,2) + (A(N,3)+A(N-1,1))*1.5 + A(N,1)*3
430 FOR J=3 TO INT((N+1)/2)
440 B(N,J)=(A(N-1,J-1)+A(N-1,J)) + (A(N,J-1)+A(N,J+1))*1.5
450 B(N,N+1-J)=B(N,J)
460 NEXT J
470 B(N,N-1)=B(N,2)
480 B(N,N)=B(N,1)
490 !
500 MAT A=B
550 NEXT X
555 ! output result
560 DIM C(N)
570 MAT C=RSUM(A)
580 DISP "PROB.="; C(N)/6^S, TIME-T0

>RUN
prob.= 9.51234350205E-6            69.54

Since most vertices can go 6 ways, scale by 6 also speed up calculations
Compare to original version, scaling speed up = 101.77/69.54 - 1 ≈ 46%

---

Trivia: if S goes infinite, eventually we reach an equilibrium

P(corners) : P(edges) : P(rest) = 1 : 2 : 3

3 + 3*(n-2) + (n-2)*(n-3)/2 = n*(n+1)/2

[3, 3*(n-2), (n-2)*(n-3)/2] * [1,2,3] = 3*n*(n-1)/2

P(bottom row, s=inf) = [2, n-2, 0] * [1,2,3] / (3*n*(n-1)/2) = 4/(3*n)

>N=5 @ S=100 @ RUN 40
prob.= .266666666668            2.68

>4/(3*N)
 .266666666667

[J-F Garnier]

(10-13-2022 09:44 PM)C.Ret Wrote:  

(10-13-2022 08:14 AM)J-F Garnier Wrote:  I updated my code above accordingly.

I get exactly the same value. But it's not surprising, even if our codes don't look alike, I actually do the calculation in the same direction as Jean-François.

It takes 5708.4 seconds (nearly 2 hours) to display the result using this compact code that is more readable by using line breaks and indents to show all nested loops:
[...]

Your code is remarkably compact, at the expense of lower speed due to the two inner loops to explore the neighbouring points.
Thanks for formatting your code in a legible way, multi-statement line codes are usually very hard to read!

Quote:Note the power of the HP-71B, while other systems do not have enough registers to memorize all the points of the triangle, the HP-71B allows, for simplicity, to use half of two matrices.
Moreover, I use a third W array to store the weights of each point. So much memory!

Yes, the BASIC language is really much conformable to use for these kinds of problems.

Quote:Besides, I discover with horror that my 1A version does not have an RSUM command or any other means of obtaining the sum of the probabilities of the last row in one instruction.

You may consider using the enhanced MATH 2B

J-F

[Vincent Weber]

(10-13-2022 11:27 PM)Valentin Albillo Wrote:  

Vincent Weber Wrote:Before posting my code in another thread, let me try to remember how to convert a recursive algorithm into an iterative one

[...]
As a general remark, I'm somewhat mystified that no RPL solutions have been posted or even discussed so far. I know that writing RPL code adds an enormous layer of sheer incomprehensibleness to the task but still ...

Hi Valentin,

Unfortunately I have to renounce making an interative solution for the vanilla 17B/27S. I realize that this would require full indirect addressing of registers. By full I mean read and write. But the HP solver only allows reading values of a list indirectly (with ITEM[I]), not writing them (the L(ITEM[I] :X) construct, while straightforward, is illegal). Why HP didn't allow this, which would have brought their powerful and elegant solver language to Turing-compatibility level to solve any kind of problems, is beyond me...

Of course Plus42 allows this construct, but is out of the competition, so I renounce this.

But I might try to implement my recursive algorithm in RPL and see how it performs

Best regards

[J-F Garnier]

During my search for another (better) solution, I considered counting the different paths and attempting to deduce the probability from the path tree. For the case 5 rows/4 steps, I found the same numbers of paths as Peter, with 16 paths ending at the bottom row over a total of 144. And here I was puzzled because the probability 16/144=0.111111... is significantly different from Valentin's hint 23/288=0.07986... confirmed by the results of Fernando's program or mine.

So to remove any doubt, I attempted to use the Monte-Carlo statistical approach to confirm the probability in this particular case, since 5 rows/4 steps is still manageable in this way.
Here is my program:

 10 OPTION BASE 0
 20 DIM G(15,6)
 30 ! the nodes are coded from P=1 to 15:
 40 !        1
 50 !       2 3
 60 !      4 5 6
 70 !     7 8 9 10
 80 !   11 12 13 14 15
 90 ! matrix G(,) codes the graph:
100 ! G(P,0) codes the number of neighbours of node P
110 ! G(P,1..6) code the neighbouring nodes, up to 6
120 ! the zeros are just fillers.
130 DATA 0,0,0,0,0,0,0
140 DATA 2,2,3,0,0,0,0
150 DATA 4,1,3,4,5,0,0
160 DATA 4,1,2,5,6,0,0
170 DATA 4,2,5,7,8,0,0
180 DATA 6,2,3,4,6,8,9
190 DATA 4,3,5,9,10,0,0
200 DATA 4,4,8,11,12,0,0
210 DATA 6,4,5,7,9,12,13
220 DATA 6,5,6,8,10,13,14
230 DATA 4,6,9,14,15,0,0
240 DATA 2,7,12,0,0,0,0
250 DATA 4,7,8,11,13,0,0
260 DATA 4,8,9,12,14,0,0
270 DATA 4,9,10,13,15,0,0
280 DATA 2,10,14,0,0,0,0
290 READ G(,)
300 !
310 K=0
320 M=10000 ! number of trials
330 FOR I=1 TO M
340   P=1 ! start at node 1
350   FOR J=1 TO 4
360      N=G(P,0) ! number of neightbours
370      X=INT(RND*N+1) ! select one randomly 1..N
380      P=G(P,X) ! move to this node
390   NEXT J
400   IF P>=11 THEN K=K+1 ! count if at last row
410 NEXT I
430 DISP K/M

and the result is ... indeed around 0.08 and not 0.11, confirming the value 23/288 given by Valentin and the validity of the published programs above :-)

Now, I think I understood why 16/144 is NOT the probability of ending in the last row: there are indeed 16 paths over 144 that end in the last row, but all paths are not equally probable, we can't just do 16/144.
However, I thought that my little program was worth to be published here.

J-F

[Albert Chan]

(10-14-2022 10:24 AM)J-F Garnier Wrote:   30 ! the nodes are coded from P=1 to 15:
 40 !        1
 50 !       2 3
 60 !      4 5 6
 70 !     7 8 9 10
 80 !   11 12 13 14 15

there are indeed 16 paths over 144 that end in the last row, but all paths are not equally probable,
we can't just do 16/144.

Yes. all paths are not equally probable.
Above triangle, from 1 to bottom row in 4 steps, all steps must go down a row.

In 2 steps:
P(1 to 4) = 1/2*1/4 = 1/8
P(1 to 5) = 1/4

In 2 steps:
P(4 to bottom) = 1/4*1/2 + 1/4*1/3 = 5/24
P(5 to bottom) = 1/3*1/3 = 1/9

In 4 steps:
P(1 to 4 to bottom) = 1/8*5/24 = 5/192 ≈ 2.604% = P(1 to 6 to bottom)
P(1 to 5 to bottom) = 1/4*1/9   = 1/36   ≈ 2.778% , slightly more probable.

  →  P(1 to bottom) = 2*5/192 + 1/36 = 23/288 ≈ 7.986%

Or, we can simply list all paths, then sum the probabilities.

Paths            1/Probability
1 2 4  7 11      2*4*4*4 = 128
1 2 4  7 12      2*4*4*4 = 128
1 2 4  8 12      2*4*4*6 = 192
1 2 4  8 13      2*4*4*6 = 192
1 2 5  8 12      2*4*6*6 = 288
1 2 5  8 13      2*4*6*6 = 288
1 2 5  9 13      2*4*6*6 = 288
1 2 5  9 14      2*4*6*6 = 288
1 3 5  8 12      2*4*6*6 = 288
1 3 5  8 13      2*4*6*6 = 288
1 3 5  9 13      2*4*6*6 = 288
1 3 5  9 14      2*4*6*6 = 288
1 3 6  9 13      2*4*4*6 = 192
1 3 6  9 14      2*4*4*6 = 192
1 3 6 10 14      2*4*4*4 = 128
1 3 6 10 15      2*4*4*4 = 128

P(1 to bottom) = 4/128 + 4/192 + 8/288 = 23/288

[Albert Chan]

(10-13-2022 09:44 PM)C.Ret Wrote:  IF A<>I OR J<>B THEN P(I,J) = P(I,J) + Q(A,B)/W(A,B)

1. Most vertices can go 6 ways, same Q(A,B)/W(A,B) calculated upto 6 times.

It is more efficient to calculate weighted Q first.
Bonus: MAT Q=P line is not needed anymore.

2. we can remove IF THEN statement, NOT(A<>I OR J<>B) ≡ (A==I AND B==J)

3. it is faster to use regular variable, sum it, then assign to P(I,J)

Combined 1,2,3: quoted line is simply T = T+Q(A,B), with T initially set to -Q(I,J)

4. it may be more accurate to go for scaled probability (by 6^S), then unscale it.

With above optimizations, we have:

10 DESTROY ALL @ OPTION BASE 1 @ INPUT "[VA]SRC012A R,S= ";R,S @ T0=TIME
20 REAL T,W(R,R),P(R,R),Q(R,R) @ INTEGER A,B,I,J,K,M
30 FOR I=1 TO R @ FOR J=1 TO I @ W(I,J)=3/(3-(J=1)-(I=J)-(I=R)) @ NEXT J @ NEXT I
40 MAT P=ZER @ P(1,1)=1 @ M=2
50 FOR K=1 TO S
60 FOR I=1 TO M @ FOR J=1 TO I @ Q(I,J)=P(I,J)*W(I,J) @ NEXT J @ NEXT I
80 FOR I=1 TO M
90 FOR J=1 TO CEIL(I/2) @ T=-Q(I,J)
100 FOR A=MAX(1,I-1) TO MIN(I+1,M)
110 FOR B=MAX(1,J-(A<=I)) TO MIN(J+(I<=A),A)
120 T=T+Q(A,B)
130 NEXT B
140 NEXT A
150 P(I,J)=T @ P(I,1+I-J)=T
160 NEXT J
170 NEXT I
180 IF M<R THEN M=M+1
190 DISP K;TIME-T0
200 NEXT K
210 T=0 @ FOR J=1 TO R @ T=T+P(R,J) @ NEXT J
220 DISP TIME-T0;R;S;T/6^S


>RUN
[VA]SRC012A R,S= 5,4
1      .04
2      .08
3      .19
4      .36
.4            5 4         7.98611111111E-2
>T, 6^S                 ! P = 23/288      
103.5      1296

>RUN
[VA]SRC012A R,S= 30,60
...
173.47     30 60      9.51234350207E-6

Compare to original version, speed up = 245/173.47 - 1 ≈ 41%