#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_poly.h>

int main(void){
	double landa1=0;
	double landa2=0;
	double landa3=0;
	int i,numero;
	
	double a[4] = {-1.331,3.63,-3.3,1.0};
	double z[6];
	gsl_complex ax,ay,az;
	
	gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(4);
	
	/* (x-1)^3 */
	numero = gsl_poly_solve_cubic ( -3,3, -1, &landa1, &landa2, &landa3);
	printf("--------------------\n");
	printf("Numero = \%d\t(x-1.0)^3\tgsl_poly_solve_cubic\n",numero);
	printf("Root[1] = %5.10lf\n",landa1);
	printf("Root[2] = %5.10lf\n",landa2);
	printf("Root[3] = %5.10lf\n",landa3);
	
	/* (x-1.1)^3 */
	landa1 = 0;
	landa2 = 0;
	landa3 = 0;
	numero = gsl_poly_solve_cubic( -3.3000,3.63000,-1.331000, &landa1, &landa2,&landa3);
	printf("--------------------\n");
	printf("Numero = \%d\t(x-1.1)^3\tgsl_poly_solve_cubic\n",numero);
	printf("Root[1] = %5.10lf\n",landa1);
	printf("Root[2] = %5.10lf\n",landa2);
	printf("Root[3] = %5.10lf\n",landa3);
	
	/* (x-1.1)^3 */
	numero = gsl_poly_complex_solve_cubic (-3.3000,3.63000,-1.331000, &ax, &ay, &az);
	printf("--------------------\n");
	printf("Numero = \%d\t(x-1.1)^3\tgsl_poly_complex_solve_cubic\n",numero);
	printf("Root[1-Real] = %+.10lf\tRoot[1-Im] = %+.10lf\n",ax.dat[0],ax.dat[1]);
	printf("Root[2-Real] = %+.10lf\tRoot[2-Im] = %+.10lf\n",ay.dat[0],ay.dat[1]);
	printf("Root[3-Real] = %+.10lf\tRoot[3-Im] = %+.10lf\n",az.dat[0],az.dat[1]);
	
	/* coefficients of P(x) =  (x-1.1)^3  */
	gsl_poly_complex_solve (a, 4, w, z);
	gsl_poly_complex_workspace_free (w);
	printf("--------------------\n");
	printf("Numero = 3\t(x-1.1)^3\tgsl_poly_complex_solve\n",numero);
	for (i = 0; i < 3; i++)
		printf ("Root[%d-Real] = %+.10lf\tRoot[%d-Im] = %+.10lf\n",i+1,z[2*i],i+1,z[2*i+1]);
	return 0;
}