#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++) {
// Normalize pivot row
double pivot = A[i][i];
for (int j = 0; j < 3; j++) A[i][j] /= pivot;
b[i] /= pivot;

// Eliminate other rows
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("Model 2 coefficients:\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("Squared error E = %.10f\n", E);

printf("\nModel 2 predicted values:\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("%f\n", yhat);
}

return 0;
}
