import math
import numpy as np


class NumericalAnalysis:
    @staticmethod
    def c2_bisection(f, a, b, tol=1e-10, max_iter=100):
        fa = f(a)
        fb = f(b)
        if fa * fb > 0:
            raise ValueError("f(a) and f(b) must have opposite signs.")
        history = []
        c = a
        for i in range(1, max_iter + 1):
            c = (a + b) / 2.0
            fc = f(c)
            history.append(
                {
                    "iter": i,
                    "a": a,
                    "b": b,
                    "c": c,
                    "f(c)": fc,
                    "error": abs(b - a) / 2.0,
                }
            )
            if abs(fc) < tol or abs(b - a) / 2.0 < tol:
                return c, history
            if fa * fc < 0:
                b = c
                fb = fc
            else:
                a = c
                fa = fc
        return c, history

    @staticmethod
    def c2_false_position(f, a, b, tol=1e-10, max_iter=100):
        fa = f(a)
        fb = f(b)
        if fa * fb > 0:
            raise ValueError("f(a) and f(b) must have opposite signs.")
        history = []
        c = a
        for i in range(1, max_iter + 1):
            denom = fb - fa
            if abs(denom) < 1e-15:
                raise ZeroDivisionError("Division by near-zero in false position.")
            c = (a * fb - b * fa) / denom
            fc = f(c)
            history.append(
                {
                    "iter": i,
                    "a": a,
                    "b": b,
                    "c": c,
                    "f(c)": fc,
                    "error": abs(fc),
                }
            )
            if abs(fc) < tol:
                return c, history
            if fa * fc < 0:
                b = c
                fb = fc
            else:
                a = c
                fa = fc
        return c, history

    @staticmethod
    def c2_fixed_point(g, x0, tol=1e-10, max_iter=100):
        history = []
        x_old = x0
        x_new = x0
        for i in range(1, max_iter + 1):
            x_new = g(x_old)
            err = abs(x_new - x_old)
            history.append(
                {
                    "iter": i,
                    "x_old": x_old,
                    "x_new": x_new,
                    "error": err,
                }
            )
            if err < tol:
                return x_new, history
            x_old = x_new
        return x_new, history

    @staticmethod
    def c2_newton_raphson(f, df, x0, tol=1e-10, max_iter=100):
        history = []
        x = x0
        for i in range(1, max_iter + 1):
            fx = f(x)
            dfx = df(x)
            if abs(dfx) < 1e-15:
                raise ZeroDivisionError("Derivative is too small.")
            x_new = x - fx / dfx
            err = abs(x_new - x)
            history.append(
                {
                    "iter": i,
                    "x": x,
                    "f(x)": fx,
                    "df(x)": dfx,
                    "x_new": x_new,
                    "error": err,
                }
            )
            if err < tol or abs(fx) < tol:
                return x_new, history
            x = x_new
        return x, history

    @staticmethod
    def c2_secant(f, x0, x1, tol=1e-10, max_iter=100):
        history = []
        for i in range(1, max_iter + 1):
            f0 = f(x0)
            f1 = f(x1)
            denom = f1 - f0
            if abs(denom) < 1e-15:
                raise ZeroDivisionError("Division by near-zero in secant method.")
            x2 = x1 - f1 * (x1 - x0) / denom
            err = abs(x2 - x1)
            history.append(
                {
                    "iter": i,
                    "x0": x0,
                    "x1": x1,
                    "x2": x2,
                    "f(x2)": f(x2),
                    "error": err,
                }
            )
            if err < tol or abs(f(x2)) < tol:
                return x2, history
            x0, x1 = x1, x2
        return x2, history

    @staticmethod
    def c3_lagrange(x_points, y_points, x):
        n = len(x_points)
        total = 0.0
        for i in range(n):
            term = y_points[i]
            xi = x_points[i]
            for j in range(n):
                if i != j:
                    term *= (x - x_points[j]) / (xi - x_points[j])
            total += term
        return total

    @staticmethod
    def c3_divided_differences(x_points, y_points):
        n = len(x_points)
        table = np.zeros((n, n), dtype=float)
        table[:, 0] = np.array(y_points, dtype=float)
        for j in range(1, n):
            for i in range(n - j):
                table[i, j] = (table[i + 1, j - 1] - table[i, j - 1]) / (x_points[i + j] - x_points[i])
        return table

    @staticmethod
    def c3_newton_interpolation(x_points, y_points, x):
        table = NumericalAnalysis.c3_divided_differences(x_points, y_points)
        n = len(x_points)
        result = table[0, 0]
        product = 1.0
        for i in range(1, n):
            product *= x - x_points[i - 1]
            result += table[0, i] * product
        return result, table

    @staticmethod
    def c3_forward_difference_table(y_points):
        n = len(y_points)
        table = np.zeros((n, n), dtype=float)
        table[:, 0] = np.array(y_points, dtype=float)
        for j in range(1, n):
            for i in range(n - j):
                table[i, j] = table[i + 1, j - 1] - table[i, j - 1]
        return table

    @staticmethod
    def c3_backward_difference_table(y_points):
        n = len(y_points)
        table = np.zeros((n, n), dtype=float)
        table[:, 0] = np.array(y_points, dtype=float)
        for j in range(1, n):
            for i in range(j, n):
                table[i, j] = table[i, j - 1] - table[i - 1, j - 1]
        return table

    @staticmethod
    def c3_newton_forward(x_points, y_points, x):
        n = len(x_points)
        if n < 2:
            raise ValueError("At least two points are required.")
        h = x_points[1] - x_points[0]
        for i in range(1, n - 1):
            if abs((x_points[i + 1] - x_points[i]) - h) > 1e-12:
                raise ValueError("x points must be equally spaced.")
        table = NumericalAnalysis.c3_forward_difference_table(y_points)
        p = (x - x_points[0]) / h
        result = y_points[0]
        term = 1.0
        for i in range(1, n):
            term *= (p - (i - 1))
            result += term * table[0, i] / math.factorial(i)
        return result, table

    @staticmethod
    def c3_newton_backward(x_points, y_points, x):
        n = len(x_points)
        if n < 2:
            raise ValueError("At least two points are required.")
        h = x_points[1] - x_points[0]
        for i in range(1, n - 1):
            if abs((x_points[i + 1] - x_points[i]) - h) > 1e-12:
                raise ValueError("x points must be equally spaced.")
        table = NumericalAnalysis.c3_backward_difference_table(y_points)
        p = (x - x_points[-1]) / h
        result = y_points[-1]
        term = 1.0
        for i in range(1, n):
            term *= (p + (i - 1))
            result += term * table[-1, i] / math.factorial(i)
        return result, table

    @staticmethod
    def c3_inverse_interpolation(x_points, y_points, y_target):
        return NumericalAnalysis.c3_lagrange(y_points, x_points, y_target)

    @staticmethod
    def c3_chebyshev_polynomial(n, x):
        if n == 0:
            return 1.0
        if n == 1:
            return float(x)
        t0 = 1.0
        t1 = float(x)
        for _ in range(2, n + 1):
            t0, t1 = t1, 2.0 * x * t1 - t0
        return t1

    @staticmethod
    def c3_chebyshev_nodes(a, b, n):
        nodes = []
        for k in range(1, n + 1):
            xk = 0.5 * (a + b) + 0.5 * (b - a) * math.cos((2 * k - 1) * math.pi / (2 * n))
            nodes.append(xk)
        return sorted(nodes)

    @staticmethod
    def c3_piecewise_linear(x_points, y_points, x):
        if x < x_points[0] or x > x_points[-1]:
            raise ValueError("x is outside the data range.")
            n = len(x_points)
        for i in range(n - 1):
            if x_points[i] <= x <= x_points[i + 1]:
                x0, x1 = x_points[i], x_points[i + 1]
                y0, y1 = y_points[i], y_points[i + 1]
                return y0 + (y1 - y0) * (x - x0) / (x1 - x0)
        return y_points[-1]

    @staticmethod
    def c3_piecewise_quadratic(x_points, y_points, x):
        if len(x_points) < 3:
            raise ValueError("At least three points are required.")
        if x < x_points[0] or x > x_points[-1]:
            raise ValueError("x is outside the data range.")
        n = len(x_points)
        for i in range(n - 1):
            if x_points[i] <= x <= x_points[i + 1]:
                if i == 0:
                    idx = [0, 1, 2]
                elif i == n - 2:
                    idx = [n - 3, n - 2, n - 1]
                else:
                    idx = [i - 1, i, i + 1]
                xp = [x_points[j] for j in idx]
                yp = [y_points[j] for j in idx]
                return NumericalAnalysis.c3_lagrange(xp, yp, x)
        return y_points[-1]

    @staticmethod
    def c3_cubic_spline_coefficients(x, y):
        x = np.array(x, dtype=float)
        y = np.array(y, dtype=float)
        n = len(x)
        if n < 3:
            raise ValueError("At least three points are required.")
        h = np.diff(x)
        if np.any(h <= 0):
            raise ValueError("x must be strictly increasing.")

        alpha = np.zeros(n, dtype=float)
        for i in range(1, n - 1):
            alpha[i] = (3.0 / h[i]) * (y[i + 1] - y[i]) - (3.0 / h[i - 1]) * (y[i] - y[i - 1])

        l = np.ones(n, dtype=float)
        mu = np.zeros(n, dtype=float)
        z = np.zeros(n, dtype=float)

        for i in range(1, n - 1):
            l[i] = 2.0 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1]
            if abs(l[i]) < 1e-15:
                raise ZeroDivisionError("Spline system breakdown.")
            mu[i] = h[i] / l[i]
            z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]

        a = y.copy()
        b = np.zeros(n - 1, dtype=float)
        c = np.zeros(n, dtype=float)
        d = np.zeros(n - 1, dtype=float)

        for j in range(n - 2, -1, -1):
            c[j] = z[j] - mu[j] * c[j + 1]
            b[j] = (a[j + 1] - a[j]) / h[j] - h[j] * (c[j + 1] + 2.0 * c[j]) / 3.0
            d[j] = (c[j + 1] - c[j]) / (3.0 * h[j])

        return a[:-1], b, c[:-1], d

    @staticmethod
    def c3_cubic_spline(x_points, y_points, x):
        if x < x_points[0] or x > x_points[-1]:
            raise ValueError("x is outside the data range.")
        a, b, c, d = NumericalAnalysis.c3_cubic_spline_coefficients(x_points, y_points)
        for i in range(len(x_points) - 1):
            if x_points[i] <= x <= x_points[i + 1]:
                dx = x - x_points[i]
                return a[i] + b[i] * dx + c[i] * dx**2 + d[i] * dx**3
        return y_points[-1]

    @staticmethod
    def c4_trapezoidal(f, a, b, n):
        h = (b - a) / n
        s = 0.5 * (f(a) + f(b))
        for i in range(1, n):
            s += f(a + i * h)
        return h * s

    @staticmethod
    def c4_romberg(f, a, b, max_level=6, tol=1e-10):
        R = np.zeros((max_level, max_level), dtype=float)
        for i in range(max_level):
            n = 2**i
            R[i, 0] = NumericalAnalysis.c4_trapezoidal(f, a, b, n)
            for j in range(1, i + 1):
                R[i, j] = R[i, j - 1] + (R[i, j - 1] - R[i - 1, j - 1]) / (4**j - 1)
            if i > 0 and abs(R[i, i] - R[i - 1, i - 1]) < tol:
                return R[i, i], R[: i + 1, : i + 1]
        return R[max_level - 1, max_level - 1], R


