/* myFindRoots.cp */

#include "myTypedef.h"
#include "myFindRoots.h"
#include "ExtComplex.h"
#include "LaguerreOld.h"
#include "fpConv.h"
#include "fpMath.h"
#include "stdlib.h"

// x array has dimension n+1, z array has dimension n since it will hold the roots
bool myFindRoots(Complex *z, Complex *x, int n)
{
	ExtComplex		*zE, *xE;
	int             i;
	
	zE = (ExtComplex *)malloc(n*sizeof(ExtComplex));
	for(i=0;i<n;i++)
		init(zE[i]);
	xE = (ExtComplex *)malloc((n+1)*sizeof(ExtComplex));
	for(i=0;i<=n;i++)
		init(xE[i]);
	
	for(i=0;i<=n;i++)
	{
		equate(xE[i].Re, x[i].Re);
		equate(xE[i].Im, x[i].Im);
	}
	
	if(!myFindRoots(zE, xE, n))
	{
		for(i=0;i<n;i++)
		{
			if(zE[i].Re.i.n && zE[i].Re.i.b)
				free(zE[i].Re.i.b);
			if(zE[i].Im.i.n && zE[i].Im.i.b)
				free(zE[i].Im.i.b);
		}
		free(zE);
		
		for(i=0;i<=n;i++)
		{
			if(xE[i].Re.i.n && xE[i].Re.i.b)
				free(xE[i].Re.i.b);
			if(xE[i].Im.i.n && xE[i].Im.i.b)
				free(xE[i].Im.i.b);
		}
		free(xE);
		
		return false;
	}
	
	
	for(i=0;i<n;i++)
	{
		if(!equate(z[i].Re, zE[i].Re))
			return false;
		if(!equate(z[i].Im, zE[i].Im))
			return false;
	}
	
	for(i=0;i<n;i++)
	{
		if(zE[i].Re.i.n && zE[i].Re.i.b)
			free(zE[i].Re.i.b);
		if(zE[i].Im.i.n && zE[i].Im.i.b)
			free(zE[i].Im.i.b);
	}
	free(zE);
	
	for(i=0;i<=n;i++)
	{
		if(xE[i].Re.i.n && xE[i].Re.i.b)
			free(xE[i].Re.i.b);
		if(xE[i].Im.i.n && xE[i].Im.i.b)
			free(xE[i].Im.i.b);
	}
	free(xE);
	
	return true;
	
}/* myFindRoots */


// x array has dimension n+1, z array has dimension n since it will hold the roots
bool myFindRoots(ExtComplex *z, ExtComplex *x, int n)
{
	static ExtComplex		root;
	INT32					i;
	static fp				dxAbs, fudgeFactor;
	static bool				arrayGood=false;
	
	if(!arrayGood)
	{
		init(root);
		init(dxAbs);
		init(fudgeFactor);
		arrayGood = true;
	}
	
	// get the roots
	i=0;
	while(n>0)
	{
		if(x[0].Re.i.n==0 && x[0].Im.i.n==0)
		{
			//root.Re.i.n = 0;
			//root.Im.i.n = 0;
			equate(root.Re, 0.);
			equate(root.Im, 0.);
		}
		else
		{
			equate(root.Re, .271828);
			equate(root.Im, .031416);
		}
		
		
		equate(fudgeFactor, 1.);
		if(!LaguerreOld(x, n, root, dxAbs, fudgeFactor))
		{
			equate(fudgeFactor, 2.);
			if(!LaguerreOld(x, n, root, dxAbs, fudgeFactor))
			{
				equate(fudgeFactor, 1.);
				equate(root.Re, -.0271828);
				equate(root.Im, .31416);
				if(!LaguerreOld(x, n, root, dxAbs, fudgeFactor))
				{
					equate(fudgeFactor, 2.);
					if(!LaguerreOld(x, n, root, dxAbs, fudgeFactor))
						return false;
				}
			}
		}
			
		EquateC_C(z[i++], root);
		
		forwardDeflation(x, n, root);
	}
	
	return true;
	
}/* myFindRoots */
