import numpy as np
import pandas as pd
from fractions import Fraction

# ================= Helper Functions =================

def parse_number(s):
    """Convert input into float, supports decimals and fractions like 1/2."""
    s = s.strip()
    try:
        if "/" in s:  # fraction input
            return float(Fraction(s))
        else:  # decimal or integer input
            return float(s)
    except Exception:
        raise ValueError(f"Invalid number: {s}")

def is_diagonally_dominant(A):
    """Check if matrix A is diagonally dominant."""
    n = len(A)
    for i in range(n):
        diag = abs(A[i][i])
        others = sum(abs(A[i][j]) for j in range(n) if j != i)
        if diag <= others:
            return False
    return True

def gauss_seidel(A, b, tol):
    """Gauss-Seidel iterative solver with difference vector and norm tracking."""
    n = len(b)
    x = np.zeros(n)  # initial guess = 0
    iteration_data = []
    k = 0

    # Initial guess row
    iteration_data.append([k] + list(x) + [[0.0] * n] + [np.nan])

    while True:
        x_new = x.copy()
        for i in range(n):
            s1 = sum(A[i][j] * x_new[j] for j in range(i))      # use updated values
            s2 = sum(A[i][j] * x[j] for j in range(i+1, n))     # use old values
            x_new[i] = (b[i] - s1 - s2) / A[i][i]

        # difference vector and norm
        diff_vec = (x_new - x).tolist()
        norm_val = np.linalg.norm(x_new - x, ord=np.inf)

        iteration_data.append([k + 1] + list(x_new) + [diff_vec] + [norm_val])

        if norm_val <= tol:
            return x_new, k + 1, pd.DataFrame(
                iteration_data,
                columns=["k"] + [f"x{i+1}" for i in range(n)] + ["x^(k+1) - x^(k)"] + ["‖x^(k+1) - x^(k)‖"]
            )

        x = x_new
        k += 1

# ================= Main Program =================

def run_program():
    print("==========================")
    print("    GAUSS-SEIDEL METHOD   ")
    print("==========================\n")

    print("DESCRIPTION:")
    print("The Gauss-Seidel Method is an iterative algorithm used to solve Ax = b.")
    print("It improves on the Jacobi Method by using the latest updated values of x")
    print("as soon as they are available within each iteration.")
    print("Stopping criterion: if ‖x^(k+1) − x^(k)‖ < ε, then accept x^(k+1).\n")

    print("INSTRUCTIONS:")
    print("1. Enter the number of equations/variables (n).")
    print("2. Enter the coefficients of matrix A row by row. Put numbers in a row using spaces.")
    print("   Example (for 3x3 system):")
    print("   Row 1: 10 -1 2")
    print("   Row 2: -1 11 -1")
    print("   Row 3: 2 -1 10")
    print("3. Enter the constants of vector b one by one.")
    print("   Example: b[1] = 6, b[2] = 25, b[3] = -11")
    print("4. Enter an error tolerance ε (e.g., 0.001).")
    print("5. The program will display all iterations with the difference vector")
    print("   and the norm until the stopping criterion is satisfied.\n")

    try:
        n = int(input("Enter the number of equations/variables: "))

        print("\nEnter the coefficients of matrix A (row by row):")
        A = []
        for i in range(n):
            row = [parse_number(x) for x in input(f"Row {i+1}: ").split()]
            if len(row) != n:
                raise ValueError("Each row must have exactly n coefficients.")
            A.append(row)
        A = np.array(A, dtype=float)

        print("\nEnter the constants of vector b:")
        b = []
        for i in range(n):
            val = parse_number(input(f"b[{i+1}]: "))
            b.append(val)
        b = np.array(b, dtype=float)

        tol = float(input("\nEnter the error tolerance ε: "))

        # Check diagonal dominance
        if not is_diagonally_dominant(A):
            print("\n❌ Gauss-Seidel method stopped: The matrix is not diagonally dominant, convergence is not guaranteed.\n")
            return

    except Exception as e:
        print("Invalid input:", e)
        return

    solution, iterations, table = gauss_seidel(A, b, tol)

    # ====== Iteration Table ======
    print("\nGauss-Seidel Method Iterations:")

    nvars = len(b)
    col_x_width = 12
    header = f"{'k':>3}  " + "  ".join([f"{'x'+str(i+1):>{col_x_width}}" for i in range(nvars)]) \
             + f"  {'x^(k+1) - x^(k)':>45}  {'‖x^(k+1) - x^(k)‖':>20}"
    print(header)
    print("-" * len(header))

    for row in table.itertuples(index=False):
        k = f"{int(row[0]):>3}"
        xs = "  ".join([f"{float(row[i+1]):>{col_x_width}.9f}" for i in range(nvars)])
        diff_list = row[-2]
        diff_fmt = "[" + ", ".join(f"{val: .9f}" for val in diff_list) + "]"
        diff_fmt = f"{diff_fmt:>45}"
        norm = "-" if pd.isna(row[-1]) else f"{float(row[-1]):>20.9f}"
        print(f"{k}  {xs}  {diff_fmt}  {norm}")

    # ====== Final Result ======
    print("\nFinal Solution:")
    for i, val in enumerate(solution, start=1):
        print(f"x{i} = {val:.9f}")

    print(f"\nStopped after {iterations} iterations because ‖x^(k+1) − x^(k)‖ < {tol}.\n")

# ================= Program Loop =================
if __name__ == "__main__":
    while True:
        run_program()
        choice = input("Do you want to solve another system? (y/n): ").strip().lower()
        if choice != "y":
            break
    input("\nPress Enter to exit...")