def print_table(rows, max_rows=10):
    if not rows:
        print("No data")
        return
    keys = list(rows[0].keys())
    widths = []
    for k in keys:
        width = len(str(k))
        for row in rows[:max_rows]:
            width = max(width, len(f"{row[k]:.10f}" if isinstance(row[k], float) else str(row[k])))
        widths.append(width + 2)

    header = "".join(str(k).ljust(w) for k, w in zip(keys, widths))
    print(header)
    print("-" * len(header))

    for row in rows[:max_rows]:
        line = []
        for k, w in zip(keys, widths):
            v = row[k]
            if isinstance(v, float):
                s = f"{v:.10f}"
            else:
                s = str(v)
            line.append(s.ljust(w))
        print("".join(line))

    if len(rows) > max_rows:
        print("...")


def read_float_list(prompt):
    return [float(x) for x in input(prompt).strip().split()]


def read_function():
    expr = input("Enter f(x) as a Python expression using x: ").strip()
    return lambda x: eval(expr, {"__builtins__": {}}, {"x": x, "math": math, "np": np})


def read_function_one():
    expr = input("Enter g(x) as a Python expression using x: ").strip()
    return lambda x: eval(expr, {"__builtins__": {}}, {"x": x, "math": math, "np": np})


def solve_chapter_2():
    print("\nChapter 2")
    print("1. Bisection")
    print("2. False Position")
    print("3. Fixed Point")
    print("4. Newton-Raphson")
    print("5. Secant")
    choice = input("Select: ").strip()

    if choice == "1":
        f = read_function()
        a = float(input("a = "))
        b = float(input("b = "))
        root, hist = NumericalAnalysis.c2_bisection(f, a, b)
        print("Root:", root)
        print_table(hist)
    elif choice == "2":
        f = read_function()
        a = float(input("a = "))
        b = float(input("b = "))
        root, hist = NumericalAnalysis.c2_false_position(f, a, b)
        print("Root:", root)
        print_table(hist)
    elif choice == "3":
        g = read_function_one()
        x0 = float(input("x0 = "))
        root, hist = NumericalAnalysis.c2_fixed_point(g, x0)
        print("Root:", root)
        print_table(hist)
    elif choice == "4":
        f = read_function()
        df = read_function_one()
        x0 = float(input("x0 = "))
        root, hist = NumericalAnalysis.c2_newton_raphson(f, df, x0)
        print("Root:", root)
        print_table(hist)
    elif choice == "5":
        f = read_function()
        x0 = float(input("x0 = "))
        x1 = float(input("x1 = "))
        root, hist = NumericalAnalysis.c2_secant(f, x0, x1)
        print("Root:", root)
        print_table(hist)


