#include <stdio.h>
#include <math.h>
/* lsq.c : 1999 3 13 by Oguni */
#define N 12
#define M 3
main(){
  double X[N][M]={.98,.77,1.,.76,1.39,1.,.58,1.79,1.,.94,.95,1.,.73,1.42,1.,.57,1.86,1.,.97,.82,1.,.91,1.00,1.,.79,1.25,1.,.92,.96,1.,.92,.99,1.,.54,1.90,1.};
  double y[N]={1.52,1.41,1.16,1.45,1.24,1.21,1.63,1.38,1.37,1.36,1.40,1.03};
  double alpha[M], ybar[N];
  int n=12, m=3, i;
  lsq(X,y,alpha,n,m);
  for (i=0; i<m; i++){
    printf("%f", alpha[i]); }
  printf("\n");
  printf("original and estimated\n");
  mvsub2(X,alpha,ybar,n,m);
  for (i=0; i<n; i++){
    printf("%f %f\n", y[i],ybar[i]); }
}
/* function */
lsq(X,y,alpha,n,m)
double X[N][M], y[N], alpha[M];
int n,m;
{
  double A[M][M],b[M],Q[N][N],Xt[M][N];
  double tol=0.00000001, s;
  int i,j,k;
  for (i=0; i<n; i++){
    for (j=0; j<m; j++){
      Xt[j][i]=X[i][j]; }
  }
  matmul1(Xt,X,A,n,m);
  mvsub1(Xt,y,b,n,m);
  printf("Matrix A\n");
  for (i=0; i<m; i++){
    printf("%f %f %f %f\n", A[i][0],A[i][1],A[i][2],b[i]); }
  gauss(A,b,alpha,tol,m);
  return;
}
/* function */
matmul1(A,B,C,n,m)
double A[M][N], B[N][M], C[M][M];
int n,m;
{
  int i, j, k;
  double s;
  for (k=0; k<m; k++){
    for (i=0; i<m; i++){
      s=0.0;
      for (j=0; j<n; j++){ s+=A[i][j]*B[j][k]; }
      C[i][k]=s;
    }
  }
  return;
}
/* function */
mvsub1(A,b,c,n,m)
double A[M][N], b[N], c[M];
int n,m;
{
  int i, j;
  double s;
  for (i=0; i<m; i++){
    s=0.0;
    for (j=0; j<n; j++){ s+=A[i][j]*b[j]; }
    c[i]=s;
  }
  return;
}
/* function */
mvsub2(A,b,c,n,m)
double A[N][M], b[M], c[N];
int n,m;
{
  int i,j;
  double s;
  for (i=0; i<n; i++){
    s=0.0;
    for (j=0; j<m; j++){ s+=A[i][j]*b[j]; }
    c[i]=s;
  }
  return;
}
/* function */
gauss(A, b, x, tol, n)
double A[M][M], b[M], x[M], tol;
int n;
{
  int i, j, k;
  double akk, aik, s;
  for (k=0; k<n-1; k++){
    if (fabs(A[k][k])<tol){
      printf("The pivot is too small at step # %d", k); return; }
    akk=1.0/A[k][k];
    for (i=k+1; i<n; i++){
      aik=-A[i][k]*akk;
      for (j=k+1; j<n; j++){
        A[i][j]+=aik*A[k][j]; }
      b[i]+=aik*b[k];
    }
    for (j=k+1; j<n; j++){
      A[k][j]*=akk; }
    b[k]*=akk;
  }
  x[n-1]=b[n-1]/A[n-1][n-1];
  for (k=n-2; k>=0; k--){
    s=0.0;
    for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; }
    x[k]=b[k]-s;
  }
  return;
}