/*
 *  findPolyRoots.cpp
 *  fp
 *
 *  Created by Robert Delaney on 1/16/11.
 *  Copyright 2011 Bob Delaney's Science Software. All rights reserved.
 *
 */

#include "fpComplex.h"
#include "findPolyRoots.h"
#include "Laguerre.h"
#include "stringUtils.h"
#include <string.h>
#include "charPoly.hpp"
#include "matcConv.hpp"

/*
 Give this inputStr and it looks for an 'x'. If it does not find 'x' it
 returns false.
 
 If it does find an x, it finds the highest power of x and returns that
 in the n variable
 */

bool xPresence(long& n, const char* inputStr)
{
	long		i, j, length, thePower;
	char		*dataString, ch;
	bool		xExists;
	
	n = 0;
	
	length = strlen(inputStr);
	if(length==0)
		return false;
	
	// search for any x
	xExists = false;
	for(i=0;i<length;i++)
		if(inputStr[i]=='x')
		{
			xExists = true;
			break;
		}
	
	if(!xExists)
		return false;
	
	// now we have x's
	
	dataString = (char *)malloc((length+1)*sizeof(char));
	
	i = 0;
getTheNext:
	// search for next x
	if(i>=length)
	{
		goto myExitTrue;
	}
	ch = inputStr[i];
	while(i<length && ch!='x')
		ch = inputStr[++i]; // i points to current value of ch
	
	// if ch='x' on first try, i=0 and points to 'x'
	// on later tries i still points to 'x' or else i=length and ch=0
	
	if(i>=length)
	{
		goto myExitTrue;
	}
	
	// now i points to 'x'
	// get rid of white space, if there is any
	ch = inputStr[++i];
	while(i<length && ch<=32)
		ch = inputStr[++i]; // i points to current value of ch
	
	if(i>=length)
	{
		thePower = 1;
		if(thePower>n)
			n = thePower;
		goto myExitTrue;
	}
	
	// 
	
	// is ch a '^' ?
	
	if(ch!='^')
	{
		thePower = 1;
		if(thePower>n)
			n = thePower;
		i++;
		goto getTheNext;
	}
	
	// now ch='^' and i points to it
	
	if(i>=(length-1))
	{
		goto myExitFalse;
	}
	
	// get rid of white space, if there is any
	ch = inputStr[++i];
	while(i<length && ch<=32)
		ch = inputStr[++i]; // i points to current value of ch
	
	if(i>=length)
	{
		goto myExitFalse;
	}
	
	// get the power
	j = 0;
	if(ch=='+' || ch=='-')
	{
		dataString[j++] = ch;
		ch = inputStr[++i];
	}
	
	while(i<length && ch!='+' && ch!='-')
	{
		dataString[j++] = ch;
		ch = inputStr[++i]; // i points to current value of ch
	}
	dataString[j] = 0;
	
	if(sscanf(dataString, "%d", &thePower)!=1)
		goto myExitFalse;
	
	if(thePower>n)
		n = thePower;
	
	goto getTheNext;
	
myExitTrue:
	free(dataString);
	return true;
	
myExitFalse:
	n = 0;
	free(dataString);
	return false;
	
}/* xPresence */


bool xNoPresence(long& n, char *sep, const char* inputStr)
{
	long		numFields;
	
	numFields = CountFields(inputStr, sep);
	// worry about last field being null; note that a null field fails to compare to ""!
	if(strlen(NthField(inputStr, sep, numFields))==0)
		numFields--;
	
	n = numFields - 1;
	
	if(!n)
		return false;
	
	return true;
	
}/* xNoPresence */


// given theField string of maximum form num/den,num/den this returns fpComplex theNumber
static bool getTheNumber(fpComplex& theNumber, const char* theField)
{
	char	*reString, *imString, *numString, *denString;
	fp		num, den;
	long		i, len;
	char	ch;
	
	len = strlen(theField);
	
	// for inputPoly_x
	if(!len)
	{
		theNumber.re = 1.;
		theNumber.im = 0.;
		return true;
	}
	
	if(len==1)
	{
		ch = theField[0];
		if(ch=='+')
		{
			theNumber.re = 1.;
			theNumber.im = 0.;
			return true;
		}
		
		if(ch=='-')
		{
			theNumber.re = -1.;
			theNumber.im = 0.;
			return true;
		}
	}
	
	// check no unwanted characters
	for(i=0;i<len;i++)
	{
		ch = theField[i];
		if(!(ch==' ' || ch=='(' || ch==')' || ch==',' || ch=='/' || ch=='^' || ch=='+' || ch=='-' || ch=='.' || ch=='e' || ch=='E' || (ch>='0' && ch<='9')))
		{
			return false;
		}
	}
	
	if(!VerifyNumComplexFracString2(reString, imString, theField))
		return false;
	if(!VerifyNumFracString2(numString, denString, reString))
		return false;
	if(!equate(num, numString))
		return false;
	if(!equate(den, denString))
		return false;
	if(den.i.n==0)
		return false;
	
	
	
	theNumber.re = num/den;
	free(numString);
	free(denString);
	
	if(!VerifyNumFracString2(numString, denString, imString))
		return false;
	if(!equate(num, numString))
		return false;
	if(!equate(den, denString))
		return false;
	if(den.i.n==0)
		return false;
	theNumber.im = num/den;
	
	free(reString);
	free(imString);
	free(numString);
	free(denString);
	
	return true;
	
}/* getTheNumber */


