545 lines
15 KiB
Plaintext
545 lines
15 KiB
Plaintext
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:
|
||
if x > t.CDouble(709.0):
|
||
return t.CDouble(1.7976931348623157e308)
|
||
if x < t.CDouble(-709.0):
|
||
return t.CDouble(0.0)
|
||
result: t.CDouble = t.CDouble(1.0)
|
||
term: t.CDouble = t.CDouble(1.0)
|
||
term = term * x / t.CDouble(1.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(2.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(3.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(4.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(5.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(6.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(7.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(8.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(9.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(10.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(11.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(12.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(13.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(14.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(15.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(16.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(17.0)
|
||
result = result + term
|
||
term = term * x / t.CDouble(18.0)
|
||
result = result + term
|
||
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:
|
||
u: t.CUInt64T = t.CUInt64T(x)
|
||
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:
|
||
u: t.CUInt64T = t.CUInt64T(x)
|
||
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
|