def solve_chapter_3():
    print("\nChapter 3")
    print("1. Lagrange")
    print("2. Newton Divided Differences")
    print("3. Newton Forward")
    print("4. Newton Backward")
    print("5. Inverse Interpolation")
    print("6. Chebyshev Polynomial")
    print("7. Chebyshev Nodes")
    print("8. Piecewise Linear")
    print("9. Piecewise Quadratic")
    print("10. Cubic Spline")
    choice = input("Select: ").strip()

    if choice in {"1", "2", "3", "4", "5", "8", "9", "10"}:
        x_points = read_float_list("x points separated by spaces: ")
        y_points = read_float_list("y points separated by spaces: ")
        if len(x_points) != len(y_points):
            raise ValueError("x and y lengths must match.")

    if choice == "1":
        x = float(input("x = "))
        val = NumericalAnalysis.c3_lagrange(x_points, y_points, x)
        print("Value:", val)
    elif choice == "2":
        x = float(input("x = "))
        val, table = NumericalAnalysis.c3_newton_interpolation(x_points, y_points, x)
        print("Value:", val)
        print(table)
    elif choice == "3":
        x = float(input("x = "))
        val, table = NumericalAnalysis.c3_newton_forward(x_points, y_points, x)
        print("Value:", val)
        print(table)
    elif choice == "4":
        x = float(input("x = "))
        val, table = NumericalAnalysis.c3_newton_backward(x_points, y_points, x)
        print("Value:", val)
        print(table)
    elif choice == "5":
        y_target = float(input("y target = "))
        val = NumericalAnalysis.c3_inverse_interpolation(x_points, y_points, y_target)
        print("x value:", val)
    elif choice == "6":
        n = int(input("n = "))
        x = float(input("x = "))
        print("Value:", NumericalAnalysis.c3_chebyshev_polynomial(n, x))
    elif choice == "7":
        a = float(input("a = "))
        b = float(input("b = "))
        n = int(input("n = "))
        print(NumericalAnalysis.c3_chebyshev_nodes(a, b, n))
    elif choice == "8":
        x = float(input("x = "))
        print("Value:", NumericalAnalysis.c3_piecewise_linear(x_points, y_points, x))
    elif choice == "9":
        x = float(input("x = "))
        print("Value:", NumericalAnalysis.c3_piecewise_quadratic(x_points, y_points, x))
    elif choice == "10":
        x = float(input("x = "))
        print("Value:", NumericalAnalysis.c3_cubic_spline(x_points, y_points, x))


