#!/usr/local/bin/tcc -run
// 数値解析課題3　5S17　西井智紀　2017May23

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

#define N 15 //データ数
#define M 1 //回帰曲線の次数
#define EPS .0001　//許容誤差

double a[M+1][M+2];

double x[N]; 
double X1[N];
double y[N];
double Y1[N];
double Answer[N];
double b;

int jordan(void);
void array(void);

int main (int argc, char **argv){
    array();
    x[N] = log(X1[N]);
    y[N] = log(Y1[N]);
    int i, j, k;

    for(i=0; i<=M; i++ ){
        for(j=0; j<=M; j++ ){
            for(k=0; k<=N; k++ ){
                a[j][i] += pow(x[k], (double)(i+j));//式(3.23)右辺
            }
        }
    }
        
    
    for(j=0; j<=M; j++ ){
            for(k=0; k<=N; k++ ){
                a[j][M+1] += y[k] * pow(x[k], j);
        }
    }

    if(jordan() == 1) return 1;
                    //解答打ち出し

    for(i=0; i<=M; i++){
        Answer[i] = a[i][M+1];
    }
    b = exp(Answer[1]);
    printf("y = %7.3lf × x^%7.3lf \n ",b, Answer[0]);
     
        return 0;
}
void array(void){
    X1[N] = {1.5, 2.0, 3.5,3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5};
    Y1[N] = {15.2,57.1,161.8,380.5,778.2,1451.2,2510.1,4092.3,6374.1,9555.6,13860.3,19563.8,26962.8,36390.8,48256.4};
}
//ガウスジョルダン法による連立方程式の計算
int jordan(void){
    double pivot, delta;
    int i, j, k;
    for(i=0; i<=M; i++){
        pivot = a[i][i];
        if(fabs(pivot) < EPS){
            printf("ピボットが許容差以下\n");
            return 1;
        }
        for(j=0; j<=M+1; j++)
            a[i][j] /=pivot;
        for(k=0; k<=M; k++){
            delta = a[k][i];
            for(j=0; j<=M+1; j++)
                if(k != i)
                    a[k][j] -= delta * a[i][j];
        }
    }
    return 0;
}