50g: an interesting RAND anomaly

[Joe Horn]

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

[DavidM]

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.