50g: an interesting RAND anomaly - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: 50g: an interesting RAND anomaly (/thread-10333.html) Pages: 1 2 RE: 50g: an interesting RAND anomaly - Joe Horn - 03-20-2018 05:47 PM And here is HP's original annotated source code for the HP 48's RAND and RANDOMIZE functions, which I think uses the same algorithm as Cyrille's code in the previous posting in this thread. Code: ```        STITLE  %RAN ***************************************************************** ** ** Name: %RAN ** ** Abstract: Returns the next element of a pseudo-random **           number sequence. Updates RNSEED. ** ** Stack: --> % ** ** Exceptions: none ** ** Author: All inquiries may be addressed to the late great **         Homer Russel, wherever he is. ** ***************************************************************** =%RAN   CON(5)  (*)+5         GOSBVL  =SAVPTR         SETDEC                  Mode for rest         D0=(5)  =RNSEED         A=DAT0  15         A=0     S         ?A#0    W                GOYES   RND1            Has been set after coldstart         LCHEX   0999500333083533        C[W] = initial seed         A=C     W               A[W] = initial seed RND1    B=0     W         ABEX    W               B[W] = X(n-1)         LCHEX   0002851130928467        C[W] = C         GOTO    RND3 RND2    A=A+B   W                RND3    C=C-1   P         GONC    RND2         BSL     W         P=P+1         GONC    RND3         DAT0=A  15     [RNSEED] = X(n) = C*X(n-1) mod 10^15 > 0         ABEX    W               B[W] = X(n), A[A] = 0, A[S] = 0         P=      14         A=A-1   A               Exponent = -1 RNDSHF  ?B#0    P               Normalize X(n)         GOYES   TRUNC         A=A-1   A         BSL     WP         GONC    RNDSHF          (BET) TRUNC   A=B     M               Pack in mantissa         GOTO    push%L *****************************************************************         STITLE  %RANDOMIZE ***************************************************************** ** ** Name: %RANDOMIZE ** ** Abstract: Updates [RNSEED] by user input (% <> 0) or by **           the system clock (% = 0). ** ** Stack: % --> ** ** Exceptions: none ** ** Author: All inquiries may be addressed to the late great **         Homer Russel, wherever he is. ** ** Algorithm: ** **  x(n+1) = a * x(n) + c  (mod m) ** **  here: a = 2851130928467  [= 67 (mod 200)] **        c = 0 **        m = 10^15 ** **  If x(0) # 0  (mod 2) and x(0) # 0 (mod 5) then the period **  is 5*10^13 (Theorem D, pg. 19 of Knuth Vol. 2). ** ***Date Changed      Programmer            Reason ***------------      ----------      ------------------- **  1/15/86         P. McClellan           packed ** 08/24/89         SB             BugFix: Use SysTime (Not TIME) ***************************************************************** =%RANDOMIZE RPL :: DUP %0<> case DORANDOMIZE    DROP SysTime HXS>#    UNCOERCE     DORANDOMIZE ; ASSEMBLE DORANDOMIZE         CON(5)  (*)+5         GOSUB   pop1%           A[W] = %         A=0     S               A[W] = |%| RNDMF   B=0     A         ABEX    X               B[A] = exponent(%) RNDM2   B=B+1   A         ?A=0    W         GOYES   RNDM3         C=B     A         C=C+C   X         GONC    RNDM3         ?C=0    X         GOYES   RNDM3         ASR     W               Normalize         GONC    RNDM2           (BET) RNDM3   A=B     X         ASL     X         A=A+1   A               insure odd but not multiple of 5         D0=(5)  =RNSEED         DAT0=A  15              Write new seed         P=      0         SETHEX getptrloop         GOVLNG  =GETPTRLOOP``` RE: 50g: an interesting RAND anomaly - DavidM - 03-21-2018 02:10 AM Cyrille and Joe, thanks for the source listings. For what it's worth, I already knew the algorithm (I had previously disassembled it with Nosy). I've actually been using it in the last few versions of the list randomizer in the ListExt library (LSHUF). While I'm sure I should have already recognized why the anomaly occurred, I'm afraid it wasn't (and still isn't) abundantly obvious to this casual observer. But the identified anomaly was simply evidence of the bigger issue for me, which is that using MOD on an unaltered seed is definitely not an appropriate way to reduce the seed to an integer in a specified range. I get the sense from Cyrille's and others' responses that I should simply do the extra multiplication (seed range * CEIL) instead of trying to find a faster integer-based method of converting the seed. I was hoping that there would be a more efficient way, given that the RAND function itself had already done most of the heavy lifting. I'm already using the Saturn+ opcodes in other parts of this command (when available, of course), so I'll do the same here. That should mitigate some of the performance degradation from the extra multiplication. I had thought of doing that with the RAND function in the past (specifically the section where the previous seed is multiplied by 2851130928467), but thought it would be safer not to mess with it. I may revisit that now.