#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL x0 AS DOUBLE
    LOCAL x1 AS DOUBLE
    LOCAL x2 AS DOUBLE
    LOCAL fx0 AS DOUBLE
    LOCAL fx1 AS DOUBLE
    LOCAL tolerance AS DOUBLE
    LOCAL iter AS LONG
    LOCAL max_iter AS LONG
    LOCAL temp_str AS STRING

    ' Secant method: x_{n+1} = x_n - f(x_n)*(x_n - x_{n-1})/(f(x_n) - f(x_{n-1}))
    ' No derivative needed, like Newton with finite difference approximation

    PRINT "Secant Method Root Finding"
    PRINT "========================="
    PRINT ""

    ' f(x) = x^2 - 4, roots at +/-2
    PRINT "f(x) = x^2 - 4:"
    x0 = 1.0
    x1 = 4.0
    tolerance = 0.0001
    max_iter = 20
    iter = 0

    fx0 = x0 * x0 - 4.0
    fx1 = x1 * x1 - 4.0

    WHILE iter < max_iter
        IF fx1 - fx0 = 0 THEN
            iter = max_iter
        ELSE
            x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0)
            IF ABS(x2 - x1) < tolerance THEN
                iter = max_iter
            ELSE
                x0 = x1
                fx0 = fx1
                x1 = x2
                fx1 = x1 * x1 - 4.0
            END IF
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(x2)
    PRINT "  Root between 1 and 4: " + temp_str
    temp_str = STR$(x2 * x2 - 4.0)
    PRINT "  f(root) = " + temp_str
    temp_str = STR$(iter)
    PRINT "  Iterations: " + temp_str

    ' f(x) = x^3 - 2x - 5
    PRINT ""
    PRINT "f(x) = x^3 - 2x - 5:"
    x0 = 2.0
    x1 = 3.0
    iter = 0

    fx0 = x0 * x0 * x0 - 2.0 * x0 - 5.0
    fx1 = x1 * x1 * x1 - 2.0 * x1 - 5.0

    WHILE iter < max_iter
        IF fx1 - fx0 = 0 THEN
            iter = max_iter
        ELSE
            x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0)
            IF ABS(x2 - x1) < tolerance THEN
                iter = max_iter
            ELSE
                x0 = x1
                fx0 = fx1
                x1 = x2
                fx1 = x1 * x1 * x1 - 2.0 * x1 - 5.0
            END IF
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(x2)
    PRINT "  Root between 2 and 3: " + temp_str
    temp_str = STR$(x2 * x2 * x2 - 2.0 * x2 - 5.0)
    PRINT "  f(root) = " + temp_str

    ' f(x) = sin(x) - 0.5
    PRINT ""
    PRINT "f(x) = sin(x) - 0.5:"
    x0 = 0.0
    x1 = 1.0
    iter = 0

    fx0 = SIN(x0) - 0.5
    fx1 = SIN(x1) - 0.5

    WHILE iter < max_iter
        IF fx1 - fx0 = 0 THEN
            iter = max_iter
        ELSE
            x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0)
            IF ABS(x2 - x1) < tolerance THEN
                iter = max_iter
            ELSE
                x0 = x1
                fx0 = fx1
                x1 = x2
                fx1 = SIN(x1) - 0.5
            END IF
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(x2)
    PRINT "  Root between 0 and 1: " + temp_str
    temp_str = STR$(SIN(x2))
    PRINT "  sin(root) = " + temp_str
    temp_str = STR$(iter)
    PRINT "  Iterations: " + temp_str

    ' Comparison with Newton: secant needs 2 initial points, no derivative
    PRINT ""
    PRINT "Secant vs Newton comparison:"
    PRINT "  Secant: no derivative needed, ~1.618 order convergence"
    PRINT "  Newton: derivative needed, ~2.0 order convergence"
    PRINT "  Secant: 2 function evaluations per 1 step"
    PRINT "  Newton: 1 function + 1 derivative per step"

    ' Test with different starting points
    PRINT ""
    PRINT "Sensitivity to initial points (x^2-4=0):"
    PRINT "  x0=0, x1=5:"
    x0 = 0.0
    x1 = 5.0
    iter = 0
    fx0 = x0 * x0 - 4.0
    fx1 = x1 * x1 - 4.0

    WHILE iter < 8
        IF fx1 - fx0 = 0 THEN iter = 8
        x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0)
        temp_str = STR$(x2)
        PRINT "    iter " + STR$(iter) + ": x=" + temp_str
        x0 = x1
        fx0 = fx1
        x1 = x2
        fx1 = x1 * x1 - 4.0
        iter = iter + 1
    WEND

    PRINT ""
    PRINT "Secant method demo complete."

    FUNCTION = 0
END FUNCTION
