import math
import numpy as np


# =========================================================
# فصل 2 : حل معادلات غیرخطی
# =========================================================

def bisection(f, a, b, tol=1e-10, max_iter=100):
    """
    روش دوبخشی برای یافتن ریشه f(x)=0 در بازه [a, b]
    شرط لازم: f(a)*f(b) < 0
    """
    fa = f(a)
    fb = f(b)

    if fa * fb > 0:
        raise ValueError("شرط تغییر علامت در بازه برقرار نیست.")

    results = []

    for i in range(1, max_iter + 1):
        c = (a + b) / 2
        fc = f(c)

        results.append({
            "iter": i,
            "a": a,
            "b": b,
            "c": c,
            "f(c)": fc,
            "error": abs(b - a) / 2
        })

        if abs(fc) < tol or abs(b - a) / 2 < tol:
            return c, results

        if fa * fc < 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc

    return c, results


def false_position(f, a, b, tol=1e-10, max_iter=100):
    """
    روش نابجایی / رگولا فالسی
    """
    fa = f(a)
    fb = f(b)

    if fa * fb > 0:
        raise ValueError("شرط تغییر علامت در بازه برقرار نیست.")

    results = []

    for i in range(1, max_iter + 1):
        c = (a * fb - b * fa) / (fb - fa)
        fc = f(c)

        results.append({
            "iter": i,
            "a": a,
            "b": b,
            "c": c,
            "f(c)": fc,
            "error": abs(fc)
        })

        if abs(fc) < tol:
            return c, results

        if fa * fc < 0:
            b = c
            fb = fc
        else:
            a = c
            fa = fc

    return c, results


def fixed_point(g, x0, tol=1e-10, max_iter=100):
    """
    روش تکرار نقطه ثابت: x_{n+1} = g(x_n)
    """
    results = []
    x_old = x0

    for i in range(1, max_iter + 1):
        x_new = g(x_old)
        err = abs(x_new - x_old)

        results.append({
            "iter": i,
            "x_old": x_old,
            "x_new": x_new,
            "error": err
        })

        if err < tol:
            return x_new, results

        x_old = x_new

    return x_new, results


def newton_raphson(f, df, x0, tol=1e-10, max_iter=100):
    """
    روش نیوتون-رافسون
    x_{n+1} = x_n - f(x_n)/f'(x_n)
    """
    results = []
    x = x0

    for i in range(1, max_iter + 1):
        fx = f(x)
        dfx = df(x)

        if abs(dfx) < 1e-15:
            raise ZeroDivisionError("مشتق در نقطه بسیار کوچک یا صفر است.")

        x_new = x - fx / dfx
        err = abs(x_new - x)

        results.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, results

        x = x_new

    return x, results


def secant(f, x0, x1, tol=1e-10, max_iter=100):
    """
    روش وتری (Secant)
    """
    results = []

    for i in range(1, max_iter + 1):
        f0 = f(x0)
        f1 = f(x1)

        if abs(f1 - f0) < 1e-15:
            raise ZeroDivisionError("اختلاف مقادیر تابع بسیار کوچک است.")

        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        err = abs(x2 - x1)

        results.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, results

        x0, x1 = x1, x2

    return x2, results


# =========================================================
# فصل 3 : درون‌یابی
# =========================================================

def lagrange_interpolation(x_points, y_points, x):
    """
    محاسبه مقدار درون‌یابی لاگرانژ در x
    """
    n = len(x_points)
    result = 0.0

    for i in range(n):
        term = y_points[i]
        for j in range(n):
            if i != j:
                term *= (x - x_points[j]) / (x_points[i] - x_points[j])
        result += term

    return result


def lagrange_polynomial_string(x_points, y_points):
    """
    نمایش تقریبی چندجمله‌ای لاگرانژ به صورت رشته
    """
    n = len(x_points)
    terms = []

    for i in range(n):
        numerator_parts = []
        denom = 1
        for j in range(n):
            if i != j:
                numerator_parts.append(f"(x - {x_points[j]})")
                denom *= (x_points[i] - x_points[j])

        term_str = f"{y_points[i]/denom:.6f}*" + "*".join(numerator_parts)
        terms.append(term_str)

    return " + ".join(terms)


def divided_differences(x_points, y_points):
    """
    جدول تفاضلات تقسیم‌شده نیوتون
    """
    n = len(x_points)
    table = np.zeros((n, n))
    table[:, 0] = y_points

    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


