/***********************************************/ /* T.Kouya's GSL sample program collection */ /* Solving Algebraic Equations */ /* Written by Tomonori Kouya */ /* */ /* Version 0.01: 2007-10-01 */ /***********************************************/ #include #include #include #include /* 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; }