/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*                 Solving Algebraic Equations */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-01                    */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_poly.h>

/* not supported but very needed functions! */
gsl_complex _gsl_poly_complex_eval(double coef[], const int len, const gsl_complex x)
{
	gsl_complex ret;
	int i;

	/* Horner's Method */
	GSL_SET_COMPLEX(&ret, coef[len - 1], 0.0);

	/* ret = ret * x + coef[i] */
	for(i = len - 2; i >= 0; i--)
		ret = gsl_complex_add_real(gsl_complex_mul(ret, x), coef[i]);

	return ret;
}

#define DEGREE 10
#define DEGREEP1 11 // degree + 1

int main(void)
{
	int i;
	double coef[DEGREEP1];
	gsl_complex csol[DEGREE], eval;
	gsl_poly_complex_workspace *work;

/* Solve 10th degree Equations */
	printf("Solve Algebraic Equations\n");

	/* Allocate workspace */
	work = gsl_poly_complex_workspace_alloc(DEGREEP1);

	/* 10 * x^10 + 9 * x^9 + ... + x + 0 = 0 */
	for(i = 0; i < DEGREEP1; i++)
		coef[i] = (double)i;

	/* Solve equation */
	gsl_poly_complex_solve(coef, DEGREEP1, work, csol);

	/* print solutions */
	for(i = 0; i < DEGREE; i++)
	{
		eval = _gsl_poly_complex_eval(coef, DEGREEP1, csol[i]);
		printf("re(x%d), im(x%d) = %g, %g ( p(x%d) = %7.1e + %7.1e * i )\n", i, i, GSL_REAL(csol[i]), GSL_IMAG(csol[i]), i, GSL_REAL(eval), GSL_IMAG(eval));
	}

	/* Free workspace */
	gsl_poly_complex_workspace_free(work);

	return 0;
}