// calculus.sy -- Differential, integral, and vector calculus as native Scyla.
//
// Numerical methods (finite differences, quadrature) and analytical
// derivatives/integrals for common functions. Every formula cited.
//
// Sources:
//   Stewart J. Calculus: Early Transcendentals, 9th ed, Cengage, 2020.
//   Burden RL, Faires JD. Numerical Analysis, 10th ed, Cengage, 2016.
//
// Marjorie McCubbins + Aether, April 2026
// Patent pending: 64/034,536

// =========================================================
// DIFFERENTIAL CALCULUS -- Numerical Derivatives
// =========================================================

// Forward difference: f'(x) ~ (f(x+h) - f(x)) / h
// First-order accurate, O(h). Burden 10th ed, Sec 4.1.
fn derivative(f_at_x: Float, f_at_x_plus_h: Float, h: Float) -> Float {
    return (f_at_x_plus_h - f_at_x) / h
}

// Central difference: f'(x) ~ (f(x+h) - f(x-h)) / (2h)
// Second-order accurate, O(h^2). Burden 10th ed, Sec 4.1.
fn derivative_central(f_at_x_minus_h: Float, f_at_x_plus_h: Float, h: Float) -> Float {
    return (f_at_x_plus_h - f_at_x_minus_h) / (2.0 * h)
}

// Second derivative: f''(x) ~ (f(x+h) - 2f(x) + f(x-h)) / h^2
// Second-order accurate, O(h^2). Burden 10th ed, Sec 4.1.
fn derivative_second(f_at_x_minus_h: Float, f_at_x: Float, f_at_x_plus_h: Float, h: Float) -> Float {
    return (f_at_x_plus_h - 2.0 * f_at_x + f_at_x_minus_h) / (h * h)
}

// Partial derivative with respect to x (holding y constant)
fn partial_x(f_xy: Float, f_x_plus_h_y: Float, h: Float) -> Float {
    return (f_x_plus_h_y - f_xy) / h
}

// Partial derivative with respect to y (holding x constant)
fn partial_y(f_xy: Float, f_x_y_plus_h: Float, h: Float) -> Float {
    return (f_x_y_plus_h - f_xy) / h
}

// Gradient x-component via central difference
fn gradient_2d_x(f_x_minus: Float, f_x_plus: Float, h: Float) -> Float {
    return (f_x_plus - f_x_minus) / (2.0 * h)
}

// Gradient y-component via central difference
fn gradient_2d_y(f_y_minus: Float, f_y_plus: Float, h: Float) -> Float {
    return (f_y_plus - f_y_minus) / (2.0 * h)
}

// Five-point stencil for first derivative: O(h^4) accuracy
// f'(x) ~ (-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)) / (12h)
// Burden 10th ed, Sec 4.1.
fn derivative_five_point(f_m2h: Float, f_mh: Float, f_ph: Float, f_p2h: Float, h: Float) -> Float {
    return (0.0 - f_p2h + 8.0 * f_ph - 8.0 * f_mh + f_m2h) / (12.0 * h)
}

// =========================================================
// INTEGRAL CALCULUS -- Numerical Quadrature
// =========================================================

// Trapezoid rule (single panel): integral ~ h * (f(a) + f(b)) / 2
// Error: -h^3/12 * f''(xi). Burden 10th ed, Sec 4.3.
fn trapezoid(f_a: Float, f_b: Float, h: Float) -> Float {
    return h * (f_a + f_b) / 2.0
}

// Simpson's 1/3 rule: integral ~ (h/3) * (f(a) + 4*f(mid) + f(b))
// Error: -h^5/90 * f^(4)(xi). Burden 10th ed, Sec 4.3.
fn simpson(f_a: Float, f_mid: Float, f_b: Float, h: Float) -> Float {
    return (h / 3.0) * (f_a + 4.0 * f_mid + f_b)
}

// Simpson's 3/8 rule: integral ~ (3h/8) * (f0 + 3*f1 + 3*f2 + f3)
// Uses cubic interpolation. Burden 10th ed, Sec 4.3.
fn simpson_38(f0: Float, f1: Float, f2: Float, f3: Float, h: Float) -> Float {
    return (3.0 * h / 8.0) * (f0 + 3.0 * f1 + 3.0 * f2 + f3)
}

