[VA] SRC #012a - Then and Now: Probability
|
10-21-2022, 05:52 PM
Post: #81
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
Here, we do P(R, S=R-1) by hand, to discover its patterns.
Even rows are previous row probabilities, scaling to 1 way, edge /4, inside /6. 1/2 1/2 --> P(2,1) = 1 1/8 1/8 1/8 1/4 1/8 --> P(3,2) = 1/2 1/32 1/24 1/32 1/32 7/96 7/96 1/32 --> P(4,3) = 5/24 1/128 7/576 7/576 1/128 1/128 23/1152 7/288 23/1152 1/128 --> P(5,4) = 23/288 We can also get P(5,4), directly from scaled P(4,3): P(5,4) = 2*(1/128 + 7/576 + 7/576 + 1/128) = 4*(1/128 + 7/576) = 23/288 (10-21-2022 04:10 PM)Albert Chan Wrote: P(R, S=R-1) = 3^(3-R) - 2^(5-2*R) We are now ready to proof above, by induction Assume formula is correct, we split it to two types, to get P(R+1, S=R): P(R, S=R-1) = P(edges) + P(inside) 3^(3-R) - 2^(5-2*R) = 2^(4-2*R) + (3^(-R+3) - 3*2^(4-2*R)) P(R+1, S=R) = 2 * sum(scaled to 1 way of P(R, S=R-1) = 2 * (2^(4-2*R)/4 + (3^(-R+3) - 3*2^(4-2*R))/6) = 9*3^(-R) - 8*2^(-2*R) = 3^(3-(R+1)) - 2^(5-2*(R+1)) QED |
|||
10-21-2022, 10:00 PM
Post: #82
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
(10-21-2022 04:10 PM)Albert Chan Wrote:(10-12-2022 12:10 PM)PeterP Wrote: 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). Very neat Albert! It converts the summation into a formula for the sum as its a geometric progression (which I did not recognize). You clearly did not need a computer and could have proven the result to be correct on an airplane with just a simple calculator, your pen and pencil :-) Thank you for sharing. Cheers, PeterP |
|||
10-27-2022, 02:07 PM
Post: #83
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
(10-05-2022 08:38 PM)Valentin Albillo Wrote: using EXCLUSIVELY VINTAGE HP CALCULATORS (physical or virtual,) coding in either RPN, RPL or HP-71B language AND NOTHING ELSE Well, to make Valentin's wishes come true, here's an entry that will solve the 30/60 problem on a real 42S, the only vintage RPN calculator able to do it. To make it fit the 42S' memory, I have taken Albert Chan's flattened code to the extreme: you don't need a full P *and* Q, they can largely overlap, all you need is an extra buffer row at the end. The memory requirements are then (R+4)/2 x (R+3), and I define REGS as such. When P is calculated it is shifted down a full row with regard to Q, in rows 2..(R+4)/2, and we move it one row up by deleting the first row and adding a new empty row at the end (which, incidentally, you can't do with INSR). here's the code. Not much time has been spent in trying to improve it, just to make it work ;-) Estimate of real 42S running time: 3h05m I use VARMENU "TRW" to set R and S, EXIT the menu and do XEQ "TRW" 00 { 325-Byte Prgm } 01▸LBL "TRW" 02 MVAR "R" 03 MVAR "S" 04 4 05 RCL+ "R" 06 2 07 STO "M" 08 ÷ 09 3 10 RCL+ "R" 11 CLV "REGS" 12 DIM "REGS" 13 1 14 STO 02 15 RCL "S" 16 STO "K" 17 EDITN "REGS" 18 GROW 19▸LBL 20 @ --------------------------------- @ P-> Q, adjust corners and edges @ --------------------------------- 20 3 21 STO× 02 @ top 22 RCL "M" 23 RCL "R" 24 X=Y? 25 DSE ST Y 26 SIGN 27 - 28 1ᴇ3 29 STO+ ST Y 30 ÷ @ I=1..M-1-(M=R) 31 2 32▸LBL 02 @ left and right edges 33 2 34 + 35 RCL+ ST Y 36 1.5 37 STO× IND ST Y 38 STO× IND ST L 39 R↓ 40 IP 41 ISG ST Y 42 GTO 02 43 RCL "M" 44 RCL "R" 45 X>Y? 46 GTO 00 47 RCL ST Z 48 ENTER 49 ENTER 50 RCL+ "R" 51 1ᴇ3 52 ÷ 53 + 54 3 55 + 56 1.5 57▸LBL 03 @ bottom edge 58 STO× IND ST Y 59 ISG ST Y 60 GTO 03 61 R^ 62 2 63 + 64 3 65 STO× IND ST Y 66 STO× IND ST T 67▸LBL 00 @ --------------------------------------------------------- @ Q->P @ P(X) := Q(X-1)+Q(X+1)+Q(X-I-1)+Q(X-I)+Q(X+I+1)+Q(X+I+2) @ and P(X) is just Q(X+R+3) @ --------------------------------------------------------- @ find I,J of P(M,M) in the (R+4)/2 x (R+3) matrix @ qmm = Reg(M*(M+1)/2 + M) @ pmm = Reg(qmm + R+3) @ J = pmm MOD (R+3) + 1 @ I = (pmm + 1 - J)/(R+3) + 1 68 RCL "M" 69 STO "I" 70 ENTER 71 XEQ 99 @ qmm 72 RCL ST X 73 3 74 RCL+ "R" 75 + 76 RCL ST X 77 LASTX @ R+3 pmm+1 pmm+1 qmm 78 MOD 79 STO- ST Y 80 X<>Y 81 LASTX 82 STO+ ST Y 83 ÷ 84 X<>Y 85 1 86 + 87 STOIJ 88 R^ 89 RCL- "M" 90 LASTX 91 2 92 + 93 RCL+ "M" 94 LASTX 95▸LBL 04 96 RCL "I" 97 STO "J" 98 DSE ST Y 99▸LBL 05 100 CLX 101 RCL IND ST T 102 RCL+ IND ST Z 103 RCL+ IND ST Y 104 DSE ST T 105 DSE ST Z 106 DSE ST Y 107 DSE ST Y 108 RCL+ IND ST T 109 RCL+ IND ST Z 110 RCL+ IND ST Y 111 ISG ST Y 112▸LBL 00 113 ← 114 DSE "J" 115 GTO 05 116 DSE ST Z 117 DSE ST Z 118 CLX 119 ← 120 R↓ 121 DSE "I" 122 GTO 04 123 I- 124 DELR @ we are at 1,1 now 125 CLX 126 ← 127 → @ GROW mode causes an extra row now 128 RCL "M" 129 RCL "R" 130 X>Y? 131 ISG "M" 132▸LBL 00 133 DSE "K" 134 GTO 20 135 RCLEL 136 EXITALL 137 RCL "R" 138 ENTER 139 ENTER 140 XEQ 99 141 0 142▸LBL 06 143 RCL+ IND ST Y 144 DSE ST Y 145 DSE ST Z 146 GTO 06 147 6 148 RCL "S" 149 Y^X 150 ÷ 151 RTN 152▸LBL 99 153 ENTER 154 X^2 155 + 156 2 157 ÷ 158 + 159 END Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
11-02-2022, 10:56 PM
Post: #84
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
Hi, all, After 4 weeks to the day, it seems this SRC #012a has run its course, so time for a few additional comments and a few new results. The aditional comments Gjermund Skailand Wrote:This has been a very interesting thread. This is a sys-rpl version for HP50g of CReth SRC12a. On an actual HP50g the calculation time for the 30 60 problem is 5 min 21sec. p= 9.51234350207E-6 Thank you for your appreciation. Your SysRPL version looks amazing, kinda assembly language, producing the correct 12-digit result at least 10x faster than a physical HP-71B, which is truly awesome. Can someone please confirm that the listing is correct and will produce the stated result in the stated time ? Albert Chan Wrote:For S = R-1, we can treat triangle as without bottom edge (and the 2 corners) [...] It simplified to: Very nice exact symbolic result for that particular case, Albert Chan, congratulations ! Normally this could be construed as going against my stated rules but as you previously posted tons of actual HP-71B code, I'm not complaining. Werner Wrote:Well, to make Valentin's wishes come true, here's an entry that will solve the 30/60 problem on a real 42S, the only vintage RPN calculator able to do it. [...] Estimate of real 42S running time: 3h05m. Thank you very much, Werner, for taking my wishes into consideration and producing such a fine HP-42S solution. Your running time estimation on a physical HP-42S seems to be 3x slower than my solution running on a physical HP-71B but, as you say, optimization would probably reduce the timing considerably and some compromises had to be made to fit it into the available RAM. Nevertheless, running your program in Free42 on my Samsung tablet takes just 4.5 seconds, while still within the rules. The new results As I've stated oftentimes, one of my main goals is to get people who like vintage HP calculators to not consider them as obsolete gadgets only fit for collecting or nostalgia, with no real place in the real world, but as still useful devices which can indeed be used to solve modern problems and best of all, to improve one's sleuthing and programming abilities while attempting the solution, in the way of "Experimental Mathematics" (EM for short); quoting from Wikipedia:
--------------------------------------------------- 5 4 7.986111e-2 20 20 5.006369e-8 40 2.666558e-1 200 5.204890e-2 50 2.666660e-1 1000 6.666566e-2 60 2.666666e-1 2000 6.666667e-2 70 2.666667e-1 10 10 1.317953e-3 30 30 1.289121e-12 100 1.317596e-1 60 9.512344e-6 1000 1.333333e-1 120 1.694782e-3 240 1.531651e-2 300 2.225664e-2 480 3.539453e-2 960 4.368177e-2 P_lim(N) = 4/(3*N) Now we can check additional cases (say N = 7, 13, 22, ... rows) to see if the conjectured formula holds, and if it does we can attempt to find a symbolical proof for it, A. Chan-style ! So far this applies to the probability of being in the bottom row at the end of the walk, but what about the probability in the limit of being in a particular, single location as the number of steps grows indefinitely ? Adding this line to my program will display the resulting probability matrix which gives the probability for each and every grid point, using the FRAC$ keyword from the JPC ROM to output exact rational results:
Rows,Steps=5,100 -> 2.666667e-1
1/15 1/15 1/15 1/10 1/15 1/15 1/10 1/10 1/15 1/30 1/15 1/15 1/15 1/30
Now you may be wondering if there's some way to automatically get the exact probability matrix in the limit for a given number N of rows (for reasonable N and running times,) and indeed there is a simple procedure we might try out. First create a copy of my program and edit these five lines to be as follows:
15 MAT A=(2/(M*(M+1))) @ W=M-1 @ K=1E-6 @ FOR I=1 TO INF @ MAT B=ZER 70 NEXT Y @ NEXT X @ MAT A=A-B @ DIM A(M*M) @ DISP I;CNORM(A) @ IF RES<K THEN 80 75 DIM A(M,M) @ MAT A=B @ NEXT I 80 FOR I=1 TO M @ FOR J=1 TO I @ DISP FRAC$(B(I,J),6);" "; @ NEXT J @ DISP @ NEXT I Once the initialization is over the program then computes the probability matrix for steps 1, 2, 3, ...., comparing each matrix with the previous one. When the difference is less than a hardcoded tolerance (K=1E-6 at line 15) the process is over and the limit probability matrix is output in exact rational form.
Step Difference ---------------------- 1.000000 0.991667 2.000000 0.043056 3.000000 0.020833 4.000000 0.012770 ... 128.000000 0.000001 129.000000 0.000001 Limit Probability matrix: 1/315 2/315 2/315 2/315 1/105 2/135 ... 2/315 1/105 ... 1/105 2/315 1/315 2/315 ... ... ... 2/315 1/315 Now that we have a working program, we can create a simplified version which just checks the difference of the probability for only the top corner location (1,1) at successive steps, simply comparing A(1,1) vs. B(1,1) instead of considering the whole matrices, and once the tolerance is met we simply output B(1,1), the top corner's probability. If we run this simpler and faster program for various numbers of rows, we get:
------------------------------------------------------------------------------ P(1,1) 1/30 1/63 1/108 1/135 1/165 1/198 1/234 1/273 1/315
Well, I hope this has provided a good example of how you can use your vintage HP calc to do some sleuthing and get nice symbolic results in the spirit of Experimental Mathematics. Once you get the numeric results, conjecturing the symbolic ones and afterwards attempting to prove them is that much easier. Regards. V. All My Articles & other Materials here: Valentin Albillo's HP Collection |
|||
11-03-2022, 11:19 AM
Post: #85
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
(11-02-2022 10:56 PM)Valentin Albillo Wrote: Running the program for 5 rows and a sufficiently large number of steps (S=100), we get in SCI 6: (10-27-2022 02:07 PM)Werner Wrote: @ P-> Q, adjust corners and edges Asymptotic P ratios can be easily derived from the code. If Q entries are all the same, nextP = P Say, all entries in Q = k (for above example, k = 1/60), Q -> P: P(corners) : P(edges) : P(inner) = 2k : 4k : 6k = 1 : 2 : 3 |
|||
11-04-2022, 03:31 PM
(This post was last modified: 11-04-2022 03:38 PM by DavidM.)
Post: #86
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
(11-02-2022 10:56 PM)Valentin Albillo Wrote:Gjermund Skailand Wrote:This has been a very interesting thread. This is a sys-rpl version for HP50g of CReth SRC12a. On an actual HP50g the calculation time for the 30 60 problem is 5 min 21sec. p= 9.51234350207E-6 There were a few typos in Gjermund's posted SysRPL code. I've made some corrections to his version below to produce a program which appears to give the proper results, so hopefully I haven't introduced any new bugs in the process. I can confirm that his code produces the correct real result for 30-60 input in about 317 seconds on my physical 50g. His version is really a hybrid mix of SysRPL, Saturn+ (the plus is important here), and even embedded UserRPL. The main "get" and "put" commands for arrays in SysRPL require that the index given for the element in question is expressed as a single binary integer. So "x,y" must be converted to a single integer representing the linear index of the element as if the array were actually a vector. In particular, PULLREALEL and PUTREALEL are the commands which do the dirty work of recalling and storing the array elements. Given that this program (along with the others, of course) depends very heavily on array indexing, anything that can speed up the conversion of x-y coordinates to a single vector index will have a substantial impact on the runtime. Gjermund has a single subroutine for this conversion ("Tind"), which I suspect he originally wrote in standard SysRPL. Given that it is a critical routine, re-coding that in Saturn assembly made good sense. Not only is it running at Saturn speed, though, he also used a "Saturn+" opcode for squaring a 5-nibble integer. This allows squaring the value in 1 assembly step instead of having to call a subroutine for that purpose. That opcode is only available on the ARM-based RPL calcs, but speeds up the conversion of x,y coordinates into a single index noticeably. Combine the above with the usual speed increase of SysRPL over UserRPL and you've got a significantly faster program than the original UserRPL version that C.Ret posted. Here's Gjermund's program with the typos corrected: !RPL !NO CODE !JAZZ :: CK2NOLASTWD CK&DISPATCH2 #11 :: COERCE2 ' CODE GOSBVL POP2# RSTK=C SAVE B=A.A A*A.A A-B.A ASRB.A C=RSTK A+C.A GOSBVL PUSH#ALOOP ENDCODE 3PICK DUP 3PICK EVAL SWAP 2 SWAPOVER {{ r s d m Tind ii }} r #1+_ONE_DO (i) * INDEX @ #1+_ONE_DO (j) INDEX@ #1+_ONE_DO (j) ( ***CORRECTION*** ) 3 DUP INDEX@ #1<> ?SKIP #1- JINDEX@ INDEX@ #<> ?SKIP #1- JINDEX@ r #<> ?SKIP #1- * UNCOERCE %/ UNCOERCE2 %/ ( ***CORRECTION*** ) LOOP (j) LOOP (i) d UNCOERCE ONE{}N FPTR2 ^XEQ>ARRY DUP %0 xCON %1 BINT1 PUTREALEL s #1+_ONE_DO (k) INDEX@ #>$ BIGDISPROW1 2DUP m m 2GETEVAL DUP4UNROLL TOTEMPOB (n) #1+_ONE_DO (q) SWAP INDEX@ PULLREALEL ROT INDEX@ PULLREALEL ROT %* 4UNROLL LOOP (q) 2DROP UNCOERCE ONE{}N x>ARRY m #1+_ONE_DO (i) INDEX@ #2 #/ #+ #1+_ONE_DO (j) JINDEX@ INDEX@ 2GETEVAL PULLREALEL %CHS JINDEX@ TOTEMPOB !ii 1 JINDEX@ #1- #MAX JINDEX@ #1+ m #MIN #1+ SWAP DO (a) * 1 JINDEX@ INDEX@ ii #> ?SKIP #1- MAX 1 JINDEX@ INDEX@ ii #> ?SKIP #1- #MAX ( ***CORRECTION*** ) * JINDEX@ ii INDEX@ #> ?SKIP #1+ INDEX@ MIN JINDEX@ ii INDEX@ #> ?SKIP #1+ INDEX@ #MIN ( ***CORRECTION*** ) #1+SWAP DO (b) SWAP JINDEX@ INDEX@ 2GETEVAL PULLREALEL ROT %+ LOOP (b) LOOP (a) ROT JINDEX@ INDEX@ 2GETEVAL 3PICKSWAP PUTREALEL * JINDEX@ #1+ INDEX@ #- JINDEX@ DUP #1+ INDEX@ #- ( ***CORRECTION*** ) 2GETEVAL ROTSWAP PUTREALEL SWAP LOOP (j) LOOP (i) DROP m DUP r #>=_ ?SKIP #1+ !m LOOP (k) SWAP %0 xCON r 1 2GETEVAL UNCOERCE xDO %1 xPUTI xUNTIL % -64 xFS? xENDDO xDROP xDOT %6 s UNCOERCE x^ x/ ABND ; ; @ |
|||
11-04-2022, 10:51 PM
Post: #87
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
.
Hi, DavidM, DavidM Wrote:Valentin Albillo Wrote:Can someone please confirm that the listing is correct and will produce the stated result in the stated time ? As I said, truly impressive, about 10x faster than a physical HP-71B running my original BASIC program. I guess that only an assembly language version of it would be able to approach the 50g's performance, perhaps J-F Garnier would obligue, it's certainly within his outstanding capabilities ... Quote:Combine the above with the usual speed increase of SysRPL over UserRPL and you've got a significantly faster program than the original UserRPL version that C.Ret posted. Indeed ! Quote:Here's Gjermund's program with the typos corrected: [...] Fantastic ! It exceeds my best expectations of having a "certified" version of the SysRPL program, both for correct listing and accurate running time, a real gem of reliability and efficiency. Thank you very, very much, DavidM, for taking the trouble to fulfill my request, and for your really detailed explanations of its inner workings, never mind the corrections and improvements. Posts like yours is what adds the most value to these humble threads of mine and makes them everlasting contributions to the art of programming these wonderful vintage HP calcs, for the benefit of all of us. Much, much appreciated ! Best regards. V. All My Articles & other Materials here: Valentin Albillo's HP Collection |
|||
02-17-2023, 05:49 AM
(This post was last modified: 02-17-2023 05:50 AM by kostrse.)
Post: #88
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
(11-04-2022 03:31 PM)DavidM Wrote: Here's Gjermund's program with the typos corrected: Sorry, I'm trying to compile it on my HP-50G and it complains at the GOSBVL (Invalid Syntax). Is it a SysRPL or Assembly instruction or what? What should I do to make it work on a stock calculator? |
|||
02-18-2023, 04:25 PM
Post: #89
|
|||
|
|||
RE: [VA] SRC #012a - Then and Now: Probability
I responded to your questions in a separate post in order to avoid hijacking Valentin's original thread with implementation-specific details.
|
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)