/* Laguerre.cp */

#include "Complex.h"
#include "myTypedef.h"
#include "ExtComplex.h"
#include "LaguerreOld.h"
#include <stdlib.h>
#include "fpMath.h"
#include "fpConv.h"
#include "fpFuncs.h"

extern INT32	blockPrec;

// see Numerical Recipes, p. 139
// I'll store coefficients of a polynomial from low degree (constant term),c[0] to high, c[n}
// where n is the degree of the polynomial

// in NR the degree of polynomial os N-1, and coefficients are C[1] (constant term) to C[N]
// coefficient of x^N-1, so N = n+1 and my indices are one less than NR's

// divides polynomial of degree n by (x-a), giving new polynomial of degree n-1 (a becomes x here)
// roots should be divided out from smallest to largest if forward deflation is used

void forwardDeflation(ExtComplex *c, INT32& n, ExtComplex& x)
{
	static ExtComplex	temp, temp1;
	static bool			arrayGood=false;
	long				i;
	
	if(!arrayGood)
	{
		init(temp);
		init(temp1);
		arrayGood = true;
	}
	
	if(n<=0)
		return;
		
	EquateC_C(temp, c[n]);
	assignZeroC(c[n]);
	
	for(i=n-1;i>=0;i--)
	{
		EquateC_C(temp1, c[i]);
		EquateC_C(c[i], temp);
		Mul(temp, temp, x);
		Add(temp, temp1, temp);
	}
	
	n--;
	
}/* forwardDeflation */


// for input x calculates value of polynomial
void Evaluate(Complex *a, INT32 n, Complex& x, Complex& P)
{
	long			i;
	
	P.Re = a[n].Re;
	P.Im = a[n].Im;
	
	for(i=n-1;i>=0;i--)
	{
		Mul(P, P, x);
		Add(P, P, a[i]);
	}
	
}/* Evaluate */


// for input x calculates value of polynomial
void Evaluate(ExtComplex *a, INT32 n, ExtComplex& x, ExtComplex& P)
{
	long			i;
	
	EquateC_C(P, a[n]);
	
	for(i=n-1;i>=0;i--)
	{
		Mul(P, P, x);
		Add(P, P, a[i]);
	}
	
}/* Evaluate */


// for input x calculates value of polynomial, P, its 1st derivative, P1, and 2nd derivative, P2
void derivPoly(ExtComplex *c, INT32 n, ExtComplex& x,
							ExtComplex& P, ExtComplex& P1, ExtComplex& P2)
{
	long				i;
	static fp			two;
	static bool			arrayGood=false;
	
	if(!arrayGood)
	{
		two.i.n = two.i.nn = 1;
		two.i.b = (UINT32*)malloc(two.i.nn*sizeof(UINT32));
		two.i.b[0] = 2;
		two.e = 0;
		arrayGood = true;
	}
	
	EquateC_C(P, c[n]);
	assignZeroC(P1);
	assignZeroC(P2);
	
	for(i=n-1;i>=0;i--)
	{
		Mul(P2, P2, x);
		Add(P2, P2, P1);
		
		Mul(P1, P1, x);
		Add(P1, P1, P);
		
		Mul(P, P, x);
		Add(P, P, c[i]);
	}
	
	
	mul(P2.Re, two, P2.Re);
	
	mul(P2.Im, two, P2.Im);
	
}/* derivPoly */


