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

[J-F Garnier]

(10-15-2022 06:28 PM)Albert Chan Wrote:  With above optimizations, we have:
[...]
20 REAL T,W(R,R),P(R,R),Q(R,R) @ INTEGER A,B,I,J,K,M

In Series 70 BASIC, there is no speed gain to use INTEGER variables; on the contrary it adds some overhead to convert INTEGER to REAL before evaluating the expression and convert the result back to INTEGER. This also applies to FOR .. NEXT loops and matrix indexes.

J-F

[Albert Chan]

Here is optimized C.Ret version, without "neighbor hunting" 2 inner loops.
I also remove INTEGER declaration, for some speedup (J-F Garnier, thanks for the tip!)

10 DESTROY ALL @ INPUT "[VA]SRC012A R,S= ";R,S @ T0=TIME
20 OPTION BASE 1 @ REAL W(R,R),P(R,R)
30 OPTION BASE 0 @ REAL Q(R+1,CEIL(R/2)+1)
40 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
50 P(1,1)=1 @ M=2
60 FOR K=1 TO S
70 FOR I=1 TO M @ FOR J=1 TO CEIL(I/2)+1 @ 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)
100 P(I,J) = Q(I-1,J-1)+Q(I-1,J) + Q(I,J-1)+Q(I,J+1) + Q(I+1,J)+Q(I+1,J+1)
110 P(I,1+I-J)=P(I,J)
120 NEXT J
130 NEXT I
140 IF M<R THEN M=M+1
150 DISP K;TIME-T0
160 NEXT K
170 T=0 @ FOR J=1 TO R @ T=T+P(R,J) @ NEXT J
180 DISP TIME-T0;R;S;T/6^S


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

From this thread so far, this is the most elegant code, and the fastest!

[Vincent Weber]

I'm very impressed by the HP-71B (which I discovered too late) elegance and capabilities.

And all this without even user-defined functions or recursion. The only advanced concept in use is 2-dimensional arrays, which makes me think that this could be ported to any SHARP pocket above the 1211, the HP-41 with advantage module (or even without, emulating matrices with single arrays), or the HP-42S.

Except that... Memory requirements are huge. The needed registers seem in the range of 2.5*R^2 at bare minimum. For R=30 you need something like 18Kb if each register takes 8 bytes, plus extra variables, plus program stack... So you need something like 20Kb RAM. This rules out the 41, the 42, the SHARP pc-1261... You need a 32K machine basically!

Best regards

[C.Ret]

(10-16-2022 10:25 AM)Vincent Weber Wrote:  Except that... Memory requirements are huge. The needed registers seem in the range of 2.5*R^2 at bare minimum. For R=30 you need something like 18Kb if each register takes 8 bytes, plus extra variables, plus program stack... So you need something like 20Kb RAM.

For the HP-41, I am currently thinking about a version where the among of registers would depend only on the number of steps.

It is based on the idea given by Jean-François and illustrated by Chan for the case (R,S)=(5,4).
However, exploring all the paths of S STEPS is likely to take an infinite amount of time...
Especially where S is greater than R !
On the other hand, this solution requires only S=60 well-used registers (sign, integer part and fractional cleverly exploited) and a few others for counters, coordinates, multiplicative factor and final probability...

So far, my buggy prototypes have only worked on reduced grids and for ridiculously small path’s lengths. I am not very confident in the practical feasibility.
The real & practical solution would be to find a more direct and efficient calculation than listing all the paths in depth.

(10-16-2022 08:32 AM)J-F Garnier Wrote:  In Series 70 BASIC, there is no speed gain to use INTEGER variables; on the contrary it adds some overhead to convert INTEGER to REAL before evaluating the expression and convert the result back to INTEGER. This also applies to FOR .. NEXT loops and matrix indexes.

This was one of the questions I asked myself while writing my code. I didn't take the time to measure the gain or loss of time with all REAL variables.
I imagine that variables of the SHORT type do not offer any gain in speed either.

Thanks to Jean-François who answers and confirms my doubts without me having to time anything.

(10-15-2022 06:28 PM)Albert Chan Wrote:  Compare to original version, speed up = 245/173.47 - 1 ≈ 41%

Thanks to Chan who optimize my code. The time savings are far from trivial or insignificant. I like the new versions

By reading his post, I had already modified the code in my HP-71B in order to make the rounding errors disappear. But I hadn't thought of merging weight and probability in the matrix Q(,) = W(,) .* P(,). Too bad that the MATH module has no instruction making the dot-product. Is there a pretext for a third version of the MATH module?

I love the last version that smashes! thank you Chan!


This made me want to transcribe this algorithm for my HP-28S.

I get for the problem (R,S)=(30,60) exactly P =9.51234350207E-6 after a very long time of 3 hours 12 min 48sec. It must be said that this version has the two internal loops FOR a and FOR b.

Now that Chan gave the idea, I will remake my RPN code to avoid the hunting for neighbors. The next version will therefore use matrices of dimension (R+1)x(R+1), the HP-28S having enough memory (32 kb).

I transcribed this version below. Will it inspired someone or will anyone, as CHAN already does, found some elegant and effective way to improve it?

