380 lines
12 KiB
Python
380 lines
12 KiB
Python
# vrandom - 硬件随机数生成库
|
||
# 基于 x86 RDRAND / RDSEED 指令，使用内嵌汇编实现
|
||
import t, c
|
||
import vipermath
|
||
|
||
|
||
# 重试次数限制（RDRAND 通常不需要重试，RDSEED 可能需要）
|
||
RDRAND_RETRY_LIMIT: t.CDefine = 10
|
||
RDSEED_RETRY_LIMIT: t.CDefine = 100
|
||
|
||
|
||
# ============================================================
|
||
# RDRAND - 硬件随机数（伪随机数发生器输出）
|
||
# ============================================================
|
||
|
||
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 - 硬件随机种子（熵源直接输出）
|
||
# ============================================================
|
||
|
||
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 版本（带成功标志，匹配 GCC builtin 接口）
|
||
# 返回 1 表示成功，0 表示失败；随机值写入 *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)
|
||
|
||
|
||
# ============================================================
|
||
# 带重试的便捷接口（自动重试直到成功或达到上限）
|
||
# ============================================================
|
||
|
||
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 标准库函数实现
|
||
# 默认类型：t.CUInt64T (u64), t.CDouble (f64), t.CInt (i32)
|
||
# ============================================================
|
||
|
||
# ---------- 一、随机整数 ----------
|
||
|
||
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
|
||
|
||
|
||
# ---------- 二、随机浮点数 ----------
|
||
|
||
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)
|