def solve_chapter_4():
    print("\nChapter 4")
    print("1. Romberg Integration")
    choice = input("Select: ").strip()

    if choice == "1":
        f = read_function()
        a = float(input("a = "))
        b = float(input("b = "))
        max_level = int(input("max level = "))
        val, table = NumericalAnalysis.c4_romberg(f, a, b, max_level=max_level)
        print("Integral:", val)
        print(table)


def demo():
    f = lambda x: x**3 - x - 2
    df = lambda x: 3 * x**2 - 1
    g = lambda x: (x + 2) ** (1 / 3)

    print("Chapter 2 Demo")
    print(NumericalAnalysis.c2_bisection(f, 1, 2)[0])
    print(NumericalAnalysis.c2_false_position(f, 1, 2)[0])
    print(NumericalAnalysis.c2_fixed_point(g, 1.5)[0])
    print(NumericalAnalysis.c2_newton_raphson(f, df, 1.5)[0])
    print(NumericalAnalysis.c2_secant(f, 1, 2)[0])

    x_data = [1, 2, 3, 4]
    y_data = [1, 4, 9, 16]

    print("Chapter 3 Demo")
    print(NumericalAnalysis.c3_lagrange(x_data, y_data, 2.5))
    print(NumericalAnalysis.c3_newton_interpolation(x_data, y_data, 2.5)[0])
    print(NumericalAnalysis.c3_inverse_interpolation(x_data, y_data, 6.25))
    print(NumericalAnalysis.c3_chebyshev_polynomial(4, 0.5))
    print(NumericalAnalysis.c3_chebyshev_nodes(0, 1, 5))
    print(NumericalAnalysis.c3_piecewise_linear(x_data, y_data, 2.5))
    print(NumericalAnalysis.c3_piecewise_quadratic(x_data, y_data, 2.5))
    print(NumericalAnalysis.c3_cubic_spline(x_data, y_data, 2.5))

    print("Chapter 4 Demo")
    val, table = NumericalAnalysis.c4_romberg(math.sin, 0, math.pi)
    print(val)
    print(table)


if name == "__main__":
    print("Numerical Analysis Project")
    print("1. Chapter 2")
    print("2. Chapter 3")
    print("3. Chapter 4")
    print("4. Demo")
    print("0. Exit")
    while True:
        try:
            choice = input("Select: ").strip()
            if choice == "1":
                solve_chapter_2()
            elif choice == "2":
                solve_chapter_3()
            elif choice == "3":
                solve_chapter_4()
            elif choice == "4":
                demo()
            elif choice == "0":
                break
        except Exception as e:
            print("Error:", e)