The equilateral triangle of probabilities P is stored in the form of a one-dimension vector of length \( d = \frac{r^²+r}{2} \). It is not memorized in a register but remains throughout the calculation in the stack.
Similarly, the weight coefficients, which are calculated once at the start of the computation, stay in the stack in the form of a vector of identical length.
On the other hand, the size of the vector Q is reduced (at least at the beginning) and depends on the progress variable m in order to save some effort.

When designing, I intended to calculate vector Q using the DPrdct instruction. But I realized that this instruction is not a native instruction of the HP-28S but is part of the personal programs that I usually use on this machine. So, I had to add the code of this program within the code. I took the opportunity to limit the size of Q and therefore reduce the calculations somewhat. But not the execution time, the HP-28S not being particularly quick to execute my codes.

Similarly, a Tind instruction is used to calculate the index of the triangular position (I,J) with J<=I. This is an instruction that I use elsewhere. It's very short, but I leave it in the code below which makes it 'a much more easy to read'.
As it's RPL, you have to understand 'a little less unreadable'.

SRC12a:
« OVER DUP Tind 2 → r s d m
   « 1 r FOR i
          1 i FOR j
               3 DUP j 1 == - i j == - i r == - /
          NEXT
     NEXT
     { d } →ARRY                                                   %% That's W
     DUP 0 CON 1 1 PUT                                             %% That's P = MAT ZER & P(1,1)=1
     1 s FOR k
          k 1 DISP
          DUP2 m m Tind → W P n                                    %% Embedded version of personal Dot.Product code
             « 1 n FOR q                                           %% size limited to m row since next rows all zero
                    W q GET P q GET *
               NEXT
               { n } →ARRY »                                       %% That's Q = W .* P
          1 m FOR i
               1 i 2 / CEIL FOR j
                    DUP i j Tind GET NEG                           %% Init t =- Q(i,j)
                    1 i 1 - MAX 1 i + m MIN FOR a
                         1 j a i <= - MAX j i a <= + a MIN FOR b
                              OVER a b Tind GET +                  %% Add  t += Q(a,b)
                         NEXT
                    NEXT
                    ROT                                            %% stack W t Q P order
                    i j         Tind 3 PICK PUT
                    i 1 i + j - Tind    ROT PUT
                    SWAP                                           %% back initial stack order
               NEXT
          NEXT
          DROP                                                     %% Drop Q 
          m DUP r < + 'm' STO                                      %% Increase m ( STO+ doesn't work with local variables )
     NEXT
     SWAP 0 CON
     r 1 Tind                                                      %% W = 0 except for last row where W = 1  
     DO
          1 PUTI
     UNTIL 46 FS? END                                              %% Flag 46 set when PUTI reach end of vector 
     DOT 6 s ^ / 
     CLMF 1430 .4 BEEP »                                           %% Restore stack display and bip    
»

 
Tind:                                                              %% Personal Triangular Indice
« OVER SQ ROT - 2 / + »                                            %%   Tind(i,j) = j+(i²-i)/2



Feel free to comment or ask any questions that you deem useful.

Best regards.
C.Ret

[J-F Garnier]

(10-16-2022 09:20 AM)Albert Chan Wrote:  Here is optimized C.Ret version, without "neighbor hunting" 2 inner loops.

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

From this thread so far, this is the most elegant code, and the fastest!

And perfectly exact too, down to the last place! Well maybe also by a little bit of chance.

Also thanks to the slightly optimized memory required for Q(R+1,CEIL(R/2)+1), it fits in a 24k-RAM HP-75 !

Here are the few simple changes to run your code on the HP-75:
7 OPTION BASE 0 @ REAL W(30,30),P(30,30),Q(31,16)
10 INPUT "[VA]SRC012A R,S= "; R,S @ T0=TIME
20 REDIM W(R,R),P(R,R)
30 REDIM Q(R+1,CEIL(R/2)+1)
35 MAT P=ZER @ MAT Q=ZER @ MAT W=ZER ! vars must be initialized
The rest is unchanged.

HP-75 result is: 9.51234350244E-6

Slightly OT: why is the HP-75 result different, and less accurate (I like to investigate these kind of questions) ?
After all, it's all about additions, multiplications and divisions, and the HP-75 and HP-71 share the same algorithms and 12-digit (15 internally) accuracy.

Well, almost the same algorithms.
The Saturn machine algorithms (from the HP-71) introduce a very small improvement, known as the "round-to-even" or "banker's rounding" rule.
When an internal 15-digit result is rounded to the user 12-digit form, and it ends exactly with ...500, then it is rounded to the closest even value, instead of being rounded upward.
It can be demonstrated for instance with 1/(2^18) = 3.81469726562(500) e-6
HP-71: 3.81469726562e-6 (round to even)
HP-75: 3.81469726563e-6 (round up)

This difference occurs in the first steps of the calculation, and introduces a small bias that is enough to bring a different and less accurate answer at the end.

J-F

[Albert Chan]

(10-16-2022 02:10 PM)J-F Garnier Wrote:  HP-75 result is: 9.51234350244E-6
...
The Saturn (from the HP-71) introduces a very small improvement, known as the "round-to-even"

My code of scaling by P by 6^S take account of round-to-even rule.
For machine that round-away-from-zero, I would scale by gcd(2,4,6) = 12 instead.
Or, scale by 1.2 for decimal machine, to keep denominator (1.2^S) from overflow, even with big S

< 40 FOR I=1 TO R @ FOR J=1 TO C @ W(I,J)=3/(3-(J=1)-(I=J)-(I=R)) @ NEXT J @ NEXT I
> 40 FOR I=1 TO R @ FOR J=1 TO C @ W(I,J)=0.6/(3-(J=1)-(I=J)-(I=R)) @ NEXT J @ NEXT I

< 180 DISP TIME-T0;R;S;T/6^S
> 180 DISP TIME-T0;R;S;T/1.2^S

With above changes, on emu71, I get: (don't worry about error of few ULP's)

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

[J-F Garnier]

(10-16-2022 03:07 PM)Albert Chan Wrote:  

(10-16-2022 02:10 PM)J-F Garnier Wrote:  HP-75 result is: 9.51234350244E-6

For machine that round-away-from-zero, I would scale by gcd(2,4,6) = 12 instead.
Or, scale by 1.2 for decimal machine, to keep denominator (1.2^S) from overflow, even with big S

< 40 FOR I=1 TO R @ FOR J=1 TO C @ W(I,J)=3/(3-(J=1)-(I=J)-(I=R)) @ NEXT J @ NEXT I
> 40 FOR I=1 TO R @ FOR J=1 TO C @ W(I,J)=0.6/(3-(J=1)-(I=J)-(I=R)) @ NEXT J @ NEXT I

< 180 DISP TIME-T0;R;S;T/6^S
> 180 DISP TIME-T0;R;S;T/1.2^S

With above changes, on emu71, I get: [..]
9.51234350204E-6

With these changes, I now get on the HP-75 (actually emu75):
9.51234350227E-6
which is indeed better, still not at the level of the best HP-71 performance but closer to my first HP-71 result (9.51234350213E-6).

J-F

[C.Ret]

The latest version for HP-71B proposed by Chan knocks the beast out.

The probability P( r=30 s=60 ) = 9.51234350207E-6 is displayed in just 31’27.4” with the last version.

As explained by Chan, since the same values P(i,j)/W(i,j) were used six times for the vast majority of points; the execution time is divided accordingly now that they are precomputed in the new matrix Q.

Moreover, as Jean-françois pointed out, a significant amount of memory is saved by using only the columns located in the right part of the triangle.

In this new version, the BASIC program is 402 bytes.
For (R,S)=(30,60), it uses 18855 bytes of data (18.4 kB) mainly for the three matrices W, P and Q.

I can't resist the pleasure of sharing this new version with you. Do not hesitate to dissect it:


10 DEF FNL(X)=1+CEIL(X/2)                                  @@ User function for easy column limit indice computation
20 DESTROY ALL
 @ DELAY 0
 @ INPUT"[VA]SRC012a R,S=";R,S @ T0=TIME                        
30 OPTION BASE 1 @ DIM W(R,R),P(R,R)                          @@ MAT W and P ranging ( 1..R   , 1..  R   )  2x7.04 ko
 @ OPTION BASE 0 @ DIM Q(1+R,FNL(R))                          @@ MAT    Q   ranging ( 0..R+1   , 1..1+R/2 )  4.12 ko
 @ P(1,1)=1 @ M=2
40 FOR I=1 TO R
 @   FOR J=1 TO I
 @     W(I,J)=3/(3-(J=1)-(I=J)-(I=R))                         @@ W(I,J) is either 0 or 1 or 1.5 or 3 
 @   NEXT J
 @ NEXT I
50 FOR K=1 TO S
 @   FOR I=1 TO M
 @     FOR J=1 TO FNL(I)
 @       Q(I,J)=P(I,j)*W(I,J)                             @@ no division, spare time, preserve precision
 @     NEXT J
 @   NEXT I
60   FOR I=1 TO M
 @     FOR J=1 TO CEIL(I/2)
70       P(I,J)=Q(I-1,J-1)+Q(I-1,J)+Q(I,J-1)+Q(I,J+1)+Q(1+I,J)+Q(1+I,1+J)   @@ no neighbors hunting, all I±1 or J±1 in Q’s subscripts range 
 @       P(I,1+I-J)=RES                                @@ RES is previous computed arithmetic sum stored in P(I,J)
80     NEXT J
 @   NEXT I
 @   M=M+(M<R)                                       @@ slow but faster than the IF THEN statment which need an extra line!
 @   DISP K;TIME-T0
 @ NEXT K
90 T=0 @ FOR J=1 TO R @ T=T+P(R,J) @ NEXT J
 @ DISP TIME-J;R;S;T/6^S
 @ BEEP

[Albert Chan]

We can cut memory required more than half, by using symmetry.
Removed W, we reduced memory to 2 arrays, P and Q:

Let C = ceil(R/2)

P array elements =     R   *   C
Q array elements = (R+2)*(C+2)

Below code had P and Q with same dimensions, to take advantage of fast MAT Q=P (*)

Total memory (elements) = (R+2)*(C+2) * 2 < (R+3.5)^2


10 DESTROY ALL @ INPUT "[VA]SRC012A R,S= ";R,S @ T0=TIME
20 C=CEIL(R/2) @ OPTION BASE 0 @ REAL P(R+1,C+1),Q(R+1,C+1)
30 P(1,1)=1 @ M=2
40 FOR K=1 TO S
50 MAT Q=P @ Q(1,1)=Q(1,1)*3 ! BUILD Q
60 FOR I=2 TO M-(M=R) @ Q(I,1)=Q(I,1)*1.5 @ J=CEIL(I/2) @ Q(I,J+1)=Q(I,I-J) @ NEXT I
70 IF M<>R THEN 100
80 FOR J=2 TO C @ Q(R,J)=Q(R,J)*1.5 @ NEXT J
90 Q(R,C+1)=Q(R,R-C) @ Q(R,1)=Q(R,1)*3
100 FOR I=1 TO M @ FOR J=1 TO CEIL(I/2) ! BUILD P
110 P(I,J)=Q(I+1,J)+Q(I+1,J+1)+Q(I,J-1)+Q(I,J+1)+Q(I-1,J-1)+Q(I-1,J)
120 NEXT J @ NEXT I
130 IF M<R THEN M=M+1
140 DISP K;TIME-T0
150 NEXT K
160 T=0 @ FOR J=1 TO R-C @ T=T+P(R,J) @ NEXT J @ T=2*T+(2*C-R)*P(R,C)
170 DISP TIME-T0;R;S;T/6^S

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

Code run even faster! 5× speed, against C.Ret orignial code (245 sec)


(*) I don't really need a copy, swapping name of array (P,Q) is enough.
Too bad ... VARSWAP don't work for arrays.

>VARSWAP P, Q
ERR:Data Type

[Valentin Albillo]

  
Hi, all,

Thank you very much for your interest in my SRC #012a, and in particular thanks to Vincent Weber, Fernando del Rey, J-F Garnier, C.Ret, Werner, Albert Chan, PeterP, ijabbot, ttw, rprosperi, rawi, and pier4r for your solutions and/or comments, much appreciated and very fine efforts throughout.

In my OP I mentioned "takes random steps" with the idea of misleading people into thinking that Monte Carlo simulations would be the way to go, but as some of you discovered upon attempting the feat, even a large numer of simulations would get just one correct digit for the toy case, let alone the (30, 60) case, as demonstrated by J-F Garnier's awesome little piece of code.

Getting the 10-12 digits asked for required a fully deterministic method and there are at least two ways to proceed: (a) for each point, identify its neighbors and take the due part from them, or (b) for each location, identify its neighbors and distribute the point's probability among them. It's a zero-sum game and the total probability is the same at the beginning (P=1 for the top location, 0 for all the rest), after every step and thus at the end.

I tried both approaches (a) and (b) above and both worked alright but approach (b), distributing the probability among the neighbors, proved to be faster so that's what my original solution uses. Though tempting and perfectly adequate for the (30, 60) case as originally stated, I don't use symmetry for the reasons given in the fourth Note below, mainly because it would lose generality, which I wanted to preserve so that I could also solve arbitrary starting positions (non-symmetric, one or more as long as their initial peobabilities added up to exactly 1).

So, this is my original solution, a 13-line, 683-byte program. For convenience, it uses two matrix keywords from the Math ROM, namely MAT..ZER and MAT=, easily replaced by trivial loops if the Math ROM's not available):

PROBLEM1 (683 bytes)

10 DESTROY ALL @ OPTION BASE 1 @ INPUT "Rows,Steps=";M,N @ DIM A(M,M),B(M,M) @ SETTIME 0
15 A(1,1)=1 @ W=M-1 @ FOR I=1 TO N @ MAT B=ZER
20 P=A(1,1)/2 @ IF P THEN B(2,1)=B(2,1)+P @ B(2,2)=B(2,2)+P
25 P=A(M,1)/2 @ IF P THEN B(W,1)=B(W,1)+P @ B(M,2)=B(M,2)+P
30 P=A(M,M)/2 @ IF P THEN B(W,W)=B(W,W)+P @ B(M,W)=B(M,W)+P
35 FOR X=2 TO W @ U=X-1 @ V=X+1 @ P=A(X,1)/4 @ Q=A(X,X)/4 @ R=A(M,X)/4
40 IF P THEN B(U,1)=B(U,1)+P @ B(V,1)=B(V,1)+P @ B(X,2)=B(X,2)+P @ B(V,2)=B(V,2)+P
45 IF Q THEN B(U,U)=B(U,U)+Q @ B(V,V)=B(V,V)+Q @ B(X,U)=B(X,U)+Q @ B(V,X)=B(V,X)+Q
50 IF R THEN B(M,U)=B(M,U)+R @ B(M,V)=B(M,V)+R @ B(W,U)=B(W,U)+R @ B(W,X)=B(W,X)+R
55 NEXT X @ FOR X=3 TO W @ U=X-1 @ V=X+1 @ FOR Y=2 TO U @ R=Y-1 @ S=Y+1 @ P=A(X,Y)/6
60 IF P THEN B(U,R)=B(U,R)+P @ B(U,Y)=B(U,Y)+P @ B(X,R)=B(X,R)+P
65 IF P THEN B(X,S)=B(X,S)+P @ B(V,Y)=B(V,Y)+P @ B(V,S)=B(V,S)+P
70 NEXT Y @ NEXT X @ MAT A=B @ NEXT I @ P=0 @ FOR Y=1 TO M @ P=P+A(M,Y) @ NEXT Y @ DISP P;TIME

>STD @ RUN
Rows,Steps= 5,4 [ENDLINE] -> 7.98611111112E-2 (=23/288) .03"

After running the program, we can display the just-built probability matrix (which contains the probability that the man is in each of the 5*(5+1)/2 = 15 points in the grid after he's walked the 4 steps in this toy case, which is small enough that we can get it in exact rational form (requires the JPC ROM) right from the command line:

>FOR I=1 TO 5 @ FOR J=1 TO I @ FRAC$(A(I,J));" "; @ NEXT J @@ NEXT I


                         11/96
                   49/384     49/384
           113/1152     101/576     113/1152
      35/1152       17/288     17/288       35/1152
1/128       23/1152      7/288      23/1152       1/128

Now for the big (30, 60) case:

>RUN
Rows,Steps= 30,60 [ENDLINE] -> 9.51234350207E-6   27.46"

We can now use the probability matrix right from the command line to answer all kinds of questions, e.g.:

(a) check the matrix correctness by adding up the probabilities for all points:

>S=0 @ FOR I=1 TO M @ FOR J=1 TO I @ S=S+A(I,J) @ NEXT J @ NEXT I @ S

      1.00000000002

which gives 1 as it should, after all the man must be in one of the 465 points.

(b) probability that he ends in the left border:

>L=0 @ FOR I=1 TO M @ L=L+A(I,1) @ NEXT I @ L

      .117372130231

which is 12338.9289091 times more probable than ending in the bottom row.

(c) probability that he ends in any border as compared to ending inside the grid:

>D=0 @ FOR I=1 TO M @ D=D+A(M,I) @ NEXT I @ D

      9.51234350207E-6 (down edge = bottom row)

>R=L @ L+R+D

      .234753772806

so he ends in a border 23.48% of the time and thus 76.52% of the time he ends inside the grid, so he's 3.26 times more likely to end inside the grid than in a border.

(d) probability that he ends in each of the corners or in any of them:

>A(1,1);A(M,1);A(M,M) @ A(1,1)+A(M,1)+A(M,M)

      8.30321536116E-3   2.21000524869E-8   2.21000524869E-8 (symmetric)

      8.30325956126E-3

(e) probability that he ends up in any of the first N rows:

>FOR N=1 TO M @ P=0 @ FOR I=1 TO N @ FOR J=1 TO I @ P=P+A(I,J) @ NEXT J @ NEXT I @ N;P @ NEXT N


 N      P
---------------------------------------------
 1   0.00830      top corner
 2   0.04100      top two rows
       ...
 5   0.25643      top 5 rows
       ...
 7   0.44905      top 7 rows, less than 50%
 8   0.54414      top 8 rows, more than 50%
       ...
10   0.71182      top 10 rows, 71%
15   0.94283      top 15 rows, 94%
20   0.99430      top 20 rows, 99%
25   0.99972      top 25 rows, almost sure
29   0.99999      top 29 rows, 99.999% certain
30   1.00000      all 30 rows, 100% certain.

Notes:

• All runtimes are for a virtual HP-71B (go71b) running on a mid-range Samsung Galaxy Tab A 6 tablet (Android). A physical HP-71B should take ~58 min for the (30, 60) case.
    
• The (30, 60) case needs 14,533 bytes of RAM to run. If the matrices are dimensioned as SHORT then it requires only 8,245 bytes of RAM but it runs slower and gives only 5 correct digits save 2 ulp (P=0.0000095121), instead of 12 correct digits.
    
• The program checks whether a location's probability is currently zero, in which case it skips it when distributing the probability among its neighbors. This saves about 20% running time for the (30,60) case.
    
• Another pretty obvious optimization would be to take advantage of the symmetry, as the man's starting location is at the top corner; this would cut running time and required RAM roughly in half but it loses generality, i.e.: arbitrary non-symmetrically placed starting positions (one or more) would not run at all, and it would complicate the coding somewhat; also, for a one-time program which already runs suitably fast it's really a moot point.
    
• It's a real pity that the HP-71B BASIC dialect doesn't include the += , etc., operators, as it would make my program that much shorter and faster, i.e. line
  
  
  50 IF R THEN B(M,U)=B(M,U)+R @ B(M,V)=B(M,V)+R @ B(W,U)=B(W,U)+R @ B(W,X)=B(W,X)+R
  
  would become
  
  
  50 IF R THEN B(M,U)+=R @ B(M,V)+=R @ B(W,U)+=R @ B(W,X)+=R
  
  and so on. Perhaps J-F Garnier could do something about it.
    
• The exact value is (465 points updated per step * 60 steps = ~ 200,000 arithmetic operations in all):
  
  
     3722200777884626618385530906788866022689096963173522895529  
  391302183008102676318141068027642364938466415279908207464546304
  
  = 9.5123435020743320973531347375383316350767310071737e−6 (50 digits)

Frankly, I thought that some accomplished RPL programmers would use the multiprecision capabilities of their advanced RPL models to compute the exact rational solution above but no such luck ...


Thanks to all of you for your interest, solutions and comments. A number of you produced correct solutions for this Problem 1 but it's the easiest of the lot and I wonder how many of you will achieve a perfect 6 for 6 score by the time all six parts of SRC #012 are over ...

See you in Problem 2 featuring soon.
V.

[Albert Chan]

(10-16-2022 05:52 PM)Valentin Albillo Wrote:  Another pretty obvious optimization would be to take advantage of the symmetry, as the man's starting location is at the top corner; this would cut running time and required RAM roughly in half but it loses generality, i.e.: arbitrary non-symmetrically placed starting positions (one or more) would not run at all, and it would complicate the coding somewhat; also, for a one-time program which already runs suitably fast it's really a moot point.

Here is my code, removed symmetry for generality.

For other starting positions, edit LINE 30, with M = first empty row.
Tips: we can *still* use symmetry, by rotation for the smallest starting M

10 DESTROY ALL @ INPUT "[VA]SRC012A R,S= ";R,S @ SETTIME 0
20 OPTION BASE 0 @ T=R+1 @ REAL P(T,T),Q(T,T)
30 P(1,1)=1 @ M=2
40 FOR K=1 TO S
50 MAT Q=P @ Q(1,1)=Q(1,1)*3 ! BUILD Q
60 FOR I=2 TO M-(M=R) @ Q(I,1)=Q(I,1)*1.5 @ Q(I,I)=Q(I,I)*1.5 @ NEXT I
70 IF M<>R THEN 100
80 FOR J=2 TO R-1 @ Q(R,J)=Q(R,J)*1.5 @ NEXT J
90 Q(R,1)=Q(R,1)*3 @ Q(R,R)=Q(R,R)*3
100 FOR I=1 TO M @ FOR J=1 TO I ! BUILD P
110 P(I,J)=Q(I+1,J)+Q(I+1,J+1)+Q(I,J-1)+Q(I,J+1)+Q(I-1,J-1)+Q(I-1,J)
120 NEXT J @ NEXT I
130 M=M+(M<R)
140 DISP K;TIME
150 NEXT K
160 T=0 @ FOR J=1 TO R @ T=T+P(R,J) @ NEXT J
170 DISP TIME;R;S;T/6^S


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

As expected, without symmetry, speed almost cut in half.

[Xorand]

(10-16-2022 10:25 AM)Vincent Weber Wrote:  The only advanced concept in use is 2-dimensional arrays, which makes me think that this could be ported to any SHARP pocket above the 1211...

Except that... Memory requirements are huge. The needed registers seem in the range of 2.5*R^2 at bare minimum. For R=30 you need something like 18Kb if each register takes 8 bytes, plus extra variables, plus program stack... So you need something like 20Kb RAM. This rules out the 41, the 42, the SHARP pc-1261... You need a 32K machine basically!

I ported Fernando del Rey's BASIC program to my Sharp PC-1500 with 8K memory module (10k total storage). It ran the size 5, step 6 case easily. I ran a 10/20 case (took just a couple minutes - didn't time it). When I try the suggested size 30, 60 steps case, the computer unsurprisingly gives an out of memory error. Trial and error shows that I could go up to a grid size of 23.

[Vincent Weber]

(10-17-2022 01:15 AM)Xorand Wrote:  

(10-16-2022 10:25 AM)Vincent Weber Wrote:  The only advanced concept in use is 2-dimensional arrays, which makes me think that this could be ported to any SHARP pocket above the 1211...

Except that... Memory requirements are huge. The needed registers seem in the range of 2.5*R^2 at bare minimum. For R=30 you need something like 18Kb if each register takes 8 bytes, plus extra variables, plus program stack... So you need something like 20Kb RAM. This rules out the 41, the 42, the SHARP pc-1261... You need a 32K machine basically!

I ported Fernando del Rey's BASIC program to my Sharp PC-1500 with 8K memory module (10k total storage). It ran the size 5, step 6 case easily. I ran a 10/20 case (took just a couple minutes - didn't time it). When I try the suggested size 30, 60 steps case, the computer unsurprisingly gives an out of memory error. Trial and error shows that I could go up to a grid size of 23.

