codon/stdlib/math.codon

e = 2.718281828459045
pi = 3.141592653589793
tau = 6.283185307179586

inf = 1.0 / 0.0
nan = 0.0 / 0.0

def factorial(x: int) -> int:
    _F = (1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880,  3628800,  39916800,  479001600,  6227020800,  87178291200,  1307674368000,  20922789888000,  355687428096000,  6402373705728000,  121645100408832000,  2432902008176640000)
    if not (0 <= x <= 20):
        raise ValueError("factorial is only supported for 0 <= x <= 20")
    return _F[x]

def isnan(x: float) -> bool:
    """
    isnan(float) -> bool

    Return True if float arg is a NaN, else False.
    """
    test = x == x
    # if it is true then it is a number
    if test:
        return False
    return True

def isinf(x: float) -> bool:
    """
    isinf(float) -> bool:

    Return True if float arg is an INF, else False.
    """
    return x == inf or x == -inf

def isfinite(x: float) -> bool:
    """
    isfinite(float) -> bool

    Return True if x is neither an infinity nor a NaN,
    and False otherwise.
    """
    if isnan(x) or isinf(x):
        return False
    return True

def ceil(x: float) -> float:
    """
    ceil(float) -> float

    Return the ceiling of x as an Integral.
    This is the smallest integer >= x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.ceil.f64(double)
        %y = call double @llvm.ceil.f64(double %x)
        ret double %y

    return f(x)

def floor(x: float) -> float:
    """
    floor(float) -> float

    Return the floor of x as an Integral.
    This is the largest integer <= x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.floor.f64(double)
        %y = call double @llvm.floor.f64(double %x)
        ret double %y

    return f(x)

def fabs(x: float) -> float:
    """
    fabs(float) -> float

    Returns the absolute value of a floating point number.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.fabs.f64(double)
        %y = call double @llvm.fabs.f64(double %x)
        ret double %y

    return f(x)

def fmod(x: float, y: float) -> float:
    """
    fmod(float, float) -> float

    Returns the remainder of x divided by y.
    """
    @pure
    @llvm
    def f(x: float, y: float) -> float:
        %z = frem double %x, %y
        ret double %z

    return f(x, y)

def exp(x: float) -> float:
    """
    exp(float) -> float

    Returns the value of e raised to the xth power.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.exp.f64(double)
        %y = call double @llvm.exp.f64(double %x)
        ret double %y

    return f(x)

def expm1(x: float) -> float:
    """
    expm1(float) -> float

    Return e raised to the power x, minus 1. expm1 provides
    a way to compute this quantity to full precision.
    """
    return _C.expm1(x)

def ldexp(x: float, i: int) -> float:
    """
    ldexp(float, int) -> float

    Returns x multiplied by 2 raised to the power of exponent.
    """
    return _C.ldexp(x, i)

def log(x: float, base: float = e) -> float:
    """
    log(float) -> float

    Returns the natural logarithm (base-e logarithm) of x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.log.f64(double)
        %y = call double @llvm.log.f64(double %x)
        ret double %y

    if base == e:
        return f(x)
    else:
        return f(x)/f(base)

def log2(x: float) -> float:
    """
    log2(float) -> float

    Return the base-2 logarithm of x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.log2.f64(double)
        %y = call double @llvm.log2.f64(double %x)
        ret double %y

    return f(x)

def log10(x: float) -> float:
    """
    log10(float) -> float

    Returns the common logarithm (base-10 logarithm) of x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.log10.f64(double)
        %y = call double @llvm.log10.f64(double %x)
        ret double %y

    return f(x)

def degrees(x: float) -> float:
    """
    degrees(float) -> float

    Convert angle x from radians to degrees.
    """
    radToDeg = 180.0/pi
    return x * radToDeg

def radians(x: float) -> float:
    """
    radians(float) -> float

    Convert angle x from degrees to radians.
    """
    degToRad = pi/180.0
    return x * degToRad

def sqrt(x: float) -> float:
    """
    sqrt(float) -> float

    Returns the square root of x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.sqrt.f64(double)
        %y = call double @llvm.sqrt.f64(double %x)
        ret double %y

    return f(x)

def pow(x: float, y: float) -> float:
    """
    pow(float, float) -> float

    Returns x raised to the power of y.
    """
    @pure
    @llvm
    def f(x: float, y: float) -> float:
        declare double @llvm.pow.f64(double, double)
        %z = call double @llvm.pow.f64(double %x, double %y)
        ret double %z

    return f(x, y)

def acos(x: float) -> float:
    """
    acos(float) -> float

    Returns the arc cosine of x in radians.
    """
    return _C.acos(x)

def asin(x: float) -> float:
    """
    asin(float) -> float

    Returns the arc sine of x in radians.
    """
    return _C.asin(x)

def atan(x: float) -> float:
    """
    atan(float) -> float

    Returns the arc tangent of x in radians.
    """
    return _C.atan(x)

def atan2(y: float, x: float) -> float:
    """
    atan2(float, float) -> float

    Returns the arc tangent in radians of y/x based
    on the signs of both values to determine the
    correct quadrant.
    """
    return _C.atan2(y, x)

