543 lines
15 KiB
Python
543 lines
15 KiB
Python
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(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
|