/*
 Allows input of polynomial in forms:
 x^2-4x+4
 4-4x+x^2
 -4x+x^2+4
 x^2-(2,3)x+(4,-5)
 x^5+1
 etc.
 We presume that white spaces have been removed, and that the order n of the polynomial
 has been determined. coeffArray has already been created, but not yet initialized to zero, except for index n.
 This sub fills this array and returns true or false.
 */
bool inputPoly_x(fpComplex*& coeffArray, long n, const char *inputStr)
{
	long				i, j, length, thePower;
	char				ch;
	char				*dataString, *singleString;
	fpComplex			theNumber;
	bool				isNegative;
	char				*newInputStr;
	
	length = strlen(inputStr);
	
	if(n==0 || length==0)
		return false;
	
	// check no unwanted characters
	for(i=0;i<length;i++)
	{
		ch = inputStr[i];
		if(!(ch==0xA || ch==0xD || ch=='x' || ch==' ' || ch=='(' || ch==')' || ch==',' || ch=='/' || ch=='^' || ch=='+' || ch=='-' || ch=='.' || ch=='e' || ch=='E' || (ch>='0' && ch<='9')))
		{
			return false;
		}
	}
	
	if(inputStr[0]=='-' || inputStr[0]=='+')
	{
		newInputStr = (char *)malloc((length+1)*sizeof(char));
		for(i=0;i<=length;i++)
			newInputStr[i] = inputStr[i];
	}
	else
	{
		newInputStr = (char *)malloc((length+2)*sizeof(char));
		newInputStr[0] = '+';
		for(i=0;i<=length;i++)
			newInputStr[i+1] = inputStr[i];
	}
		
	length = strlen(newInputStr);
	
	dataString = (char *)malloc((length+1)*sizeof(char));
	singleString = (char *)malloc((length+1)*sizeof(char));
	
	for(i=0;i<n;i++)
	{
		coeffArray[i] = 0.;
	}
	
	i = 0;
getNext:
	j = 0;
	ch = newInputStr[i];
	if(ch=='+' || ch=='-')
	{
		dataString[j++] = ch;
		ch = newInputStr[++i];
	}
	if(ch=='(')
	{	
		ch = newInputStr[++i];  // added to not use '('
		while(ch!=')' && i<length)
		{
			dataString[j++] = ch;
			ch = newInputStr[++i];
		}
		if(i>=length)
			goto exitFalse;
		//dataString[j++] = ch;  // don't use ')'
		ch = newInputStr[++i];
	}
	else
	{
		while(ch!='x' && ch!=0xA && ch!=0xD && i<length)
		{
			dataString[j++] = ch;
			ch = newInputStr[++i]; // i points to current value of ch
		}
	}
	dataString[j] = 0;
	
	
	if(dataString[0]=='-')
	{
		isNegative = true;
		dataString[0] ='+';
	}
	else {
		isNegative=false;
	}
	
	// get theNumber
	if(!getTheNumber(theNumber, dataString))
		goto exitFalse;
	
	if(isNegative)
	{
		theNumber.re = -theNumber.re;
		theNumber.im = -theNumber.im;
	}
	
	//cout << "getTheNumber " << theNumber << endl;
	
	// get thePower
	
	if(ch!='x' &&(i>=(length-1) || ch=='+' || ch=='-'))
	{
		thePower = 0;
		coeffArray[thePower] = coeffArray[thePower] + theNumber;
	}
	
	if(ch!='x' && i>=(length-1))
		goto exitTrue;
	
	if(ch=='+' || ch=='-')
		goto getNext;
	
	// ch should be an x
	if(ch!='x')
		goto exitFalse;
	
	// is next a ^ ?
	
	ch = newInputStr[++i];
	if(ch!='^')
	{
		thePower = 1;
		coeffArray[thePower] = coeffArray[thePower] + theNumber;
		if(i>=(length-1))
			goto exitTrue;
		goto getNext;
	}
	
	// now we must have a power
	if(i>=(length-1))
		goto exitFalse;
	
	// ch is now ^ and i points to it
	
	// get the power
	j = 0;
	ch = newInputStr[++i];
	if(ch=='+' || ch=='-')
	{
		singleString[j++] = ch;
		ch = newInputStr[++i];  // got '='
	}
	while(i<length && ch!='+' && ch!='-')
	{
		singleString[j++] = ch;
		ch = newInputStr[++i]; // i points to current value of ch
	}
	singleString[j] = 0;
	
	if(sscanf(singleString, "%ld", &thePower)!=1)
		goto exitFalse;
	
	if(thePower<0)
		goto exitFalse;
	
	coeffArray[thePower] = coeffArray[thePower] + theNumber;
	
	if(i>=(length-1))
		goto exitTrue;
	
	goto getNext;
	
exitTrue:
	while(coeffArray[n]==0. && n>0)
		n--;  // do realloc?
	
    free(dataString);
    free(singleString);
    free(newInputStr);
    
	if(!n || n==-1)
    {
        
        return false;
    }
	
	return true;
	
exitFalse:
	free(dataString);
	free(singleString);
	free(newInputStr);
	return false;
	
}/* inputPoly_x */


