/********************************************************************
* Program 2.1: Bairstow Method for the roots of polynomials.		*
*				Coded by Dilip Datta								*
*																	*
*********************************************************************/

# include <stdio.h>
# include <math.h>
# include <stdlib.h>
# include <malloc.h>

int rcnt = 0;				/* Counters for roots & points of infinity */
double *rx, *ix;				/* Pointers to store roots (both real & imaginary parts) */
/********************************************************************/

int main()
{
	int i, n, *nit;
	double *a, r, s, eps;
	double *aa;

	void bairstow(int n, int *nit, double *a, double r, double s, double eps);

	printf("\nDegree of the polynomial                     : ");				/* Step-1 */
	scanf("%d", &n);

	nit = (int *)malloc(n*sizeof(int));
	a = (double *)malloc( (n+1+1)*sizeof(double) );
	aa = (double *)malloc( (n+1+1)*sizeof(double) );
	rx = (double *)malloc( (n+1)*sizeof(double) );
	ix = (double *)malloc( (n+1)*sizeof(double) );

	printf("\nCo-efficients of the polynomial ");
	printf("(a_%d.x^%d", n, n);
	for(i = n-1; i >= 2; i--)  printf(" + a_%d.x^%d", i, i);
	printf(" + a_1.x + a_0 = 0) :\n");

	for(i = n; i >= 0; i--)
	{
		printf("a_%d = ", i);
		scanf("%lf", &a[i]);
	}

	printf("\nTwo initial guesses to the roots             : ");
	scanf("%lf %lf", &r, &s);
	printf("\nA small positive number as converging factor : ");
	scanf("%lf", &eps);

	for(i = n; i >= 0; i--)  aa[i] = a[i]; 		/* Copy the initial co-efficients */

	bairstow(n, nit, aa, r, s, eps);

// 	for(i = n; i >= 0; i--)  a[i] = aa[i];
    
	printf("\n\n");
	printf("Results from Bairstow Method for roots of polynomials\n");
	printf("-----------------------------------------------------\n");
	printf("Given polynomial                 : ");

	if(a[n] < 0)		printf("-");
	if(fabs(a[n]) == 1)	printf("x^%d", n);
	else				printf("%lfx^%d", fabs(a[n]), n);

	for(i = n-1; i >= 1; i--)
	{
		if(a[i] == 0)	continue;

		if(a[i] > 0)		printf(" + ");
		else				printf(" - ");

		if(fabs(a[i]) != 1)	printf("%lf", fabs(a[i]));

		if(i > 1)			printf("x^%d", i);
		else				printf("x");
	}

	if(a[0] > 0)	printf(" + %lf = 0\n", a[0]);
	else			printf(" - %lf = 0\n", a[0]);

	printf("Two initial guesses to the roots : %lf, %lf\n", r, s);
	printf("Converging factor                : %lf\n", eps);

	printf("\n			Roots of the polynomial");
	printf("\n			-----------------------\n");
	printf("SN		Real Part	Imaginary Part		Iteration\n");
	printf("--		---------	--------------		---------\n");

	for(i = n; i >= 1; i--)  printf("%2d  %20.4f %17.4f %18d\n", i, rx[i], ix[i], nit[i]);

	return(0);
}
/********************************************************************/

void bairstow(int n, int *nit, double *a, double ri, double si, double eps)
/* Sub-routine for Bairstow Method for roots of polynomial */
{
	int i, it, nrt = 0;
	double error, denom, *b, *c, r, s, dr, ds, rl[3], im[3];

	void qdroot(double aa, double bb, double cc, double rl[3], double im[3]);

	b = (double *)malloc( (n+1+1)*sizeof(double) );
	c = (double *)malloc( (n+1)*sizeof(double) );

	while(n >= 3)																	/* Step-2 */
	{
// 		for(i = n-1; i >= 0; i--)  a[i] = a[i]/a[n];								/* Step-3 */

		it = 0;
		r = ri;
		s = si;

		do																			/* Step-4 */
		{
			b[n] = a[n];
			b[n-1] = a[n-1] + r*b[n];
			for(i = n-2; i >= 0; i--)  b[i] = a[i] + r*b[i+1] + s*b[i+2];

			c[n] = b[n];															/* Step-5 */
			c[n-1] = b[n-1] + r*c[n];
			for(i = n-2; i >= 1; i--)  c[i] = b[i] + r*c[i+1] + s*c[i+2];
	     
			denom = c[1]*c[3] - c[2]*c[2];											/* Step-6 */

			if(denom == 0)															/* Step-7 */
			{
				printf("\nWrong initial approximation to roots. Try again...\n");
				exit(1);
			}

			dr = (b[1]*c[2]-b[0]*c[3])/denom;										/* Step-8 */
			ds = (b[0]*c[2]-b[1]*c[1])/denom;

			r += dr;																/* Step-9 */
			s += ds;

// 			error = fabs(dr)+fabs(ds);												/* Step-10 */
			error = fabs(b[1])+fabs(b[0])+fabs(dr)+fabs(ds);
			it += 1;
		}while(error > eps);														/* Step-11 */
	  
		qdroot(1, -r, -s, rl, im);													/* Step-12 */

		rx[++nrt] = rl[1];
		ix[nrt] = im[1];
		nit[nrt] = it;
		rx[++nrt] = rl[2];
		ix[nrt] = im[2];
		nit[nrt] = it;


		n -= 2;																		/* Step-13 */
		for(i = n; i >= 0; i--)  a[i] = b[i+2];										/* Step-14 */
	}																				/* Step-15 */

	if(n == 2)																		/* Step-16 */
	{
		it = 1;

		qdroot(a[2], a[1], a[0], rl, im);											/* Step-17 */

		rx[++nrt] = rl[1];
		ix[nrt] = im[1];
		nit[nrt] = it;
		rx[++nrt] = rl[2];
		ix[nrt] = im[2];
		nit[nrt] = it;
	}
       
	if(n == 1)																		/* Step-19 */
	{
		it = 1;
		rx[++nrt] = -a[0]/a[1];
		ix[nrt] = 0;
		nit[nrt] = it;
	}

	return;
}
/********************************************************************/

void qdroot(double a, double b, double c, double real[3], double imag[3])
/* Sub-routine for quadratic roots of a.x^2+b.x+c=0 */
{
	double disc;

	disc = b*b-4*a*c;																/* Step-1 */
	if(disc >= 0)																	/* Step-2 */
	{	
		real[1] = (-b+sqrt(disc))/(2*a);
		real[2] = (-b-sqrt(disc))/(2*a);
		imag[1] = 0;
		imag[2] = 0;
	}
	else																			/* Step-3 */
	{
		real[1] = -b/(2*a);
		real[2] = real[1];
		imag[1] = sqrt(-disc)/(2*a);
		imag[2] = -imag[1];
	}

	return;
}
/********************************************************************/