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

[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