// ============================================================================
// Numerical integration — trapezoidal and Simpson's rule
// ============================================================================
//
// Approximate the area under f over [a,b] with n panels. The trapezoid rule
// averages the endpoints of each panel; Simpson's rule fits parabolas and is
// exact for cubics. Here f(x) = x², whose integral over [0,1] is 1/3.
//
//   soma run integration.cell run
// ============================================================================

cell Integration {
    face {
        signal f(x: Float) -> Float
        signal trapezoid(a: Float, b: Float, n: Int) -> Float
        signal simpson(a: Float, b: Float, n: Int) -> Float
        signal run() -> Map
    }

    on f(x: Float) [native] {
        return x * x
    }

    on trapezoid(a: Float, b: Float, n: Int) [native] {
        let h = (b - a) / to_float(n)
        let s = (f(a) + f(b)) / 2.0
        let i = 1
        while i < n {
            s = s + f(a + to_float(i) * h)
            i = i + 1
        }
        return s * h
    }

    on simpson(a: Float, b: Float, n: Int) [native] {
        // n must be even
        let h = (b - a) / to_float(n)
        let s = f(a) + f(b)
        let i = 1
        while i < n {
            let x = a + to_float(i) * h
            if i % 2 == 1 { s = s + 4.0 * f(x) } else { s = s + 2.0 * f(x) }
            i = i + 1
        }
        return s * h / 3.0
    }

    on run() {
        return map(
            "trapezoid", trapezoid(0.0, 1.0, 1000),
            "simpson", simpson(0.0, 1.0, 100)
        )
    }
}

cell test IntegrationTests {
    rules {
        // integral of x^2 over [0,1] is 1/3
        assert abs(trapezoid(0.0, 1.0, 1000) - 0.3333333333333) < 0.0001
        // Simpson is exact for polynomials up to cubic
        assert abs(simpson(0.0, 1.0, 100) - 0.3333333333333) < 0.0000001
        // integral over [0,3] of x^2 = 9
        assert abs(simpson(0.0, 3.0, 100) - 9.0) < 0.0000001
    }
}