module math

// gamma function computed by Stirling's formula.
// The pair of results must be multiplied together to get the actual answer.
// The multiplication is left to the caller so that, if careful, the caller can avoid
// infinity for 172 <= x <= 180.
// The polynomial is valid for 33 <= x <= 172 larger values are only used
// in reciprocal and produce denormalized floats. The lower precision there
// masks any imprecision in the polynomial.
fn stirling(x f64) (f64, f64) {
    if x > 200 {
        return inf(1), 1.0
    }
    sqrt_two_pi := 2.506628274631000502417
    max_stirling := 143.01608
    mut w := 1.0 / x
    w = 1.0 + w * ((((gamma_s[0] * w + gamma_s[1]) * w + gamma_s[2]) * w + gamma_s[3]) * w +
        gamma_s[4])
    mut y1 := exp(x)
    mut y2 := 1.0
    if x > max_stirling { // avoid Pow() overflow
        v := pow(x, 0.5 * x - 0.25)
        y1_ := y1
        y1 = v
        y2 = v / y1_
    } else {
        y1 = pow(x, x - 0.5) / y1
    }
    return y1, f64(sqrt_two_pi) * w * y2
}

// gamma returns the gamma function of x.
//
// special ifs are:
// gamma(+inf) = +inf
// gamma(+0) = +inf
// gamma(-0) = -inf
// gamma(x) = nan for integer x < 0
// gamma(-inf) = nan
// gamma(nan) = nan
pub fn gamma(a f64) f64 {
    mut x := a
    euler := 0.57721566490153286060651209008240243104215933593992 // A001620
    if is_neg_int(x) || is_inf(x, -1) || is_nan(x) {
        return nan()
    }
    if is_inf(x, 1) {
        return inf(1)
    }
    if x == 0.0 {
        return copysign(inf(1), x)
    }
    mut q := abs(x)
    mut p := floor(q)
    if q > 33 {
        if x >= 0 {
            y1, y2 := stirling(x)
            return y1 * y2
        }
        // Note: x is negative but (checked above) not a negative integer,
        // so x must be small enough to be in range for conversion to i64.
        // If |x| were >= 2⁶³ it would have to be an integer.
        mut signgam := 1
        ip := i64(p)
        if (ip & 1) == 0 {
            signgam = -1
        }
        mut z := q - p
        if z > 0.5 {
            p = p + 1
            z = q - p
        }
        z = q * sin(pi * z)
        if z == 0 {
            return inf(signgam)
        }
        sq1, sq2 := stirling(q)
        absz := abs(z)
        d := absz * sq1 * sq2
        if is_inf(d, 0) {
            z = pi / absz / sq1 / sq2
        } else {
            z = pi / d
        }
        return f64(signgam) * z
    }
    // Reduce argument
    mut z := 1.0
    for x >= 3 {
        x = x - 1
        z = z * x
    }
    for x < 0 {
        if x > -1e-09 {
            unsafe {
                goto small
            }
        }
        z = z / x
        x = x + 1
    }
    for x < 2 {
        if x < 1e-09 {
            unsafe {
                goto small
            }
        }
        z = z / x
        x = x + 1
    }
    if x == 2 {
        return z
    }
    x = x - 2
    p = (((((x * gamma_p[0] + gamma_p[1]) * x + gamma_p[2]) * x + gamma_p[3]) * x +
        gamma_p[4]) * x + gamma_p[5]) * x + gamma_p[6]
    q = ((((((x * gamma_q[0] + gamma_q[1]) * x + gamma_q[2]) * x + gamma_q[3]) * x +
        gamma_q[4]) * x + gamma_q[5]) * x + gamma_q[6]) * x + gamma_q[7]
    if true {
        return z * p / q
    }
    small:
    if x == 0 {
        return inf(1)
    }
    return z / ((1.0 + euler * x) * x)
}

