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

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

	・マラー法を用いて3次代数方程式の実根を求めます.

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

#define delta pow(10,-8)	//δ(textに対応)
#define ix0  0			//初期の近似値x0
#define ix1 -1			//初期の近似値x1
#define ix2 -2			//初期の近似値x2
#define kmax 1000		//収束判定条件との比較回数の上限

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

int main(void){
	int k;		
	double x0, x1, x2, x3, y0, y1, y2, y3, h1, g, D;
	double denomi1, denomi2, denomi, lambda1;
	x0 = ix0;
	x1 = ix1;
	x2 = ix2;
	y0 = f(x0);
	y1 = f(x1);
	y2 = f(x2);
	h1 = x2-x1;			//step1
	lambda1 = h1/(x1-x0);		//step2

	for(k=3; k<kmax+1; k++){

		g = (1+2*lambda1)*(y2-y1)-pow(lambda1,2)*(y1-y0);
					//step3

		D = pow(g,2)-4*y2*(1+lambda1)*lambda1
			*(y2-y1-lambda1*(y1-y0));

		if(D <= 0){		//複素数のとき

			denomi = (g*g-D)/g;

		}else{

			denomi1 = g+sqrt(D);
					//g+rootD

			denomi2 = g-sqrt(D);
					//g-rootD

			if(fabs(denomi1)>fabs(denomi2)){
					denomi = denomi1;
			}else{
					denomi = denomi2;
			}

		}

		lambda1 = -2*y2*(1+lambda1)/denomi;
					//step4

		x3 = x2+h1*lambda1;	//step5
		y3 = f(x3);

		printf("k=%2d x=%12.9lf f(x)=%15.7e\n", k, x3, y3);
		if(fabs(y3) < delta){
					return 0;
		}else{

					//step5
			x0 = x1;
			x1 = x2;
			x2 = x3;
			y0 = f(x0);
			y1 = f(x1);
			y2 = f(x2);
			h1 = h1*lambda1;
		}
	}
	return 1;
}