028_trapezoidal_rule.sio

//@ run-pass
// HumanEval 028: Numerical Integration via Trapezoidal Rule
//
// Approximate the integral of f(x) over [a,b] using the trapezoidal rule
// with n intervals. Uses integer arithmetic scaled by 10000 to avoid f64.
// f is passed as a function reference. a, b are scaled integers (x * 1000).
// Returns integral * 1000000 (product of two 1000-scales).

fn identity(x: i64) -> i64 with Mut, Panic, Div { x }

fn square(x: i64) -> i64 with Mut, Panic, Div { x * x / 1000 }

fn cube(x: i64) -> i64 with Mut, Panic, Div { x * x / 1000 * x / 1000 }

fn trapezoid(f: fn(i64) -> i64, a: i64, b: i64, n: i64) -> i64 with Mut, Panic, Div {
    // a, b are x * 1000; f takes x*1000 and returns y*1000
    // h = (b - a) / n
    let h = (b - a) / n

    // sum = f(a)/2 + f(b)/2
    var sum = f(a) + f(b)

    // Add interior points
    var i: i64 = 1
    while i < n {
        let xi = a + i * h
        sum = sum + 2 * f(xi)
        i = i + 1
    }

    // integral = h * sum / 2 (already in 1000*1000 scale)
    h * sum / 2
}

fn main() -> i64 with IO, Mut, Panic, Div {
    // Test 1: integral of identity f(x)=x over [0, 2] = 2.0
    // a=0, b=2000, n=100
    // Expected: 2.0 -> 2000000 in scaled units
    let r1 = trapezoid(identity, 0, 2000, 100)
    // Should be close to 2000000
    let d1 = r1 - 2000000
    let ad1 = if d1 < 0 { 0 - d1 } else { d1 }
    assert(ad1 < 1000)

    // Test 2: integral of f(x)=x over [0, 1] = 0.5
    // Expected: 500000
    let r2 = trapezoid(identity, 0, 1000, 100)
    let d2 = r2 - 500000
    let ad2 = if d2 < 0 { 0 - d2 } else { d2 }
    assert(ad2 < 1000)

    // Test 3: integral of f(x)=x^2 over [0, 1] = 1/3 ~ 333333
    let r3 = trapezoid(square, 0, 1000, 100)
    let d3 = r3 - 333333
    let ad3 = if d3 < 0 { 0 - d3 } else { d3 }
    assert(ad3 < 5000)

    // Test 4: integral of f(x)=x^2 over [0, 2] = 8/3 ~ 2666666
    let r4 = trapezoid(square, 0, 2000, 100)
    let d4 = r4 - 2666666
    let ad4 = if d4 < 0 { 0 - d4 } else { d4 }
    assert(ad4 < 10000)

    // Test 5: integral over zero-width interval [3, 3] = 0
    let r5 = trapezoid(identity, 3000, 3000, 1)
    assert(r5 == 0)

    // Test 6: integral of f(x)=x^3 over [0, 1] = 1/4 = 250000
    let r6 = trapezoid(cube, 0, 1000, 200)
    let d6 = r6 - 250000
    let ad6 = if d6 < 0 { 0 - d6 } else { d6 }
    assert(ad6 < 5000)

    println("028_trapezoidal_rule: ALL TESTS PASSED")
    0
}