#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() {
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];
}
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("(2): y = b1 + b2*x^2 + b3*x^4 計算結果\n"); printf(" b1 = %.6f\n", sol
[0]); printf(" b2 = %.6f\n", sol
[1]); printf(" b3 = %.6f\n", sol
[2]);
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");
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;
}
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+Cgp2b2lkIGdhdXNzMyhkb3VibGUgQVszXVszXSwgZG91YmxlIGJbM10sIGRvdWJsZSB4WzNdKSB7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IDM7IGkrKykgewogICAgICAgIGRvdWJsZSBwaXZvdCA9IEFbaV1baV07CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCAzOyBqKyspIEFbaV1bal0gLz0gcGl2b3Q7CiAgICAgICAgYltpXSAvPSBwaXZvdDsKCiAgICAgICAgZm9yIChpbnQgayA9IDA7IGsgPCAzOyBrKyspIHsKICAgICAgICAgICAgaWYgKGsgPT0gaSkgY29udGludWU7CiAgICAgICAgICAgIGRvdWJsZSBmYWN0b3IgPSBBW2tdW2ldOwogICAgICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IDM7IGorKykgQVtrXVtqXSAtPSBmYWN0b3IgKiBBW2ldW2pdOwogICAgICAgICAgICBiW2tdIC09IGZhY3RvciAqIGJbaV07CiAgICAgICAgfQogICAgfQoKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgMzsgaSsrKSB4W2ldID0gYltpXTsKfQoKaW50IG1haW4oKSB7CiAgICBkb3VibGUgeHZhbFs4XSA9IHswLjE1NzA4LCAwLjIzOTgyLCAwLjM3NDAwLCAwLjU3MTIwLCAwLjgyNjc0LCAxLjA0NzIwLCAxLjIzMjAwLCAxLjQzNDUyfTsKICAgIGRvdWJsZSB5dmFsWzhdID0gezAuOTg3NjksIDAuOTcxMzgsIDAuOTMwODcsIDAuODQxMjUsIDAuNjc3MjgsIDAuNTAwMDAsIDAuMzMyMzYsIDAuMTM1ODZ9OwoKICAgIGRvdWJsZSBTMiA9IDAsIFM0ID0gMCwgUzYgPSAwLCBTOCA9IDA7CiAgICBkb3VibGUgU3kgPSAwLCBTMnkgPSAwLCBTNHkgPSAwOwoKCiAgICBmb3IgKGludCBpID0gMDsgaSA8IDg7IGkrKykgewogICAgICAgIGRvdWJsZSB4MiA9IHh2YWxbaV0gKiB4dmFsW2ldOwogICAgICAgIGRvdWJsZSB4NCA9IHgyICogeDI7CiAgICAgICAgZG91YmxlIHg2ID0geDQgKiB4MjsKICAgICAgICBkb3VibGUgeDggPSB4NCAqIHg0OwoKICAgICAgICBTMiArPSB4MjsKICAgICAgICBTNCArPSB4NDsKICAgICAgICBTNiArPSB4NjsKICAgICAgICBTOCArPSB4ODsKCiAgICAgICAgU3kgKz0geXZhbFtpXTsKICAgICAgICBTMnkgKz0geDIgKiB5dmFsW2ldOwogICAgICAgIFM0eSArPSB4NCAqIHl2YWxbaV07CiAgICB9CgogICAKICAgIGRvdWJsZSBBWzNdWzNdID0gewogICAgICAgIHs4LCAgUzIsIFM0fSwKICAgICAgICB7UzIsIFM0LCBTNn0sCiAgICAgICAge1M0LCBTNiwgUzh9CiAgICB9OwoKICAgIGRvdWJsZSBiWzNdID0ge1N5LCBTMnksIFM0eX07CiAgICBkb3VibGUgc29sWzNdOwoKICAgIGdhdXNzMyhBLCBiLCBzb2wpOwogICAgcHJpbnRmKCIoMik6IHkgPSBiMSArIGIyKnheMiArIGIzKnheNCDoqIjnrpfntZDmnpxcbiIpOwogICAgcHJpbnRmKCLmsYLjgoHjgonjgozjgZ/kv4LmlbA6XG4iKTsKICAgIHByaW50ZigiICBiMSA9ICUuNmZcbiIsIHNvbFswXSk7CiAgICBwcmludGYoIiAgYjIgPSAlLjZmXG4iLCBzb2xbMV0pOwogICAgcHJpbnRmKCIgIGIzID0gJS42ZlxuIiwgc29sWzJdKTsKCiAgICBkb3VibGUgRSA9IDA7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IDg7IGkrKykgewogICAgICAgIGRvdWJsZSB4MiA9IHh2YWxbaV0gKiB4dmFsW2ldOwogICAgICAgIGRvdWJsZSB4NCA9IHgyICogeDI7CiAgICAgICAgZG91YmxlIHloYXQgPSBzb2xbMF0gKyBzb2xbMV0gKiB4MiArIHNvbFsyXSAqIHg0OwogICAgICAgIGRvdWJsZSBkaWZmID0geXZhbFtpXSAtIHloYXQ7CiAgICAgICAgRSArPSBkaWZmICogZGlmZjsKICAgIH0KCiAgICBwcmludGYoIlxu6L+R5Ly85byP44Gr44KI44KL5YGP5beu44GuMuS5l+WSjCAo6Kqk5beuIEUpOlxuIik7CiAgICBwcmludGYoIiAgRSA9ICUuMTBmXG4iLCBFKTsKICAgIAogICAgcHJpbnRmKCJcbui/keS8vOW8j+OBq+OCiOOCiyB5IOOBruS6iOa4rOWApDpcbiIpOwogICAgZm9yIChpbnQgaSA9IDA7IGkgPCA4OyBpKyspIHsKICAgICAgICBkb3VibGUgeDIgPSB4dmFsW2ldICogeHZhbFtpXTsKICAgICAgICBkb3VibGUgeDQgPSB4MiAqIHgyOwogICAgICAgIGRvdWJsZSB5aGF0ID0gc29sWzBdICsgc29sWzFdICogeDIgKyBzb2xbMl0gKiB4NDsKICAgICAgICBwcmludGYoIiVmXG4iLCB5aGF0KTsKICAgIH0KCiAgICByZXR1cm4gMDsKfQ==