/*
 *  orderRoots.cpp
 *  PolyRoots
 *
 *  Created by Robert Delaney on 4/24/06.
 *  Copyright 2006 __MyCompanyName__. All rights reserved.
 *
 */

#include "myTypedef.h"
#include "ExtComplex.h"
#include "orderRoots.h"
#include "fpConv.h"
#include <stdlib.h>
#include "fpMath.h"

/*
	We take the array root. Remember bitPrecision, then change to
	bitPrecision=4*outputPrecision. 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. Restore
	original bitPrecision. Precision stuff didn't work as expected.
*/

void orderRoots(ExtComplex *root, int n)
{
	int				i, j, iMax, n1;
	ExtComplex		*rootCopy;
	static fp		temp, theMax;
	static bool		arrayGood=false;
	
	if(!arrayGood)
	{
		init(temp);
		init(theMax);
		arrayGood = true;
	}
	
	//bitPrecision = RR::precision();
	//RR::SetPrecision(4*RR::OutputPrecision());
	
	rootCopy = (ExtComplex *)malloc(n*sizeof(ExtComplex));
	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
		EquateC_C(rootCopy[j++], root[iMax]);
		for(i=iMax;i<(n1-1);i++)
			EquateC_C(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++)
		EquateC_C(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);
	
	//RR::SetPrecision(bitPrecision);
	
}/* orderRoots */