Excellent, thanks !

[Dave Britten]

(10-16-2022 05:52 PM)Valentin Albillo Wrote:  So, this is my original solution, a 13-line, 683-byte program. For convenience, it uses two matrix keywords from the Math ROM, namely MAT..ZER and MAT=, easily replaced by trivial loops if the Math ROM's not available):[font=Courier]

Darn, I was hoping your mention of "which is a short program capable of quickly solving the generic problem" meant you had a clever multinomial approach that could solve the problem in mere seconds on an HP 67.

Still, it's a nice elegant solution, even if matrix iteration is the only feasible way to do it. I have a direct adaptation of your program running on my Sharp PC-1403H now. We'll see how long it takes for the large case - it's been going for about 40 minutes...

[PeterP]

Here is my code for the HP41, which does not have enough memory for the matrix approach. It only works for S = N-1 and uses that one only has to worry about the paths that increase R in every step, as all other paths do not reach the last row in time. It then uses that the number of ways of reaching the last row is the binomial tree (from each node there are two ways to reach the next row, with the exception of the edges).
As such, its rather trivial.

Inspired by Albert and others I was then trying to figure out if there is a way to scale the result for S=N-1 to S=N, S=N+1, etc given that it trends asymptotically to a fixed value. However, I was not able to do so and now it has closed, so apologies for posting my code only now.

