#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL a AS DOUBLE
    LOCAL b AS DOUBLE
    LOCAL n AS LONG
    LOCAL h AS DOUBLE
    LOCAL x AS DOUBLE
    LOCAL sum_val AS DOUBLE
    LOCAL total AS DOUBLE
    LOCAL prev_total AS DOUBLE
    LOCAL tolerance AS DOUBLE
    LOCAL iter AS LONG
    LOCAL max_iter AS LONG
    LOCAL temp_str AS STRING
    LOCAL i AS LONG

    ' Simpson's rule: integral of f(x) from a to b
    ' Approx = h/3 * (f(a) + f(b) + 4*sum_odd + 2*sum_even)
    ' Adaptive: double segments until convergence

    PRINT "Simpson's Rule Adaptive Integration"
    PRINT "==================================="
    PRINT ""

    ' Integral of x^2 from 0 to 2 = 8/3 = 2.6667
    PRINT "Integral of x^2 from 0 to 2 (exact = 2.6667):"
    a = 0.0
    b = 2.0
    tolerance = 0.0001
    max_iter = 20
    iter = 0
    n = 2
    prev_total = 0

    WHILE iter < max_iter
        h = (b - a) / n
        sum_val = (a * a) + (b * b)
        x = a + h
        i = 1
        WHILE i < n
            IF i MOD 2 = 1 THEN
                sum_val = sum_val + 4.0 * (x * x)
            ELSE
                sum_val = sum_val + 2.0 * (x * x)
            END IF
            x = x + h
            i = i + 1
        WEND
        total = sum_val * h / 3.0

        IF iter > 0 AND ABS(total - prev_total) < tolerance THEN
            iter = max_iter
        ELSE
            prev_total = total
            n = n * 2
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(total)
    PRINT "  Result: " + temp_str
    temp_str = STR$(n / 2)
    PRINT "  Segments: " + temp_str

    ' Integral of sin(x) from 0 to pi = 2
    PRINT ""
    PRINT "Integral of sin(x) from 0 to pi (exact = 2.0):"
    a = 0.0
    b = 3.14159265
    n = 4
    prev_total = 0
    iter = 0

    WHILE iter < max_iter
        h = (b - a) / n
        sum_val = SIN(a) + SIN(b)
        x = a + h
        i = 1
        WHILE i < n
            IF i MOD 2 = 1 THEN
                sum_val = sum_val + 4.0 * SIN(x)
            ELSE
                sum_val = sum_val + 2.0 * SIN(x)
            END IF
            x = x + h
            i = i + 1
        WEND
        total = sum_val * h / 3.0

        IF iter > 0 AND ABS(total - prev_total) < tolerance THEN
            iter = max_iter
        ELSE
            prev_total = total
            n = n * 2
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(total)
    PRINT "  Result: " + temp_str

    ' Integral of 1/x from 1 to e = 1
    PRINT ""
    PRINT "Integral of 1/x from 1 to e (exact = 1.0):"
    a = 1.0
    b = 2.71828183
    n = 4
    prev_total = 0
    iter = 0

    WHILE iter < max_iter
        h = (b - a) / n
        sum_val = (1.0 / a) + (1.0 / b)
        x = a + h
        i = 1
        WHILE i < n
            IF i MOD 2 = 1 THEN
                sum_val = sum_val + 4.0 * (1.0 / x)
            ELSE
                sum_val = sum_val + 2.0 * (1.0 / x)
            END IF
            x = x + h
            i = i + 1
        WEND
        total = sum_val * h / 3.0

        IF iter > 0 AND ABS(total - prev_total) < tolerance THEN
            iter = max_iter
        ELSE
            prev_total = total
            n = n * 2
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(total)
    PRINT "  Result: " + temp_str

    ' Integral of cos(x) from 0 to pi/2 = 1
    PRINT ""
    PRINT "Integral of cos(x) from 0 to pi/2 (exact = 1.0):"
    a = 0.0
    b = 1.57079633
    n = 4
    prev_total = 0
    iter = 0

    WHILE iter < max_iter
        h = (b - a) / n
        sum_val = COS(a) + COS(b)
        x = a + h
        i = 1
        WHILE i < n
            IF i MOD 2 = 1 THEN
                sum_val = sum_val + 4.0 * COS(x)
            ELSE
                sum_val = sum_val + 2.0 * COS(x)
            END IF
            x = x + h
            i = i + 1
        WEND
        total = sum_val * h / 3.0

        IF iter > 0 AND ABS(total - prev_total) < tolerance THEN
            iter = max_iter
        ELSE
            prev_total = total
            n = n * 2
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(total)
    PRINT "  Result: " + temp_str

    ' Convergence table
    PRINT ""
    PRINT "Convergence for x^2 from 0 to 2:"
    a = 0.0
    b = 2.0
    n = 2
    iter = 0

    WHILE iter < 8
        h = (b - a) / n
        sum_val = (a * a) + (b * b)
        x = a + h
        i = 1
        WHILE i < n
            IF i MOD 2 = 1 THEN
                sum_val = sum_val + 4.0 * (x * x)
            ELSE
                sum_val = sum_val + 2.0 * (x * x)
            END IF
            x = x + h
            i = i + 1
        WEND
        total = sum_val * h / 3.0
        temp_str = STR$(total)
        PRINT "  n=" + STR$(n) + ": " + temp_str
        n = n * 2
        iter = iter + 1
    WEND

    PRINT ""
    PRINT "Simpson integration demo complete."

    FUNCTION = 0
END FUNCTION
