#include <stdio.h>
#include <math.h>
 
int main(void) {
    // データセットの定義 
    double x[] = {
        0.15708, 0.23982, 0.37400, 0.57120,
        0.82674, 1.04720, 1.23200, 1.43452
    };
    double y[] = {
        0.98769, 0.97138, 0.93087, 0.84125,
        0.67728, 0.50000, 0.33236, 0.13586
    };
 
    int n = sizeof(x) / sizeof(x[0]); // データ数
 
    // 正規方程式の各項の和を格納する変数
    // 行列 A = [[s_x4, s_x6], [s_x6, s_x8]]
    // ベクトル b = [s_yx2, s_yx4]
    double s_x4 = 0.0, s_x6 = 0.0, s_x8 = 0.0;
    double s_yx2 = 0.0, s_yx4 = 0.0;
 
    int i;
    for (i = 0; i < n; i++) {
        double x2 = x[i] * x[i];       // x^2
        double x4 = x2 * x2;           // x^4
        double x6 = x4 * x2;           // x^6
        double x8 = x4 * x4;           // x^8
        double Y = y[i] - 1.0;         // Y = y - 1 (定数項を移項)
 
        s_x4 += x4;
        s_x6 += x6;
        s_x8 += x8;
 
        s_yx2 += Y * x2;               // (y-1) * x^2
        s_yx4 += Y * x4;               // (y-1) * x^4
    }
 
    // 行列式 (Determinant) の計算 (2x2行列なので直接計算)
    double det = s_x4 * s_x8 - s_x6 * s_x6;
 
    if (fabs(det) < 1.0e-10) {
        printf("解けません（行列式が0に近いため）\n");
        return 1;
    }
 
    // クラメルの公式、または逆行列を用いて a1, a2 を求める
    double a1 = (s_yx2 * s_x8 - s_yx4 * s_x6) / det;
    double a2 = (s_x4 * s_yx4 - s_yx2 * s_x6) / det;
 
    // 結果の表示
    printf("=== 実行結果 ===\n");
    printf("データ数 n = %d\n", n);
    printf("係数 a1 = %f\n", a1);
    printf("係数 a2 = %f\n", a2);
    printf("\n求められた近似式:\n");
    printf("y = 1 + (%.5f)x^2 + (%.5f)x^4\n", a1, a2);
 
    // Excel等でのグラフ作成用データの出力（CSV形式）
    printf("\n=== グラフ作成用データ (CSV) ===\n");
    printf("x,y_measured,y_approx\n");
    for (i = 0; i < n; i++) {
        double y_approx = 1.0 + a1 * pow(x[i], 2) + a2 * pow(x[i], 4);
        printf("%f,%f,%f\n", x[i], y[i], y_approx);
    }
 
    return 0;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgppbnQgbWFpbih2b2lkKSB7CiAgICAvLyDjg4fjg7zjgr/jgrvjg4Pjg4jjga7lrprnvqkgCiAgICBkb3VibGUgeFtdID0gewogICAgICAgIDAuMTU3MDgsIDAuMjM5ODIsIDAuMzc0MDAsIDAuNTcxMjAsCiAgICAgICAgMC44MjY3NCwgMS4wNDcyMCwgMS4yMzIwMCwgMS40MzQ1MgogICAgfTsKICAgIGRvdWJsZSB5W10gPSB7CiAgICAgICAgMC45ODc2OSwgMC45NzEzOCwgMC45MzA4NywgMC44NDEyNSwKICAgICAgICAwLjY3NzI4LCAwLjUwMDAwLCAwLjMzMjM2LCAwLjEzNTg2CiAgICB9OwoKICAgIGludCBuID0gc2l6ZW9mKHgpIC8gc2l6ZW9mKHhbMF0pOyAvLyDjg4fjg7zjgr/mlbAKICAgIAogICAgLy8g5q2j6KaP5pa556iL5byP44Gu5ZCE6aCF44Gu5ZKM44KS5qC857SN44GZ44KL5aSJ5pWwCiAgICAvLyDooYzliJcgQSA9IFtbc194NCwgc194Nl0sIFtzX3g2LCBzX3g4XV0KICAgIC8vIOODmeOCr+ODiOODqyBiID0gW3NfeXgyLCBzX3l4NF0KICAgIGRvdWJsZSBzX3g0ID0gMC4wLCBzX3g2ID0gMC4wLCBzX3g4ID0gMC4wOwogICAgZG91YmxlIHNfeXgyID0gMC4wLCBzX3l4NCA9IDAuMDsKCiAgICBpbnQgaTsKICAgIGZvciAoaSA9IDA7IGkgPCBuOyBpKyspIHsKICAgICAgICBkb3VibGUgeDIgPSB4W2ldICogeFtpXTsgICAgICAgLy8geF4yCiAgICAgICAgZG91YmxlIHg0ID0geDIgKiB4MjsgICAgICAgICAgIC8vIHheNAogICAgICAgIGRvdWJsZSB4NiA9IHg0ICogeDI7ICAgICAgICAgICAvLyB4XjYKICAgICAgICBkb3VibGUgeDggPSB4NCAqIHg0OyAgICAgICAgICAgLy8geF44CiAgICAgICAgZG91YmxlIFkgPSB5W2ldIC0gMS4wOyAgICAgICAgIC8vIFkgPSB5IC0gMSAo5a6a5pWw6aCF44KS56e76aCFKQoKICAgICAgICBzX3g0ICs9IHg0OwogICAgICAgIHNfeDYgKz0geDY7CiAgICAgICAgc194OCArPSB4ODsKICAgICAgICAKICAgICAgICBzX3l4MiArPSBZICogeDI7ICAgICAgICAgICAgICAgLy8gKHktMSkgKiB4XjIKICAgICAgICBzX3l4NCArPSBZICogeDQ7ICAgICAgICAgICAgICAgLy8gKHktMSkgKiB4XjQKICAgIH0KCiAgICAvLyDooYzliJflvI8gKERldGVybWluYW50KSDjga7oqIjnrpcgKDJ4MuihjOWIl+OBquOBruOBp+ebtOaOpeioiOeulykKICAgIGRvdWJsZSBkZXQgPSBzX3g0ICogc194OCAtIHNfeDYgKiBzX3g2OwoKICAgIGlmIChmYWJzKGRldCkgPCAxLjBlLTEwKSB7CiAgICAgICAgcHJpbnRmKCLop6PjgZHjgb7jgZvjgpPvvIjooYzliJflvI/jgYww44Gr6L+R44GE44Gf44KB77yJXG4iKTsKICAgICAgICByZXR1cm4gMTsKICAgIH0KCiAgICAvLyDjgq/jg6njg6Hjg6vjga7lhazlvI/jgIHjgb7jgZ/jga/pgIbooYzliJfjgpLnlKjjgYTjgaYgYTEsIGEyIOOCkuaxguOCgeOCiwogICAgZG91YmxlIGExID0gKHNfeXgyICogc194OCAtIHNfeXg0ICogc194NikgLyBkZXQ7CiAgICBkb3VibGUgYTIgPSAoc194NCAqIHNfeXg0IC0gc195eDIgKiBzX3g2KSAvIGRldDsKCiAgICAvLyDntZDmnpzjga7ooajnpLoKICAgIHByaW50ZigiPT09IOWun+ihjOe1kOaenCA9PT1cbiIpOwogICAgcHJpbnRmKCLjg4fjg7zjgr/mlbAgbiA9ICVkXG4iLCBuKTsKICAgIHByaW50Zigi5L+C5pWwIGExID0gJWZcbiIsIGExKTsKICAgIHByaW50Zigi5L+C5pWwIGEyID0gJWZcbiIsIGEyKTsKICAgIHByaW50ZigiXG7msYLjgoHjgonjgozjgZ/ov5HkvLzlvI86XG4iKTsKICAgIHByaW50ZigieSA9IDEgKyAoJS41Zil4XjIgKyAoJS41Zil4XjRcbiIsIGExLCBhMik7CgogICAgLy8gRXhjZWznrYnjgafjga7jgrDjg6njg5XkvZzmiJDnlKjjg4fjg7zjgr/jga7lh7rlipvvvIhDU1blvaLlvI/vvIkKICAgIHByaW50ZigiXG49PT0g44Kw44Op44OV5L2c5oiQ55So44OH44O844K/IChDU1YpID09PVxuIik7CiAgICBwcmludGYoIngseV9tZWFzdXJlZCx5X2FwcHJveFxuIik7CiAgICBmb3IgKGkgPSAwOyBpIDwgbjsgaSsrKSB7CiAgICAgICAgZG91YmxlIHlfYXBwcm94ID0gMS4wICsgYTEgKiBwb3coeFtpXSwgMikgKyBhMiAqIHBvdyh4W2ldLCA0KTsKICAgICAgICBwcmludGYoIiVmLCVmLCVmXG4iLCB4W2ldLCB5W2ldLCB5X2FwcHJveCk7CiAgICB9CgogICAgcmV0dXJuIDA7Cn0=