erand48.py

text/x-python

Download

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)