//#include"stdafx.h"
#include<iostream>
#include<algorithm>
#include<math.h>
#define eps 0.000001
#define INF 999999

using namespace std;

double min(double a, double b)
{
	return (a < b)? a:b;
}

double max(double a, double b)

{
	return (a > b)? a :b;
}

double f(double x)
{
	return (sin(3*x) * exp(2*cos(x)) - 5*x);
}

double f_prime(double x)
{
	return( exp(2*cos(x)) * (3*cos(3*x)+cos(4*x)-cos(2*x)) - 5 );
}

double f_2prime(double x)
{
	return( exp(2*cos(x)) * ( (-2*sin(x))*(3*cos(3*x)+cos(4*x)-cos(2*x)) + (2*sin(2*x)-4*sin(4*x)-9*sin(3*x)) ) );
}

double l = -7, r = 7, n = 103, fx1, fa, f2a, fx0, f1x0, x1, x0, m1, M1;
double f1p, f2p, h = (r-l)/n;
int cnt;

void newton(double a, double b)
{
	m1 = INF;
	M1 = 0;
	for(double i = 0; i <= n; i++)
	{
		f1p = fabs(f_prime(a+i*h));
		f2p = fabs(f_2prime(a+i*h));
		m1 = min(m1, f1p);
		M1 = max(M1, f2p);
	}

	fa = f(a);
	f2a = f_2prime(a);
	if(fa*f2a > 0)
		x0 = a;
	else
		x0 = b;

	fx0 = f(x0);
	f1x0 = f_prime(x0);
	x1 = x0 - fx0/f1x0;
	fx1 = f(x1);
	cnt++;

	while((M1*fabs(x1-x0)*fabs(x1-x0)/2/m1 >= eps) || (fabs(fx1) >= eps))
	{
		cnt++;
		x0 = x1;
		fx0 = fx1;
		f1x0 = f_prime(x0);
		x1 = x0 - fx0/f1x0;
		fx1 = f(x1);
	}
}

int main(int argc, char* argv[])
{
	printf("Hello World!\n");

	for(double i = 1; i <= n; i++)
	{
		if(f(l+(i-1)*h)*f(l+i*h) < 0)
		{
			cnt = 0;
			newton(l+(i-1)*h, l+i*h);
			printf("cnt=%d x=%.10f fx=%.10f \n\n",cnt,x1,fx1);
		}
	}

	return 0;
}