import numpy as np
import matplotlib.pyplot as plt
def simpson_38_plot(f, a, b, n):
    if n % 3 !=0:
        raise ValueError(" To n πρέπει να ειναι πολλαπλάσιο του 3")

    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    fx = f(x)
    total = fx[0] + fx[-1]
    for i in range(1,n):
        if i%3 == 0:
            total += 2* fx[i]
        else:
            total += 3* fx[i]
    result = (3 * h / 8) * total
    print(f"h = {h}\n")
    print("i\t x_i\t f(x_i)\t coefficient\t term")
    print("-" *50)
    
    for i in range(n+1):
        if i == 0 or i == n:
            coeff = 1
        elif i % 3 == 0:
            coeff = 2
        else:
            coeff = 3
        term = coeff * fx[i]
        print(f"{i}\t {x[i]:.2f}\t {fx[i]:.2f}\t\t {coeff}\t\t {term:.2f}")
    
    print("\n Telikos ypologismos:")
    print(f"(3h/8) * sum = (3 * {h} / 8) * {total} = {result}")
    x_fine = np.linspace(a, b, 500)
    plt.figure(figsize=(10,5))
    plt.plot(x_fine, f(x_fine), 'b', label='f(x)')
    plt.scatter(x, fx, color='red', zorder=5, label='sample of points')

    for i in range(0, n, 3):
        xi = x[i:i+4]
        yi = fx[i:i+4]
        poly = np.polyfit(xi, yi, 3)
        x_segment = np.linspace(xi[0], xi[-1], 100)
        plt.plot(x_segment, np.polyval(poly, x_segment), 'g', alpha=0.5)

    plt.title(f"Simpson 3/8 approximation (Itegral  {result:.4f})")
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.legend
    plt.grid(True)
    plt.show()
    return result

def f(x):
     return x**2
a = 0
b = 3
n = 6

result = simpson_38_plot(f, a, b, n)
print ("\n Προσεγγιστικά το ολοκλήρωμα θα ειναι:", result)