The code uses that only steps that increase the row can work. It also uses the symmetry between left and right.

When walking down the left edge, starting at some point in row y, the first step then always has a probability of 1/4 to go towards the middle. Afterwards all probabilities to go downwards are 2 * 1/6 (one to the left and one to the right) or 1/3. This goes on for the remaining rows x to the last row.
The first y = R-3, with R being the total number of rows. (The step from row 1 to 2 is just 1/2 to go to the left, then you have 1/4 to go right down towards the middle, row 3, from which the 1/6 probabilities to go down left and right to row 4 start)
x = (R-3) - y
For each of these possible vertices the probabilities reach final row are (1/4)^x * (1/3)^y
At the end we have to take care of the symmetry and the first step.

STO 00......R-3
STO 01......y
STO 02...... probability of ending in last row
STO 03...... 1/3 (for speed, number entry in 41 is very slow)
STO 04...... 1/4 (for speed, number entry in 41 is very slow)

Enter into X the number of rows, press R/S

LBL "SRC12A"
3
-
STO 00 ! R-3
STO 01 ! y
LastX
1/x
STO 03 ! (1/3)
X<>Y
y^x
STO 02 ! prob P
4
1/x
STO 04 ! (1/4)
LBL 00 ! loop over y = (R-3) to 0
RCL 03
RCL 01
DECx
x<0?
GTO 01 ! -> we are done
STO 01
y^x
RCL 04
RCL 00
RCL 01
-
y^x
*
ST+ 02
GTO 00
LBL 01
RCL 04
RCL 00
y^x
ST+ 02
RCL 04
ST* 02
View 02
STOP
end