// log_gamma returns the natural logarithm and sign (-1 or +1) of Gamma(x).
//
// special ifs are:
// log_gamma(+inf) = +inf
// log_gamma(0) = +inf
// log_gamma(-integer) = +inf
// log_gamma(-inf) = -inf
// log_gamma(nan) = nan
pub fn log_gamma(x f64) f64 {
    y, _ := log_gamma_sign(x)
    return y
}

pub fn log_gamma_sign(a f64) (f64, int) {
    mut x := a
    ymin := 1.461632144968362245
    tiny := exp2(-70)
    two52 := exp2(52) // 0x4330000000000000 ~4.5036e+15
    two58 := exp2(58) // 0x4390000000000000 ~2.8823e+17
    tc := 1.46163214496836224576e+00 // 0x3FF762D86356BE3F
    tf := -1.21486290535849611461e-01 // 0xBFBF19B9BCC38A42
    // tt := -(tail of tf)
    tt := -3.63867699703950536541e-18 // 0xBC50C7CAA48A971F
    mut sign := 1
    if is_nan(x) {
        return x, sign
    }
    if is_inf(x, 1) {
        return x, sign
    }
    if x == 0.0 {
        return inf(1), sign
    }
    mut neg := false
    if x < 0 {
        x = -x
        neg = true
    }
    if x < tiny { // if |x| < 2**-70, return -log(|x|)
        if neg {
            sign = -1
        }
        return -log(x), sign
    }
    mut nadj := 0.0
    if neg {
        if x >= two52 {
            // x| >= 2**52, must be -integer
            return inf(1), sign
        }
        t := sin_pi(x)
        if t == 0 {
            return inf(1), sign
        }
        nadj = log(pi / abs(t * x))
        if t < 0 {
            sign = -1
        }
    }
    mut lgamma := 0.0
    if x == 1 || x == 2 { // purge off 1 and 2
        return 0.0, sign
    } else if x < 2 { // use lgamma(x) = lgamma(x+1) - log(x)
        mut y := 0.0
        mut i := 0
        if x <= 0.9 {
            lgamma = -log(x)
            if x >= (ymin - 1 + 0.27) { // 0.7316 <= x <=  0.9
                y = 1.0 - x
                i = 0
            } else if x >= (ymin - 1 - 0.27) { // 0.2316 <= x < 0.7316
                y = x - (tc - 1)
                i = 1
            } else { // 0 < x < 0.2316
                y = x
                i = 2
            }
        } else {
            lgamma = 0
            if x >= (ymin + 0.27) { // 1.7316 <= x < 2
                y = f64(2) - x
                i = 0
            } else if x >= (ymin - 0.27) { // 1.2316 <= x < 1.7316
                y = x - tc
                i = 1
            } else { // 0.9 < x < 1.2316
                y = x - 1
                i = 2
            }
        }
        if i == 0 {
            z := y * y
            gamma_p1 := lgamma_a[0] + z * (lgamma_a[2] + z * (lgamma_a[4] + z * (lgamma_a[6] +
                z * (lgamma_a[8] + z * lgamma_a[10]))))
            gamma_p2 := z * (lgamma_a[1] + z * (lgamma_a[3] + z * (lgamma_a[5] + z * (lgamma_a[7] +
                z * (lgamma_a[9] + z * lgamma_a[11])))))
            p := y * gamma_p1 + gamma_p2
            lgamma += (p - 0.5 * y)
        } else if i == 1 {
            z := y * y
            w := z * y
            gamma_p1 := lgamma_t[0] + w * (lgamma_t[3] + w * (lgamma_t[6] + w * (lgamma_t[9] +
                w * lgamma_t[12]))) // parallel comp
            gamma_p2 := lgamma_t[1] + w * (lgamma_t[4] + w * (lgamma_t[7] + w * (lgamma_t[10] +
                w * lgamma_t[13])))
            gamma_p3 := lgamma_t[2] + w * (lgamma_t[5] + w * (lgamma_t[8] + w * (lgamma_t[11] +
                w * lgamma_t[14])))
            p := z * gamma_p1 - (tt - w * (gamma_p2 + y * gamma_p3))
            lgamma += (tf + p)
        } else if i == 2 {
            gamma_p1 := y * (lgamma_u[0] + y * (lgamma_u[1] + y * (lgamma_u[2] + y * (lgamma_u[3] +
                y * (lgamma_u[4] + y * lgamma_u[5])))))
            gamma_p2 := 1.0 + y * (lgamma_v[1] + y * (lgamma_v[2] + y * (lgamma_v[3] +
                y * (lgamma_v[4] + y * lgamma_v[5]))))
            lgamma += (-0.5 * y + gamma_p1 / gamma_p2)
        }
    } else if x < 8 { // 2 <= x < 8
        i := int(x)
        y := x - f64(i)
        p := y * (lgamma_s[0] + y * (lgamma_s[1] + y * (lgamma_s[2] + y * (lgamma_s[3] +
            y * (lgamma_s[4] + y * (lgamma_s[5] + y * lgamma_s[6]))))))
        q := 1.0 + y * (lgamma_r[1] + y * (lgamma_r[2] + y * (lgamma_r[3] + y * (lgamma_r[4] +
            y * (lgamma_r[5] + y * lgamma_r[6])))))
        lgamma = 0.5 * y + p / q
        mut z := 1.0 // lgamma(1+s) = log(s) + lgamma(s)
        if i == 7 {
            z *= (y + 6)
            z *= (y + 5)
            z *= (y + 4)
            z *= (y + 3)
            z *= (y + 2)
            lgamma += log(z)
        } else if i == 6 {
            z *= (y + 5)
            z *= (y + 4)
            z *= (y + 3)
            z *= (y + 2)
            lgamma += log(z)
        } else if i == 5 {
            z *= (y + 4)
            z *= (y + 3)
            z *= (y + 2)
            lgamma += log(z)
        } else if i == 4 {
            z *= (y + 3)
            z *= (y + 2)
            lgamma += log(z)
        } else if i == 3 {
            z *= (y + 2)
            lgamma += log(z)
        }
    } else if x < two58 { // 8 <= x < 2**58
        t := log(x)
        z := 1.0 / x
        y := z * z
        w := lgamma_w[0] + z * (lgamma_w[1] + y * (lgamma_w[2] + y * (lgamma_w[3] +
            y * (lgamma_w[4] + y * (lgamma_w[5] + y * lgamma_w[6])))))
        lgamma = (x - 0.5) * (t - 1.0) + w
    } else { // 2**58 <= x <= Inf
        lgamma = x * (log(x) - 1.0)
    }
    if neg {
        lgamma = nadj - lgamma
    }
    return lgamma, sign
}

// sin_pi(x) is a helper function for negative x
fn sin_pi(x_ f64) f64 {
    mut x := x_
    two52 := exp2(52) // 0x4330000000000000 ~4.5036e+15
    two53 := exp2(53) // 0x4340000000000000 ~9.0072e+15
    if x < 0.25 {
        return -sin(pi * x)
    }
    // argument reduction
    mut z := floor(x)
    mut n := 0
    if z != x { // inexact
        x = mod(x, 2)
        n = int(x * 4)
    } else {
        if x >= two53 { // x must be even
            x = 0
            n = 0
        } else {
            if x < two52 {
                z = x + two52 // exact
            }
            n = 1 & int(f64_bits(z))
            x = f64(n)
            n <<= 2
        }
    }
    if n == 0 {
        x = sin(pi * x)
    } else if n == 1 || n == 2 {
        x = cos(pi * (0.5 - x))
    } else if n == 3 || n == 4 {
        x = sin(pi * (1.0 - x))
    } else if n == 5 || n == 6 {
        x = -cos(pi * (x - 1.5))
    } else {
        x = sin(pi * (x - 2))
    }
    return -x
}