#include <iostream.h>
#include <stdlib.h>
#include <math.h>
/////////////////////////////////////////////////////////////
// Define the inner product function.
/////////////////////////////////////////////////////////////

float inner(float a[], float b[], int dim)
{
  int n;
  float prod=0;
  for(n=0;n<dim;n++)
  {
    prod+=a[n]*b[n];
  }
   return prod;
}

/////////////////////////////////////////////////////////////
// Define the tabphi vectors.
/////////////////////////////////////////////////////////////

float * tabphi(int n, int M, float x[])
{
 float * array = new float[M];
 int k;
 for(k=0;k<M;k++)
 { if(n==0){array[k]=1;}
   if(n==1){array[k]=sin(2*3.14*x[k]/12);}
   if(n==2){array[k]=cos(2*3.14*x[k]/12);}
 }
 return array;
}


//////////////////////////////////////////////////////
// Function which finds the inverse of a matrix:
//////////////////////////////////////////////////////

// This function takes the flattened version of a square matrix A as an
// argument, i.e., an NxN matrix written as an array of N^2 elements. It also
// takes the dimension N as an argument.

float * FlatInv(int N, float FlatA[])
{
  float A[N][N];
  int k; int l;
  for(k=0;k<N;k++)
  {
    for(l=0; l<N;l++)
    {
      A[k][l]=FlatA[k*N+l];
    }
  }


//////////////////////////////////////////////////////
// Form the augmented matrix
//////////////////////////////////////////////////////

   float Aug[N][2*N];

   for(k=0;k<N;k++)
   {
     for(l=0;l<N;l++)
     {
       Aug[k][l]=A[k][l];
       if(k==l){Aug[k][N+l]=1;}
       else{Aug[k][N+l]=0;}
     }
   }

//////////////////////////////////////////////////////
// Now start the Gaussian algorithm:
//////////////////////////////////////////////////////

//a) First row echelon form:

   int d;

   float M[2*N];

   for(d=0;d<N;d++)               // Start at the d th row.
   {                               //
     while(Aug[d][d]==0)           // If A[d][d]==0 then put that row
     {                             // at the bottom. This way one eventually
       for(l=0;l<2*N;l++)         // gets a nonzero pivot element if the d th
       {                           // column is not zero altogether.
         M[l]=Aug[d][l];
       }
        for(k=d;k<N;k++)
        {
          for(l=0;l<2*N;l++)
          {
            Aug[k][l]=Aug[k+1][l];
          }
        }
        for(l=0;l<2*N;l++)
        {
          Aug[N-1][l]=M[l];
        }
     }                              // End of resorting.

     for(l=2*N-1; l>=d; l--)             // Divide the d th row by the pivot element.
     {
        Aug[d][l]=Aug[d][l]/Aug[d][d];
     }

     for(k=d+1;k<N;k++)               // Use the pivot at the d th
     {                                 // row to cancel the terms
        for(l=2*N-1; l>=d; l--)          // under it.
        {
         Aug[k][l]=Aug[k][l]-Aug[k][d]*Aug[d][l];
        }
     }                                 // End of business with d th row.
                                       // Now go back and do it for d+1 st row.
   }


// b) Now the reduced row echelon form:

 for(d=N-1; d>=0; d--)
 {
   for(k=d-1;k>=0;k--)           // Use the pivot at row d to cancel the terms
   {                             // above it.
     for(l=2*N-1; l>=0; l--)
     {
       Aug[k][l]=Aug[k][l]-Aug[k][d]*Aug[d][l];
     }
   }
 }

//////////////////////////////////////////////////////
// Register the inverse of the matrix:
//////////////////////////////////////////////////////

   float Inv[N][N];

   for(k=0;k<N;k++)
   {
     for(l=0;l<N;l++)
     {
       Inv[k][l]=Aug[k][N+l];
     }
   }

//////////////////////////////////////////////////////
// Form the flat inverse matrix:
//////////////////////////////////////////////////////


   float * FlatInv=new float[N*N];

   for(k=0;k<N;k++)
   {
     for(l=0;l<N;l++)
     {
       FlatInv[k*N+l]=Inv[k][l];
     }
   }

return FlatInv;

}