// Laguerre refines guess x to a root
bool LaguerreOld(ExtComplex *c, INT32 n, ExtComplex& x, fp& dxAbs, fp& fudgeFactor)
{
	static ExtComplex		x1, dx, G, G2, GP, GM, H, P, P1, P2; // P-poly value, P1-1st derivative, P2- 2nd derivative/2
	static fp				epss;
	static fp				err, xOld, xAbs, tempX, two, ten, tempn, tempn1, tempGP, tempGM, myBitPow;
	INT32					i, iter, bitPrecision;
	static ExtComplex		temp;
	INT32					maxNumIters;
	static bool				arrayGood=false;
	
	if(!arrayGood)
	{
		init(x1);
		init(dx);
		init(G);
		init(G2);
		init(GM);
		init(H);
		init(P);
		init(P1);
		init(P2);
		init(temp);
		
		init(epss);
		init(err);
		init(xOld);
		init(xAbs);
		init(tempX);
		init(tempn);
		init(tempn1);
		init(tempGP);
		init(tempGM);
		init(myBitPow);

		two.i.n = two.i.nn = 1;
		two.i.b = (UINT32*)malloc(two.i.nn*sizeof(UINT32));
		two.i.b[0] = 2;
		two.e = 0;
		
		ten.i.n = ten.i.nn = 1;
		ten.i.b = (UINT32*)malloc(ten.i.nn*sizeof(UINT32));
		ten.i.b[0] = 10;
		ten.e = 0;
		arrayGood = true;
	}
	
	maxNumIters = 400 + 5*n + n*n/2; // from investigation of x^n-1
	
	bitPrecision = 32*blockPrec;
	
	equate(myBitPow, (double)-bitPrecision/4); // sensitive!
	
	pow(epss, ten, myBitPow);
	div(epss, epss, ten);
	div(epss, epss, ten);
	
	//conv(epss, 1.e-16); // to be adjusted was 6e-32, now 6.e-34
	// 1.e-18 gives problem for truncated version
	// 1.e-16 gives about 8 sig. figs.
	
	if(n<1)
		return false;
		
	if(isEqualZero(c[0]))
	{
		assignZeroC(x);
		return true;
	}
	
	
	init(tempn, n);
	init(tempn1, n-1);
	AbsC(xOld, x);
	//dxOldAbs = 0;
	dxAbs.i.n = 0;
	assignZeroC(dx);
	
	for(iter=0;iter<maxNumIters;iter++)
	{
		EquateC_C(P, c[n]);
		AbsC(err, P);
		assignZeroC(P1);
		assignZeroC(P2);
		AbsC(xAbs, x);
		
		for(i=n-1;i>=0;i--)
		{
			Mul(P2, P2, x);
			Add(P2, P2, P1);
			
			Mul(P1, P1, x);
			Add(P1, P1, P);
			
			Mul(P, P, x);
			Add(P, P, c[i]);
			mul(err, xAbs, err);
			
			AbsC(tempX, P);
			add(err, err, tempX);
		}
		
		if(n==1)
			assignZeroC(P2);
		
		mul(err, err, epss);
		AbsC(tempX, P);

		if(compare(tempX, err)<=0)
			return true; //we have the root, This is where all my successful examples exit.
			
		// now we use Laguerre's formula
		
		Div(G, P1, P);
		Mul(G2, G, G);
		Div(H, P2, P);
		// without the two factor the truncated example works, but the exact does not!
		// with the two factor the exact works, but the truncated does not
		// start here was root.Re = .271828; root.Im = .031416;
		// I wonder if we fail should we try the other alternative? Make the two
		// (or fudgeFactor) an input with the value one or two
		
		mul(H.Re, H.Re, fudgeFactor);
		mul(H.Im, H.Im, fudgeFactor);
		Sub(H, G2, H);
		
		EquateC_C(temp, H);
		
		mul(temp.Re, temp.Re, tempn);
		mul(temp.Im, temp.Im, tempn);
		Sub(temp, temp, G2);
		
		mul(temp.Re, temp.Re, tempn1);
		mul(temp.Im, temp.Im, tempn1);
		mySqrtCC(temp, temp);
		Add(GP, G, temp);
		Sub(GM, G, temp);
		
		AbsC(tempGP, GP);
		AbsC(tempGM, GM);
		//if(compareE_E(tempGP, tempGM)==-1)
		
		if(compare(tempGP, tempGM)<0)
			EquateC_C(GP, GM);
		
		equate(dx.Re, tempn);
		dx.Im.i.n = 0;
		Div(dx, dx, GP);
		
		Sub(x1, x, dx);
		
		if(!compare(x.Re, x1.Re) && !compare(x.Im, x1.Im))
			return true;
			
		EquateC_C(x, x1);
		AbsC(dxAbs, dx);
		/*
		if(iter>6 && dxAbs>=dxOldAbs)
			return true;
		*/
		
		//if(iter>40 && compareE_E(dxAbs, dxOldAbs)==1)
		/*
		if(iter>40 && dxAbs>dxOldAbs)
			return true;
		 */
		
		//dxOldAbs = dxAbs;
	}
	
	return false;

}/* Laguerre */