// Midpoint rule: integral ~ f(midpoint) * h
fn midpoint_rule(f_mid: Float, h: Float) -> Float {
    return f_mid * h
}

// Gauss-Legendre 2-point quadrature on [a,b]
// Nodes: x = (a+b)/2 +/- (b-a)/(2*sqrt(3))
// Weights: w1 = w2 = 1. Exact for polynomials up to degree 3.
// f_neg = f at left node, f_pos = f at right node
// Burden 10th ed, Sec 4.7.
fn gauss_legendre_2pt(f_neg: Float, f_pos: Float, a: Float, b: Float) -> Float {
    return (b - a) / 2.0 * (f_neg + f_pos)
}

// Gauss-Legendre node positions on [a,b]
fn gauss_legendre_2pt_node_left(a: Float, b: Float) -> Float {
    return (a + b) / 2.0 - (b - a) / (2.0 * sqrt(3.0))
}

fn gauss_legendre_2pt_node_right(a: Float, b: Float) -> Float {
    return (a + b) / 2.0 + (b - a) / (2.0 * sqrt(3.0))
}

// =========================================================
// VECTOR CALCULUS
// =========================================================

// Divergence in 2D: div F = dFx/dx + dFy/dy
// Stewart 9th ed, Sec 16.5.
fn divergence_2d(df_dx: Float, df_dy: Float) -> Float {
    return df_dx + df_dy
}

// Divergence in 3D: div F = dFx/dx + dFy/dy + dFz/dz
fn divergence_3d(df_dx: Float, df_dy: Float, df_dz: Float) -> Float {
    return df_dx + df_dy + df_dz
}

// Scalar curl in 2D: curl F = dFy/dx - dFx/dy
// Stewart 9th ed, Sec 16.5.
fn curl_2d(dfy_dx: Float, dfx_dy: Float) -> Float {
    return dfy_dx - dfx_dy
}

// Laplacian in 2D: nabla^2 f = d2f/dx2 + d2f/dy2
// Stewart 9th ed, Sec 16.5.
fn laplacian_2d(d2f_dx2: Float, d2f_dy2: Float) -> Float {
    return d2f_dx2 + d2f_dy2
}

// Laplacian in 3D: nabla^2 f = d2f/dx2 + d2f/dy2 + d2f/dz2
fn laplacian_3d(d2f_dx2: Float, d2f_dy2: Float, d2f_dz2: Float) -> Float {
    return d2f_dx2 + d2f_dy2 + d2f_dz2
}

// Directional derivative: D_u f = grad(f) dot u_hat
fn directional_derivative_2d(grad_x: Float, grad_y: Float, dir_x: Float, dir_y: Float) -> Float {
    let mag = sqrt(dir_x * dir_x + dir_y * dir_y)
    return grad_x * (dir_x / mag) + grad_y * (dir_y / mag)
}

// Gradient magnitude in 2D: |grad f| = sqrt((df/dx)^2 + (df/dy)^2)
fn gradient_magnitude_2d(grad_x: Float, grad_y: Float) -> Float {
    return sqrt(grad_x * grad_x + grad_y * grad_y)
}

// =========================================================
// ANALYTICAL DERIVATIVES (exact, common functions)
// Stewart 9th ed, differentiation rules.
// =========================================================

// d/dx[x^n] = n * x^(n-1)
fn d_power(x: Float, n: Float) -> Float {
    return n * pow(x, n - 1.0)
}

// d/dx[e^x] = e^x
fn d_exp(x: Float) -> Float {
    return exp(x)
}

// d/dx[a^x] = a^x * ln(a)
fn d_exp_base(x: Float, a: Float) -> Float {
    return pow(a, x) * ln(a)
}

// d/dx[ln(x)] = 1/x
fn d_ln(x: Float) -> Float {
    return 1.0 / x
}

// d/dx[log_a(x)] = 1/(x * ln(a))
fn d_log_base(x: Float, a: Float) -> Float {
    return 1.0 / (x * ln(a))
}

// d/dx[sin(x)] = cos(x)
fn d_sin(x: Float) -> Float {
    return cos(x)
}

// d/dx[cos(x)] = -sin(x)
fn d_cos(x: Float) -> Float {
    return 0.0 - sin(x)
}

