GSL sample: Solving Algebraic Equations

Last Update: Oct. 4, 2007

Sample code

algebraic_eq.c

/***********************************************/
/* 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;
}

Example execution

$ ./algebraic_eq
Solve Algebraic Equations
re(x0), im(x0) = 0.621203, 0.536707 ( p(x0) = -7.8e-15 + -1.0e-15 * i )
re(x1), im(x1) = 0.621203, -0.536707 ( p(x1) = -7.8e-15 + 1.0e-15 * i )
re(x2), im(x2) = 0.189254, 0.757699 ( p(x2) = 3.6e-15 + -8.9e-17 * i )
re(x3), im(x3) = 0.189254, -0.757699 ( p(x3) = 3.6e-15 + 8.9e-17 * i )
re(x4), im(x4) = 3.21965e-15, 0 ( p(x4) = 3.2e-15 + 0.0e+00 * i )
re(x5), im(x5) = -0.747054, 0 ( p(x5) = 3.5e-15 + -0.0e+00 * i )
re(x6), im(x6) = -0.617518, 0.426103 ( p(x6) = 2.1e-15 + -5.8e-15 * i )
re(x7), im(x7) = -0.617518, -0.426103 ( p(x7) = 2.1e-15 + 5.8e-15 * i )
re(x8), im(x8) = -0.269412, 0.711294 ( p(x8) = 4.1e-15 + -7.5e-17 * i )
re(x9), im(x9) = -0.269412, -0.711294 ( p(x9) = 4.1e-15 + 7.5e-17 * i )

<- Go back