Files
Viper_VRTE/vipermath.vp
T
TermiNexus 52af24e85d fix(includes): asm 内联汇编输出占位符/clobber、vipermath 位重解释、string.split 容量参数
- asm.vp: inb/inw/inl 模板引用输出占位符 \ 并补 ~{rax}(原取硬编码 eax 且静默破坏 EAX,对齐 mm.vp cmos_read 模式)

- vipermath.vp: isnan/isinf 改 c.Load 位重解释(原 fptosi 位模式丢失导致恒返回 0)

- string.vp: split 增加 cap 容量参数,token 数 <= cap-1,防止越界写
2026-09-11 14:26:47 +08:00

552 lines
16 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
import t, c
U_M_PI: t.CDefine = 3.14159265358979323846
U_M_E: t.CDefine = 2.71828182845904523536
U_M_PI_2: t.CDefine = 1.57079632679489661923
U_M_PI_4: t.CDefine = 0.78539816339744830962
U_M_1_PI: t.CDefine = 0.31830988618379067154
U_M_2_PI: t.CDefine = 0.63661977236758134308
U_M_LN2: t.CDefine = 0.69314718055994530942
U_M_LN10: t.CDefine = 2.30258509299404568402
U_M_2_SQRT_PI: t.CDefine = 1.77245385090551602730
def radians(degrees: t.CDouble) -> t.CDouble:
return degrees * (U_M_PI / t.CDouble(180.0))
def degrees(radians: t.CDouble) -> t.CDouble:
return radians * (t.CDouble(180.0) / U_M_PI)
def sin(x: t.CDouble) -> t.CDouble:
while x > U_M_PI:
x -= t.CDouble(2.0) * U_M_PI
while x < -U_M_PI:
x += t.CDouble(2.0) * U_M_PI
x2: t.CDouble = x * x
result: t.CDouble = x
term: t.CDouble = x
term = term * x2 / t.CDouble(6.0)
result = result - term
term = term * x2 / t.CDouble(20.0)
result = result + term
term = term * x2 / t.CDouble(42.0)
result = result - term
term = term * x2 / t.CDouble(72.0)
result = result + term
term = term * x2 / t.CDouble(110.0)
result = result - term
term = term * x2 / t.CDouble(156.0)
result = result + term
term = term * x2 / t.CDouble(210.0)
result = result - term
term = term * x2 / t.CDouble(272.0)
result = result + term
return result
def cos(x: t.CDouble) -> t.CDouble:
while x > U_M_PI:
x -= t.CDouble(2.0) * U_M_PI
while x < -U_M_PI:
x += t.CDouble(2.0) * U_M_PI
x2: t.CDouble = x * x
result: t.CDouble = t.CDouble(1.0)
term: t.CDouble = t.CDouble(1.0)
term = term * x2 / t.CDouble(2.0)
result = result - term
term = term * x2 / t.CDouble(12.0)
result = result + term
term = term * x2 / t.CDouble(30.0)
result = result - term
term = term * x2 / t.CDouble(56.0)
result = result + term
term = term * x2 / t.CDouble(90.0)
result = result - term
term = term * x2 / t.CDouble(132.0)
result = result + term
term = term * x2 / t.CDouble(182.0)
result = result - term
term = term * x2 / t.CDouble(240.0)
result = result + term
return result
def tan(x: t.CDouble) -> t.CDouble:
return sin(x) / cos(x)
def asin(x: t.CDouble) -> t.CDouble:
if x > t.CDouble(1.0):
x = t.CDouble(1.0)
if x < t.CDouble(-1.0):
x = t.CDouble(-1.0)
x2: t.CDouble = x * x
result: t.CDouble = x
result += x * x2 / t.CDouble(6.0)
result += t.CDouble(3.0) * x * x2 * x2 / t.CDouble(40.0)
result += t.CDouble(5.0) * x * x2 * x2 * x2 / t.CDouble(112.0)
result += t.CDouble(35.0) * x * x2 * x2 * x2 * x2 / t.CDouble(1152.0)
result += t.CDouble(63.0) * x * x2 * x2 * x2 * x2 * x2 / t.CDouble(2816.0)
return result
def acos(x: t.CDouble) -> t.CDouble:
return U_M_PI / t.CDouble(2.0) - asin(x)
def _atan_core(x: t.CDouble) -> t.CDouble:
x2: t.CDouble = x * x
result: t.CDouble = x
term: t.CDouble = x
term = term * x2 / t.CDouble(3.0)
result = result - term
term = term * x2 / t.CDouble(5.0)
result = result + term
term = term * x2 / t.CDouble(7.0)
result = result - term
term = term * x2 / t.CDouble(9.0)
result = result + term
term = term * x2 / t.CDouble(11.0)
result = result - term
term = term * x2 / t.CDouble(13.0)
result = result + term
term = term * x2 / t.CDouble(15.0)
result = result - term
term = term * x2 / t.CDouble(17.0)
result = result + term
term = term * x2 / t.CDouble(19.0)
result = result - term
term = term * x2 / t.CDouble(21.0)
result = result + term
return result
def atan(x: t.CDouble) -> t.CDouble:
ax: t.CDouble = x
if ax < t.CDouble(0.0):
ax = -ax
if ax > t.CDouble(1.0):
return t.CDouble(0.0) - _atan_core(t.CDouble(1.0) / ax) + U_M_PI_2 if x > t.CDouble(0.0) else _atan_core(t.CDouble(1.0) / ax) - U_M_PI_2
return _atan_core(x) if x >= t.CDouble(0.0) else t.CDouble(0.0) - _atan_core(ax)
def atan2(y: t.CDouble, x: t.CDouble) -> t.CDouble:
if x > t.CDouble(0.0):
return atan(y / x)
elif x < t.CDouble(0.0):
if y >= t.CDouble(0.0):
return atan(y / x) + U_M_PI
else:
return atan(y / x) - U_M_PI
else:
if y > t.CDouble(0.0):
return U_M_PI_2
elif y < t.CDouble(0.0):
return t.CDouble(0.0) - U_M_PI_2
else:
return t.CDouble(0.0)
def sinh(x: t.CDouble) -> t.CDouble:
ex: t.CDouble = exp(x)
e_negx: t.CDouble = exp(-x)
return (ex - e_negx) / t.CDouble(2.0)
def cosh(x: t.CDouble) -> t.CDouble:
ex: t.CDouble = exp(x)
e_negx: t.CDouble = exp(-x)
return (ex + e_negx) / t.CDouble(2.0)
def tanh(x: t.CDouble) -> t.CDouble:
ex: t.CDouble = exp(x)
e_negx: t.CDouble = exp(-x)
return (ex - e_negx) / (ex + e_negx)
def exp(x: t.CDouble) -> t.CDouble:
"""e^x
范围归约 + 短泰勒: 直接对 x 展开 18 项泰勒在 |x| 稍大时就会发散
(x=-10 得 54.5, 真值 4.5e-5; x=-30 得 3.7e10), 负数端更是给出巨大正值,
会让 softmax 的分母被撑爆 → 所有概率 ≈ 0。改为先做
x = k*ln2 + r (|r| <= ln2/2), 再算 exp(r) 后乘 2^k。
"""
# 先用"有限且落在范围里"的正向判断过一遍: NaN 与任何数比较都为假, 会同时漏过
# 下面两个判断, 而 t.CInt(NaN) 在 x86 上是 INT_MIN -> 后面 2^k 缩放循环要跑
# 21 亿次, 表现为"卡死"。故 NaN/越界必须在这里挡掉。
if x <= t.CDouble(709.0) and x >= t.CDouble(-745.0):
pass
else:
if x > t.CDouble(709.0):
return t.CDouble(1.7976931348623157e308)
if x < t.CDouble(-745.0):
return t.CDouble(0.0)
return x # NaN: 原样返回, 让 NaN 传播而不是卡死
# 归约: 把 x 折到 [-ln2/2, ln2/2], 整数部分累积到 k (内联取整, 避免前向引用)
ln2: t.CDouble = U_M_LN2
q: t.CDouble = x / ln2
k: t.CInt = t.CInt(q)
frac: t.CDouble = q - t.CDouble(k)
if frac > t.CDouble(0.5):
k = k + 1
elif frac < t.CDouble(-0.5):
k = k - 1
r: t.CDouble = x - t.CDouble(k) * ln2
# exp(r), |r| <= 0.3466: 12 项后截断误差 < 1e-17
result: t.CDouble = t.CDouble(1.0)
term: t.CDouble = t.CDouble(1.0)
n: t.CInt = 1
while n <= 12:
term = term * r / t.CDouble(n)
result = result + term
n += 1
while k > 0:
result = result * t.CDouble(2.0)
k = k - 1
while k < 0:
result = result * t.CDouble(0.5)
k = k + 1
return result
def log(x: t.CDouble) -> t.CDouble:
if x <= t.CDouble(0.0):
return t.CDouble(-1e308)
exponent: t.CInt = 0
while x >= t.CDouble(2.0):
x /= t.CDouble(2.0)
exponent += 1
while x < t.CDouble(1.0):
x *= t.CDouble(2.0)
exponent -= 1
y: t.CDouble = (x - t.CDouble(1.0)) / (x + t.CDouble(1.0))
y2: t.CDouble = y * y
result: t.CDouble = t.CDouble(2.0) * y
term: t.CDouble = t.CDouble(2.0) * y
term = term * y2
result = result + term / t.CDouble(3.0)
term = term * y2
result = result + term / t.CDouble(5.0)
term = term * y2
result = result + term / t.CDouble(7.0)
term = term * y2
result = result + term / t.CDouble(9.0)
term = term * y2
result = result + term / t.CDouble(11.0)
term = term * y2
result = result + term / t.CDouble(13.0)
term = term * y2
result = result + term / t.CDouble(15.0)
term = term * y2
result = result + term / t.CDouble(17.0)
term = term * y2
result = result + term / t.CDouble(19.0)
term = term * y2
result = result + term / t.CDouble(21.0)
result += t.CDouble(exponent) * U_M_LN2
return result
def sqrt(x: t.CDouble) -> t.CDouble:
if x < t.CDouble(0.0):
return t.CDouble(0.0)
if x == t.CDouble(0.0):
return t.CDouble(0.0)
guess: t.CDouble = x / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
guess = (guess + x / guess) / t.CDouble(2.0)
return guess
def abs(x: t.CInt) -> t.CInt:
return x if x >= 0 else -x
def fabs(x: t.CDouble) -> t.CDouble:
return x if x >= t.CDouble(0.0) else -x
def labs(x: t.CLong) -> t.CLong:
return x if x >= 0 else -x
def factorial(n: t.CInt) -> t.CDouble:
if n < 0: return t.CDouble(0.0)
if n == 0 or n == 1: return t.CDouble(1.0)
result: t.CDouble = t.CDouble(1.0)
for i in range(2, n + 1):
result *= i
return result
def combination(n: t.CInt, k: t.CInt) -> t.CDouble:
if k < 0 or k > n: return t.CDouble(0.0)
if k == 0 or k == n: return t.CDouble(1.0)
k = k if k < n - k else n - k
result: t.CDouble = t.CDouble(1.0)
for i in range(1, k + 1):
result = result * (n - k + i) / i
return result
def permutation(n: t.CInt, k: t.CInt) -> t.CDouble:
if k < 0 or k > n: return t.CDouble(0.0)
result: t.CDouble = t.CDouble(1.0)
for i in range(0, k):
result *= (n - i)
return result
def pow(x: t.CDouble, y: t.CDouble) -> t.CDouble:
if x == t.CDouble(0.0):
if y > t.CDouble(0.0):
return t.CDouble(0.0)
return t.CDouble(1.0)
if x < t.CDouble(0.0):
iy: t.CDouble = y * log(-x)
return exp(iy)
return exp(y * log(x))
def powf(x: t.CDouble, y: t.CDouble) -> t.CDouble:
return pow(x, y)
def cbrt(x: t.CDouble) -> t.CDouble:
return pow(x, t.CDouble(1.0) / t.CDouble(3.0))
def hypot(x: t.CDouble, y: t.CDouble) -> t.CDouble:
return sqrt(x * x + y * y)
def floor(x: t.CDouble) -> t.CDouble:
if x >= t.CDouble(0.0):
return t.CDouble(t.CInt(x))
i: t.CInt = t.CInt(x)
if x < t.CDouble(i):
return t.CDouble(i - 1)
else:
return t.CDouble(i)
def ceil(x: t.CDouble) -> t.CDouble:
if x < t.CDouble(0.0):
return t.CDouble(t.CInt(x))
i: t.CInt = t.CInt(x)
if x > t.CDouble(i):
return t.CDouble(i + 1)
else:
return t.CDouble(i)
def round(x: t.CDouble) -> t.CDouble:
return floor(x + t.CDouble(0.5))
def trunc(x: t.CDouble) -> t.CDouble:
if x >= t.CDouble(0.0):
return floor(x)
else:
return ceil(x)
def fmod(x: t.CDouble, y: t.CDouble) -> t.CDouble:
if y == t.CDouble(0.0):
return t.CDouble(0.0)
result: t.CDouble = x - y * floor(x / y)
if result < t.CDouble(0.0):
result += y
return result
def fmodf(x: float, y: float) -> float:
if y == float(0.0):
return float(0.0)
result: float = x - y * floorf(x / y)
if result < float(0.0):
result += y
return result
def modf(x: t.CDouble, iptr: t.CDouble | t.CPtr) -> t.CDouble:
int_part: t.CDouble = floor(x)
# c.Set(a, b) = a = b:a 为 c.Deref(iptr) 左值表达式 → *iptr = int_part
# TPV 的 c.Deref 立即 load 返回右值,故 c.Set 需识别 Deref 左值形态(translate_c_set 形态 1)
c.Set(c.Deref(iptr), int_part)
return x - int_part
def log10(x: t.CDouble) -> t.CDouble:
return log(x) / U_M_LN10
def log2(x: t.CDouble) -> t.CDouble:
return log(x) / U_M_LN2
def exp2(x: t.CDouble) -> t.CDouble:
return pow(t.CDouble(2.0), x)
def expm1(x: t.CDouble) -> t.CDouble:
return exp(x) - t.CDouble(1.0)
def log1p(x: t.CDouble) -> t.CDouble:
return log(t.CDouble(1.0) + x)
def asinh(x: t.CDouble) -> t.CDouble:
return log(x + sqrt(x * x + t.CDouble(1.0)))
def acosh(x: t.CDouble) -> t.CDouble:
return log(x + sqrt(x * x - t.CDouble(1.0)))
def atanh(x: t.CDouble) -> t.CDouble:
return (log(t.CDouble(1.0) + x) - log(t.CDouble(1.0) - x)) / t.CDouble(2.0)
def gamma(x: t.CDouble) -> t.CDouble:
if x <= t.CDouble(0.0): return t.CDouble(0.0)
if x == t.CDouble(1.0) or x == t.CDouble(2.0): return t.CDouble(1.0)
ln_gamma: t.CDouble = t.CDouble(0.5) * log(t.CDouble(2.0) * U_M_PI / x) + x * log(x / U_M_E)
return exp(ln_gamma)
def erf(x: t.CDouble) -> t.CDouble:
sign: t.CDouble = t.CDouble(1.0) if x >= t.CDouble(0.0) else t.CDouble(-1.0)
abs_x: t.CDouble = fabs(x)
t0: t.CDouble = t.CDouble(4.0) * abs_x * abs_x / U_M_PI
exp_term: t.CDouble = exp(-t0)
return sign * sqrt(t.CDouble(1.0) - exp_term)
def erfc(x: t.CDouble) -> t.CDouble:
return t.CDouble(1.0) - erf(x)
def sqrtf(x: t.CFloat) -> t.CFloat:
if x < t.CFloat(0.0):
return t.CFloat(0.0)
if x == t.CFloat(0.0):
return t.CFloat(0.0)
guess: t.CFloat = x / t.CFloat(2.0)
for i in range(10):
guess = (guess + x / guess) / t.CFloat(2.0)
return guess
def sinf(x: t.CFloat) -> t.CFloat:
while x > t.CFloat(3.14159265):
x -= t.CFloat(6.28318530)
while x < t.CFloat(-3.14159265):
x += t.CFloat(6.28318530)
x2: t.CFloat = x * x
result: t.CFloat = x
term: t.CFloat = x
term = term * x2 / t.CFloat(6.0)
result = result - term
term = term * x2 / t.CFloat(20.0)
result = result + term
term = term * x2 / t.CFloat(42.0)
result = result - term
term = term * x2 / t.CFloat(72.0)
result = result + term
term = term * x2 / t.CFloat(110.0)
result = result - term
term = term * x2 / t.CFloat(156.0)
result = result + term
term = term * x2 / t.CFloat(210.0)
result = result - term
term = term * x2 / t.CFloat(272.0)
result = result + term
return result
def cosf(x: t.CFloat) -> t.CFloat:
while x > t.CFloat(3.14159265):
x -= t.CFloat(6.28318530)
while x < t.CFloat(-3.14159265):
x += t.CFloat(6.28318530)
x2: t.CFloat = x * x
result: t.CFloat = t.CFloat(1.0)
term: t.CFloat = t.CFloat(1.0)
term = term * x2 / t.CFloat(2.0)
result = result - term
term = term * x2 / t.CFloat(12.0)
result = result + term
term = term * x2 / t.CFloat(30.0)
result = result - term
term = term * x2 / t.CFloat(56.0)
result = result + term
term = term * x2 / t.CFloat(90.0)
result = result - term
term = term * x2 / t.CFloat(132.0)
result = result + term
term = term * x2 / t.CFloat(182.0)
result = result - term
term = term * x2 / t.CFloat(240.0)
result = result + term
return result
def tanf(x: t.CFloat) -> t.CFloat:
c: t.CFloat = cosf(x)
if c == t.CFloat(0.0):
return t.CFloat(0.0)
return sinf(x) / c
def fabsf(x: t.CFloat) -> t.CFloat:
return x if x >= t.CFloat(0.0) else -x
def floorf(x: t.CFloat) -> t.CFloat:
if x >= t.CFloat(0.0):
return t.CFloat(t.CInt(x))
i: t.CInt = t.CInt(x)
if x < t.CFloat(i):
return t.CFloat(i - 1)
else:
return t.CFloat(i)
def ceilf(x: t.CFloat) -> t.CFloat:
if x < t.CFloat(0.0):
return t.CFloat(t.CInt(x))
i: t.CInt = t.CInt(x)
if x > t.CFloat(i):
return t.CFloat(i + 1)
else:
return t.CFloat(i)
class _U(t.CUnion):
d: t.CDouble
u: t.CUInt64T
def isnan(x: t.CDouble) -> t.CStatic | t.CInt:
# IEEE-754 位重解释(不能用 t.CUInt64T(x) 数值转换:那是 fptosi,位模式全丢
# → exp_part 恒 0 → 永远返回 0)。用 c.Load 按位拷贝(viperlib
# _check_special_float 同款惯用法)。
xv: t.CDouble = x
u: t.CUInt64T = 0
c.Load(c.Addr(u), c.Addr(xv))
exp_part: t.CUInt64T = (u >> 52) & 0x7FF
mantissa_part: t.CUInt64T = u & 0xFFFFFFFFFFFFF
if exp_part == 0x7FF:
if mantissa_part != 0:
return 1
return 0
def isinf(x: t.CDouble) -> t.CStatic | t.CInt:
xv2: t.CDouble = x
u2: t.CUInt64T = 0
c.Load(c.Addr(u2), c.Addr(xv2))
exp_part2: t.CUInt64T = (u2 >> 52) & 0x7FF
mantissa_part2: t.CUInt64T = u2 & 0xFFFFFFFFFFFFF
if exp_part2 == 0x7FF:
if mantissa_part2 == 0:
return 1
return 0