tanya4

import sympy as sym
import numpy as np

from fractions import Fraction

def rank(a):
    return len([1 for i in a if any(i)])

def gauss_method(a):
    n, m = len(a), len(a[0])
    for k in range(n):
        t = k
        if k < n and t < m and not a[k][t] or not any(a[k]):
            a = np.append(a, [a[k]], axis=0)
            a = np.delete(a, k, axis=0)
        for i in range(1+k, n):
            while t < m-1 and not a[k][t]:
                if a[k][t]:
                    break
                t += 1
            if t == m or not a[k][t]:
                continue

            K = -a[i][t]/a[k][t]
            for j in range(m):
                a[i][j] = a[k][j]*K + a[i][j]

    for x in range(n):
        if not any(a[x]):
            a = np.append(a, [a[x]], axis=0)
            a = np.delete(a, x, axis=0)
    return a

def find_free(a):
    n, m = len(a), len(a[0])
    res = set()
    for i in range(n):
        for j in range(m):
            if a[i][j] != 0:
                res.add(j)
                break
    return set(range(m)).difference(res)

def identity_matrix(a):
    n, m = len(a), len(a[0])
    for i in range(n):
        a[i] /= a[i][i]

    for i in range(n-1, -1, -1):
        for j in range(i-1, -1, -1):
            a[j] += -a[j][i]*a[i]

    return a

def fundamental_solve(a, general_variables, free_variables):
    n, m, k = len(a), len(a[0]), len(general_variables)
    ans = np.array([m*[Fraction(0.)] for x in range(len(free_variables))])
    for ind, x in enumerate(general_variables):
        for i, elem in enumerate(a[ind][k:]):
            if elem != 0:
                ans[i][x] = elem
    for i, x in enumerate(free_variables):
        ans[i][x] = 1
    return ans

def solve_slau(a):
    n, m = len(a), len(a[0])
    a = gauss_method(a)
    a = a[:rank(a)]
    free_variables = list(find_free(a))
    general_variables = list(set(range(m)).difference(set(free_variables)))
    temp = [[] for i in range(len(a))]
    for i in sorted(free_variables, reverse=True):
        temp = np.append(temp, -1*a[:, i:i+1], axis=1)
        a = np.delete(a, i, axis=1)
    a = np.append(a, temp[:, ::-1], axis=1)
    return fundamental_solve(identity_matrix(a), general_variables, free_variables)

def gram(a):
    n, m = len(a), len(a[0])
    b = np.array([[Fraction(0)]*m for x in range(n)])
    b[0] = a[0]
    for i in range(1, n):
        b[i] = a[i] + sum([(-np.dot(a[i], b[j])/np.dot(b[j], b[j]))*b[j] for j in range(i)])
    return b

def main():
    print("Введите размерность: ")
    n = int(input())
    print("Введите элементы матрицы: ")
    a = sym.Matrix([[Fraction(x) for x in input().split()] for i in range(n)])
    b = sym.matrices.eye(n, n)*sym.Symbol('λ')
    print('\nСоставляем характеристическое уравнение:')
    sym.pprint(a-b)
    print()
    equation = sym.Eq(sym.det(a-b), 0)
    roots = sym.roots(equation, sym.Symbol('λ'))
    print('Матрица линейного оператора a в базисе, состоящем из собственных векторов:')
    sym.pprint(sym.matrices.diag(*[x for x, i in sorted(roots.items()) for t in range(i)]))
    print()

    answer = sym.Matrix([])
    for step, root in enumerate(sorted(roots.keys()), start=1):
        # Находим собственные вектора
        temp = solve_slau(np.array(a - sym.matrices.eye(n, n)*root))
        # Ортонормируем собственные вектора
        temp = [sym.Matrix(i)/sym.Matrix(i).norm() for i in gram(temp[::-1])]
        for i in temp:
            answer = answer.col_insert(step, sym.Matrix(i))

    print('Ортнормированный базис')
    sym.pprint(answer)

if __name__ == '__main__':
    main()