#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.1527901621; 
    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;
}
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgovLyDlhoblkajnjofjga7lrprnvqkKI2lmbmRlZiBNX1BJCiNkZWZpbmUgTV9QSSAzLjE0MTU5MjY1MzU4OTc5MzIzODQ2CiNlbmRpZgoKLy8gLS0tIOmWouaVsOWumue+qSAtLS0KLy8g5a++6LGh6Zai5pWwOiB5ID0gc2luKHgpIC8geApkb3VibGUgZihkb3VibGUgeCkgewogICAgaWYgKGZhYnMoeCkgPCAxZS05KSByZXR1cm4gMS4wOyAvLyB4PTDku5jov5Hjga7mpbXpmZDlh6bnkIYKICAgIHJldHVybiBzaW4oeCkgLyB4Owp9CgovLyDlsI7plqLmlbDjga7nnJ/lgKQ6IHknID0gKHgqY29zKHgpIC0gc2luKHgpKSAvIHheMgpkb3VibGUgZGZfdHJ1ZShkb3VibGUgeCkgewogICAgaWYgKGZhYnMoeCkgPCAxZS05KSByZXR1cm4gMC4wOwogICAgcmV0dXJuICh4ICogY29zKHgpIC0gc2luKHgpKSAvICh4ICogeCk7Cn0KCi8vIC0tLSAoMSkg5b6u5YiG44Gu5a6f6KGM6Zai5pWwICg1562J5YiGLCAxMOetieWIhueUqCkgLS0tCnZvaWQgcnVuX2RpZmZlcmVudGlhdGlvbihpbnQgTikgewogICAgZG91YmxlIHhfc3RhcnQgPSAyLjA7CiAgICBkb3VibGUgeF9lbmQgPSA0LjA7CiAgICBkb3VibGUgaCA9ICh4X2VuZCAtIHhfc3RhcnQpIC8gTjsKICAgIGludCBudW1fcG9pbnRzID0gTiArIDE7CgogICAgZG91YmxlIHhbMTAwXTsgCiAgICBkb3VibGUgeVsxMDBdOwogICAgZG91YmxlIGR5XzNwdFsxMDBdOwogICAgZG91YmxlIGR5XzVwdFsxMDBdOwoKICAgIHByaW50ZigiXG49PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwogICAgcHJpbnRmKCIg4piFIOW+ruWIhuiqsumhjDogJWQg562J5YiGIChoID0gJS40ZilcbiIsIE4sIGgpOwogICAgcHJpbnRmKCI9PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT1cbiIpOwoKICAgIC8vIOODh+ODvOOCv+OBruS9nOaIkAogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBudW1fcG9pbnRzOyBpKyspIHsKICAgICAgICB4W2ldID0geF9zdGFydCArIGkgKiBoOwogICAgICAgIHlbaV0gPSBmKHhbaV0pOwogICAgfQoKICAgIC8vIOW+ruS/guaVsOOBruioiOeulwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCBudW1fcG9pbnRzOyBpKyspIHsKICAgICAgICAvLyAtLS0gM+eCueW+ruWIhuWFrOW8jyAtLS0KICAgICAgICBpZiAoaSA9PSAwKSB7CiAgICAgICAgICAgIC8vIOWJjemAsuW3ruWIhgogICAgICAgICAgICBkeV8zcHRbaV0gPSAoLTMqeVtpXSArIDQqeVtpKzFdIC0geVtpKzJdKSAvICgyKmgpOwogICAgICAgIH0gZWxzZSBpZiAoaSA9PSBudW1fcG9pbnRzIC0gMSkgewogICAgICAgICAgICAvLyDlvozpgIDlt67liIYKICAgICAgICAgICAgZHlfM3B0W2ldID0gKHlbaS0yXSAtIDQqeVtpLTFdICsgMyp5W2ldKSAvICgyKmgpOwogICAgICAgIH0gZWxzZSB7CiAgICAgICAgICAgIC8vIOS4reWkruW3ruWIhgogICAgICAgICAgICBkeV8zcHRbaV0gPSAoeVtpKzFdIC0geVtpLTFdKSAvICgyKmgpOwogICAgICAgIH0KCiAgICAgICAgLy8gLS0tIDXngrnlvq7liIblhazlvI8gLS0tCiAgICAgICAgaWYgKGkgPT0gMCkgewogICAgICAgICAgICAvLyDnq6/ngrkgKHgxKQogICAgICAgICAgICBkeV81cHRbaV0gPSAoLTExKnlbaV0gKyAxOCp5W2krMV0gLSA5KnlbaSsyXSArIDIqeVtpKzNdKSAvICg2KmgpOwogICAgICAgIH0gZWxzZSBpZiAoaSA9PSAxKSB7CiAgICAgICAgICAgIC8vIOerr+eCueOBrumaoyAoeDIpCiAgICAgICAgICAgIGR5XzVwdFtpXSA9ICgtMip5W2ktMV0gLSAzKnlbaV0gKyA2KnlbaSsxXSAtIHlbaSsyXSkgLyAoNipoKTsKICAgICAgICB9IGVsc2UgaWYgKGkgPT0gbnVtX3BvaW50cyAtIDIpIHsKICAgICAgICAgICAgLy8g57WC54K544Gu6ZqjICh4bi0xKQogICAgICAgICAgICBkeV81cHRbaV0gPSAoeVtpLTJdIC0gNip5W2ktMV0gKyAzKnlbaV0gKyAyKnlbaSsxXSkgLyAoNipoKTsKICAgICAgICB9IGVsc2UgaWYgKGkgPT0gbnVtX3BvaW50cyAtIDEpIHsKICAgICAgICAgICAgLy8g57WC54K5ICh4bikKICAgICAgICAgICAgZHlfNXB0W2ldID0gKC0yKnlbaS0zXSArIDkqeVtpLTJdIC0gMTgqeVtpLTFdICsgMTEqeVtpXSkgLyAoNipoKTsKICAgICAgICB9IGVsc2UgewogICAgICAgICAgICAvLyDkuK3lpK4KICAgICAgICAgICAgZHlfNXB0W2ldID0gKHlbaS0yXSAtIDgqeVtpLTFdICsgOCp5W2krMV0gLSB5W2krMl0pIC8gKDEyKmgpOwogICAgICAgIH0KICAgIH0KCiAgICAvLyDntZDmnpzooajnpLoKICAgIHByaW50ZigiICAgeCAgICAgIHkoeCkgICAgIDPngrnlhazlvI8gICAgNeeCueWFrOW8jyAgICAg55yf5YCkICAgICAgIDPngrnoqqTlt64gICAgIDXngrnoqqTlt65cbiIpOwogICAgcHJpbnRmKCItLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tXG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbnVtX3BvaW50czsgaSsrKSB7CiAgICAgICAgZG91YmxlIHRydWVfdmFsID0gZGZfdHJ1ZSh4W2ldKTsKICAgICAgICBwcmludGYoIiU2LjJmICAlOC41ZiAgJTkuNWYgICU5LjVmICAlOS41ZiAgJTEwLjJlICAlMTAuMmVcbiIsIAogICAgICAgICAgICAgICB4W2ldLCB5W2ldLCBkeV8zcHRbaV0sIGR5XzVwdFtpXSwgdHJ1ZV92YWwsIAogICAgICAgICAgICAgICBmYWJzKGR5XzNwdFtpXSAtIHRydWVfdmFsKSwgZmFicyhkeV81cHRbaV0gLSB0cnVlX3ZhbCkpOwogICAgfQp9CgovLyAtLS0gKDIpIOepjeWIhuOBruWun+ihjOmWouaVsCAoNOetieWIhiwgMTDnrYnliIbnlKgpIC0tLQp2b2lkIHJ1bl9pbnRlZ3JhdGlvbihpbnQgTikgewogICAgZG91YmxlIHhfc3RhcnQgPSAyLjA7CiAgICBkb3VibGUgeF9lbmQgPSA0LjA7CiAgICBkb3VibGUgaCA9ICh4X2VuZCAtIHhfc3RhcnQpIC8gTjsKICAgIGludCBudW1fcG9pbnRzID0gTiArIDE7CgogICAgZG91YmxlIHhbMTAwXTsKICAgIGRvdWJsZSB5WzEwMF07CgogICAgcHJpbnRmKCJcblxuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKICAgIHByaW50ZigiIOKYhSDnqY3liIboqrLpoYw6ICVkIOetieWIhiAoaCA9ICUuNGYpXG4iLCBOLCBoKTsKICAgIHByaW50ZigiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09XG4iKTsKCiAgICAvLyDjg4fjg7zjgr/jga7kvZzmiJDjgajooajnpLogKOmbouaVo+ODh+ODvOOCv+OCkueiuuiqjeOBmeOCi+OBn+OCgSkKICAgIHByaW50ZigiIFvkvb/nlKjjgZnjgovpm6LmlaPjg4fjg7zjgr9dXG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgbnVtX3BvaW50czsgaSsrKSB7CiAgICAgICAgeFtpXSA9IHhfc3RhcnQgKyBpICogaDsKICAgICAgICB5W2ldID0gZih4W2ldKTsKICAgICAgICBwcmludGYoIiAgeFslZF09JS4yZiwgeVslZF09JS41ZlxuIiwgaSwgeFtpXSwgaSwgeVtpXSk7CiAgICB9CgogICAgLy8gMS4g6ZW35pa55b2i5YWs5byPICjlt6bnq6/ln7rmupbjgajku67lrpo6IHkwICsgLi4uICsgeW4tMSkKICAgIGRvdWJsZSBTX3JlY3QgPSAwLjA7CiAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgIFNfcmVjdCArPSBoICogeVtqXTsKICAgIH0KCiAgICAvLyAyLiDlj7DlvaLlhazlvI8KICAgIGRvdWJsZSBTX3RyYXAgPSAwLjA7CiAgICBmb3IgKGludCBqID0gMDsgaiA8IE47IGorKykgewogICAgICAgIFNfdHJhcCArPSBoICogKHlbal0gKyB5W2orMV0pIC8gMi4wOwogICAgfQoKICAgIC8vIDMuIOOCt+ODs+ODl+OCveODs+OBruWFrOW8jwogICAgZG91YmxlIFNfc2ltcCA9IDAuMDsKICAgIGlmIChOICUgMiA9PSAwKSB7IC8vIOWBtuaVsOWIhuWJsuOBruOBv+mBqeeUqOWPr+iDvQogICAgICAgIGRvdWJsZSBzdW1fb2RkID0gMC4wOwogICAgICAgIGRvdWJsZSBzdW1fZXZlbiA9IDAuMDsKCiAgICAgICAgZm9yIChpbnQgaSA9IDE7IGkgPCBOOyBpKyspIHsKICAgICAgICAgICAgaWYgKGkgJSAyICE9IDApIHN1bV9vZGQgKz0geVtpXTsgICAvLyB5MSwgeTMuLi4KICAgICAgICAgICAgZWxzZSAgICAgICAgICAgIHN1bV9ldmVuICs9IHlbaV07ICAvLyB5MiwgeTQuLi4KICAgICAgICB9CiAgICAgICAgU19zaW1wID0gKGggLyAzLjApICogKHlbMF0gKyB5W05dICsgNC4wICogc3VtX29kZCArIDIuMCAqIHN1bV9ldmVuKTsKICAgIH0gZWxzZSB7CiAgICAgICAgcHJpbnRmKCIgKOKAu+WIhuWJsuaVsE49JWQg44Gv5YG25pWw44Gn44Gv44Gq44GE44Gf44KB44CB44K344Oz44OX44K944Oz44Gu5YWs5byP44Gv6YGp55So44Gn44GN44G+44Gb44KTKVxuIiwgTik7CiAgICB9CgogICAgLy8g57WQ5p6c6KGo56S6CiAgICBwcmludGYoIlxuIFvnqY3liIblgKTjga7mr5TovINdXG4iKTsKICAgIHByaW50ZigiICDplbfmlrnlvaLlhazlvI8gOiAlLjhmXG4iLCBTX3JlY3QpOyAgCiAgICBwcmludGYoIiAg5Y+w5b2i5YWs5byPICAgOiAlLjhmXG4iLCBTX3RyYXApOwogICAgaWYoTiAlIDIgPT0gMCkgewogICAgICAgIHByaW50ZigiICDjgrfjg7Pjg5fjgr3jg7Pms5UgOiAlLjhmXG4iLCBTX3NpbXApOyAKICAgIH0KICAgIAogICAgZG91YmxlIHRydWVfaW50ZWdyYWwgPSAwLjE1Mjc5MDE2MjE7IAogICAgcHJpbnRmKCIgICjlj4LogIPnnJ/lgKQpIDogJS44ZlxuIiwgdHJ1ZV9pbnRlZ3JhbCk7Cn0KCmludCBtYWluKCkgewogICAgLy8gMS4g5b6u5YiG44Gu6Kqy6aGMICg1562J5YiGLCAxMOetieWIhikKICAgIHJ1bl9kaWZmZXJlbnRpYXRpb24oNSk7CiAgICBydW5fZGlmZmVyZW50aWF0aW9uKDEwKTsKCiAgICAvLyAyLiDnqY3liIbjga7oqrLpoYwgKDTnrYnliIYsIDEw562J5YiGKQogICAgcnVuX2ludGVncmF0aW9uKDQpOwogICAgcnVuX2ludGVncmF0aW9uKDEwKTsKCiAgICByZXR1cm4gMDsKfQo=