/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*                  Solving Nonlinear Equation */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-04                    */
/***********************************************/
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

/* Limit of Iterative Times */
#define MAXTIMES 100

/* func(x) = 0 */
/* x^2 - 2 = 0 */
double func(double x, void *param)
{
	return x * x - 2;
}

/* f'(x) */
double dfunc(double x, void *param)
{
	return 2.0 * x;
}

/* simultaneous eval of f and df */
void fdfunc(double x, void *param, double *ret_f, double *ret_df)
{
	*ret_f  = func(x, param);
	*ret_df = dfunc(x, param);

	return;
}

int main(void)
{
	int i, times, status;
	gsl_function f;
	gsl_root_fsolver *workspace_f;
	gsl_function_fdf fdf;
	gsl_root_fdfsolver *workspace;
	double x, x_l, x_r;

/* Fsolver */

	/* Define Solver */
	workspace_f = gsl_root_fsolver_alloc(gsl_root_fsolver_bisection);

	printf("F solver: %s\n", gsl_root_fsolver_name(workspace_f));

	f.function = &func;
	f.params = 0;

	/* set initial interval */
	x_l = 0.0;
	x_r = 2 * M_PI;

	/* set solver */
	gsl_root_fsolver_set(workspace_f, &f, x_l, x_r);

	/* main loop */
	for(times = 0; times < MAXTIMES; times++)
	{
		status = gsl_root_fsolver_iterate(workspace_f);

		x_l = gsl_root_fsolver_x_lower(workspace_f);
		x_r = gsl_root_fsolver_x_upper(workspace_f);
		printf("%d times: [%10.3e, %10.3e]\n", times, x_l, x_r);

		status = gsl_root_test_interval(x_l, x_r, 1.0e-13, 1.0e-20);
		if(status != GSL_CONTINUE)
		{
			printf("Status: %s\n", gsl_strerror(status));
			printf("\n Root = [%25.17e, %25.17e]\n\n", x_l, x_r);
			break;
		}
	}

	/* free */
	gsl_root_fsolver_free(workspace_f);

/* FDFsolver */

	/* Define Solver */
	workspace = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);

	printf("FDF solver: %s\n", gsl_root_fdfsolver_name(workspace));

	fdf.f = &func;
	fdf.df = &dfunc;
	fdf.fdf = &fdfunc;
	fdf.params = 0;

	/* set initial value */
	x = 2.0;

	/* set solver */
	gsl_root_fdfsolver_set(workspace, &fdf, x);

	/* main loop */
	for(times = 0; times < MAXTIMES; times++)
	{
		status = gsl_root_fdfsolver_iterate(workspace);
		x = gsl_root_fdfsolver_root(workspace);
		status = gsl_root_test_residual(func(x, NULL), 1.0e-5);

		printf("%d times: %10.3e\n", times, x);

		if((status != GSL_CONTINUE) || (status == GSL_EZERODIV) || (status == GSL_EZERODIV))
		{
			printf("Status: %s\n", gsl_strerror(status));
			printf("\n Root = %25.17e\n\n", x);
			break;
		}
	}

	/* free */
	gsl_root_fdfsolver_free(workspace);

	return 0;
}