def calc_integral(ys, h, a, b):
    func = lambda x: calc_polynominal(ys, x, h)
    #n = 7
    #k = 17280
    #H0 = 751
    #H1 = 3577
    #H2 = 1323
    #H3 = 2989
    #coeff_cotes = (Decimal(H0/k), Decimal(H1/k),
    #               Decimal(H2/k), Decimal(H3/k),
    #               Decimal(H3/k), Decimal(H2/k),
    #               Decimal(H1/k), Decimal(H0/k))
    n = 8
    k = 28350
    H0 = 939
    H1 = 5888
    H2 = -928
    H3 = 10496
    H4 = -4540
    coeff_cotes = (Decimal(H0/k), Decimal(H1/k),
               Decimal(H2/k), Decimal(H3/k), Decimal(H4/k),
               Decimal(H3/k), Decimal(H2/k),
               Decimal(H1/k), Decimal(H0/k))
    result = 0
    h_int = Decimal((b-a)/Decimal(n))
    x = a
    i = 0
    while x <= b:
        result += func(x) * coeff_cotes[i]
        print(i, x, result)
        x += h_int
        i += 1
    return Decimal(b-a)*result


def calc_rn(a, b, h):
    y_7d = lambda x_f: Decimal(-(92160*(-1+84*x_f**2-560*x_f**4+448*x_f**6))/(1+4*x_f**2)**7)
    pn = Decimal(1)
    n = 6
    h_int = Decimal((b-a)/Decimal(n))
    x_n = a
    y1 = y_7d(a)
    y2 = y_7d(b)
    if y1 > y2:
        x = a
    else:
        x = b
    for i in range(n):
        pn *= x - x_n
        x_n += h_int
    if y1 > y2:
        return (y1*pn)/Decimal(math.factorial(n))
    else:
        return (y2*pn)/Decimal(math.factorial(n))