#include<stdio.h>
#include<math.h>

/*************** 表4.3の数値解を求めるプログラム ****************

	・(4.27)，(4.28)の式を用いて
		3次代数方程式の虚根を求めます．

*****************************************************************/

#define eps pow(10,-10)
#define x0 -1				//xの初期値
#define y0 -2				//yの初期値
#define kmax 1000			//収束判定条件との比較回数の上限

double f(double x, double y){		//f(x,y)
	double z;
	z = pow(x,3)-3*x*pow(y,2)+6*pow(x,2)-6*pow(y,2)+21*x+32;
	return z;
}

double dfx(double x, double y){		//∂f(x,y)/∂x
	double z;
	z = 3*pow(x,2)-3*pow(y,2)+12*x+21;
	return z;
}

double dfy(double x, double y){		//∂f(x,y)/∂y
	double z;
	z = -6*x*y-12*y;
	return z;
}

double g(double x, double y){		//g(x,y)
	double z;
	z = 3*pow(x,2)*y-pow(y,3)+12*x*y+21*y;
	return z;
}

double dgx(double x, double y){		//∂g(x,y)/∂x
	double z;
	z = 6*x*y+12*y;
	return z;
}

double dgy(double x, double y){		//∂g(x,y)/∂y
	double z;
	z = 3*pow(x,2)-3*pow(y,2)+12*x+21;
	return z;
}

double J(double x, double y){		//J(x,y)
	double j;
	j = dfx(x,y)*dgy(x,y)-dfy(x,y)*dgx(x,y);
	return j;
}

int main(void){
	int k=0;
	double x1, x2, y1, y2;
	x1 = x0;
	y1 = y0;
	while(pow(f(x1,y1),2)+pow(g(x1,y1),2) >= pow(eps,2)){

		//式(4.27)
		x2 = x1-(f(x1,y1)*dgy(x1,y1)-g(x1,y1)*dfy(x1,y1))/J(x1,y1);
		y2 = y1-(g(x1,y1)*dfx(x1,y1)-f(x1,y1)*dgx(x1,y1))/J(x1,y1);
		x1 = x2;
		y1 = y2;
		k++;
		printf("k=%2d, xk=%13.7lf, yk=%13.7lf,							f(xk,yk)=%16.9e, g(xk,yk)=%16.9e\n",						k, x1, y1, f(x1,y1), g(x1,y1));
		if(k > kmax){
			printf("error!!\n");
			return 1;
		}
	}
	return 0;
}