from fractions import Fraction
import sympy 
import math

def isZero(x):
    return math.isclose(x, 0, abs_tol=1e-09)

def isZeroVector(v):
    return all(isZero(vi) for vi in v)

def addVector(v, w):
    return [vi + wi for vi, wi in zip(v, w)]

def subVector(v, w):
    return [vi - wi for vi, wi in zip(v, w)]

def mulScalarVector(alpha, v):
    return [alpha*vi for vi in v]

def equalVector(v, w):
    return isZeroVector(subVector(v, w))

def switchRow(A, i, k):    
    A[i], A[k] = A[k], A[i] #"di <-> dk"
    
def mulRow(A, i, num):  
    A[i] = mulScalarVector(num, A[i]) #"di = num*di"

def addRow(A, i, k, num):  
    A[i] = addVector(A[i], mulScalarVector(num, A[k])) #"di = di + num*dk"

def convertToFraction(A):
    if isinstance(A, list):
        return [convertToFraction(i) for i in A]
    else:
        return Fraction(A)
def printMatrix(A, sep=" "):
    if isinstance(A, list) and A:
        if isinstance(A[0], list):                  #list of list
            m, n = len(A), len(A[0])                #size of matrix
            widths = [max(len(str(di[j])) for di in A) for j in range(n)]
            rows = [sep.join(format(str(di[j]),  f">{widths[j]}") for j in range(n)) for di in A]

            print("[" + "\n".join((" [" if i > 0 else "[") + rows[i] + "]" for i in range(m)) + "]")
        else: #list
            print("[" + sep.join(str(q) for q in A) + "]")
    else:
        print(A)   

def Gauss_elimination(A):
    R = A.copy()             #Result
    m, n = len(R), len(R[0]) #Size of matrix
    row = col = 0
    while row < m:
        #Step 1: Determines the leftmost column that does not contain all zeroes
        while col < n and all(isZero(R[i][col]) for i in range(row, m)):
            col += 1
        if col == n: #Echelon matrix
            break       
        #Step 2: Select the first line with a non-zero term. Switch two rows
        pivot_row = row + [not isZero(R[i][col]) for i in range(row, m)].index(True)
        switchRow(R, row, pivot_row)  
        #Step 3: With the leading term from Step 2 being a≠0 , multiply the row containing it by 1/a
        mulRow(R, row, 1/R[row][col])    
        #Step 4: Add an appropriate multiple of the top line to each bottom line to turn the terms below the leading number into 0
        for i in range(row + 1, m):
            multiplier = R[i][col] / R[row][col]
            addRow(R, i, row, -multiplier)
        #Step 5: Repeat the steps until you get the echolon matrix.
        row += 1
    print("Ma trận bậc thang của A là: ")
    printMatrix(R)
    print("\n")
    return R
 

def additionalMatrix(A, b):    
    return [ai + [bi] for ai, bi in zip(A, b)] #Create additional matrix [A|b]


def back_substitution(A): #R is the echolon matrix of the complement matrix of the system of equations Ax = b
    m, n = len(A), len(A[0]) #Size of matrix 
    solution = [None] * (n - 1)  #vector solution (hidden order list)
    #Find the bottom line #0
    row = m - 1
    while row >= 0 and all(isZero(A[row][j]) for j in range(n)):
        row -= 1
    if row >= 0 and [not isZero(A[row][j]) for j in range(n)].index(True) == n - 1:
        print("Hệ PT vô nghiệm\n")
        return 
    lastposCol = n - 1
    while row >= 0:
        posCol = [not isZero(A[row][j]) for j in range(n)].index(True)
        for i in range(posCol, lastposCol): #Hide freedom
            solution[i] = sympy.symbols(f"x{i + 1}")
        Sum = sum(A[row][j]*solution[j] for j in range(posCol + 1, n - 1))
        solution[posCol] = (A[row][n - 1] - Sum) / A[row][posCol]
        lastposCol = posCol
        row -= 1
    check = True
    for i in range(len(solution)):
        if not isinstance(solution[i], Fraction):
            x = sympy.symbols('x')
            check = x in solution[i].free_symbols
        if check == False:
                break      
    if check == True:
        print("Hệ PT có nghiệm duy nhất là: ",)
    else:
        print("Hệ PT có nghiệm tổng quát là: ",)
    printMatrix(solution, sep=", ")

def createMatrixUnit(n):
    return [[1 if i == j else 0 for j in range(n)] for i in range(n)]

def mulScalarMatrix(scalar, matrix):
    return [[scalar * a for a in matrix_row] for matrix_row in matrix]





#Example in excercise chapter 1   
A1 = [[1, 2, -1, -1],       #4, -3 ,-1
      [2, 2, 1, 1],
      [3, 5, -2, -1]]

A2 =  [[1, -2, -1, 1],      # x1 = 9 - 5a, x2 = 4 - 3a, x3 = a
      [2, -3, 1, 6],
      [3, -5, 0, 7],
      [1, 0, 5, 9]] 

A3 = [[1, 2, 0, 2, 6],      # 2, 3, -2, -1
      [3, 5, -1, 6, 17],
      [2, 4, 1, 2, 12],
      [2, 0, -7, 11, 7]] 

A4 = [[2, -4,  -1, 1],      #No solution 
      [1, -3, 1, 1],
      [3, -5, -3, 2]]

A5 = [[1, 2,  -2, 3],       #x1= 5/7, x2 = 8/7 + a,  x3 = a
      [3, -1, 1, 1],
      [-1, 5, -5, 5]]

A8 = [[1, -2, 3, -3],      # -5, 5, 4
      [2, 2, 0, 0],
      [0, -3, 4, 1],
      [1, 0, 1, -1]]

A10 = [[1, -1, 1, -3, 0],   # x1 = -3a - B, x2 = -2a - 4B, x3 = a, x4 =B
       [2, -1, 4, -2, 0]]

A12 = [[1, 6, 4, 0],        # x1 = 11/4a,  x2 = -9/8a, x3 =a
       [2, 4, -1, 0],
       [-1, 2, 5, 0]]


#Additional matrix constructor from column b and matrix A. (additionalMatrix)
A = [[2, -4, 6],
     [1, -1, 1],
     [1, -3, 4]]
b = [8, -1, 0]

A6 = additionalMatrix(A, b)     # 3, 13, 9 

temp = Gauss_elimination(convertToFraction(A2))  

back_substitution(temp)