' Polynomial: evaluation (Horner), root finding, interpolation, arithmetic
#COMPILE EXE
#DIM ALL

FUNCTION PolyEval(BYVAL a0 AS DOUBLE, BYVAL a1 AS DOUBLE, BYVAL a2 AS DOUBLE, BYVAL a3 AS DOUBLE, BYVAL x AS DOUBLE) AS DOUBLE
    FUNCTION = ((a3 * x + a2) * x + a1) * x + a0
END FUNCTION

FUNCTION HornerEval(coeffs() AS DOUBLE, BYVAL n AS LONG, BYVAL x AS DOUBLE) AS DOUBLE
    LOCAL i AS LONG
    LOCAL result AS DOUBLE
    result = coeffs(n)
    FOR i = n - 1 TO 1 STEP -1
        result = result * x + coeffs(i)
    NEXT i
    FUNCTION = result
END FUNCTION

FUNCTION BisectionRoot(BYVAL a0 AS DOUBLE, BYVAL a1 AS DOUBLE, BYVAL a2 AS DOUBLE, BYVAL a3 AS DOUBLE, BYVAL a AS DOUBLE, BYVAL b AS DOUBLE) AS DOUBLE
    LOCAL fa  AS DOUBLE
    LOCAL fb  AS DOUBLE
    LOCAL fc  AS DOUBLE
    LOCAL c  AS DOUBLE
    LOCAL iter AS LONG
    fa = PolyEval(a0, a1, a2, a3, a)
    fb = PolyEval(a0, a1, a2, a3, b)
    IF fa * fb > 0 THEN FUNCTION = 0: EXIT FUNCTION
    DO WHILE (b - a) > 0.0000000001
        c = (a + b) / 2.0
        fc = PolyEval(a0, a1, a2, a3, c)
        IF fc = 0.0 THEN FUNCTION = c: EXIT FUNCTION
        IF PolyEval(a0, a1, a2, a3, a) * fc < 0 THEN b = c ELSE a = c
        iter = iter + 1
        IF iter > 200 THEN EXIT LOOP
    LOOP
    FUNCTION = (a + b) / 2.0
END FUNCTION

FUNCTION LagrangeInterp(BYVAL x1 AS DOUBLE, BYVAL y1 AS DOUBLE, BYVAL x2 AS DOUBLE, BYVAL y2 AS DOUBLE, BYVAL x3 AS DOUBLE, BYVAL y3 AS DOUBLE, BYVAL x AS DOUBLE) AS DOUBLE
    LOCAL term1  AS DOUBLE
    LOCAL term2  AS DOUBLE
    LOCAL term3  AS DOUBLE
    term1 = y1 * (x - x2) * (x - x3) / ((x1 - x2) * (x1 - x3))
    term2 = y2 * (x - x1) * (x - x3) / ((x2 - x1) * (x2 - x3))
    term3 = y3 * (x - x1) * (x - x2) / ((x3 - x1) * (x3 - x2))
    FUNCTION = term1 + term2 + term3
END FUNCTION

FUNCTION PolyQuadratic(BYVAL a AS DOUBLE, BYVAL b AS DOUBLE, BYVAL c AS DOUBLE, root1 AS DOUBLE, root2 AS DOUBLE) AS LONG
    LOCAL disc AS DOUBLE
    disc = b * b - 4.0 * a * c
    IF disc < 0 THEN FUNCTION = 0: EXIT FUNCTION
    root1 = (-b + SQR(disc)) / (2.0 * a)
    root2 = (-b - SQR(disc)) / (2.0 * a)
    FUNCTION = 1
END FUNCTION

FUNCTION PBMAIN() AS LONG
    LOCAL c(1 TO 4) AS DOUBLE
    LOCAL r1  AS DOUBLE
    LOCAL r2  AS DOUBLE
    PRINT "=== Polynomial Library ==="
    PRINT
    PRINT "f(x) = x^3 - 6x^2 + 11x - 6"
    PRINT "f(0) = " + FORMAT$(PolyEval(-6, 11, -6, 1, 0.0), "#0.000")
    PRINT "f(1) = " + FORMAT$(PolyEval(-6, 11, -6, 1, 1.0), "#0.000")
    PRINT "f(2) = " + FORMAT$(PolyEval(-6, 11, -6, 1, 2.0), "#0.000")
    PRINT "f(3) = " + FORMAT$(PolyEval(-6, 11, -6, 1, 3.0), "#0.000")
    PRINT
    PRINT "f(x) = x^3 - 1"
    PRINT "Bisection root in [0,2]: " + FORMAT$(BisectionRoot(-1, 0, 0, 1, 0.0, 2.0), "#0.000000")
    PRINT
    PRINT "Quadratic x^2 - 5x + 6 = 0:"
    IF PolyQuadratic(1.0, -5.0, 6.0, r1, r2) THEN
        PRINT "  Roots: " + FORMAT$(r1, "#0.000") + ", " + FORMAT$(r2, "#0.000")
    END IF
    PRINT
    PRINT "Lagrange interpolation (0,1),(1,2),(2,5):"
    PRINT "  f(0.5) = " + FORMAT$(LagrangeInterp(0, 1, 1, 2, 2, 5, 0.5), "#0.000")
    PRINT "  f(1.5) = " + FORMAT$(LagrangeInterp(0, 1, 1, 2, 2, 5, 1.5), "#0.000")
    PRINT
    c(1) = -6: c(2) = 11: c(3) = -6: c(4) = 1
    PRINT "Horner eval f(10) = " + FORMAT$(HornerEval(c(), 4, 10.0), "#0.000")
    PRINT
    PRINT "Polynomial tests complete."
    FUNCTION = 0
END FUNCTION
