#COMPILE EXE
#DIM ALL

FUNCTION PBMAIN() AS LONG
    LOCAL x0 AS DOUBLE
    LOCAL fx AS DOUBLE
    LOCAL fpx AS DOUBLE
    LOCAL x1 AS DOUBLE
    LOCAL tolerance AS DOUBLE
    LOCAL iter AS LONG
    LOCAL max_iter AS LONG
    LOCAL temp_str AS STRING

    ' Newton-Raphson root finding: x_{n+1} = x_n - f(x_n)/f'(x_n)

    PRINT "Newton-Raphson Root Finding"
    PRINT "==========================="
    PRINT ""

    ' Find root of f(x) = x^2 - 4 (roots at x=2, x=-2)
    PRINT "f(x) = x^2 - 4 (roots at +/-2):"
    x0 = 3.0
    tolerance = 0.0001
    max_iter = 20
    iter = 0

    WHILE iter < max_iter
        fx = x0 * x0 - 4.0
        fpx = 2.0 * x0
        IF fpx = 0 THEN
            iter = max_iter
        ELSE
            x1 = x0 - fx / fpx
            IF ABS(x1 - x0) < tolerance THEN
                iter = max_iter
            ELSE
                x0 = x1
            END IF
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(x1)
    PRINT "  Root near 3: " + temp_str
    temp_str = STR$(x1 * x1 - 4.0)
    PRINT "  f(root) = " + temp_str

    ' Find root of f(x) = x^3 - 2x - 5 (Babbage's equation)
    PRINT ""
    PRINT "f(x) = x^3 - 2x - 5:"
    x0 = 2.0
    iter = 0

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

    temp_str = STR$(x1)
    PRINT "  Root: " + temp_str
    temp_str = STR$(x1 * x1 * x1 - 2.0 * x1 - 5.0)
    PRINT "  f(root) = " + temp_str

    ' Find root of f(x) = sin(x) - 0.5 (root at pi/6 = 0.5236)
    PRINT ""
    PRINT "f(x) = sin(x) - 0.5 (root at 0.5236):"
    x0 = 0.5
    iter = 0

    WHILE iter < max_iter
        fx = SIN(x0) - 0.5
        fpx = COS(x0)
        IF fpx = 0 THEN
            iter = max_iter
        ELSE
            x1 = x0 - fx / fpx
            IF ABS(x1 - x0) < tolerance THEN
                iter = max_iter
            ELSE
                x0 = x1
            END IF
        END IF
        iter = iter + 1
    WEND

    temp_str = STR$(x1)
    PRINT "  Root: " + temp_str
    temp_str = STR$(SIN(x1))
    PRINT "  sin(root) = " + temp_str

    ' Convergence behavior
    PRINT ""
    PRINT "Convergence table for x^2 - 4 = 0 starting at x0=4:"
    x0 = 4.0
    iter = 0

    WHILE iter < 5
        temp_str = STR$(x0)
        PRINT "  iter " + STR$(iter) + ": x = " + temp_str
        fx = x0 * x0 - 4.0
        fpx = 2.0 * x0
        IF fpx <> 0 THEN
            x1 = x0 - fx / fpx
            x0 = x1
        END IF
        iter = iter + 1
    WEND

    ' Quadratic convergence demonstration
    PRINT ""
    PRINT "Quadratic convergence: error roughly squares each step"
    PRINT "Starting from x=8, iter 5 should reach near 2"

    x0 = 8.0
    iter = 0
    WHILE iter < 6
        temp_str = STR$(x0)
        PRINT "  iter " + STR$(iter) + ": x = " + temp_str
        fx = x0 * x0 - 4.0
        fpx = 2.0 * x0
        IF fpx <> 0 THEN
            x1 = x0 - fx / fpx
            x0 = x1
        END IF
        iter = iter + 1
    WEND

    PRINT ""
    PRINT "Newton-Raphson demo complete."

    FUNCTION = 0
END FUNCTION
