Least Square Method (최소 자승법)

// Least Square Method (최소 자승법)

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

int  lstsq(double x[], double y[], int n, int number, double c[]);

const int n = 1; // 구하려는 함수의 차수
const int number = 10; // 데이터의 갯수

double x[number] = { 4.0, 4.2, 4.5, 4.7, 5.1, 5.5, 5.9, 6.3, 6.8, 7.1 },
       y[number] = { 102.56, 113.18, 130.11, 142.05, 167.53, 195.14,
                     224.87, 256.73,299.50,326.72};
double p;

void main()
{
   double c[number], error=0.0;
   int chkerr, i;

   printf(" < Least Squre Method > \n\n");

   chkerr = lstsq(x, y, n, number, c);

   // 결과 출력
   for(i=0; i<=n; i++) printf("a%d = %10.8f\n",i, c[i]);
   if(chkerr==0) printf("\nNo Error");
   else printf("\nError");
}

// 0 출력-> 에러없음, 999 출력-> 계산에러
int  lstsq(double x[], double y[], int n, int number, double c[])
{
   int  i, j, k, l;
   double  w1, w2, w3, pivot, aik, a[21][22], w[42];

   if(n >= number || n < 1 || n > 20) return(999);
   for(i = 0; i < n*2; i++) {
      w1 = 0.0;
      for(j = 0; j < number; j++) {
         w2 = w3 = x[j];
	 for(k = 0; k < i; k++)	w2 *= w3;
  	 w1 += w2;
      }
      w[i] = w1;
   }

   // matrix 입력
   for(i = 0; i < n+1; i++) {
      for(j = 0; j < n+1; j++) {
         l = i + j - 1;
         a[i][j] = w[l];
      }
   }
   a[0][0] = number;
   w1 = 0.0;
   for(i = 0; i < number; i++) w1 += y[i];
   a[0][n+1] = w1;

   // sigma(Yi Xi) 계산해서 대입
   for(i = 0; i < n; i++) {
      w1 = 0.0;
      for(j = 0; j < number; j++) {
         w2 = w3 = x[j];
	 for(k = 0; k < i; k++) w2 *= w3;
  	 w1 += y[j] * w2;
      }
      a[i+1][n+1] = w1;
   }

   // matrix 계산
   for(k = 0; k < n+1; k++) {
      pivot = a[k][k];
      for(j = k; j < n+2; j++) a[k][j] /= pivot;
      for(i = 0; i < n+1; i++) {
         if(i != k) {
            aik = a[i][k];
	    for(j = k; j < n+2; j++) a[i][j] -= aik * a[k][j];
         }
      }
   }

   // 다항식의 계수를 최종적으로 전달 (pass by reference)
   for(i = 0; i < n+1; i++) c[i] = a[i][n+1];
   return(0);
}