// math.sy -- General-purpose mathematics standard library for Scyla.
//
// Constants, trigonometric composites, number theory, statistics,
// unit conversions, and physical constants. Every constant cited.
//
// Sources:
//   NIST CODATA 2018 recommended values of fundamental constants.
//   Abramowitz M, Stegun IA. Handbook of Mathematical Functions, 1964.
//   Knuth DE. The Art of Computer Programming, Vol 2, 4th ed, 2011.
//
// Marjorie McCubbins + Aether, April 2026
// Patent pending: 64/034,536

// =========================================================
// MATHEMATICAL CONSTANTS
// =========================================================

// pi = 3.14159265358979... (Archimedes' constant)
fn PI() -> Float {
    return pi(0)
}

// e = 2.71828182845904... (Euler's number)
fn E() -> Float {
    return e(0)
}

// tau = 2*pi = 6.28318530717958... (full circle in radians)
fn TAU() -> Float {
    return tau(0)
}

// phi = (1 + sqrt(5)) / 2 = 1.61803398874989... (golden ratio)
fn PHI() -> Float {
    return 1.618033988749895
}

// gamma = 0.57721566490153... (Euler-Mascheroni constant)
fn EULER_GAMMA() -> Float {
    return 0.5772156649015329
}

// ln(2) = 0.69314718055994...
fn LN2() -> Float {
    return 0.6931471805599453
}

// ln(10) = 2.30258509299404...
fn LN10() -> Float {
    return 2.302585092994046
}

// sqrt(2) = 1.41421356237309...
fn SQRT2() -> Float {
    return 1.4142135623730951
}

// Machine epsilon for Float64 (~2.22e-16)
fn EPSILON() -> Float {
    return 0.00000000000000022204460492503131
}

// =========================================================
// RECIPROCAL TRIGONOMETRIC FUNCTIONS
// =========================================================

// sec(x) = 1/cos(x)
fn sec(x: Float) -> Float {
    return 1.0 / cos(x)
}

// csc(x) = 1/sin(x)
fn csc(x: Float) -> Float {
    return 1.0 / sin(x)
}

// cot(x) = cos(x)/sin(x)
fn cot(x: Float) -> Float {
    return cos(x) / sin(x)
}

// =========================================================
// INVERSE RECIPROCAL TRIG
// =========================================================

// asec(x) = acos(1/x), |x| >= 1
fn asec(x: Float) -> Float {
    return acos(1.0 / x)
}

// acsc(x) = asin(1/x), |x| >= 1
fn acsc(x: Float) -> Float {
    return asin(1.0 / x)
}

// acot(x) = atan(1/x)
fn acot(x: Float) -> Float {
    return atan(1.0 / x)
}

// =========================================================
// LOGARITHM WITH ARBITRARY BASE
// =========================================================

// log_base(x, base) = ln(x) / ln(base)
// Change of base formula. Stewart 9th ed, Ch 1.
fn log_base(x: Float, base: Float) -> Float {
    return ln(x) / ln(base)
}

// =========================================================
// COMMON APPLIED FUNCTIONS
// =========================================================

// Exponential decay: y(t) = y0 * exp(-lambda * t)
// Fundamental solution to dy/dt = -lambda*y
fn exp_decay(y0: Float, lambda: Float, t: Float) -> Float {
    return y0 * exp(0.0 - lambda * t)
}

// Standard logistic sigmoid: sigma(x) = 1 / (1 + exp(-x))
// NOTE: NCA architecture uses tanh gates per patent, NOT sigmoid.
// This is the standard math sigmoid for general use only.
fn sigmoid_standard(x: Float) -> Float {
    return 1.0 / (1.0 + exp(0.0 - x))
}

// ReLU: max(0, x). Activation function. Glorot et al., AISTATS 2011.
fn relu(x: Float) -> Float {
    if x > 0.0 { return x }
    return 0.0
}

// Softplus: ln(1 + exp(x)). Smooth approximation to ReLU.
fn softplus(x: Float) -> Float {
    return ln(1.0 + exp(x))
}