def newton_interpolation(x_points, y_points, x):
    """
    درون‌یابی نیوتون با تفاضلات تقسیم‌شده
    """
    n = len(x_points)
    table = divided_differences(x_points, y_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


def newton_forward_difference_table(y):
    """
    جدول تفاضلات پیشرو
    """
    n = len(y)
    table = np.zeros((n, n))
    table[:, 0] = y

    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


def factorial(n):
    return math.factorial(n)


def newton_forward(x_points, y_points, x):
    """
    درون‌یابی نیوتون پیشرو
    شرط: نقاط با فاصله مساوی باشند
    """
    n = len(x_points)
    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("فاصله نقاط مساوی نیست.")

    table = newton_forward_difference_table(y_points)
    p = (x - x_points[0]) / h

    result = y_points[0]
    p_term = 1.0

    for i in range(1, n):
        p_term *= (p - (i - 1))
        result += (p_term * table[0, i]) / factorial(i)

    return result, table


def newton_backward_difference_table(y):
    """
    جدول تفاضلات پسرو
    """
    n = len(y)
    table = np.zeros((n, n))
    table[:, 0] = y

    for j in range(1, n):
        for i in range(n - 1, j - 1, -1):
            table[i, j] = table[i, j - 1] - table[i - 1, j - 1]

    return table


def newton_backward(x_points, y_points, x):
    """
    درون‌یابی نیوتون پسرو
    شرط: نقاط با فاصله مساوی باشند
    """
    n = len(x_points)
    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("فاصله نقاط مساوی نیست.")

    table = newton_backward_difference_table(y_points)
    p = (x - x_points[-1]) / h

    result = y_points[-1]
    p_term = 1.0

    for i in range(1, n):
        p_term *= (p + (i - 1))
        result += (p_term * table[-1, i]) / factorial(i)

    return result, table


def inverse_interpolation(x_points, y_points, y_target):
    """
    درون‌یابی معکوس:
    با جابجایی نقش x و y، مقدار x متناظر با y_target را تقریب می‌زند
    """
    return lagrange_interpolation(y_points, x_points, y_target)


# =========================================================
# چندجمله‌ای‌های چبیشف
# =========================================================

def chebyshev_polynomial(n, x):
    """
    محاسبه T_n(x) به روش بازگشتی
    """
    if n == 0:
        return 1
    elif n == 1:
        return x

    T0, T1 = 1, x
    for _ in range(2, n + 1):
        T0, T1 = T1, 2 * x * T1 - T0

    return T1


def chebyshev_nodes(a, b, n):
    """
    تولید n گره چبیشف در بازه [a, b]
    """
    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)


# =========================================================
# درون‌یابی تکه‌ای
# =========================================================

def piecewise_linear_interpolation(x_points, y_points, x):
    """
    درون‌یابی خطی تکه‌ای
    """
    n = len(x_points)

    if x < x_points[0] or x > x_points[-1]:
        raise ValueError("x خارج از بازه داده‌هاست.")

    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]


def piecewise_quadratic_interpolation(x_points, y_points, x):
    """
    درون‌یابی سهمی تکه‌ای
    برای هر بازه از 3 نقطه نزدیک استفاده می‌شود
    """
    n = len(x_points)

    if n < 3:
        raise ValueError("برای درون‌یابی سهمی حداقل 3 نقطه لازم است.")

    if x < x_points[0] or x > x_points[-1]:
        raise ValueError("x خارج از بازه داده‌هاست.")

    for i in range(n - 1):
        if x_points[i] <= x <= x_points[i + 1]:
            if i == 0:
                indices = [0, 1, 2]
            elif i == n - 2:
                indices = [n - 3, n - 2, n - 1]
            else:
                indices = [i - 1, i, i + 1]

            xp = [x_points[j] for j in indices]
            yp = [y_points[j] for j in indices]
            return lagrange_interpolation(xp, yp, x)

    return y_points[-1]


