import math
from decimal import Decimal


# Расчет значений исходной функции на промежутке [a; b] с шагом h
# h должен быть такой, что (b - a)/h - целое число
def calc_Ys_Xs(func, a, b, h):
    # Список со всеми значениями
    l_values = [[], []]
    # Кол-во значений X и Y
    k = int((b - a)/h) + 1
    x = Decimal(a)
    for i in range(k):
        l_values[0].append(x)
        l_values[1].append(func(x))
        x += h
    for i in range(2, k + 1):
        if len(l_values[i-1]) != 1:
            l_values.append([])
            for j in range(1, len(l_values[i - 1])):
                l_values[i].append(l_values[i - 1][j] - l_values[i - 1][j - 1])
        else:
            break
    return l_values


# Найти y(x) при помощи полинома (выбор полинома зависит от fl)
def calc_polyominal(l_values, x, h, fl = 0):
    # Грубо говоря, это значение произведения q(q-1)(q-2)...(q-n+1) для текущей
    # итерации
    n_q = Decimal(1)
    if fl == 0:
        q = Decimal((x - l_values[0][0])/h)
        y = Decimal(l_values[1][0])
    else:
        q = Decimal((x - l_values[0][len(l_values[0]) - 1])/h)
        y = Decimal(l_values[1][len(l_values[1]) - 1])
    # Цикл начинается с 2, т.к. 0 и 1 списки это Х и Y, соответственно
    for i in range(2, len(l_values)):
        if fl == 0:
            n_q *= q - i + 2
            y += (n_q/Decimal(math.factorial(i - 1)))*l_values[i][0]
        else:
            n_q *= q + i - 2
            y += (n_q/Decimal(math.factorial(i - 1)))*l_values[i][len(l_values[i]) - 1]
    return Decimal(y)


def calc_diff(func, x, eps):
    x = Decimal(x)
    return Decimal((func(x+eps) - func(x))/eps)


def first_diff(q):
    return Decimal(Decimal(2.5)*(Decimal(-0.157617)+Decimal(0.09025)*q-Decimal(0.027825)*(q**2)+Decimal(0.0045)*(q**3)-Decimal(0.000291667)*(q**4)))


def main():
    func = lambda x: Decimal(math.atan(1/(2*x)))
    a = Decimal(1.0)
    b = Decimal(5.0)
    n = Decimal(10)
    eps = Decimal(10**(-10))
    h = Decimal((b-a)/n)
    ys = calc_Ys_Xs(func, a, b, h)
    func_fd = lambda x: Decimal(calc_diff(func, x, Decimal(10**(-7))))
    print("|%-10s|%-10s|%-10s|%-10s|%-10s|%-10s|%-10s|%-10s|%-10s|%-10s|%-10s|"
          % ('X', 'Y', 'd^1 Y', 'd^2 Y', 'd^3 Y', 'd^4 Y', 'd^5 Y', 'd^6 Y',
              'd^7 Y', 'd^8 Y', 'd^9 Y'), "%-10s|" % 'd^10 Y', sep="")
    print("-" * 133)
    for i in range(len(ys[0])):
        print("|%-10.6f|%-10.6f" % (ys[0][i], ys[1][i]), sep="", end="")
        for j in range(2, len(ys)):
            try:
                print("|%-10.6f" % (ys[j][i]), end="")
            except IndexError:
                print("|     -    ", end="")
        print('|')
    print()
    print("|%-10s|%-10s|%-10s|%-10s|%-12s|%-10s|%-10s|%-13s|"
          % ('X', 'Y', 'Y\'', 'F(x)H\'', 'Погрешность\'', 'Y\'\'', 'F(x)H\'\'',
             'Погрешность\'\''))
    print("-" * 94)
    for x in (Decimal(1.0), Decimal(1.2), Decimal(4.8), Decimal(5.0)):
        y = func(x)
        q = Decimal((x - ys[0][0])/h)
        print("|%-10.6f|%-10.6f" % (x, y), sep="", end="")
        if x <= (a + b)/2:
            y_polynominal = lambda x: Decimal(calc_polyominal(ys, x, h))
            y_polynominal_fd = lambda x: Decimal(calc_diff(y_polynominal, x, eps))
        else:
            y_polynominal = lambda x: Decimal(calc_polyominal(ys, x, h, 1))
            y_polynominal_fd = lambda x: Decimal(calc_diff(y_polynominal, x, eps))
        y_fd = Decimal(func_fd(x))
        print('|%-10f' % y_fd, sep='', end='')
        y_poly_fd = Decimal(y_polynominal_fd(x))
        print('|%-10f' % y_poly_fd, sep='', end='')
        print('|%-12f' % abs(y_poly_fd-y_fd), sep='', end='')
        y_sd = Decimal(calc_diff(func_fd, x, Decimal(10**(-7))))
        print('|%-10f' % y_sd, sep='', end='')
        y_poly_sd = Decimal(calc_diff(y_polynominal_fd, x, eps))
        print('|%-10f' % y_poly_sd, sep='', end='')
        print('|%-13f|' % abs(y_poly_sd-y_sd), sep='', end='')
        if x <= (a + b)/2:
            print('|%-10f' % first_diff(q), sep='')
        else:
            print()


if __name__ == "__main__":
          main()