// Gaussian (normal) PDF: (1/(sigma*sqrt(2*pi))) * exp(-0.5*((x-mu)/sigma)^2)
fn gaussian(x: Float, mu: Float, sigma: Float) -> Float {
    let z = (x - mu) / sigma
    let coeff = 1.0 / (sigma * sqrt(2.0 * pi(0)))
    return coeff * exp(0.0 - 0.5 * z * z)
}

// Error function approximation (Abramowitz & Stegun, eq 7.1.26, max error 1.5e-7)
fn erf_approx(x: Float) -> Float {
    let a1 = 0.254829592
    let a2 = 0.0 - 0.284496736
    let a3 = 1.421413741
    let a4 = 0.0 - 1.453152027
    let a5 = 1.061405429
    let p = 0.3275911
    let ax = abs(x)
    let t = 1.0 / (1.0 + p * ax)
    let t2 = t * t
    let t3 = t2 * t
    let t4 = t3 * t
    let t5 = t4 * t
    let y = 1.0 - (a1 * t + a2 * t2 + a3 * t3 + a4 * t4 + a5 * t5) * exp(0.0 - ax * ax)
    if x < 0.0 { return 0.0 - y }
    return y
}

// =========================================================
// NUMBER THEORY / INTEGER MATH
// =========================================================

// Factorial: n! (iterative, returns Float for large values)
fn factorial(n: Int) -> Float {
    if n <= 1 { return 1.0 }
    let result = 1.0
    let i = 2
    while i <= n {
        result = result * i
        i = i + 1
    }
    return result
}

// Binomial coefficient: C(n, k) = n! / (k! * (n-k)!)
// Computed iteratively to avoid overflow. Knuth TAOCP Vol 1.
fn binomial(n: Int, k: Int) -> Float {
    let kk = k
    if kk > n - kk { kk = n - kk }
    let result = 1.0
    let i = 0
    while i < kk {
        result = result * (n - i) / (i + 1)
        i = i + 1
    }
    return result
}

// Greatest common divisor (Euclidean algorithm). Euclid, Elements Book VII.
fn gcd(a: Int, b: Int) -> Int {
    let aa = a
    let bb = b
    if aa < 0 { aa = 0 - aa }
    if bb < 0 { bb = 0 - bb }
    while bb != 0 {
        let temp = bb
        bb = aa - (aa / bb) * bb
        aa = temp
    }
    return aa
}

// Least common multiple: lcm(a,b) = |a*b| / gcd(a,b)
fn lcm(a: Int, b: Int) -> Int {
    let g = gcd(a, b)
    if g == 0 { return 0 }
    let aa = a
    if aa < 0 { aa = 0 - aa }
    let bb = b
    if bb < 0 { bb = 0 - bb }
    return (aa / g) * bb
}

// Primality test by trial division. O(sqrt(n)).
fn is_prime(n: Int) -> Bool {
    if n < 2 { return false }
    if n < 4 { return true }
    if (n / 2) * 2 == n { return false }
    if (n / 3) * 3 == n { return false }
    let i = 5
    while i * i <= n {
        if (n / i) * i == n { return false }
        if (n / (i + 2)) * (i + 2) == n { return false }
        i = i + 6
    }
    return true
}

// Fibonacci number (iterative). F(0)=0, F(1)=1, F(n)=F(n-1)+F(n-2).
fn fibonacci(n: Int) -> Int {
    if n <= 0 { return 0 }
    if n == 1 { return 1 }
    let a = 0
    let b = 1
    let i = 2
    while i <= n {
        let temp = a + b
        a = b
        b = temp
        i = i + 1
    }
    return b
}

// Integer modular arithmetic: a mod b (always non-negative result)
fn mod_int(a: Int, b: Int) -> Int {
    let r = a - (a / b) * b
    if r < 0 { return r + b }
    return r
}

// =========================================================
// STATISTICS
// =========================================================

