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