#include <stdio.h>
#include <math.h>
/* nsq.c : 1999 3 27 by T.Oguni */
#define N 5
#define M 2
double alpha[M]={1.3,-1.};
/* function */  
jacobnsq(A,x,y,b,n,m)
double A[M][M],x[N],y[N],b[M];
int n,m;
{
  double f0[N],f1[N];
  double s1,s2,s3,b1,b2,b3,r;
  int i;
  for (i=0; i<n; i++){
    f0[i]=exp(alpha[1]*x[i]);
    f1[i]=alpha[0]*x[i]*exp(alpha[1]*x[i]);
  }
  s1=0.; s2=0.; s3=0.;
  for (i=0; i<n; i++){
    s1+=f0[i]*f0[i]; s2+=f0[i]*f1[i]; s3+=f1[i]*f1[i]; }
  A[0][0]=s1; A[0][1]=s2; A[1][0]=s2; A[1][1]=s3;
  b1=0.; b2=0.;
  for (i=0; i<n; i++){
    r=alpha[0]*exp(alpha[1]*x[i])-y[i];
    b1-=f0[i]*r;  b2-=f1[i]*r;
  }
  b[0]=b1; b[1]=b2;
  return;
}
/* function */
nsq(x,y,delta,tol,A,b,n,m,itr)
double x[N],y[N],delta[M],tol,A[M][M],b[M];
int n,m,itr[2];
{
  double s, err;
  int k, i, j;
  k=0;
  while (k<itr[0]){
    k++; itr[1]=k;
    jacobnsq(A,x,y,b,n,m);
    gauss(A,b,delta,tol,m);
    for (i=0; i<m; i++){
      alpha[i]+=delta[i]; }
    s=0.;
    for (i=0; i<m; i++){
      s+=delta[i]*delta[i]; }
    err=sqrt(s);
    if (err<tol) return;
  }
  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;
}
/* main */
main(){
  double x[N]={0.1,0.2,0.3,0.4,0.5};
  double y[N]={1.1112,0.8232,0.6099,0.4518,0.3347};
  double A[M][M], b[M], delta[M];
  double tol=0.00001;
  int n=5, m=2, i, itr[2]={10,1};
  nsq(x,y,delta,tol,A,b,n,m,itr);
  printf("%d\n", itr[1]);
  for (i=0; i<m; i++){
    printf("%f\n", alpha[i]); }
}