/////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// The main part:
/////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////

int main()
{
    cout <<
    "This program finds the periodic function which fits best to a given set of data points."
    << endl << endl;
    ///////////////////////////////////////////////////////
    // Get the data.
    ///////////////////////////////////////////////////////
    cout << "Please enter the number of your data points: ";
    int M;
    cin >> M;
    float x[M];
    float tabf[M];

    cout << endl << "Please enter your data points:" << endl;
    int m;
    for(m=0; m<M; m++)
    {
      cout << "x["<< m << "]=";
      cin >> x[m];
      cout << "   f(x["<< m << "])=";
      cin >> tabf[m];
    }
    ///////////////////////////////////////////////////////
    // Ask about the polynomial.
    ///////////////////////////////////////////////////////

//   cout << "What is the degree of the polynomial that you would like to fit to this data? " ;
//    int N;
//    cin >> N;

     int N=2;

    ///////////////////////////////////////////////////////
    // Write down the tabphi and tabf vectors.
    ///////////////////////////////////////////////////////

    int n; int k;
    for(n=0;n<=N;n++)
    {
      cout << "tabphi("<< n << ")= ";
      for(k=0;k<M;k++)
      {
         cout << tabphi(n,M,x)[k]<<"   ";
      }
      cout << endl<< endl;
    }

    cout << "tabf= ";
    for(m=0;m<M;m++)
    {
       cout << tabf[m]<<"   ";
    }
    cout << endl<< endl;

    ///////////////////////////////////////////////////////
    // Define and write down the A Matrix and the B vector.
    ///////////////////////////////////////////////////////


    cout<< "The A matrix is" << endl <<endl;
    float A[N+1][N+1];

    for(n=0;n<N+1;n++)
    {
      for(k=0;k<N+1;k++)
      {
      A[n][k]=inner(tabphi(n,M,x),tabphi(k,M,x),M);
      cout << A[n][k]<< "   ";
      }
      cout<<endl<<endl;
    }

    cout<< "The B vector is" << endl<< endl;
    float B[M];
    cout << "B= ";
    for(n=0;n<=N;n++)
    {
       B[n]=inner(tabf,tabphi(n,M,x),M);
       cout << B[n]<<"   ";
    }

    cout << endl<< endl;

    ///////////////////////////////////////////////////////
    // Find the inverse of the matrix A:
    ///////////////////////////////////////////////////////


    float FlatA[(N+1)*(N+1)];       // Flatten the matrix
    int l;
    for(k=0;k<N+1;k++)
    {
      for(l=0;l<N+1;l++)
      {
        FlatA[k*(N+1)+l]=A[k][l];
      }
    }


   float * FInv;                  // Find the flattened inverse
   FInv=FlatInv(N+1,FlatA);         // and obtain the inverse matrix:

   float Inv[N+1][N+1];

  for(k=0;k<N+1;k++)
  {
    for(l=0; l<N+1;l++)
    {
      Inv[k][l]=FInv[k*(N+1)+l];
    }

  }

    ///////////////////////////////////////////////////////
    // Find c coefficients:
    ///////////////////////////////////////////////////////

    float C[N];

    for(n=0;n<N+1;n++)
    {
      C[n]=0;
      for(k=0;k<N+1;k++)
      {
      C[n]+=Inv[n][k]*B[k];
      }
      cout << "C["<< n << "]="<< C[n] << endl;
    }

    ///////////////////////////////////////////////////////
    // Write down the function:
    ///////////////////////////////////////////////////////

    cout << endl<<endl<<"The function which best fits the data is "
    << endl<<endl;

    cout << "("<< C[0]<<")+"<<"("<< C[1]<<")*sin(2*Pi*x/12)+"
    <<"("<< C[2]<<")*cos(2*Pi*x/12)+";

    ///////////////////////////////////////////////////////
    // End the main program.
    ///////////////////////////////////////////////////////

    cout << endl<< endl;
    system("PAUSE");
    return 0;
}