/* bairstow.c -- Bairstow's complex root finder
 *
 *	(C) 1991, C. Bond
 *
 *	Finds all real and complex roots of polynomials
 *	with real coefficients.
 *
 *	Features:
 *		o Global, self-starting method,
 *		o Does not require initial estimates,
 *		o Finds all roots and quadratic factors,
 *
 *	Caveats:
 *		All polynomial root finders have problems
 *		maintaining accuracy in the presence of
 *		repeated roots.
 *		This is no exception!
 */

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

#define MAX_TERMS 21
int precision_error_flag;
double a[MAX_TERMS],b[MAX_TERMS],c[MAX_TERMS],d[MAX_TERMS];

void find_poly_roots(int n)
{
        double r,s,dn,dr,ds,drn,dsn,eps;
	int i,iter;

	r=s=0;
	dr=1.0;
	ds=0;
	eps=1e-14;
	iter=1;

	while ((fabs(dr)+fabs(ds)) > eps) {
		if ((iter % 200) == 0) {
			r=(double)rand()/16000.;
		}
		if ((iter % 500) == 0) {
			eps*=10.0;
			precision_error_flag=1;
			printf("Loss of precision\n");
		}
		b[1]=a[1] - r;
		c[1]=b[1] - r;

		for (i=2;i<=n;i++){
			b[i]=a[i] - r * b[i-1] -s * b[i-2];
			c[i]=b[i] - r * c[i-1] -s * c[i-2];
		}
		dn=c[n-1] * c[n-3] - c[n-2] * c[n-2];
		drn=b[n] * c[n-3] - b[n-1] * c[n-2];
		dsn=b[n-1]*c[n-1] - b[n]*c[n-2];

		if (fabs(dn) < 1e-16) {
			dn=1;
			drn=1;
			dsn=1;
		}
		dr=drn/dn;
		ds=dsn/dn;

		r+=dr;
		s+=ds;
		iter++;
	}
	for (i=0;i<n-1;i++) a[i]=b[i];
	a[n]=s;
	a[n-1]=r;
}


main()
{
	int i,n,order=0;
	double tmp;

	while ((order < 2) || (order > MAX_TERMS-1)) {
		printf("Polynomial order (2-20): ");
		scanf("%d",&order);
	}
	printf("Enter coefficients, high order to low order.\n");
	for (i=0;i<=order;i++){
		printf("a(%d)=",i);
		scanf("%lf",&a[i]);
		if (i==0) {
			tmp=a[i];
			a[i]=1.0;
		}
		else a[i]/=tmp;
		d[i]=a[i];
 	}
	b[0]=c[0]=1.0;
	n=order;
	precision_error_flag=0;
	while (n > 2) {
		find_poly_roots(n);
                n-=2;
	}
	printf("The quadratic factors are:\n");
	
	for (i=order;i>=2;i-=2) /* print quadratics */
		printf("t^2  %+.15lg t  %+.15lg\n",a[i-1],a[i]);
	if ((n % 2) == 1)
		printf("The linear term is: \nt  %+.15lg\n",a[1]);
}