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