[math-fun] a few cheap linear congruential randgen cases
X := RotateLeft(X,L) + X X := X + carry (look Ma, no multiply instruction!) implements a linear congruential generator modulo the "Mersenne" modulus 2^wordsize-1, with multiplier 1+2^L, where the first add generates boolean carry, used in the second add. This code however has an off-by-1 error when the second add involves both X=2^wordsize-1 and carry=1; thw wrong result is X=0. That, when it happens, is embarrassing, since then all subsequent "random" numbers will be 0. That error happens only once per period and if desired could be corrected by adding a third statement if(X==0){ replace X by some other appropriate value } someplace, the new value being chosen to make the period be whatever we like, e.g. to get a prime period skip over a few entries in the hypothetical full period, just enough to reduce period to some large prime a bit below it. We actually can, in basically the same way, implement any of the four multipliers 1+-2^K and 2^K+-1 by using some subtractions instead of adds. "Fermat" moduli 2^wordsize+1 can also be done, but unfortunately it costs more instructions: if(X==0){ X := replacement value; } Y := ShiftLeft(X,L); Z := ShiftRight(X,L); if(Y==0 && Z!=0){ Y:=1; } X := (+-X) + (+-)(Y - Z); X := X +- carry; with three appropriate sign choices, most simply +, +, and -. If we skipped the final X+carry, so only employ the single line X := RotateLeft(X,L) + X then we are implementing a linear congruential randgen modulo 2^wordsize-1, except that "noise" of value -1 or 0 is added each step (the noise being the skipped carries) -- and this noise will tend to shorten the period drastically by sending us back to some already-used X. Since the "Fermat" number 2^64+1 = 274177*67280421310721 we can in this way achieve period arbitrarily below 274176*67280421310720=2^63.999994738; the elements of the period are the numbers relatively prime to 2^64+1, except for the ones skipped. Similarly 2^32+1=641*6700417 allowing any period below 640*6700416=2^31.99774733329, and 2^16+1 actually is the (last?) Fermat prime. -- Warren D. Smith http://RangeVoting.org <-- add your endorsement (by clicking "endorse" as 1st step)
participants (1)
-
Warren D Smith