[Albert Chan]

(10-16-2022 10:25 AM)Vincent Weber Wrote:  Memory requirements are huge. The needed registers seem in the range of 2.5*R^2 at bare minimum. For R=30 you need something like 18Kb if each register takes 8 bytes, plus extra variables, plus program stack... So you need something like 20Kb RAM. This rules out the 41, the 42, the SHARP pc-1261... You need a 32K machine basically!

I was thinking of minimum memory requirement, without code of checking corners, edges.
Also, code that would work with non-symmetrical starting conditions. (no symmetry trick)

In other words, this would still hold, collecting probabilities from hexagon vertices.
Q = weighted probability of previous P

P(I,J) = Q(I-1,J-1) + Q(I-1,J) + Q(I,J-1) + Q(I,J+1) + Q(I+1,J) + Q(I+1,J+1)

Current implementations treated (P, Q) as square matrix, with lots of 0's on the right.
But, if we flatten it, say into (p,q) of 1 dimension, all we need is 1 zero between rows.

Q(I,J) ≡ q(I*(I+1)/2 + J)

Q(1,1) = q(2)
Q(2,1), Q(2,2) = q(4), q(5)
Q(3,1), Q(3,2), Q(3,3) = q(7), q(8), q(9)
Q(4,1), Q(4,2), Q(4,3), Q(4,4) = q(11), q(12), q(13), q(14)
...

