// russian : cp1251 encoding
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "ext_value.h"

void g(double *P, double *Q);
void dg(double *P, double *dQdP);

void g(double *P, double *Q)
{
	double x,y,z,lambda1,lambda2;
	 x=P[0];
	 y=P[1];
	 z=P[2];
	 lambda1=P[3];
	 lambda2=P[4];
	Q[0]=y*z+2*lambda1*x+lambda2;
	Q[1]=x*z+2*lambda1*y+lambda2;
	Q[2]=x*y+2*lambda1*z+lambda2;
	Q[3]=x*x+y*y+z*z-1;
	Q[4]=x+y+z;
}

void dg(double *P, double *dQdP)
{
	double x,y,z,lambda1;
//	double lambda2;
// значение lambda2 в вычислении вторых производных не требуется
	 x=P[0];
	 y=P[1];
	 z=P[2];
	 lambda1=P[3];
//	 lambda2=P[4];
   dQdP[0]=2*lambda1; // dQ[0]/dx
   dQdP[1]=z; // dQ[0]/dy
   dQdP[2]=y; // dQ[0]/dz
   dQdP[3]=2*x; // dQ[0]/d1
   dQdP[4]=1.0; // dQ[0]/d2   
   
   dQdP[5]=z; // dQ[1]/dx
   dQdP[6]=2*lambda1; // dQ[1]/dy
   dQdP[7]=x; // dQ[1]/dz
   dQdP[8]=2*y; // dQ[1]/d1
   dQdP[9]=1.0; // dQ[1]/d2   
   
   dQdP[10]=y; // dQ[2]/dx
   dQdP[11]=x; // dQ[2]/dy
   dQdP[12]=2*lambda1; // dQ[2]/dz
   dQdP[13]=2*z; // dQ[2]/d1
   dQdP[14]=1.0; // dQ[2]/d2   

   dQdP[15]=2*x; // dQ[3]/dx
   dQdP[16]=2*y; // dQ[3]/dy
   dQdP[17]=2*z; // dQ[3]/dz
   dQdP[18]=0.0; // dQ[3]/d1
   dQdP[19]=0.0; // dQ[3]/d2   
   
   dQdP[20]=1.0; // dQ[4]/dx
   dQdP[21]=1.0; // dQ[4]/dy
   dQdP[22]=1.0; // dQ[4]/dz
   dQdP[23]=0.0; // dQ[4]/d1
   dQdP[24]=0.0; // dQ[4]/d2   
}

int main()
{ 
	double P[5],dQdP[5*5],Q[5];
	int i,j;
// Это массив переменных	
	P[0]=0.409248; // x
	P[1]=0.406248; // y
	P[2]=-0.812128; // z
	P[3]=0.205724; // \lambda_1
	P[4]=0.150102; // \lambda_2
// Вычисление значений вектор-функции:
	g(P,Q);
// Явное вычисление матрицы производных:
	dg(P,dQdP);

// Использование проекта ext_value для вычисления значений вектор-функции и производных
// состоит в описании переменных (вызываются "правильные" конструкторы)
	ext_value2 xe(5,P[0],0), ye(5,P[1],1), ze(5,P[2],2), le1(5,P[3],3), le2(5,P[4],4), L(5);
// и вычислении функции Лагранжа:
	L=xe*ye*ze+le1*(xe*xe+ye*ye+ze*ze-1)+le2*(xe+ye+ze);
	
// Сравнение полученных результатов:
	printf("Значения вектор-функции:\n");
	for(i=0;i<5;i++) printf("%le\t%le\t%le\n",Q[i],L.get(i),Q[i]-L.get(i));
	printf("\n");
	printf("Матрица Якоби:\n");
	for(i=0;i<5;i++) {
	for(j=0;j<5;j++) printf("%le\t%le\t%le\n",dQdP[i*5+j],L.get(i,j),dQdP[i*5+j]-L.get(i,j));
					 printf("\n");
	}

	return 0;
}