def cubic_spline_coefficients(x, y):
    """
    محاسبه ضرایب اسپلاین مکعبی طبیعی
    """
    n = len(x)
    a = np.array(y, dtype=float)
    h = np.diff(x)

    alpha = np.zeros(n)
    for i in range(1, n - 1):
        alpha[i] = (3 / h[i]) * (a[i + 1] - a[i]) - (3 / h[i - 1]) * (a[i] - a[i - 1])

    l = np.ones(n)
    mu = np.zeros(n)
    z = np.zeros(n)

    for i in range(1, n - 1):
        l[i] = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1]
        mu[i] = h[i] / l[i]
        z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]

    c = np.zeros(n)
    b = np.zeros(n - 1)
    d = np.zeros(n - 1)

    for j in range(n - 2, -1, -1):
        c[j] = z[j] - mu[j] * c[j + 1] if j < n - 1 else 0
        if j < n - 1:
            b[j] = (a[j + 1] - a[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3
            d[j] = (c[j + 1] - c[j]) / (3 * h[j])

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


def cubic_spline_interpolation(x_points, y_points, x):
    """
    درون‌یابی با اسپلاین مکعبی طبیعی
    """
    if x < x_points[0] or x > x_points[-1]:
        raise ValueError("x خارج از بازه داده‌هاست.")

    a, b, c, d = 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]


# =========================================================
# فصل 4 : انتگرال‌گیری رامبرگ
# =========================================================

def trapezoidal_rule(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


def romberg_integration(f, a, b, max_level=5, tol=1e-10):
    """
    انتگرال‌گیری رامبرگ
    """
    R = np.zeros((max_level, max_level), dtype=float)

    for i in range(max_level):
        n = 2**i
        R[i, 0] = trapezoidal_rule(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_iteration_table(results, max_rows=10):
    """
    چاپ بخشی از جدول تکرارها
    """
    if not results:
        print("نتیجه‌ای برای نمایش وجود ندارد.")
        return

    keys = list(results[0].keys())
    print("-" * 100)
    print(" | ".join(f"{k:^15}" for k in keys))
    print("-" * 100)

    for row in results[:max_rows]:
        print(" | ".join(f"{str(round(v, 10)) if isinstance(v, float) else str(v):^15}" for v in row.values()))
        if len(results) > max_rows:
        print("...")
    print("-" * 100)


# =========================================================
# مثال‌های تست
# =========================================================

if name == "__main__":
    # -------------------------------
    # فصل 2 : حل معادلات غیرخطی
    # -------------------------------
    f = lambda x: x**3 - x - 2
    df = lambda x: 3*x**2 - 1
    g = lambda x: (x + 2)**(1/3)

    print("\n===== روش دوبخشی =====")
    root, res = bisection(f, 1, 2)
    print("ریشه تقریبی:", root)
    print_iteration_table(res)

    print("\n===== روش نابجایی =====")
    root, res = false_position(f, 1, 2)
    print("ریشه تقریبی:", root)
    print_iteration_table(res)

    print("\n===== روش نقطه ثابت =====")
    root, res = fixed_point(g, 1.5)
    print("ریشه تقریبی:", root)
    print_iteration_table(res)

    print("\n===== روش نیوتون-رافسون =====")
    root, res = newton_raphson(f, df, 1.5)
    print("ریشه تقریبی:", root)
    print_iteration_table(res)

    print("\n===== روش وتری =====")
    root, res = secant(f, 1, 2)
    print("ریشه تقریبی:", root)
    print_iteration_table(res)

    # -------------------------------
    # فصل 3 : درون‌یابی
    # -------------------------------
    x_data = [1, 2, 3, 4]
    y_data = [1, 4, 9, 16]
    x_val = 2.5

    print("\n===== لاگرانژ =====")
    y_interp = lagrange_interpolation(x_data, y_data, x_val)
    print(f"f({x_val}) ≈ {y_interp}")

    print("\n===== نیوتون تقسیم‌شده =====")
    y_interp, table = newton_interpolation(x_data, y_data, x_val)
    print(f"f({x_val}) ≈ {y_interp}")
    print("جدول تفاضلات تقسیم‌شده:")
    print(table)

    x_eq = [0, 1, 2, 3]
    y_eq = [1, 2, 4, 8]

    print("\n===== نیوتون پیشرو =====")
    y_interp, table = newton_forward(x_eq, y_eq, 1.5)
    print("مقدار تقریبی:", y_interp)
    print(table)

    print("\n===== نیوتون پسرو =====")
    y_interp, table = newton_backward(x_eq, y_eq, 2.5)
    print("مقدار تقریبی:", y_interp)
    print(table)

    print("\n===== درون‌یابی معکوس =====")
    x_inv = inverse_interpolation(x_data, y_data, 6.25)
    print("x تقریبی برای y=6.25 :", x_inv)

    print("\n===== چندجمله‌ای چبیشف =====")
    print("T_4(0.5) =", chebyshev_polynomial(4, 0.5))
    print("گره‌های چبیشف در [0,1] برای n=5:")
    print(chebyshev_nodes(0, 1, 5))

    print("\n===== درون‌یابی تکه‌ای خطی =====")
    print(piecewise_linear_interpolation(x_data, y_data, 2.5))

    print("\n===== درون‌یابی تکه‌ای سهمی =====")
    print(piecewise_quadratic_interpolation(x_data, y_data, 2.5))

    print("\n===== اسپلاین مکعبی =====")
    print(cubic_spline_interpolation(x_data, y_data, 2.5))

    # -------------------------------
    # فصل 4 : انتگرال‌گیری رامبرگ
    # -------------------------------
    f_int = lambda x: math.sin(x)
    integral, table = romberg_integration(f_int, 0, math.pi, max_level=6)
    print("\n===== انتگرال‌گیری رامبرگ =====")
    print("مقدار تقریبی انتگرال:", integral)
    print("جدول رامبرگ:")
    print(table)