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
