import numpy as np #библиотека для матриц, векторов и т.д.
import random as rand #ну, бля, догадался, думаю

max = 10 #максимальные и минимальные значения для рандома
min = -10
d = 1. #число d позволяет управлять кол-вом точек в исходных значениях (на вычисления не влияет)
N = 5 #кол-во элементов в массиве

#рандомим основные векторы
a = np.random.randint(min, max, N-1)/d
b = np.random.randint(min, max, N)/d
c = np.random.randint(min, max, N-1)/d

#рандомим N-1-ый и N-ый столбец, но учитываем общие точки с векторами диагоналей
p = np.random.randint(min, max, N)/d
p[N-1] = a[N-2]
p[N-2] = b[N-2]
p[N-3] = c[N-3]
q = np.random.randint(min, max, N)/d
q[N-1] = b[N-1]
q[N-2] = c[N-2]
x1 = np.ones(N) #массив из единиц размера N (будет являться решением нашей системы)

def Gen (N, mat, x): #имея на руках матрицу СЛАУ и вектор-решение (x - который будем потом искать) генерируем правый вектор столбец
    f = np.zeros(N)
    for i in range(N):
        sum = 0
        for j in range(N):
            sum = sum + mat[i, j]*x[j]
        f[i] = sum
    return f
def PrinAll (a, b, c, p, q, f): #информация о том, что у тебя сейчас есть - нахуй не надо по сути
    print("Векторы диагонали:")
    print(a)
    print(b)
    print(c)
    print("Правый вектор-столбец: ", f)
    print("Обидная строка и столюец: ")
    print(p)
    print(q)
    print()
    pass
def Mat (N, a, b, c, p, q): #генерация двумерного массива из всех этих векторов
    A = np.zeros((N, N))
    A[0, 0] = b[0]
    A[0, 1] = c[0]
    A[0, N-1] = q[0]
    A[0, N-2] = p[0]
    for i in range(1, N):
            A[i, i-1] = a[i-1]
            A[i, i] = b[i]
            if i<N-1:  A[i, i+1] = c[i]
            A[i, N-1] = q[i]
            A[i, N-2] = p[i]
    return A
f = Gen (N, Mat(N, a, b, c, p, q), x1)

def Step1 (N, a, b, c, p, q, f): #ПЕРВЫЙ ШАГ
    for i in range(N-2):
        R = 1/b[i]
        b[i] = 1
        c[i] = c[i]*R
        if i == N-3: p[N-3] = c[N-3]
        else: p[i] = p[i]*R
        q[i] = q[i]*R
        f[i] = f[i]*R
        R = a[i]
        a[i] = 0
        b[i+1] = b[i+1] - R*c[i]
        if i==N-3:
            p[N-2] = b[N-2]
        else:
            p[i+1] = p[i+1] - R*p[i]
            if i==N-4: c[N-3] = p[N-3]
        q[i+1] = q[i+1] - R*q[i]
        if i==N-3: c[N-2] = q[N-2]
        f[i+1] = f[i+1] - R*f[i]
    pass
def Step2 (N, a, b, c, p, q, f): #ВТОРОЙ ШАГ
    R = 1/b[N-1]
    b[N-1] = 1
    q[N-1] = b[N-1]
    a[N-2] = a[N-2]*R
    p[N-1] = a[N-2]
    f[N-1] = f[N-1]*R
    R = c[N-2]
    c[N-2] = 0
    q[N-2] = c[N-2]
    b[N-2] = b[N-2] - R*a[N-2]
    p[N-2] = b[N-2]
    f[N-2] = f[N-2] - R*f[N-1]
    pass
def Step3 (N, a, b, c, p, q, f): #ТРЕТИЙ ШАГ
    f[N-2] = f[N-2]/b[N-2]
    b[N-2] = 1
    p[N-2] = b[N-2]
    f[N-1] = f[N-1] - f[N-2]*p[N-1]
    p[N-1] = 0
    a[N-2] = 0
    for i in range(N-3, -1, -1):
            f[i] = f[i] - f[N-2]*p[i]
            p[i] = 0
            if i==N-3: c[N-3] = p[N-3]
            f[i] = f[i] - f[N-1]*q[i]
            q[i] = 0
    pass
def Step4 (N, a, c, f): #ХУЙ
    x = np.zeros(N)
    for i in range(N-1, N-4, -1):
        x[i] = f[i]
    for i in range(N-4, -1, -1):
        x[i] = f[i]-c[i]*x[i+1]
    return x

def solution(N, a, b, c, p, q, f): #Основная функция
    Step1 (N, a, b, c, p, q, f)
    Step2 (N, a, b, c, p, q, f)
    Step3 (N, a, b, c, p, q, f)
    return Step4 (N, a, c, f)

def ErrorsCal(n, mea, tr): #вычисляет n-мерный массив из относительных погрешностей
    er = np.zeros(n)
    for i in range(n):
        absolute = mea[i]-tr[i]
        er[i] = absolute/mea[i]
    return er

def Average(n, a): #вычисляет среднее арифметическое (в сочетании с предыдущей функцией найдет среднюю относительную погрешность)
    sum = 0
    for i in a:
        sum = sum + i
    return sum/n
print()
PrinAll (a, b, c, p, q, f)

A = Mat (N, a, b, c, p, q)
print(A) #пялимся на случайную матрицу
y = solution(N, a, b, c, p, q, f) #решаем
print()
print("Правильный ответ: ", x1) #сравниваем результат визуально (питон округляет очень маленькие погрешности при выводе)
print("Результат: ", y)
print()
print(Average(N, ErrorsCal(N, y, x1))) #получаем точную среднюю относительную погрешность