// coeffArray[n] is to contain coeff of x^n and so on down
// inputStr contains the coeffs from n on down to 0, but NthField numbers them n+1 down to 1
bool inputPoly(fpComplex*& coeffArray, long n, const char *sep, const char *inputStr)
{	
	fp		num, den;
	int		i;
	char	*theField;
	
	// fill coeffArray
	for(i=0;i<=n;i++)
	{
		theField = NthField(inputStr, sep, i+1);
		if(!getTheNumber(coeffArray[n-i], theField))
			return false;
	}
	
	while(coeffArray[n]==0.)
		n--;  // do realloc?
	
	if(!n || n==-1)
		return false;
	
	return true;
	
}/* inputPoly */


// x array has dimension n+1, z array has dimension n since it will hold the roots
bool myFindRoots(fpComplex*& z, fpComplex*& x, long n)
{
	fpComplex		root;
	long			i;
	fp				dxAbs, fudgeFactor;
	
    // init z
    
    z = (fpComplex*)malloc(n*sizeof(fpComplex));
    
    for(i=0;i<n;++i)
    {
        init(z[i].re);
        init(z[i].im);
    }
    
	// get the roots
	i=0;
	while(n>0)
	{
		if(x[0].re.i.n==0 && x[0].im.i.n==0)
		{
			root = 0.;
		}
		else
		{
			root.re = .271828;
			root.im = .031416;  // .031416;
		}
		
		fudgeFactor = 1.;
		if(!Laguerre(x, n, root, dxAbs, fudgeFactor)) // Laguerre now always returns true so below is ignored
		{
			fudgeFactor = 1.5;  // was 2
			//cout << "here 1" << endl;
			if(!Laguerre(x, n, root, dxAbs, fudgeFactor))
			{
				//cout << "here 2" << endl;
				fudgeFactor = 1.;
				root.re = -.0271828;
				root.im = 0.;  // .31416;
				if(!Laguerre(x, n, root, dxAbs, fudgeFactor))
				{
					//cout << "here 3" << endl;
					fudgeFactor = 2.;
					if(!Laguerre(x, n, root, dxAbs, fudgeFactor))
						return false;
				}
			}
		}
		
		z[i++] = root;
		
		forwardDeflation(x, n, root);
	}
	
	return true;
	
}/* myFindRoots */


/*
 We take the array root.
 Copy to new array from largest Re
 to smallest. Then test and for conjugate pairs have the positive
 Im part preceed the negative one. Then copy back to root.
 */

void orderRoots(fpComplex*& root, long n)
{
	long			i, j, iMax, n1;
	fpComplex		*rootCopy;
	fp				temp, theMax;
	
	
	rootCopy = (fpComplex *)malloc(n*sizeof(fpComplex));
	for(i=0;i<n;i++)
	{
		init(rootCopy[i]);
	}
	
	
	n1 = n;
	j = 0;
	while(n1>0)
	{
		equate(theMax, root[0].re);
		iMax = 0;
		for(i=1;i<n1;i++)
		{
			equate(temp, root[i].re);
			
			if(compare(temp, theMax)>0)
			{
				equate(theMax, temp);
				iMax = i;
			}
		}
		// transfer max root, collapse root, and dec n1
		rootCopy[j++] = root[iMax];
		for(i=iMax;i<(n1-1);i++)
			root[i] = root[i+1];
		n1--;
	}
	// now work on rootCopy and find pairs to reorder Im parts
	i = 0;
	while(i<(n-1))
	{
		
		if(!compare(rootCopy[i].re, rootCopy[i+1].re))
		{
			// found a pair
			
			if(compare(rootCopy[i].im, rootCopy[i+1].im)<0)
			{
				// switch
				rootCopy[i].im.i.n = -rootCopy[i].im.i.n;
				rootCopy[i+1].im.i.n = -rootCopy[i+1].im.i.n;
			}
			i+=2;
		}
		else
			i++;
	}
	
	// copy back to root
	for(i=0;i<n;i++)
		root[i] = rootCopy[i];
	
	
	for(i=0;i<n;i++)
	{
		if(rootCopy[i].re.i.b)
			free(rootCopy[i].re.i.b);
		if(rootCopy[i].im.i.b)
			free(rootCopy[i].im.i.b);
	}
	
	free(rootCopy);
	
	
}/* orderRoots */



