import numpy as np


def table6_3():
    # Integrate f(x) over [a, b] using simpson's rule
    #
    # Inputs:
    #       a = left endpoint of interval
    #       b = right endpoint of interval
    #       n = number of subintervals
    #

    a = 0
    b = 1
    exact = (np.exp(3) - 1) / 3

    # interv = number of subinternals
    interv = np.array([10, 20, 40, 80, 160, 320])
    nk = len(interv)
    S = np.zeros(nk)
    err = np.zeros(nk)
    for k in range(nk):
        n = int(interv[k])
        S[k] = simp(a, b, n)
        err[k] = abs(exact - S[k])

    print("\n RIE = relative iterative error \n")
    print(" k \t  n  \t\t I_S \t \t Error \t\t RIE")
    print(f" {1:d} \t {interv[0]:d}\t {S[0]:19.15f} \t {err[0]:5.1e}")
    # print('\\, %i \\, & ...')
    for k in range(1, nk):
        rie = abs((S[k] - S[k - 1]) / S[k])
        print(
            f" {k + 1:d} \t {interv[k]:d} \t {S[k]:19.15f} "
            f"\t {err[k]:5.1e}  \t {rie:5.1e}"
        )
        # print('\\, %i \\, & ...')
    print()


def f(x):
    y = np.exp(3 * x)
    return y


### composite simpson: n = number of subintervals
def simp(a, b, n):
    xd = np.linspace(a, b, n + 1)
    h = xd[1] - xd[0]
    s = f(xd[0])
    ff2 = 0.0
    for j in range(1, n, 2):
        ff2 = f(xd[j + 1])
        s = s + 4 * f(xd[j]) + 2 * ff2
    s = h * (s - ff2) / 3
    return s


if __name__ == "__main__":
    table6_3()