def cos(x: float) -> float:
    """
    cos(float) -> float

    Returns the cosine of a radian angle x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.cos.f64(double)
        %y = call double @llvm.cos.f64(double %x)
        ret double %y

    return f(x)

def sin(x: float) -> float:
    """
    sin(float) -> float

    Returns the sine of a radian angle x.
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.sin.f64(double)
        %y = call double @llvm.sin.f64(double %x)
        ret double %y

    return f(x)

def hypot(x: float, y: float) -> float:
    """
    hypot(float, float) -> float

    Return the Euclidean norm.
    This is the length of the vector from the
    origin to point (x, y).
    """
    return _C.hypot(x, y)

def tan(x: float) -> float:
    """
    tan(float) -> float

    Return the tangent of a radian angle x.
    """
    return _C.tan(x)

def cosh(x: float) -> float:
    """
    cosh(float) -> float

    Returns the hyperbolic cosine of x.
    """
    return _C.cosh(x)

def sinh(x: float) -> float:
    """
    sinh(float) -> float

    Returns the hyperbolic sine of x.
    """
    return _C.sinh(x)

def tanh(x: float) -> float:
    """
    tanh(float) -> float

    Returns the hyperbolic tangent of x.
    """
    return _C.tanh(x)

def acosh(x: float) -> float:
    """
    acosh(float) -> float

    Return the inverse hyperbolic cosine of x.
    """
    return _C.acosh(x)

def asinh(x: float) -> float:
    """
    asinh(float) -> float

    Return the inverse hyperbolic sine of x.
    """
    return _C.asinh(x)

def atanh(x: float) -> float:
    """
    atanh(float) -> float

    Return the inverse hyperbolic tangent of x.
    """
    return _C.atanh(x)

def copysign(x: float, y: float) -> float:
    """
    copysign(float, float) -> float

    Return a float with the magnitude (absolute value) of
    x but the sign of y.
    """
    @pure
    @llvm
    def f(x: float, y: float) -> float:
        declare double @llvm.copysign.f64(double, double)
        %z = call double @llvm.copysign.f64(double %x, double %y)
        ret double %z

    return f(x, y)

def log1p(x: float) -> float:
    """
    log1p(float) -> float

    Return the natural logarithm of 1+x (base e).
    """
    return _C.log1p(x)

def trunc(x: float) -> float:
    """
    trunc(float) -> float

    Return the Real value x truncated to an Integral
    (usually an integer).
    """
    @pure
    @llvm
    def f(x: float) -> float:
        declare double @llvm.trunc.f64(double)
        %y = call double @llvm.trunc.f64(double %x)
        ret double %y

    return f(x)

def erf(x: float) -> float:
    """
    erf(float) -> float

    Return the error function at x.
    """
    return _C.erf(x)

def erfc(x: float) -> float:
    """
    erfc(float) -> float

    Return the complementary error function at x.
    """
    return _C.erfc(x)

def gamma(x: float) -> float:
    """
    gamma(float) -> float

    Return the Gamma function at x.
    """
    return _C.tgamma(x)

def lgamma(x: float) -> float:
    """
    lgamma(float) -> float

    Return the natural logarithm of
    the absolute value of the Gamma function at x.
    """
    return _C.lgamma(x)

def remainder(x: float, y: float) -> float:
    """
    remainder(float, float) -> float

    Return the IEEE 754-style remainder of x with respect to y.
    For finite x and finite nonzero y, this is the difference
    x - n*y, where n is the closest integer to the exact value
    of the quotient x / y. If x / y is exactly halfway between
    two consecutive integers, the nearest even integer is used
    for n.
    """
    return _C.remainder(x, y)

def gcd(a: float, b: float) -> float:
    """
    gcd(float, float) -> float

    returns greatest common divisor of x and y.
    """
    a = abs(a)
    b = abs(b)
    while a:
        a, b = b%a, a
    return b

@pure
def frexp(x: float) -> Tuple[float, int]:
    """
    frexp(float) -> Tuple[float, int]

    The returned value is the mantissa and the integer pointed
    to by exponent is the exponent. The resultant value is
    x = mantissa * 2 ^ exponent.
    """
    tmp = i32(0)
    res = _C.frexp(float(x), __ptr__(tmp))
    return (res, int(tmp))

@pure
def modf(x: float) -> Tuple[float, float]:
    """
    modf(float) -> Tuple[float, float]

    The returned value is the fraction component (part after
    the decimal), and sets integer to the integer component.
    """
    tmp = 0.0
    res = _C.modf(float(x), __ptr__(tmp))
    return (res, tmp)

def isclose(a: float, b: float, rel_tol: float = 1e-09, abs_tol: float = 0.0) -> bool:
    """
    isclose(float, float) -> bool

    Return True if a is close in value to b, and False otherwise.
    For the values to be considered close, the difference between them
    must be smaller than at least one of the tolerances.
    """

    # short circuit exact equality -- needed to catch two
    # infinities of the same sign. And perhaps speeds things
    # up a bit sometimes.
    if a == b:
        return True

    # This catches the case of two infinities of opposite sign, or
    # one infinity and one finite number. Two infinities of opposite
    # sign would otherwise have an infinite relative tolerance.
    # Two infinities of the same sign are caught by the equality check
    # above.
    if a == inf or b == inf:
        return False

    # NAN is not close to anything, not even itself
    if a == nan or b == nan:
        return False

    # regular computation
    diff = fabs(b - a)

    return (((diff <= fabs(rel_tol * b)) or
            (diff <= fabs(rel_tol * a))) or
            (diff <= abs_tol))