#include <math.h>
#include <iomanip.h>
#include <iostream.h>
#include <process.h>

void ludcmp(double a[16][16],int n,int indx[16],double d)
{
	  int nmax,i,j,k,imax;
	  double tiny,aamax,sum,dum;
      nmax = 100;
      tiny = 1e-20;
      double vv[100];
      d = 1.0;
      for (i = 1;i<= n;i++)
	  {
          aamax = 0.0;
          for (j = 1;j<=n;j++)
              if (fabs(a[i][j]) > aamax ) aamax = fabs(a[i][j]);
         
          if (aamax == 0.0)  cout<< "singular matrix."<<endl;
          vv[i] = 1.0 / aamax;
	  }
      for (j = 1;j<=n;j++)
	  {
          if (j > 1) 
		  {
              for (i = 1;i<=j - 1;i++)
			  {
                  sum = a[i][j];
                  if (i > 1)
				  {
                      for (k = 1;k<=i - 1;k++)
                          sum = sum - a[i][k] * a[k][j];
                      
                      a[i][j] = sum;
                  }
              }
          }
          aamax = 0.0;
          for (i = j;i<= n;i++)
		  {
              sum = a[i][j];
              if (j > 1)
			  {
                  for (k = 1;k<=j - 1;k++)
                      sum = sum - a[i][k] * a[k][j];
                  
                  a[i][j] = sum;
              }
              dum = vv[i] * fabs(sum);
              if (dum >= aamax)
			  {
                  imax = i;
                  aamax = dum;
              }
		  }
          if (j!= imax)
		  {
              for (k = 1;k<=n;k++)
			  {
                  dum = a[imax][k];
                  a[imax][k] = a[j][k];
                  a[j][k] = dum;
              }
              d = -d;
              vv[imax] = vv[j];
          }
          indx[j] = imax;
          if (j != n) 
		  {
              if (a[j][j] == 0.0) a[j][j] = tiny;
              dum = 1.0 / a[j][j];
              for (i = j + 1;i<=n;i++)
			  {
                  a[i][j] = a[i][j] * dum;
              }
          }
      }
      if (a[n][n] == 0.0) a[n][n] = tiny;
}

void  lubksb(double a[16][16],int n,int indx[16],double b[16])
{
	int ll,ii,i,j;
	double sum;
    ii = 0;
    for( i = 1;i<=n;i++)
	{
        ll = indx[i];
        sum = b[ll];
        b[ll] = b[i];
        if (ii != 0)
		{
            for (j = ii;j<=i - 1;j++)
                sum = sum - a[i][j] * b[j];
        }
        else
		{
		if (sum != 0.0)
            ii = i;
        }
        b[i] = sum;
    }
    for (i = n ;i>=1 ;i--)
	{
        sum = b[i];
        if (i < n)
		{
            for (j = i + 1;j<=n;j++)
                sum = sum - a[i][j] * b[j];
        }
        b[i] = sum / a[i][i];
    }
}

void usrfun(double x[16],double alpha[16][16],double beta[16])
{
	int np;
    np = 15;
    alpha[1][1] = -2.0 * x[1];
    alpha[1][2] = -2.0 * x[2];
    alpha[1][3] = -2.0 * x[3];
    alpha[1][4] = 1.0;
    alpha[2][1] = 2.0 * x[1];
    alpha[2][2] = 2.0 * x[2];
    alpha[2][3] = 2.0 * x[3];
    alpha[2][4] = 2.0 * x[4];
    alpha[3][1] = 1.0;
    alpha[3][2] = -1.0;
    alpha[3][3] = 0.0;
    alpha[3][4] = 0.0;
    alpha[4][1] = 0.0;
    alpha[4][2] = 1.0;
    alpha[4][3] = -1.0;
    alpha[4][4] = 0.0;
    beta[1]= x[1]*x[1]+x[2]*x[2]+x[3]*x[3]-x[4];
    beta[2]=-x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-x[4]*x[4]+1.0;
    beta[3]=-x[1]+x[2];
    beta[4]=-x[2]+x[3];
}

void mnewt(int ntrial,double x[16],int n,double tolx,double tolf)
{
    double alpha[16][16],beta[16];
	int k,i,indx[16];
	double errf,d,errx;
    for (k = 1; k<=ntrial; k++)
	{
        usrfun(x,alpha,beta);
        errf = 0.0;
        for (i = 1; i<=n; i++)
            errf = errf + fabs(beta[i]);
        if (errf <= tolf)  _c_exit();
        ludcmp(alpha,n,indx,d);
        lubksb(alpha,n,indx,beta);
        errx = 0.0;
        for (i = 1; i<=n; i++)
		{
            errx = errx + fabs(beta[i]);
            x[i] = x[i] + beta[i];
        }
        if (errx <= tolx) _c_exit();
    }
}

void main()
{
	//program d10r13
    //driver for routine mnewt
	int i,j,kk,k,ntrial,np,n;
	double tolx,tolf,xx;
    ntrial = 5;
    tolx = 0.000001;
    n = 4;
    tolf = 0.000001;
    np = 15;
    double x[16],alpha[16][16],beta[16];
    for (kk = -1; kk<= 1; kk=kk+2)
	{
        for (k =1; k<=3; k++)
		{
            xx = 0.2 * k * kk;
            cout<< "Starting vector number "<<k<<endl;
            for (i = 1; i<=4; i++)
			{
                x[i] = xx + 0.2 * i;
                cout<<"x("<<i<<")= "<<x[i]<<endl; 
            }
            for (j = 1;j <=ntrial; j++)
			{
                mnewt(1,x,n,tolx,tolf);
                usrfun(x,alpha,beta);
                cout<<"  i        x(i)            f"<<endl;
                for (i = 1; i<=n; i++)
				{
                    cout<< setw(3)<< i;
                    cout<< setw(14)<< x[i];
                    cout<< setw(16)<< -beta[i]<<endl;
                }
            }
        }
    }
}