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

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 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;
                }
            }
        }
    }
}