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

// 円周率の定義
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// --- 関数定義 ---
// 対象関数: y = sin(x) / x
double f(double x) {
    if (fabs(x) < 1e-9) return 1.0; // x=0付近の極限処理
    return sin(x) / x;
}

// 導関数の真値: y' = (x*cos(x) - sin(x)) / x^2
double df_true(double x) {
    if (fabs(x) < 1e-9) return 0.0;
    return (x * cos(x) - sin(x)) / (x * x);
}

// --- (1) 微分の実行関数 (5等分, 10等分用) ---
void run_differentiation(int N) {
    double x_start = 2.0;
    double x_end = 4.0;
    double h = (x_end - x_start) / N;
    int num_points = N + 1;

    double x[100]; 
    double y[100];
    double dy_3pt[100];
    double dy_5pt[100];

    printf("\n======================================================\n");
    printf(" ★ 微分課題: %d 等分 (h = %.4f)\n", N, h);
    printf("======================================================\n");

    // データの作成
    for (int i = 0; i < num_points; i++) {
        x[i] = x_start + i * h;
        y[i] = f(x[i]);
    }

    // 微係数の計算
    for (int i = 0; i < num_points; i++) {
        // --- 3点微分公式 ---
        if (i == 0) {
            // 前進差分
            dy_3pt[i] = (-3*y[i] + 4*y[i+1] - y[i+2]) / (2*h);
        } else if (i == num_points - 1) {
            // 後退差分
            dy_3pt[i] = (y[i-2] - 4*y[i-1] + 3*y[i]) / (2*h);
        } else {
            // 中央差分
            dy_3pt[i] = (y[i+1] - y[i-1]) / (2*h);
        }

        // --- 5点微分公式 ---
        if (i == 0) {
            // 端点 (x1)
            dy_5pt[i] = (-11*y[i] + 18*y[i+1] - 9*y[i+2] + 2*y[i+3]) / (6*h);
        } else if (i == 1) {
            // 端点の隣 (x2)
            dy_5pt[i] = (-2*y[i-1] - 3*y[i] + 6*y[i+1] - y[i+2]) / (6*h);
        } else if (i == num_points - 2) {
            // 終点の隣 (xn-1)
            dy_5pt[i] = (y[i-2] - 6*y[i-1] + 3*y[i] + 2*y[i+1]) / (6*h);
        } else if (i == num_points - 1) {
            // 終点 (xn)
            dy_5pt[i] = (-2*y[i-3] + 9*y[i-2] - 18*y[i-1] + 11*y[i]) / (6*h);
        } else {
            // 中央
            dy_5pt[i] = (y[i-2] - 8*y[i-1] + 8*y[i+1] - y[i+2]) / (12*h);
        }
    }

    // 結果表示
    printf("   x      y(x)     3点公式    5点公式     真値       3点誤差     5点誤差\n");
    printf("-------------------------------------------------------------------------\n");
    for (int i = 0; i < num_points; i++) {
        double true_val = df_true(x[i]);
        printf("%6.2f  %8.5f  %9.5f  %9.5f  %9.5f  %10.2e  %10.2e\n", 
               x[i], y[i], dy_3pt[i], dy_5pt[i], true_val, 
               fabs(dy_3pt[i] - true_val), fabs(dy_5pt[i] - true_val));
    }
}

// --- (2) 積分の実行関数 (4等分, 10等分用) ---
void run_integration(int N) {
    double x_start = 2.0;
    double x_end = 4.0;
    double h = (x_end - x_start) / N;
    int num_points = N + 1;

    double x[100];
    double y[100];

    printf("\n\n======================================================\n");
    printf(" ★ 積分課題: %d 等分 (h = %.4f)\n", N, h);
    printf("======================================================\n");

    // データの作成と表示 (離散データを確認するため)
    printf(" [使用する離散データ]\n");
    for (int i = 0; i < num_points; i++) {
        x[i] = x_start + i * h;
        y[i] = f(x[i]);
        printf("  x[%d]=%.2f, y[%d]=%.5f\n", i, x[i], i, y[i]);
    }

    // 1. 長方形公式 (左端基準と仮定: y0 + ... + yn-1)
    double S_rect = 0.0;
    for (int j = 0; j < N; j++) {
        S_rect += h * y[j];
    }

    // 2. 台形公式
    double S_trap = 0.0;
    for (int j = 0; j < N; j++) {
        S_trap += h * (y[j] + y[j+1]) / 2.0;
    }

    // 3. シンプソンの公式
    double S_simp = 0.0;
    if (N % 2 == 0) { // 偶数分割のみ適用可能
        double sum_odd = 0.0;
        double sum_even = 0.0;

        for (int i = 1; i < N; i++) {
            if (i % 2 != 0) sum_odd += y[i];   // y1, y3...
            else            sum_even += y[i];  // y2, y4...
        }
        S_simp = (h / 3.0) * (y[0] + y[N] + 4.0 * sum_odd + 2.0 * sum_even);
    } else {
        printf(" (※分割数N=%d は偶数ではないため、シンプソンの公式は適用できません)\n", N);
    }

    // 結果表示
    printf("\n [積分値の比較]\n");
    printf("  長方形公式 : %.8f\n", S_rect);  
    printf("  台形公式   : %.8f\n", S_trap);
    if(N % 2 == 0) {
        printf("  シンプソン法 : %.8f\n", S_simp); 
    }
    
    double true_integral = 0.165809716; 
    printf("  (参考真値) : %.8f\n", true_integral);
}

int main() {
    // 1. 微分の課題 (5等分, 10等分)
    run_differentiation(5);
    run_differentiation(10);

    // 2. 積分の課題 (4等分, 10等分)
    run_integration(4);
    run_integration(10);

    return 0;
}