from math import acos, cos, sqrt, pi, fabs, log, ceil
import sympy as sp
import numpy as np


def cardano(a, b, c, d):
    p = c/a - b*b/(3*a*a)
    q = 2*b*b*b/(27*a*a*a) - b*c/(3*a*a) + d/a

    Q = (p/3)**3 + (q/2)**2
    if Q >= 0:
        print("WARNING : discriminant >= 0 !!!")

    phi = acos(3*q/(2*p) * sqrt(-3/p))

    y1 = 2 * sqrt(-p/3) * cos(phi/3)
    y2 = 2 * sqrt(-p/3) * cos(phi/3 + 2*pi/3)
    y3 = 2 * sqrt(-p/3) * cos(phi/3 - 2*pi/3)

    return [y1 - b/(3*a), y2 - b/(3*a), y3 - b/(3*a)]


def gauss3(func_x, a, b, alpha, _a):
    # степень полинома = 3
    n = 3
    # подсчет весовых моментов
    mu = [((b - _a) ** (i - alpha + 1) - (a - _a) ** (i - alpha + 1)) / (i - alpha + 1) for i in range(2*n)]

    # матрица и вектор для нахождения коэф. полинома
    system_matrix = [[mu[j + s] for j in range(n)] for s in range(n)]
    system_vector = [- mu[s] for s in range(n, 2*n)]

    # коэф. полинома
    pol_coef = np.linalg.solve(system_matrix, system_vector)

    # решение полинома
    sol = cardano(1, pol_coef[n-1], pol_coef[n-2], pol_coef[n-3])

    # новая система - для нахождения коэф. формулы Ai
    system_matrix = [[sol[j]**s for j in range(n)] for s in range(n)]
    system_vector = [mu[s] for s in range(n)]

    # коэф. формулы
    A_coef = np.linalg.solve(system_matrix, system_vector)

    # новая функция - от t
    func_t = lambda t: func_x(t + _a)

    # значение интеграла
    integral_sum = sum([A_coef[i] * func_t(sol[i]) for i in range(n)])
    return integral_sum


def compound_gauss(function, a, b, alpha, step):
    func_x = sp.lambdify(x, function, "math")

    intervals_number = ceil((b - a) / step)
    v = np.linspace(a, b, intervals_number + 1)

    integral_sum = sum([gauss3(func_x, v[j], v[j + 1], alpha, a) for j in range(intervals_number)])
    return integral_sum


def richardson(function_expr, a, b, alpha, start_step, eps):
    L, n = 2, 3
    step = start_step
    iterations = 0

    while True:
        iterations += 1
        h = [step, step/L, step/(L**2)]
        sh = [compound_gauss(function_expr, a, b, alpha, h[i]) for i in range(n)]
        m = - (log(fabs((sh[2] - sh[1]) / (sh[1] - sh[0]))) / log(L))

        system_matrix = [[1, - (h[i]**m), - (h[i]**(m+1))] for i in range(n)]
        sol = np.linalg.solve(system_matrix, sh)

        if (sol[1]*(h[-1]**m) + sol[2]*(h[-1]**(m + 1))) < eps:
            print("       J =", sol[0])
            print("min step =", h[-1])
            print("max step =", h[0])
            print("       m =", m)
            print("    iter =", iterations)
            return sol[0]

        step = h[-1]/L


x = sp.symbols('x')
sp.init_printing(use_unicode=True)

function1 = 6*sp.cos(1.5*x)*sp.exp(5*x/3) + 2*sp.sin(0.5*x)*sp.exp(-1.3*x) + 5.4*x

a, b = 3.5, 3.7
alpha = 2/3
integral_value = 2308.2875247384072
eps = 10e-6

step = b - a
print("С шага:", step)
richardson(function1, a, b, alpha, step, eps)