// if true it returns the roots of the polynomial, and the P(root), else it returns false

bool findPolyRoots(fpComplex*& root, fpComplex*& Parray, fpComplex*& coeffArray, long n, bool doRoundToZero)
{
	long		i;
	fpComplex	*coeffArrayCopy, P, P1, temp;
	fp			s, sInv; // scale factor
	fp			nInv, sPow;
	bool		isScaled;
	
	isScaled = false;
	
	coeffArrayCopy = (fpComplex*)malloc((n+1)*sizeof(fpComplex));
	
	for(i=0;i<=n;i++)
	{
		init(coeffArrayCopy[i]);
	}
	
	// we have a good coeffArray, copy it
	for(i=0;i<=n;i++)
		coeffArrayCopy[i] = coeffArray[i];
	
	// do we need to scale?
	
	s = abs(coeffArrayCopy[0]);
	if(s>1e20)
	{
		nInv = n;
		nInv = 1/nInv;
		s = pow(s, nInv);
		sInv = 1/s;
		sPow = sInv;
		
		for(i=n-1;i>=0;i--)
		{
			coeffArrayCopy[i] = sPow*coeffArrayCopy[i];
			sPow = sInv*sPow;
		}
		
		isScaled = true;
	}
	
	if(myFindRoots(root, coeffArrayCopy, n))
	{
		if(isScaled)
			for(i=0;i<n;i++)
				root[i] = s*root[i];
		
		if(doRoundToZero)
		{
			for(i=0;i<n;i++)
			{
				// log2(10^16) = 56.47
				if((Lg2(root[i].re)-Lg2(root[i].im))>56)
				{
					Evaluate(coeffArray, n, root[i], P);
					temp = root[i];
					temp.im = 0.;
					Evaluate(coeffArray, n, temp, P1);
					if(abs(P1.re)<=abs(P.re) || abs(Lg2(P.re)-Lg2(P1.re))<2)
						root[i].im = 0.;
				}
				
				if((Lg2(root[i].im)-Lg2(root[i].re))>56)
				{
					Evaluate(coeffArray, n, root[i], P);
					temp = root[i];
					temp.re = 0.;
					Evaluate(coeffArray, n, temp, P1);
					if(abs(P1.im)<=abs(P.im) || abs(Lg2(P.im)-Lg2(P1.im))<2)
						root[i].re = 0.;
				}
			}
		}
		
		orderRoots(root, n);
		
		// add polynomial evaluation
		
		for(i=0;i<n;i++)
		{
			Evaluate(coeffArray, n, root[i], P);
			Parray[i] = P;
		}
	}
	else
	{
		for(i=0;i<=n;i++)
		{
			free(coeffArrayCopy[i].re.i.b);
			free(coeffArrayCopy[i].im.i.b);
		}
		free(coeffArrayCopy);
		return false;
	}
	
	for(i=0;i<=n;i++)
	{
		free(coeffArrayCopy[i].re.i.b);
		free(coeffArrayCopy[i].im.i.b);
	}
	free(coeffArrayCopy);
	return true;

}/* findPolyRoots */


// this will be called in the ComplexMatrix plugin
// x is input as square matrix whose eigenvalues are to be found
// eigenvalue will be created and contain the eigenvalues as Complex(double, double)
bool myEigenvalues(Complex** x, Complex*& eigenvalues, long n)
{
    matc        M, z;
    fpComplex   *p, *roots;
    long        i, j;
    
    init(M, n, n);
    
    for(i=0;i<n;++i)
        for(j=0;j<n;++j)
            M.array[i][j] = x[i][j];
    
    charPoly(M, p); // p is of dim n+1 and returns coeffs. of char. poly. of M
    if(!myFindRoots(roots, p, n))
        return false;
    
    for(i=0;i<n;++i)
        eigenvalues = (Complex*)malloc(n*sizeof(Complex));
    
    // convert roots to eigenvalues
    for(i=0;i<n;++i)
    {
        equate(eigenvalues[i].Re, roots[i].re);
        equate(eigenvalues[i].Im, roots[i].im);
    }
    
    return true;
    
}/* myEigenvalues */