Note that q of triangular numbers are outside the triangle, with probability of 0.

Example, to get next iteration of P(3,2) = p(3*4/2+2 = 8)

p(8) = q(4) + q(5) + q(7) + q(9) + q(12) + q(13)

In general, for P(I,J) = p(x = I*(I+1)/2 + J)

p(x) = q(x-I-1) + q(x-I) + q(x-1) + q(x+1) + q(x+I+1) + q(x+I+2)

---

Q(0,0) .. Q(R+1,R+1) ≡ q(0) .. q((R+1)*(R+4)/2)

Total q elements = (R+1)*(R+4)/2 + 1= (R+2)*(R+3)/2

Technically p does not need as many spaces.
But for convenience, we like to keep both arrays with same dimensions.

Total array elements needed = (R+2)*(R+3) = floor((R+2.5)^2)

With theory out of the way, here is the flattened array version.

10 DESTROY ALL @ INPUT "[VA]SRC012A R,S= ";R,S @ SETTIME 0
20 OPTION BASE 0 @ T=(R+1)*(R+4)/2 @ REAL P(T),Q(T)
30 P(2)=1 @ M=2
40 FOR K=1 TO S
50 MAT Q=P @ Q(2)=Q(2)*3 @ T=3 ! BUILD Q
60 FOR I=2 TO M-(M=R) @ X=T+1 @ Q(X)=Q(X)*1.5 @ X=T+I @ Q(X)=Q(X)*1.5 @ T=X+1 @ NEXT I
70 IF M<>R THEN 100
80 FOR X=T+2 TO T+R-1 @ Q(X)=Q(X)*1.5 @ NEXT X
90 T=T+1 @ Q(T)=Q(T)*3 @ Q(X)=Q(X)*3
100 T=1 @ FOR I=1 TO M @ FOR X=T+1 TO T+I ! BUILD P
110 P(X)=Q(X-I-1)+Q(X-I)+Q(X-1)+Q(X+1)+Q(X+I+1)+Q(X+I+2)
120 NEXT X @ T=X @ NEXT I
130 M=M+(M<R)
140 DISP K;TIME
150 NEXT K
160 T=0 @ K=R*(R+1)/2 @ FOR X=K+1 TO K+R @ T=T+P(X) @ NEXT X
170 DISP TIME;R;S;T/6^S

