// lsfit_log.c
//               compile   gcc -o lsfit_log.c simeq.c -lm
//               run       lsfit_log > lsfit_log.out
//
//    | SUM(S*S)  SUM(S*T)  SUM(S*U)  SUM(S*V) |   | a |   | SUM(S*Y) |
//    | SUM(T*S)  SUM(T*T)  SUM(T*U)  SUM(T*V) | x | b | = | SUM(T*Y) |
//    | SUM(U*S)  SUM(U*T)  SUM(U*U)  SUM(U*V) |   | c |   | SUM(U*Y) |
//    | SUM(V*S)  SUM(V*T)  SUM(V*U)  SUM(V*V) |   | d |   | SUM(V*Y) |

#include "simeq.h" // needed to solve simultaneous equations
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define n 21 // data set observations
#define nt 8 // number of terms
             // 1, log(x), log(x^2), sin(x), cos(x), x, x^2, x^3 
#undef  abs
#define abs(x) (((x)<0.0)?(-(x)):(x))
#undef  max
#define max(x,y) (((x)>(y))?((x)):(y))
#undef  min
#define min(x,y) (((x)<(y))?((x)):(y))

int main(int argc, char* argv[])
{
  double xd[n];
  double yd[n];
  double A[nt][nt];
  double X[nt];
  double Y[nt];     // |A| * |X| = |Y|  |X| will be  a, b, c
  double trmi[nt];
  double trmj[nt];
  int i, j, k;
  double Yc;

  printf("lsfit_log.c run to make lsfit_log.out \n");
  // create data
  for(i=0; i<n; i++)
  {
    if(i<10)
    {
      xd[i]=1.0/(double)(11-i);
    }
    else
    {
      xd[i] = (double)(i-9);
    }
    yd[i] = log(xd[i]);
    printf("i=%2d, xd=%f, yd=%f\n", i, xd[i], yd[i]);
  }
  for(i=0; i<nt; i++)
  {
    Y[i] = 0.0; // must start sum with zero
    for(k=0; k<n; k++)
    {
      trmi[0] = 1.0;
      trmi[1] = log(xd[k]);
      trmi[2] = log(xd[k]*xd[k]);
      trmi[3] = sin(xd[k]);
      trmi[4] = cos(xd[k]);
      trmi[5] = xd[k];
      trmi[6] = xd[k]*xd[k];
      trmi[7] = xd[k]*xd[k]*xd[k];
      Y[i] += trmi[i]*yd[k];
    }
    for(j=0; j<nt; j++)
    {
      A[i][j] = 0.0;  // must start sum with zero
      for(k=0; k<n; k++) // run through observations
      {
        trmj[0] = 1.0;
        trmj[1] = log(xd[k]);
        trmj[2] = log(xd[k]*xd[k]);
        trmj[3] = sin(xd[k]);
        trmj[4] = cos(xd[k]);
        trmj[5] = xd[k];
        trmj[6] = xd[k]*xd[k];
        trmj[7] = xd[k]*xd[k]*xd[k];
        A[i][j] += trmj[i]*trmj[j];  // Xj * Xi
      }
    }
  }  // end i
  simeq2(nt, A, Y, X); // pick correct argument list
  printf("1=%f, log(x)=%f, log(x^2)=%f, sin(x)=%f \n",
         X[0], X[1], X[2], X[3]);
  printf("cos(x)=%f, x=%f, x^2=%f, x^3=%f \n",
         X[4], X[5], X[6], X[7]);
  // check
  for(k=0; k<n; k++)  // k is observation
  {
    Yc = X[0]+X[1]*log(xd[k])+X[2]*log(xd[k]*xd[k]) +
         X[3]*sin(xd[k]) + X[4]*cos(xd[k]) + X[5]*xd[k] +
	 X[6]*xd[k]*xd[k] + X[7]*xd[k]*xd[k]*xd[k];
    printf("Y_actual=%f, Y_approximate=%f, err=%e \n",
	   yd[k], Yc, yd[k]-Yc);
  }
} // end lsfit_log.c