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

void gauss3(double A[3][3], double b[3], double x[3]) {
    for (int i = 0; i < 3; i++) {
        // ピボット行の正規化
        double pivot = A[i][i];
        for (int j = 0; j < 3; j++) A[i][j] /= pivot;
        b[i] /= pivot;

        // 他の行の消去 (掃き出し処理)
        for (int k = 0; k < 3; k++) {
            if (k == i) continue;
            double factor = A[k][i];
            for (int j = 0; j < 3; j++) A[k][j] -= factor * A[i][j];
            b[k] -= factor * b[i];
        }
    }

    for (int i = 0; i < 3; i++) x[i] = b[i];
}

int main() {
    // 与えられたデータセット (8組)
    double xval[8] = {0.15708, 0.23982, 0.37400, 0.57120, 0.82674, 1.04720, 1.23200, 1.43452};
    double yval[8] = {0.98769, 0.97138, 0.93087, 0.84125, 0.67728, 0.50000, 0.33236, 0.13586};

    double S2 = 0, S4 = 0, S6 = 0, S8 = 0;
    double Sy = 0, S2y = 0, S4y = 0;

    // 正規方程式に必要な各合計値の計算
    for (int i = 0; i < 8; i++) {
        double x2 = xval[i] * xval[i];
        double x4 = x2 * x2;
        double x6 = x4 * x2;
        double x8 = x4 * x4;

        S2 += x2;
        S4 += x4;
        S6 += x6;
        S8 += x8;

        Sy += yval[i];
        S2y += x2 * yval[i];
        S4y += x4 * yval[i];
    }

    // 係数行列 A と右辺ベクトル b の設定
    double A[3][3] = {
        {8,  S2, S4},
        {S2, S4, S6},
        {S4, S6, S8}
    };

    double b[3] = {Sy, S2y, S4y};
    double sol[3];

    // 掃き出し法 (ガウスの消去法) で連立方程式を解く
    gauss3(A, b, sol);

    // 結果の出力
    printf("=========================================\n");
    printf("モデル (2): y = b1 + b2*x^2 + b3*x^4 計算結果\n");
    printf("=========================================\n");
    printf("求められた係数:\n");
    printf("  b1 = %.6f\n", sol[0]);
    printf("  b2 = %.6f\n", sol[1]);
    printf("  b3 = %.6f\n", sol[2]);

    // 残差平方和 (誤差の2乗和 E) の計算
    double E = 0;
    for (int i = 0; i < 8; i++) {
        double x2 = xval[i] * xval[i];
        double x4 = x2 * x2;
        double yhat = sol[0] + sol[1] * x2 + sol[2] * x4;
        double diff = yval[i] - yhat;
        E += diff * diff;
    }

    printf("\n近似式による偏差の2乗和 (誤差 E):\n");
    printf("  E = %.10f\n", E);

    // Excelプロット用の予測値出力
    printf("\n近似式による y の予測値 (Excelグラフ作成用):\n");
    printf("x 値\t\t計算された y 予測値\n");
    for (int i = 0; i < 8; i++) {
        double x2 = xval[i] * xval[i];
        double x4 = x2 * x2;
        double yhat = sol[0] + sol[1] * x2 + sol[2] * x4;
        printf("%.5f\t%.6f\n", xval[i], yhat);
    }
    printf("=========================================\n");

    return 0;
}