301 lines
9.8 KiB
Plaintext
301 lines
9.8 KiB
Plaintext
# vrandom - Hardware Random Number Generation Library
|
|
# Implemented using inline assembly based on x86 RDRAND / RDSEED instructions
|
|
import math
|
|
|
|
|
|
# Retry limit (RDRAND usually doesn't need retries, RDSEED might)
|
|
define RDRAND_RETRY_LIMIT = 10
|
|
define RDSEED_RETRY_LIMIT = 100
|
|
|
|
|
|
# ============================================================
|
|
# RDRAND - Hardware random numbers (output from a pseudo-random number generator)
|
|
# ============================================================
|
|
def rdrand16() -> uint16:
|
|
val: uint32 = 0
|
|
__asm__(f"rdrand {asmout(val, ASM_DESCR.OUTPUT_REG)}",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC])
|
|
return uint16(val)
|
|
|
|
def rdrand32() -> uint32:
|
|
val: uint32 = 0
|
|
__asm__(f"rdrand {asmout(val, ASM_DESCR.OUTPUT_REG)}",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC])
|
|
return val
|
|
|
|
def rdrand64() -> uint64:
|
|
val: uint64 = 0
|
|
__asm__(f"rdrand {asmout(val, ASM_DESCR.OUTPUT_REG)}",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC])
|
|
return val
|
|
|
|
|
|
# ============================================================
|
|
# RDSEED - Hardware random seeds (direct output from entropy source)
|
|
# ============================================================
|
|
def rdseed16() -> uint16:
|
|
val: uint32 = 0
|
|
__asm__(f"rdseed {asmout(val, ASM_DESCR.OUTPUT_REG)}",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC])
|
|
return uint16(val)
|
|
|
|
def rdseed32() -> uint32:
|
|
val: uint32 = 0
|
|
__asm__(f"rdseed {asmout(val, ASM_DESCR.OUTPUT_REG)}",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC])
|
|
return val
|
|
|
|
def rdseed64() -> uint64:
|
|
val: uint64 = 0
|
|
__asm__(f"rdseed {asmout(val, ASM_DESCR.OUTPUT_REG)}",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC])
|
|
return val
|
|
|
|
|
|
# ============================================================
|
|
# Step version (with success flag, matching GCC builtin interface)
|
|
# Returns success flag: 1 on success, 0 on failure; random value written to *p
|
|
# ============================================================
|
|
def rdrand16_step(p: uint16*) -> int:
|
|
val: uint32 = 0
|
|
ok: uint32 = 0
|
|
__asm__(f"""rdrand {asmout(val, ASM_DESCR.OUTPUT_REG)}
|
|
setc al
|
|
movzx {asmout(ok, ASM_DESCR.OUTPUT_REG)}, al""",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG), asmout(ok, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC, ASM_DESCR.CLOBBER_AL])
|
|
p[0] = uint16(val)
|
|
return int(ok)
|
|
|
|
def rdrand32_step(p: uint32*) -> int:
|
|
val: uint32 = 0
|
|
ok: uint32 = 0
|
|
__asm__(f"""rdrand {asmout(val, ASM_DESCR.OUTPUT_REG)}
|
|
setc al
|
|
movzx {asmout(ok, ASM_DESCR.OUTPUT_REG)}, al""",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG), asmout(ok, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC, ASM_DESCR.CLOBBER_AL])
|
|
p[0] = val
|
|
return int(ok)
|
|
|
|
def rdrand64_step(p: uint64*) -> int:
|
|
val: uint64 = 0
|
|
ok: uint32 = 0
|
|
__asm__(f"""rdrand {asmout(val, ASM_DESCR.OUTPUT_REG)}
|
|
setc al
|
|
movzx {asmout(ok, ASM_DESCR.OUTPUT_REG)}, al""",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG), asmout(ok, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC, ASM_DESCR.CLOBBER_AL])
|
|
p[0] = val
|
|
return int(ok)
|
|
|
|
def rdseed16_step(p: uint16*) -> int:
|
|
val: uint32 = 0
|
|
ok: uint32 = 0
|
|
__asm__(f"""rdseed {asmout(val, ASM_DESCR.OUTPUT_REG)}
|
|
setc al
|
|
movzx {asmout(ok, ASM_DESCR.OUTPUT_REG)}, al""",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG), asmout(ok, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC, ASM_DESCR.CLOBBER_AL])
|
|
p[0] = uint16(val)
|
|
return int(ok)
|
|
|
|
def rdseed32_step(p: uint32*) -> int:
|
|
val: uint32 = 0
|
|
ok: uint32 = 0
|
|
__asm__(f"""rdseed {asmout(val, ASM_DESCR.OUTPUT_REG)}
|
|
setc al
|
|
movzx {asmout(ok, ASM_DESCR.OUTPUT_REG)}, al""",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG), asmout(ok, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC, ASM_DESCR.CLOBBER_AL])
|
|
p[0] = val
|
|
return int(ok)
|
|
|
|
def rdseed64_step(p: uint64*) -> int:
|
|
val: uint64 = 0
|
|
ok: uint32 = 0
|
|
__asm__(f"""rdseed {asmout(val, ASM_DESCR.OUTPUT_REG)}
|
|
setc al
|
|
movzx {asmout(ok, ASM_DESCR.OUTPUT_REG)}, al""",
|
|
out=[asmout(val, ASM_DESCR.OUTPUT_REG), asmout(ok, ASM_DESCR.OUTPUT_REG)],
|
|
op=[ASM_DESCR.CLOBBER_CC, ASM_DESCR.CLOBBER_AL])
|
|
p[0] = val
|
|
return int(ok)
|
|
|
|
|
|
# ============================================================
|
|
# Convenient interface with retry (automatically retries until it succeeds or reaches the limit)
|
|
# ============================================================
|
|
def random_u16() -> uint16:
|
|
r: uint16 = 0
|
|
for i in range(int(RDRAND_RETRY_LIMIT)):
|
|
if rdrand16_step(&r) == 1: return r
|
|
return 0
|
|
|
|
def random_u32() -> uint32:
|
|
r: uint32 = 0
|
|
for i in range(int(RDRAND_RETRY_LIMIT)):
|
|
if rdrand32_step(&r) == 1: return r
|
|
return 0
|
|
|
|
def random_u64() -> uint64:
|
|
r: uint64 = 0
|
|
for i in range(int(RDRAND_RETRY_LIMIT)):
|
|
if rdrand64_step(&r) == 1: return r
|
|
return 0
|
|
|
|
def seed_u16() -> uint16:
|
|
r: uint16 = 0
|
|
for i in range(int(RDRAND_RETRY_LIMIT)):
|
|
if rdseed16_step(&r) == 1: return r
|
|
return 0
|
|
|
|
def seed_u32() -> uint32:
|
|
r: uint32 = 0
|
|
for i in range(int(RDRAND_RETRY_LIMIT)):
|
|
if rdseed32_step(&r) == 1: return r
|
|
return 0
|
|
|
|
def seed_u64() -> uint64:
|
|
r: uint64 = 0
|
|
for i in range(int(RDSEED_RETRY_LIMIT)):
|
|
if rdseed64_step(&r) == 1: return r
|
|
return 0
|
|
|
|
|
|
# ============================================================
|
|
# Python random standard library function implementations
|
|
# Default types: uint64 (u64), double (f64), int (i32)
|
|
# ============================================================
|
|
# ---------- ONE, Random integers ----------
|
|
def randint(a: int, b: int) -> int:
|
|
return int(a + int(random_u64() % uint64(b - a + 1)))
|
|
|
|
def randrange(start: int, stop: int, step: int = 1) -> int:
|
|
width: int = 0
|
|
if step > 0: width = (stop - start + step - 1) // step
|
|
else: width = (start - stop - step - 1) // (-step)
|
|
return int(start + int(random_u64() % uint64(width)) * step)
|
|
|
|
def getrandbits(k: int) -> uint64:
|
|
r: uint64 = rdrand64()
|
|
if k >= 64: return r
|
|
mask: uint64 = (uint64(1) << k) - uint64(1)
|
|
return r & mask
|
|
|
|
# ---------- TWO, Random floating-point numbers ----------
|
|
def random() -> double:
|
|
return double(rdrand64() >> 11) * (1.0 / 9007199254740992.0)
|
|
|
|
def uniform(a: double, b: double) -> double:
|
|
return a + (b - a) * random()
|
|
|
|
def triangular(low: double, high: double, mode: double) -> double:
|
|
u: double = random()
|
|
c: double = (mode - low) / (high - low)
|
|
tmp: double = 0.0
|
|
if u > c:
|
|
u = 1.0 - u; c = 1.0 - c
|
|
tmp = low; low = high; high = tmp
|
|
return low + (high - low) * math.sqrt(u * c)
|
|
|
|
def betavariate(alpha: double, beta: double) -> double:
|
|
y: double = gammavariate(alpha, 1.0)
|
|
if y != 0.0: return y / (y + gammavariate(beta, 1.0))
|
|
return 0.0
|
|
|
|
def expovariate(lambd: double) -> double:
|
|
u: double = random()
|
|
if u < 1e-300: u = 1e-300
|
|
return -math.log(1.0 - u) / lambd
|
|
|
|
def gammavariate(alpha: double, beta: double) -> double:
|
|
u: double = 0.0
|
|
d: double = 0.0
|
|
c: double = 0.0
|
|
x: double = 0.0
|
|
v: double = 0.0
|
|
if alpha <= 0.0 or beta <= 0.0: return 0.0
|
|
if alpha < 1.0:
|
|
u = random()
|
|
if u < 1e-300: u = 1e-300
|
|
return gammavariate(alpha + 1.0, beta) * math.pow(u, 1.0 / alpha)
|
|
d = alpha - 1.0 / 3.0
|
|
c = 1.0 / math.sqrt(9.0 * d)
|
|
while True:
|
|
x = gauss(0.0, 1.0)
|
|
v = 1.0 + c * x
|
|
if v > 0.0:
|
|
v = v * v * v
|
|
u = random()
|
|
if u < 1.0 - 0.0331 * x * x * x * x: return d * v * beta
|
|
if math.log(u) < 0.5 * x * x + d * (1.0 - v + math.log(v)): return d * v * beta
|
|
|
|
def gauss(mu: double, sigma: double) -> double:
|
|
u1: double = random()
|
|
u2: double = random()
|
|
if u1 < 1e-300: u1 = 1e-300
|
|
z0: double = math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.U_M_PI * u2)
|
|
return mu + sigma * z0
|
|
|
|
def normalvariate(mu: double, sigma: double) -> double:
|
|
NV_MAGICCONST: double = 1.7155277699214135
|
|
u1: double = 0.0
|
|
u2: double = 0.0
|
|
z: double = 0.0
|
|
zz: double = 0.0
|
|
while True:
|
|
u1 = random()
|
|
u2 = random()
|
|
if u2 < 1e-300: continue
|
|
z = NV_MAGICCONST * (u1 - 0.5) / u2
|
|
zz = z * z / 4.0
|
|
if zz <= 0.0 - math.log(u2): break
|
|
return mu + z * sigma
|
|
|
|
def lognormvariate(mu: double, sigma: double) -> double:
|
|
return math.exp(gauss(mu, sigma))
|
|
|
|
def vonmisesvariate(mu: double, kappa: double) -> double:
|
|
if kappa <= 1e-6: return math.U_M_PI * 2.0 * random()
|
|
s: double = 0.5 / kappa
|
|
r: double = s + math.sqrt(1.0 + s * s)
|
|
u1: double = 0.0
|
|
u2: double = 0.0
|
|
u3: double = 0.0
|
|
z: double = 0.0
|
|
d: double = 0.0
|
|
q: double = 0.0
|
|
f: double = 0.0
|
|
theta: double = 0.0
|
|
two_pi: double = math.U_M_PI * 2.0
|
|
while True:
|
|
u1 = random()
|
|
z = math.cos(math.U_M_PI * u1)
|
|
d = 1.0 / (r + z)
|
|
u2 = random()
|
|
if (u2 < 1.0 - d * d) or (u2 < (1.0 - z) * math.exp(kappa * (z - 1.0) - d * kappa * r)): break
|
|
q = 1.0 / r
|
|
f = (q + z) / (1.0 + q * z)
|
|
u3 = random()
|
|
theta = mu + (math.acos(f) if u3 > 0.5 else -math.acos(f))
|
|
while theta < 0.0: theta += two_pi
|
|
while theta >= two_pi: theta -= two_pi
|
|
return theta
|
|
|
|
def paretovariate(alpha: double) -> double:
|
|
u: double = random()
|
|
if u < 1e-300: u = 1e-300
|
|
return 1.0 / math.pow(1.0 - u, 1.0 / alpha)
|
|
|
|
def weibullvariate(alpha: double, beta: double) -> double:
|
|
u: double = random()
|
|
if u < 1e-300: u = 1e-300
|
|
return alpha * math.pow(-1.0 * math.log(1.0 - u), 1.0 / beta)
|