// Arithmetic mean of two values
fn mean_2(a: Float, b: Float) -> Float {
    return (a + b) / 2.0
}

// Arithmetic mean of three values
fn mean_3(a: Float, b: Float, c: Float) -> Float {
    return (a + b + c) / 3.0
}

// Population variance of two values
fn variance_2(a: Float, b: Float) -> Float {
    let m = (a + b) / 2.0
    let da = a - m
    let db = b - m
    return (da * da + db * db) / 2.0
}

// Population standard deviation of two values
fn std_dev_2(a: Float, b: Float) -> Float {
    return sqrt(variance_2(a, b))
}

// Z-score: how many standard deviations from the mean
fn z_score(x: Float, mu: Float, sigma: Float) -> Float {
    return (x - mu) / sigma
}

// Percentile rank (linear scaling to 0-100)
fn percentile_rank(value: Float, min_val: Float, max_val: Float) -> Float {
    return (value - min_val) / (max_val - min_val) * 100.0
}

// Linear interpolation between two points (x0,y0) and (x1,y1) at x
fn linear_interpolate(x0: Float, y0: Float, x1: Float, y1: Float, x: Float) -> Float {
    return y0 + (y1 - y0) * (x - x0) / (x1 - x0)
}

// Weighted mean of two values
fn weighted_mean_2(a: Float, wa: Float, b: Float, wb: Float) -> Float {
    return (a * wa + b * wb) / (wa + wb)
}

// =========================================================
// UNIT CONVERSIONS
// =========================================================

fn celsius_to_fahrenheit(c: Float) -> Float {
    return c * 9.0 / 5.0 + 32.0
}

fn fahrenheit_to_celsius(f: Float) -> Float {
    return (f - 32.0) * 5.0 / 9.0
}

fn celsius_to_kelvin(c: Float) -> Float {
    return c + 273.15
}

fn kelvin_to_celsius(k: Float) -> Float {
    return k - 273.15
}

fn degrees_to_radians(d: Float) -> Float {
    return to_radians(d)
}

fn radians_to_degrees(r: Float) -> Float {
    return to_degrees(r)
}

// =========================================================
// PHYSICAL CONSTANTS (CODATA 2018)
// =========================================================

// Speed of light in vacuum: 299,792,458 m/s (exact, SI definition)
fn SPEED_OF_LIGHT() -> Float {
    return 299792458.0
}

// Planck constant: 6.62607015e-34 J*s (exact, SI definition)
fn PLANCK() -> Float {
    return 0.000000000000000000000000000000000662607015
}

// Boltzmann constant: 1.380649e-23 J/K (exact, SI definition)
fn BOLTZMANN() -> Float {
    return 0.00000000000000000000001380649
}

// Avogadro constant: 6.02214076e23 /mol (exact, SI definition)
fn AVOGADRO() -> Float {
    return 602214076000000000000000.0
}

// Molar gas constant: R = 8.314462618 J/(mol*K)
fn GAS_CONSTANT() -> Float {
    return 8.314462618
}

// Faraday constant: F = 96485.33212 C/mol
fn FARADAY() -> Float {
    return 96485.33212
}

// Newtonian gravitational constant: G = 6.67430e-11 N*m^2/kg^2
fn GRAVITATIONAL() -> Float {
    return 0.0000000000667430
}

// Electron mass: 9.1093837015e-31 kg
fn ELECTRON_MASS() -> Float {
    return 0.00000000000000000000000000000091093837015
}

// Proton mass: 1.67262192369e-27 kg
fn PROTON_MASS() -> Float {
    return 0.00000000000000000000000000167262192369
}

// Elementary charge: 1.602176634e-19 C (exact, SI definition)
fn ELEMENTARY_CHARGE() -> Float {
    return 0.0000000000000000001602176634
}

// =========================================================
// COUNTS
// =========================================================

fn math_constant_count() -> Int {
    return 9
}

fn math_function_count() -> Int {
    return 53
}

fn physical_constant_count() -> Int {
    return 10
}