On line 90, I use the quirk the last X is outside loop limits. Q(X) == Q(T+R)
Same idea on line 120, "T=X" same as "T=T+I+1", next triangular number

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

>RUN
[VA]SRC012A R,S= 5,4
...
.31         5 4         7.98611111111E-2

>MAT P = (6^(-S)) * P ! scaled to have sum(P) = 1.0
>MAT DISP P ! note the single zero gap, between "row"
0
0
.114583333333
0
.127604166667
.127604166667
0
9.80902777778E-2
.175347222222
9.80902777778E-2
0
0
3.03819444445E-2
5.90277777778E-2
5.90277777778E-2
3.03819444445E-2
0
.0078125
1.99652777778E-2
2.43055555556E-2
1.99652777778E-2
.0078125
0
0
0
0
0
0
0

[Dave Britten]

Sharp PC-1403H version of Valentin's 71B program (I didn't really have to change much, aside from minor syntactic/keyword differences):

10 CLEAR:INPUT "ROWS?";M:INPUT "STEPS?";N:DIM A(M,M):DIM B(M,M)
15 A(1,1)=1:W=M-1:FOR I=1 TO N:GOSUB 100
20 P=A(1,1)/2:IF P LET B(2,1)=B(2,1)+P:B(2,2)=B(2,2)+P
25 P=A(M,1)/2:IF P LET B(W,1)=B(W,1)+P:B(M,2)=B(M,2)+P
30 P=A(M,M)/2:IF P LET B(W,W)=B(W,W)+P:B(M,W)=B(M,W)+P
35 FOR X=2 TO W:U=X-1:V=X+1:P=A(X,1)/4:Q=A(X,X)/4:R=A(M,X)/4
40 IF P LET B(U,1)=B(U,1)+P:B(V,1)=B(V,1)+P:B(X,2)=B(X,2)+P:B(V,2)=B(V,2)+P
45 IF Q LET B(U,U)=B(U,U)+Q:B(V,V)=B(V,V)+Q:B(X,U)=B(X,U)+Q:B(V,X)=B(V,X)+Q
50 IF R LET B(M,U)=B(M,U)+R:B(M,V)=B(M,V)+R:B(W,U)=B(W,U)+R:B(W,X)=B(W,X)+R
55 NEXT X:FOR X=3 TO W:U=X-1:V=X+1:FOR Y=2 TO U:R=Y-1:S=Y+1:P=A(X,Y)/6
60 IF P LET B(U,R)=B(U,R)+P:B(U,Y)=B(U,Y)+P:B(X,R)=B(X,R)+P
65 IF P LET B(X,S)=B(X,S)+P:B(V,Y)=B(V,Y)+P:B(V,S)=B(V,S)+P
70 NEXT Y:NEXT X:GOSUB 200:NEXT I:P=0:FOR Y=1 TO M:P=P+A(M,Y):NEXT Y:BEEP 3: PRINT P:END
100 FOR J=1 TO M:FOR K=1 TO J:B(J,K)=0:NEXT K:NEXT J:RETURN
200 FOR J=1 TO M:FOR K=1 TO J:A(J,K)=B(J,K):NEXT K:NEXT J:RETURN
300 "A"S=0:FOR I=1 TO M:FOR J=1 TO I:S=S+A(I,J):NEXT J:NEXT I:PRINT S:END

