---

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# 1. Определение функции и пределов (Вариант 20, левый столбец)
def f(x):
    return 1/(x*(np.sqrt((4*x**2)-1)))

a, b = 1/np.sqrt(2), 1
n_values = [10, 100, 1000]

# Истинное значение через scipy
true_val, _ = quad(f, a, b)

# --- Реализация методов ---

def left_rect(f, a, b, n):
    h = (b - a) / n
    return h * sum(f(a + i * h) for i in range(n))

def right_rect(f, a, b, n):
    h = (b - a) / n
    return h * sum(f(a + (i + 1) * h) for i in range(n))

def mid_rect(f, a, b, n):
    h = (b - a) / n
    return h * sum(f(a + h/2 + i * h) for i in range(n))

def trapezoid(f, a, b, n):
    h = (b - a) / n
    return h * (0.5 * (f(a) + f(b)) + sum(f(a + i * h) for i in range(1, n)))

def simpson(f, a, b, n):
    if n % 2 != 0: n += 1  # Для Симпсона n должно быть четным
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return h/3 * (y[0] + y[-1] + 4 * sum(y[1:-1:2]) + 2 * sum(y[2:-2:2]))

def monte_carlo(f, a, b, n):
    x_rand = np.random.uniform(a, b, n)
    # Находим максимум функции на отрезке для создания "коробки"
    # Для sqrt(1-x^2) на [0, 0.5] максимум в f(0) = 1
    f_max = 1
    y_rand = np.random.uniform(0, f_max, n)
    under_curve = y_rand < f(x_rand)
    return (b - a) * f_max * (np.sum(under_curve) / n)

# --- Вычисления и вывод таблицы ---

methods = {
    "ЛПР": left_rect,
    "ППР": right_rect,
    "СПР": mid_rect,
    "ТР": trapezoid,
    "СИМП": simpson,
    "ММК": monte_carlo
}

print(f"{'n':<10} | {'Метод':<10} | {'Результат':<12} | {'Погрешность'}")
print("-" * 55)

for n in n_values:
    for name, func in methods.items():
        res = func(f, a, b, n)
        error = abs(true_val - res)
        print(f"{n:<10} | {name:<10} | {res:.8f} | {error:.8e}")
    print("-" * 55)

print(f"Истинное значение (True Area): {true_val:.8f}")

# --- Построение графика ---

x_plot = np.linspace(a - 0.1, b + 0.1, 400)
y_plot = f(x_plot)
x_fill = np.linspace(a, b, 100)
y_fill = f(x_fill)

plt.figure(figsize=(10, 6))
plt.plot(x_plot, y_plot, 'r', label='f(x) = $1/(x*(np.sqrt((4*x**2)-1)))$')
plt.fill_between(x_fill, y_fill, alpha=0.3, color='blue', label='Искомая площадь')
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.title('График функции и область интегрирования (Вариант 20)')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()