import math


# Структура, описывающая сплайн на каждом сегменте сетки
class MySpline:
    def __init__(self, a, b, c, d, x):
        self.a = a
        self.b = b
        self.c = c
        self.d = d
        self.x = x

    def __str__(self):
        return "a = {0}, b = {1}, c = {2}, d = {3}".format(self.a,
                                                           self.b,
                                                           self.c,
                                                           self.d)


# Построение сплайна
# x - узлы сетки, должны быть упорядочены по возрастанию,
# кратные узлы запрещены
# y - значения функции в узлах сетки
# n - количество узлов сетки
def BuildSpline(x, y, n, my_h, my_n):
    # Инициализация массива сплайнов
    splines = [MySpline(0, 0, 0, 0, 0) for _ in range(0, n)]
    for i in range(0, n):
        splines[i].x = x[i]
        splines[i].a = y[i]

    splines[0].c = splines[n - 1].c = 0.0

    # Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки
    # для трехдиагональных матриц
    # Вычисление прогоночных коэффициентов - прямой ход метода прогонки
    alpha = [0.0 for _ in range(0, n - 1)]
    beta = [0.0 for _ in range(0, n - 1)]

    for i in range(1, n - 1):
        hi = x[i] - x[i - 1]
        hi1 = x[i + 1] - x[i]
        A = hi
        C = 2.0 * (hi + hi1)
        B = hi1
        F = 6.0 * ((y[i + 1] - y[i]) / hi1 - (y[i] - y[i - 1]) / hi)
        z = (A * alpha[i - 1] + C)
        alpha[i] = -B / z
        beta[i] = (F - A * beta[i - 1]) / z

    # Нахождение решения - обратный ход метода прогонки
    for i in range(n - 2, 0, -1):
        splines[i].c = alpha[i] * splines[i + 1].c + beta[i]

    # По известным коэффициентам c[i] находим значения b[i] и d[i]
    for i in range(n - 1, 0, -1):
        hi = x[i] - x[i - 1]
        splines[i].d = (splines[i].c - splines[i - 1].c) / hi
        splines[i].b = hi * (2.0 * splines[i].c +
                             splines[i - 1].c) / 6.0 + (y[i] - y[i - 1]) / hi

    print("Таблица коэффициентов сплайнов:")
    for i in range(0, n):
        print("i = {0}, {1}".format(i, splines[i]))

    print()

    print("Таблица значений сплайнов:")
    if my_n == 5:
        for i in range(0, n):
            print("i = {0}, {1}".format(i, get_spline_value(splines[i], my_h)))
    elif my_n == 25:
        for i in range(0, 5):
            i = 3 + 5 * i
            print("i = {0}, {1}".format(i, get_spline_value(splines[i], my_h)))
    else:
        for i in range(0, 5):
            i = 13 + 5 * i
            print("i = {0}, {1}".format(i, get_spline_value(splines[i], my_h)))

    return splines


def get_spline_value(spline, h):
    return spline.a + spline.b * h / 4 + spline.c * h ** 2 / 4 + spline.d * h ** 3 / 8


# Вычисление значения интерполированной функции в произвольной точке
def Interpolate(splines, x):
    if not splines:
        return None  # Если сплайны ещё не построены - возвращаем NaN

    n = len(splines)
    s = MySpline(0, 0, 0, 0, 0)

    # Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива
    if x <= splines[0].x:
        s = splines[0]
    # Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива
    elif x >= splines[n - 1].x:
        s = splines[n - 1]
    # Иначе x лежит между граничными точками сетки - производим
    # бинарный поиск нужного эл-та массива
    else:
        i = 0
        j = n - 1
        while i + 1 < j:
            k = i + (j - i) // 2
            if x <= splines[k].x:
                j = k
            else:
                i = k
        s = splines[j]

    dx = x - s.x
    # Вычисляем значение сплайна в заданной точке по схеме Горнера
    return s.a + (s.b + (s.c / 2.0 + s.d * dx / 6.0) * dx) * dx


def frange(start, stop, step):
    i = start
    while i <= stop:
        yield i
        i += step


def f(x):
    return math.log(2 * x - 1) - math.sin(pow(x, 2))


a = 1
b = 3
# n = 5
# n = 25
n = 5
h = (b - a) / n
x = [i for i in frange(a, b, h)]
y = [f(i) for i in x]
x2 = [i for i in frange(a+1, b+1, h)]
z = [f(i-0.5) for i in x2]
new_x = 2

func = "f(x) = ln(2x − 1) − sin x^2"

print("Интерполяция функции {0} на отрезке [{1}, {2}] с шагом {3}".format(
    func, a, b, h))
print("Количество узлов сетки: {0}".format(n))
print()
print("Узлы сетки:")
print("x = {0}".format(x))
print("Таблица значений функции f(a + ih):")
print("y = {0}".format(y))
print('Xci = {0}'.format(x2))
print("Таблица значений функции f(a + (i-0.5)*h):")
print("y = {0}".format(z))
print()
print("Интерполяция в точке x = {0}".format(new_x))
print()

spline = BuildSpline(x, y, len(x), h, n)

answer = Interpolate(spline, new_x)

print()
print("Значение функции {0} в точке x = {1} равно {2}".format(
    func, new_x, answer))
print("Табличное значение функции = {0}".format(f(new_x)))
print("Погрешность = {0}".format(abs(answer - f(new_x))))