#include <stdio.h>
#include <math.h>

// 円周率の定義
#define PI 3.14159265358979323846

// 対象の関数 f(x) = 1 / (2 + cos(pi * x))
double f(double x) {
    return 1.0 / (2.0 + cos(PI * x));
}

// ラグランジュ基底多項式 l_j(x) の計算
double lagrange_basis(int j, double x, const double nodes[], int n) {
    double l = 1.0;
    for (int k = 0; k <= n; k++) {
        if (k != j) {
            l *= (x - nodes[k]) / (nodes[j] - nodes[k]);
        }
    }
    return l;
}

// 補間多項式 P_n(x) の値の計算
double interpolant(double x, const double nodes[], const double values[], int n) {
    double p = 0.0;
    for (int j = 0; j <= n; j++) {
        p += values[j] * lagrange_basis(j, x, nodes, n);
    }
    return p;
}

int main(void) {
    int i, j, k;
    int n = 7; // 7次のチェビシェフ補間（節点数 8個）
    double nodes[8];
    double values[8];

    // 0. 積分の真値を高精度（十分な分割数による数値積分）で近似
    int ref_steps = 100000;
    double a = -1.0, b = 1.0;
    double h_ref = (b - a) / ref_steps;
    double true_integral = 0.5 * (f(a) + f(b));
    for (i = 1; i < ref_steps; i++) {
        true_integral += f(a + i * h_ref);
    }
    true_integral *= h_ref;
    printf("積分の真値 (近似): %.12f\n\n", true_integral);

    // ==========================================
    // (1) 7次のチェビシェフ補間の節点と値の計算
    // ==========================================
    for (k = 0; k <= n; k++) {
        // チェビシェフ節点: x_k = cos((2k+1)pi / (2(n+1)))
        nodes[k] = cos((2.0 * k + 1.0) * PI / (2.0 * (n + 1.0)));
        values[k] = f(nodes[k]);
    }

    printf("--- (1) 7次チェビシェフ補間の節点 ---\n");
    for (k = 0; k <= n; k++) {
        printf("x_%d = %10.6f, f(x_%d) = %10.6f\n", k, nodes[k], k, values[k]);
    }
    printf("\n");

    // ==========================================
    // (2) 補間多項式の [-1, 1] での積分と真値との誤差
    // ==========================================
    // 補間多項式を十分に細かい分割（例: 1000分割）で数値積分する
    int interp_steps = 1000;
    double h_interp = (b - a) / interp_steps;
    double poly_integral = 0.5 * (interpolant(a, nodes, values, n) + interpolant(b, nodes, values, n));
    for (i = 1; i < interp_steps; i++) {
        double xi = a + i * h_interp;
        poly_integral += interpolant(xi, nodes, values, n);
    }
    poly_integral *= h_interp;

    double error_interp = fabs(poly_integral - true_integral);

    printf("--- (2) 補間多項式の積分と誤差 ---\n");
    printf("補間多項式の積分値: %.12f\n", poly_integral);
    printf("真値との誤差      : %.2e\n\n", error_interp);

    // ==========================================
    // (3) 7分割・8分割の台形則による数値積分
    // ==========================================
    // 7分割の台形則
    int n7 = 7;
    double h7 = (b - a) / n7;
    double trap_7 = 0.5 * (f(-1.0) + f(1.0));
    for (i = 1; i < n7; i++) {
        trap_7 += f(-1.0 + i * h7);
    }
    trap_7 *= h7;
    double err_trap_7 = fabs(trap_7 - true_integral);

    // 8分割の台形則
    int n8 = 8;
    double h8 = (b - a) / n8;
    double trap_8 = 0.5 * (f(-1.0) + f(1.0));
    for (i = 1; i < n8; i++) {
        trap_8 += f(-1.0 + i * h8);
    }
    trap_8 *= h8;
    double err_trap_8 = fabs(trap_8 - true_integral);

    printf("--- (3) 台形則による数値積分と誤差 ---\n");
    printf("7分割の台形則 積分値: %.12f, 誤差: %.2e\n", trap_7, err_trap_7);
    printf("8分割の台形則 積分値: %.12f, 誤差: %.2e\n", trap_8, err_trap_8);

    return 0;
}