// d/dx[tan(x)] = sec^2(x) = 1/cos^2(x)
fn d_tan(x: Float) -> Float {
    let c = cos(x)
    return 1.0 / (c * c)
}

// d/dx[arcsin(x)] = 1/sqrt(1-x^2)
fn d_asin(x: Float) -> Float {
    return 1.0 / sqrt(1.0 - x * x)
}

// d/dx[arccos(x)] = -1/sqrt(1-x^2)
fn d_acos(x: Float) -> Float {
    return 0.0 - 1.0 / sqrt(1.0 - x * x)
}

// d/dx[arctan(x)] = 1/(1+x^2)
fn d_atan(x: Float) -> Float {
    return 1.0 / (1.0 + x * x)
}

// d/dx[sinh(x)] = cosh(x)
fn d_sinh(x: Float) -> Float {
    return cosh(x)
}

// d/dx[cosh(x)] = sinh(x)
fn d_cosh(x: Float) -> Float {
    return sinh(x)
}

// =========================================================
// ANALYTICAL INTEGRALS (antiderivatives at a point)
// Stewart 9th ed, integration tables.
// =========================================================

// Integral of x^n dx = x^(n+1)/(n+1), n != -1
fn int_power(x: Float, n: Float) -> Float {
    return pow(x, n + 1.0) / (n + 1.0)
}

// Integral of e^x dx = e^x
fn int_exp(x: Float) -> Float {
    return exp(x)
}

// Integral of sin(x) dx = -cos(x)
fn int_sin(x: Float) -> Float {
    return 0.0 - cos(x)
}

// Integral of cos(x) dx = sin(x)
fn int_cos(x: Float) -> Float {
    return sin(x)
}

// Integral of 1/x dx = ln|x|
fn int_recip(x: Float) -> Float {
    return ln(abs(x))
}

// Integral of sec^2(x) dx = tan(x)
fn int_sec2(x: Float) -> Float {
    return tan(x)
}

// Definite integral helper: F(b) - F(a)
fn definite_integral(antideriv_at_a: Float, antideriv_at_b: Float) -> Float {
    return antideriv_at_b - antideriv_at_a
}

// =========================================================
// TAYLOR SERIES
// =========================================================

// Factorial helper for series (returns Float)
fn series_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
}

// Taylor series for e^x: sum_{k=0}^{n-1} x^k / k!
// Stewart 9th ed, Sec 11.10.
fn taylor_exp(x: Float, n_terms: Int) -> Float {
    let result = 0.0
    let k = 0
    while k < n_terms {
        result = result + pow(x, k) / series_factorial(k)
        k = k + 1
    }
    return result
}

// Taylor series for sin(x): sum_{k=0}^{n-1} (-1)^k * x^(2k+1) / (2k+1)!
// Stewart 9th ed, Sec 11.10.
fn taylor_sin(x: Float, n_terms: Int) -> Float {
    let result = 0.0
    let k = 0
    while k < n_terms {
        let sign = 1.0
        if (k / 2) * 2 != k { sign = 0.0 - 1.0 }
        let exponent = 2 * k + 1
        result = result + sign * pow(x, exponent) / series_factorial(exponent)
        k = k + 1
    }
    return result
}

// Taylor series for cos(x): sum_{k=0}^{n-1} (-1)^k * x^(2k) / (2k)!
// Stewart 9th ed, Sec 11.10.
fn taylor_cos(x: Float, n_terms: Int) -> Float {
    let result = 0.0
    let k = 0
    while k < n_terms {
        let sign = 1.0
        if (k / 2) * 2 != k { sign = 0.0 - 1.0 }
        let exponent = 2 * k
        result = result + sign * pow(x, exponent) / series_factorial(exponent)
        k = k + 1
    }
    return result
}

// =========================================================
// SPECIAL CONSTANTS
// =========================================================

// Euler-Mascheroni constant: gamma = 0.5772156649...
fn euler_mascheroni() -> Float {
    return 0.5772156649015329
}

// Golden ratio: phi = (1 + sqrt(5)) / 2 = 1.6180339887...
fn golden_ratio() -> Float {
    return 1.618033988749895
}

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

fn calculus_function_count() -> Int {
    return 48
}
