/* Laguerre.cp */

#include "fpComplex.h"
#include "Laguerre.h"
#include <stdlib.h>
#include <iostream>

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(fpComplex*& c, long& n, fpComplex& x)
{
	fpComplex			temp, temp1;
	long				i;
	
	if(n<=0)
		return;
	
	temp = c[n];
	c[n] = 0.;
	
	for(i=n-1;i>=0;i--)
	{
		temp1 = c[i];
		c[i] = temp;
		temp = temp*x + temp1;
	}
	
	n--;
	
}/* forwardDeflation */


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


// for input x calculates value of polynomial, P, its 1st derivative, P1, and 2nd derivative, P2
void derivPoly(fpComplex*& c, long n, fpComplex& x,
							fpComplex& P, fpComplex& P1, fpComplex& P2)
{
	long	i;
	
	P = c[n];
	P1 = 0.;
	P2 = 0.;
	
	
	for(i=n-1;i>=0;i--)
	{
		P2 = P2*x + P1;
		P1 = P1*x + P;
		P = P*x + c[i];
	}
	
	P2 = 2*P2;
	
}/* derivPoly */


// Laguerre refines guess x to a root
bool Laguerre(fpComplex*& c, long n, fpComplex& x, fp& dxAbs, fp& fudgeFactor)
{
	fpComplex		x1, dx, G, G2, GP, GM, H, P, P1, P2; // P-poly value, P1-1st derivative, P2- 2nd derivative/2
	fp				epss;
	fp				err, xOld, xAbs, tempX, two, ten, tempn, tempn1, tempGP, tempGM, myBitPow;
	long			i, iter, bitPrecision;
	fpComplex		temp;
	long			maxNumIters;
	
	maxNumIters = 400 + 5*n + n*n/2; // from investigation of x^n-1 was 400
	
	// x^8 + x^6 + x^4 + x^2 - 1e100 shows dependence on size of coefficients; need 400 -> 800 for good result
	
	bitPrecision = 32*blockPrec;
	
	two = 2;
	ten = 10;
	
	myBitPow = (double)-bitPrecision/4; // sensitive!
	
	pow(epss, ten, myBitPow);
	epss = epss/ten;
	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(c[0]==0.)
	{
		x = 0.;
		return true;
	}
	
	tempn = n;
	tempn1 = n - 1;
	xOld = abs(x);
	dxAbs = 0.;
	dx = 0.;
	
	for(iter=0;iter<maxNumIters;iter++)
	{
		P = c[n];
		err = abs(P);
		P1 = 0.;
		P2 = 0.;
		xAbs = abs(x);
		
		for(i=n-1;i>=0;i--)
		{
			P2 = P2*x + P1;
			P1 = P1*x + P;
			P = P*x + c[i];
			err = xAbs*err;
			tempX = abs(P);
			err = err + tempX;
		}
		
		if(n==1)
			P2 = 0.;

		err = err*epss;
		tempX = abs(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
		
		G = P1/P;
		G2 = G*G;
		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
		
		H = H*fudgeFactor;
		H = G2 - H;
		
		temp = H;
		temp = temp*tempn;
		temp = temp - G2;
		
		temp = temp*tempn1;
		temp = sqrt(temp);  // is this sqrt same as old mySqrtCC?
		GP = G + temp;
		GM = G - temp;
		
		tempGP = abs(GP);
		tempGM = abs(GM);
		//if(compareE_E(tempGP, tempGM)==-1)
		
		if(compare(tempGP, tempGM)<0)
			GP = GM;
		
		dx.re = tempn;
		dx.im = 0.;
		dx = dx/GP;
		
		x1 = x - dx;
		
		if(!compare(x.re, x1.re) && !compare(x.im, x1.im))
			return true;
		
		x = x1;
		dxAbs = abs(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 true;

}/* Laguerre */

