On Thu, 11 Mar 2021 at 00:58, Bruce Momjian <br...@momjian.us> wrote:
>
> Maybe Dean Rasheed can help because of his math background --- CC'ing him.
>
Reading the thread I can see how such a function might be useful to
scatter non-uniformly random values.
The implementation looks plausible too, though it adds quite a large
amount of new code. The main thing that concerns me is justifying the
code. With this kind of thing, it's all too easy to overlook corner
cases and end up with trivial sequences in certain special cases. I'd
feel better about that if we were implementing a standard algorithm
with known pedigree.
Thinking about the use case for this, it seems that it's basically
designed to turn a set of non-uniform random numbers (produced by
random_exponential() et al.) into another set of non-uniform random
numbers, where the non-uniformity is scattered so that the more/less
common values aren't all clumped together.
I'm wondering if that's something that can't be done more simply by
passing the non-uniform random numbers through the uniform random
number generator to scatter them uniformly across some range -- e.g.,
given an integer n, return the n'th value from the sequence produced
by random(), starting from some initial seed -- i.e., implement
nth_random(lb, ub, seed, n). That would actually be pretty
straightforward to implement using O(log(n)) time to execute (see the
attached python example), though it wouldn't generate a permutation,
so it'd need a bit of thought to see if it met the requirements.
Regards,
Dean
import math
a = 0x5deece66d
c = 0xb
m = 0xffffffffffff
def erand48(seed):
x = (seed[2] << 32) | (seed[1] << 16) | seed[0]
x = (x * a + c) & m
seed[0] = x & 0xffff
seed[1] = (x >> 16) & 0xffff
seed[2] = (x >> 32) & 0xffff
return math.ldexp(x, -48)
def nth_erand48(seed, n):
x = (seed[2] << 32) | (seed[1] << 16) | seed[0]
# Binary exponentiation to compute a^n and the sum of the geometric
# series gs = 1+a+...+a^(n-1)
a_pow_n = 1
gs = 0
t = 1
e = a
while n > 0:
if n & 1 == 1:
a_pow_n = (a_pow_n * e) & m
gs = (gs * e + t) & m
t = ((e+1)*t) & m
e = (e*e) & m
n = n/2
# n'th value in erand48 sequence
x = (x * a_pow_n + gs * c) & m
return math.ldexp(x, -48)
iseed = [1234, 5678, 9012]
seed = list(iseed)
for n in range(1,21):
r = erand48(seed)
print n, seed, r, nth_erand48(iseed, n)