This model takes about 4 hours to run for the 30, 60 case, and uses about 15 KB RAM for the program + data. After the program finishes, you can press DEF A to calculate the total probability for the whole map (should be extremely close to 1). One could easily tack on some more lines with user-key labels for performing other calculations on the finished matrix (probability of ending in the starting position, on any edge, at specific coordinates, etc.).

I believe this model has some machine-code matrix routines that can be invoked from within BASIC. I'll have to see if those can be used to speed up any of the matrix copying.

[Valentin Albillo]

Dave Britten Wrote:Darn, I was hoping your mention of "which is a short program capable of quickly solving the generic problem" meant you had a clever multinomial approach that could solve the problem in mere seconds on an HP 67.

Sarcastic much, are we ?

Sorry to disillusion you, but quoting myself:

Valentin Albillo Wrote:[...] let's not overhype my abilities lest disillusionment ensues, ok ?

Too bad you didn't read it.

Dave Britten Wrote:Sharp PC-1403H version of Valentin's 71B program (I didn't really have to change much, aside from minor syntactic/keyword differences) [...] This model takes about 4 hours to run for the 30, 60 case [...]

Thanks for taking the trouble to key in my program into the PC-1403H, appreciated. That model is slow as molasses and the fact that you had to use FOR..NEXT loops to clear and copy matrices at each of the 60 steps slows down the program even further.

I would be curious to know how long does it take Mr. Chan's fastest non-symmetric (general) program to run the (30, 60) case in: (a) The 1403H, and (b) a physical HP-71B, if you'd obligue.

V.

[Xorand]

Fernando del Rey's program on my Sharp PC-1500 yielded the following various results. Originally, it would support a size of 23, but by adding a couple statements to grab the start and end times it dropped to 22 max.

Code:


Size  Step  Result            Elapsed
5     6     1.529224538E-01   00:00:28
10    12    4.841708817E-03   00:03:32
10    15    1.376503457E-02   00:04:13
10    20    3.34319009E-02    00:05:34
15    20    3.3016796293E-04  00:13:03
15    30    4.414569829E-03   00:19:33
20    30    4.948015608E-05   00:35:20
20    40    5.702504313E-04   00:47:04
22    44    2.51225541E-04    01:02:59

[Vincent Weber]

Hi Valentin,

I'm surprised to read that the PC-1403H is slow.. Is the PC-1475 faster ? Or the PC-1350/60 ?

Best regards,

Vincent