using namespace std;
#include <vector>
#include <cmath>

const double a = 1;
const double b = 27;

double integral_function(double x) {
    return pow(sin(x), 2);
}

double atiderivative(double x) {
    return -sin(2 * x) / 4 + x / 2;
}

double second_derivative(double x) {
    return 2 * (sin(2 * x) + x * cos(2 * x));
}

double get_step(double eps) {
    double max = 0;
    double cur;
    for (double i = a; i <= b; i += 0.01) {
        cur = abs(second_derivative(i));
        if (cur > max)
            max = cur;
    }

    double h = sqrt((12 * eps) / ((b - a)*max));
    return h;
}

double calculate_integral(double a, double b) {
    return atiderivative(b) - atiderivative(a);
}

double trapezoid(double a, double b, double h) {
    double psum = integral_function(a) + integral_function(b);
    int n = (b - a) / h;
    for (int i = 1; i < n; ++i)
        psum = psum + 2.0 * integral_function(a + i * h);

    psum = (h / 2) * psum;

    return psum;
}

void runge(double a, double b, double eps, double abs_error) {
    double n = 1 / sqrt(eps);
    double R;
    double val_n, val_2n;
    do {
        val_n = trapezoid(a, b, (b - a) / n);
        val_2n = trapezoid(a, b, (b - a) / (2 * n));
        R = fabs(val_n - val_2n) / 3;
        n *= 2;
    } while (R > eps);

    double ex_value = calculate_integral(a, b);
    printf("| %.12lf | %.5le | %.12lf |\n", abs_error, (b - a) / n, abs(val_2n - ex_value));
}

int main() {
    vector<double> errors;

    printf("Trapezoid method\n|  Eps  |     Step    | Exact value |  Error   |\n");
    for (double eps = 0.1; eps > 1E-10; eps *= 0.01) {
        double h = get_step(eps);
        double value = trapezoid(a, b, h);
        double integral = calculate_integral(a, b);
        double error = fabs(value - integral);
        printf("| %.0le | %.5le | %.9lf | %lf |\n", eps, h, integral, error);
        errors.push_back(error);
    }

    printf("\nRunge principle\n| Expected error |     Step    | Prodused error |\n");
    double eps = 0.1;
    for (int i = 0; i < errors.size(); ++i) {
        runge(a, b, eps, errors[i]);
        eps *= 0.01;
    }
}