/*
 *  fpConv.cpp
 *  fp
 *
 *  Created by Robert Delaney on 4/24/09.
 *  Copyright 2009 __Bob Delaney's Science Software__. All rights reserved.
 *
 */

#include "fp.h"

using namespace std;

extern INT32 blockPrec;


static INT32	Max(INT32 x, INT32 y)
{
	if(x>y)
		return x;
	
	return y;
	
}/* Max */

/*
bool isItINT32(const fp& x)
{
	mb		xmb;
	
	if(x.i.n==0)
		return true;
	
	equate(xmb, x);
	
	if(abs(xmb.n)>1)
	{
		return false;
	}
	
	return true;
	
}*/ /* isItINT32 */


// This should only be used on a local static variable in a function, never on a returned fp&
void equateWithCare(fp& x, const fp& y)
{
	equateWithCare(x.i, y.i);
	x.e = y.e;
	
}/* equateWithCare */


void equate(fp& x, const fp& y)
{
	equate(x.i, y.i);
	x.e = y.e;
	
}/* equate */


// x may be inited
bool equate(fp& x, const char *xString)
{
    
	if(x.i.nn && x.i.b)
	{
		free(x.i.b);
		x.i.n = x.i.nn = 0;
		x.i.b = 0;
	}
    
	return init(x, xString);
	
}/* equate */

// does this make .111111111... zero? Don't we want this 1 or -1?
// yes, this is truncation. Maybe we should also have roundEquate
void equate(mb& z, const fp& x)
{
	equate(z, x.i);
	if(x.e>=0)
		mbShiftLeft(z, z, x.e);
	else
		mbShiftRight(z, z, abs(x.e));
	
}/* equate */


void equate(fp& z, const mb& x)
{
	equate(z.i, x);
	z.e = 0;
	
}/* equate */


// x must be inited
void equate(fp& x, double y)
{
    
	if(x.i.nn && x.i.b)
	{
		free(x.i.b);
		x.i.n = x.i.nn = 0;
		x.i.b = 0;
	}
	
	init(x, y);
	
}/* equate */


bool equate(double& z, const fp& x)
{
	char	*xString;
    
	if(x.i.n==0)
	{
		z = 0;
		return true;
	}
	
	//xString = fpConvToString(x, 16); // may have to modify xString
	xString = fpToStr(x, 16);
	if(!sscanf(xString, "%lg", &z)) // %Lg is bad news on an Intel iMac!
	{
		z = 0;
		free(xString);
		return false;
	}
	   
	free(xString);
	
	return true;
	
}/* equate */


// returns true only if x is an INT32
// have to worry about x being an integer +/- a very small number
// especially worry about .99999999 which gives zmb=0, but want to round to 1
bool equate(INT32& z, const fp& x)
{
	mb			zmb;
	fp			zfp;
	
	z = 0;
	if(x.i.n==0)
		return true;
	
	equate(zmb, x);
	if(abs(zmb.n)>1 || zmb.n==0)
	   return false;
	   
	z = zmb.b[0];
	if(z<0)
		return false;
	
	if(zmb.n<0)
		z = -z;
	   
	equate(zfp, zmb);
	if(zfp!=x)
		return false;
	
	return true;
	
}/* equate */


void equate(char*& xString, const fp& x, INT32 decPrec)
{
	xString = fpToStr(x, decPrec);
	
}/* equate */


// rounds xString to decPrec and returns "nice" string for fp x
static char* FormRoundedfpString(const char* xString, INT32 xDecExp, bool isNegative, INT32 myDecPrec)
{
	INT32		i, len, length, index, numTrailZeroes;
	char		*myxString, *outStr, expStr[255];
	
	if(myDecPrec<=0)
	{
		myxString = (char*)malloc(2*sizeof(char));
		myxString[0] = '0';
		myxString[1] = 0;
		return myxString;
	}
	
	len = strlen(xString);  // we're guaranteed that xString length is greater than myDecPrec
	
	myxString = (char*)malloc((myDecPrec+1)*sizeof(char));
	
	for(i=0;i<myDecPrec;i++)
	{
		index = i;
		myxString[i] = xString[i];
	}
	myxString[index+1] = 0;
	
	// since we've lost len-myDecPrec digits we must increase xDecExp by that amount
	xDecExp+=(len-myDecPrec);
	
	length = strlen(myxString);
	
	// need to trim trailing zeroes, but watch out for all zeroes
	numTrailZeroes = 0;
	for(i=length-1;i>=0;i--)
	{
		index = i;
		if(myxString[i]=='0')
			numTrailZeroes++;
		else
			break;
	}
	
	if(numTrailZeroes>0)
	{
		// index points to the left of first zero
		myxString[index+1] = 0;
		xDecExp+=(length-index-1);
		if(myxString[0]=='0' && index==1)
		{
			myxString[1] = 0;
			return myxString;
		}
	}
	
	length = strlen(myxString);
	
	outStr = (char*)malloc((myDecPrec+35)*sizeof(char));
	outStr[0] = 0;
	
	if(xDecExp>3 || xDecExp<-length-3)
	{
		// use scientific notation
		
		if(isNegative)
		{
			strcpy(outStr, "-");
			outStr[1] = myxString[0];
			outStr[2] = 0;
		}
		else
		{
			outStr[0] = myxString[0];
			outStr[1] = 0;
		}
		
		strcat(outStr, ".");
		strcat(outStr, myxString+1);
		strcat(outStr, "e");
		sprintf(expStr, "%d", xDecExp+length-1);
		strcat(outStr, expStr);
	}
	else
	{
		if(xDecExp>=0)
		{
			if(isNegative)
				strcpy(outStr, "-");
			else
				outStr[0] = 0;
			
			strcat(outStr, myxString);
			for(i=0;i<xDecExp;i++)
				strcat(outStr, "0");
		}
		else
		{
			if(xDecExp<=-length)
			{
				if(isNegative)
					strcpy(outStr, "-");
				else
					outStr[0] = 0;
				
				strcat(outStr, "0.");
				for (i = 0; i < -length-xDecExp; i++)
   					strcat(outStr, "0");
				strcat(outStr, myxString);
			}
			else
			{
				if(isNegative)
					strcpy(outStr, "-");
				else
					outStr[0] = 0;
				
				expStr[1] = 0;
				for(i=0;i<length+xDecExp;i++)
				{
					expStr[0] = myxString[i];
					strcat(outStr, expStr);
				}
				strcat(outStr, ".");
				for(i=length+xDecExp;i<length;i++)
				{
					expStr[0] = myxString[i];
					strcat(outStr, expStr);
				}
			}
		}
	}
	
	free(myxString);
	
	return outStr;
	
}/* FormRoundedfpString */


// follows NTL method, gives "pretty" output
char* fpToStr(const fp& x, int myDecPrec)
{
	fp		xt, ten=10, xMan;
	char	*xString, *xManString;
	INT32	bp, logtenxInt, xDecExp;
	bool	isNegative;
	INT32	length;
	
	if(!x.i.n)
	{
		xString = (char*)malloc(2*sizeof(char));
		xString[0] = '0';
		xString[1] = 0;
		return xString;
	}
	
	xt = x;
	isNegative = false;
	if(xt.i.n<0)
	{
		xt.i.n = -xt.i.n;
		isNegative = true;
	}
	
	// only low precision is necessary to get integer part of log to base ten of x
	bp = blockPrec;
	blockPrec = 3;
	
	logtenxInt = to_long(floor(log(xt)/log(ten)));
	
	xDecExp = logtenxInt - myDecPrec - 1;
	
	blockPrec = bp + 1;
	
	mul(xMan, xt, pow(ten, -xDecExp));
	// so must round during shifting
	
	blockPrec = bp;
	
	if(xMan.e>0)
	{
		mbShiftLeft(xMan.i, xMan.i, xMan.e);
		xMan.e = 0;
	}
	else
	{
		mbShiftRight(xMan.i, xMan.i, abs(xMan.e));
		xMan.e = 0;
	}
	
	// xMan.i as a base 10 integer is one or two digits longer than decPrec, so add 5 or 50 to round
	xManString = mb2Str(xMan.i);
	length = strlen(xManString);
	if(length>myDecPrec)
	{
		if(strlen(xManString)>(myDecPrec+1))
			xMan.i+=50;
		else
			xMan.i+=5;
	}
	
	free(xManString);
	
	xManString = mb2Str(xMan.i);
	
	xString = FormRoundedfpString(xManString, xDecExp, isNegative, myDecPrec);
	
	free(xManString);
	
	return xString;
	
}/* fpToStr */

// returns in base ten the integer part of fp as string up to decPrec in digit accuracy
// returns decimal exponent in DecExp
// used in Decimal Plugin
char* fpToStr2(INT32& DecExp, const fp& x, INT32 decPrec)  // has problem
{
	static fp		xt, log_2_10, ten, pointFive;
	static mb		xtInt;
	static bool		initGood=false;
	fp				ln2, ln10, temp, approxX, tenToDecPrec;
	bool			isNegative;
	INT32			i, k, len, log_10_x, temp_p, new_p, saveBlockPrec;
	char			*outStr, *tempStr, ch;
	
	if(!initGood)
	{
		init(xt);
		init(log_2_10, 3.321928094887362);
		//init(temp);
		//init(tenToDecPrec);
		//init(approxX);
		init(xtInt);
		pointFive.i.n = pointFive.i.nn = 1;
		pointFive.i.b = (UINT32*)malloc(pointFive.i.n*sizeof(UINT32));
		if(!pointFive.i.b)
			exit(1);
		pointFive.i.b[0] = 1;
		pointFive.e = -1;
		ten.i.n = ten.i.nn = 1;
		ten.i.b = (UINT32*)malloc(ten.i.n*sizeof(UINT32));
		if(!ten.i.b)
			exit(1);
		ten.i.b[0] = 10;
		ten.e = 0;
		initGood = true;
	}
	
	if(!x.i.n || decPrec<=0)
	{
		outStr = (char*)malloc(2*sizeof(char));
		if(!outStr)
			exit(1);
		outStr[0] = '0';
		outStr[1] = 0;
		DecExp = 0;
		return outStr;
	}
	
	equate(xt, x);
	
	saveBlockPrec = blockPrec;
	
	temp_p = (Max(NumBits(decPrec), NumBits(Lg2(x))) + 10)/blockBits;
	if(!temp_p)
		temp_p = 1;
	
	blockPrec = temp_p;
	
	// approximate log to base 10 of x
	/*
	init(ln2, log(2.));
	init(ln10, log(10.));
	div(log_2_10, ln10, ln2);
	 */
	
	
	mul(temp, log_2_10, decPrec);
	div(temp, temp, blockBits);
	new_p = to_long(temp) + 2;
	
	div(temp, Lg2(x), log_2_10);
	
	log_10_x = to_long(temp);
	
	blockPrec = new_p;
	
	if(xt.i.n<0)
	{
		xt.i.n = -xt.i.n;
		isNegative = true;
	}
	else
		isNegative = false;
	
	k = decPrec - log_10_x;
	
	power(tenToDecPrec, ten, decPrec);
	power(approxX, ten, log_10_x);
	
	div(xt, xt, approxX);
	mul(xt, xt, tenToDecPrec);
	
	while(compare(xt, tenToDecPrec)==-1)
	{
		mul(xt, xt, ten);
		k++;
	}
	
	while(compare(xt, tenToDecPrec)>=0)
	{
		div(xt, xt, ten);
		k--;
	}
	
	add(xt, xt, pointFive); // round
	k = -k;
	
	equate(xtInt, xt);
	
	len = decPrec + 35;
	tempStr = (char*)malloc(len*sizeof(char));
	if(!tempStr)
		exit(1);
	outStr = (char*)malloc(len*sizeof(char));
	if(!outStr)
		exit(1);
	
	len = 0;
	while(len<=decPrec && xtInt.n)
	{
		tempStr[len++] = mbShortDiv(xtInt, xtInt, 10) + '0';
	}
	
	tempStr[len] = 0;
	
	// reverse string
	for(i=0;i<len/2;i++)
	{
		ch = tempStr[i];
		tempStr[i] = tempStr[len-1-i];
		tempStr[len-1-i] = ch;
	}
	
	i = len - 1;
	while(tempStr[i]=='0')
		i--;
	
	k+=(len-1-i);
	len = i + 1;
	
	DecExp = k;
	
	tempStr[len] = 0;  // tempStr should contain the integer part
	
	outStr[0] = 0; // forgot this; gave hard to find bug when did strcat for positive number
	
	// copy it to outStr
	if(isNegative)
		strcpy(outStr, "-");

	strcat(outStr, tempStr);
	
	blockPrec = saveBlockPrec;
	free(tempStr);
	
	return outStr;
	
}/* fpToStr2 */


bool checkInputString(const char *inString)
{
	char	*myString, *manString, *expString, ch;
	INT32	length, i, index;
	bool	isExp;
	
	length = strlen(inString);
	if(!length)
		return false;
	
	// check for obvious
	for(i=0;i<length;i++)
	{
		ch = inString[i];
		if(!(ch=='+' || ch=='-' || ch=='e' || ch=='E' || ch=='.' || (ch>='0' && ch<='9')))
			return false;
	}
	
	myString = (char*)malloc((length+1)*sizeof(char));
	
	// remove spaces
	index = 0;
	for(i=0;i<length;i++)
	{
		ch = inString[i];
		if(!(ch==' '))
		{
			myString[index] = ch;
			index++;
		}
			
	}
	myString[index] = 0;
	
	length = strlen(myString);
	if(!length)
	{
		free(myString);
		return false;
	}
	
	// is there exponent?
	isExp = false;
	index = 0;
	for(i=0;i<length;i++)
	{
		index = i;
		ch = myString[i];
		if(ch=='e' || ch=='E')
		{
			isExp = true;
			break;
		}
	}
	
	// if there is exponent, there must be a number in front of it
	if(isExp && index==0)
	{
		free(myString);
		return false;
	}
	
	if(isExp)
	{
		ch = myString[index-1];
		if(!((ch>='0' && ch<='9') || ch=='.'))
		{
			free(myString);
			return false;
		}
		if(ch=='.' && index==1)
		{
			free(myString);
			return false;
		}
		if(ch=='.')  // must be a digit before it; index must be at least 2
		{
			if(index<2)
			{
				free(myString);
				return false;
			}
			
			ch = myString[index-2];
			if(!(ch>='0' && ch<='9'))
			{
				free(myString);
				return false;
			}
		}
	}
	
	manString = (char*)malloc((length+1)*sizeof(char));
	if(isExp)
	{
		for(i=0;i<=index-1;i++)
			manString[i] = myString[i];
		manString[index] = 0;
	}
	else
	{
		for(i=0;i<length;i++)
			manString[i] = myString[i];
		manString[length] = 0;
	}

	// check manString, may have only one sign in front, so it's sign or digit
	length = strlen(manString);
	for(i=1;i<length;i++)
	{
		ch = manString[i];
		if(ch=='+' || ch=='-')
		{
			free(manString);
			free(myString);
			return false;
		}
	}
	
	// can have at most one decimal point
	INT32	count;
	count = 0;
	for(i=1;i<length;i++)
	{
		if(manString[i]=='.')
			count++;
	}
	
	if(count>1)
	{
		free(manString);
		free(myString);
		return false;
	}
	
	free(manString);
	
	if(!isExp)
	{
		free(myString);
		return true;
	}
	
	// check expString
	length = strlen(myString);
	expString = (char*)malloc((length+1)*sizeof(char));
	
	for(i=index+1;i<length;i++)
	{
		expString[i-index-1] = myString[i];
	}
	expString[length-index-1] = 0;
	
	// check expString, may have only one sign in front, so it's sign or digit
	// expString cannot have decimal point
	if(expString[0]=='.')
	{
		free(expString);
		free(myString);
		return false;
	}
	
	length = strlen(expString);
	for(i=1;i<length;i++)
	{
		ch = expString[i];
		if(ch=='+' || ch=='-' || ch=='e' || ch=='E' || ch=='.')
		{
			free(expString);
			free(myString);
			return false;
		}
	}
	
	free(expString);
	free(myString);
	return true;
	
}/* checkInputString */


bool fpConvFromString(fp& z, const char* inString)
{
	static fp	ten, p;
	static bool	initGood=false;
	INT32		len, i, j, ie, id, decExp;
	char		ch, expString[21], *fpString;
	
	if(!initGood)
	{
		init(p);
		ten.i.n = ten.i.nn = 1;
		ten.i.b = (UINT32*)malloc(ten.i.nn*sizeof(UINT32));
		if(!ten.i.b)
			exit(1);
		ten.i.b[0] = 10;
		ten.e = 0;
		initGood = true;
	}
	
	len = strlen(inString);
	if(!len)
		return false;
	
	if(!checkInputString(inString))
		return false;
	
	fpString = (char*)malloc((len+1)*sizeof(char));
	if(!fpString)
		exit(1);
	strcpy(fpString, inString);
	
	
	decExp = 0;
	
	// find e or E
	for(i=0;i<=len;i++)
	{
		ie = i;
		ch = fpString[i];
		if(ch=='e' || ch=='E')
			break;
	}
	
	i = ie; // being cautious
	
	if(ch)
	{
		// have exponent, put it in expString; i points to e or E
		j = 0;
		ch = fpString[++i];
		while(j<21 && ch)
		{
			expString[j++] = ch;
			ch = fpString[++i];
		}
		
		expString[j] = 0;
		
		
		if(!(sscanf(expString, "%d", &decExp)==1))
		{
			free(fpString);
			return false;
		}
		
		// strip off exponent part
		fpString[ie] = 0;
		len = strlen(fpString);
	}
	
	// does fpString have a decimal point?
	for(i=0;i<=len;i++)
	{
		id = i;
		ch = fpString[i];
		if(ch=='.')
			break;
	}
	
	if(ch)
	{
		// yes, there is a decimal point pointed to by id
		// for each number to its right we decrement decExp
		i = id + 1;
		ch = fpString[i];
		while(ch>='0' && ch<='9')  // only goes once through
		{
			fpString[i-1] = ch;
			ch = fpString[++i];
			decExp--;
		}
		fpString[i-1] = 0;
	}
		
	// now fpString should be without exponent or decimal point, so we get the initial integer part of z
	if(!mbConvFromStr(z.i, fpString))
	{
		free(fpString);
		return false;
	}
	z.e = 0;
	
	blockPrec++;
	// now we need to use the dec exponent
	if(!power(p, ten, decExp))
	{
		free(fpString);
		blockPrec--;
		return false;
	}
	
	mul(z, z, p); // this uses blockPrec!
	
	blockPrec--;
	
	free(fpString);
	
	return true;
	
}/* fpConvFromString */


bool init(fp& x)
{
	init(x.i);
	x.e = 0;
	return true;
	
}/* init */


bool init(fp& x, const char *inString)
{
	init(x);
	
	if(!fpConvFromString(x, inString))
		return false;
	
	fpNormalize(x);
	
	return true;
	
}/* init with string */


bool init(fp& x, INT32 y)
{
	init(x.i, y);
	x.e = 0;
	
	fpNormalize(x);
	
	return true;
	
}/* init with INT32 */

bool init(fp& x, INT64 y)
{
    long64tomb(x.i, y);
    //init(x.i, y);
    x.e = 0;
    
    fpNormalize(x);
    
    return true;
    
}/* init with INT64 */


// if this is missing conversion routine fpToStr is very slow
bool init(fp& x, double y)
{
	char		dataString[100], ch;
	INT32		i, j, len, eIndex, zeroCount;
	
	if(!y)
	{
		x.i.n = x.i.nn = 0;
		x.i.b = 0;
		x.e = 0;
		return true;
	}
	
	sprintf(dataString, "%.17e", y); // will allow up to 17 sig figs
	
	// for say 2, we'll see 2.00000000000000e+00 so we need to remove trailing zeroes
	//find 'e'
	len = strlen(dataString);
	
	for(i=0;i<len;i++)
	{
		if(dataString[i]=='e')
			break;
	}
	eIndex = i;
	
	// count continuous zeroes to left of e
	
	zeroCount = 0;
	ch = dataString[--i];
	while(ch=='0' && i>0)
	{
		zeroCount++;
		ch = dataString[--i];
	}
	// i now points to first non '0' as come in from right
	
	// move the exp stuff down to just after where i points, that's a move of eIndex-(i+1)
	
	for(j=eIndex;j<=len;j++)
		dataString[i+1+j-eIndex] = dataString[j];
	
	return init(x, dataString);
	
}/* init */

void long32tofp(fp&z, long x)
{
    long32tomb(z.i, x);
    z.e = 0;
    
}/* long32tofp */

void long64tofp(fp&z, long x)
{
    long64tomb(z.i, x);
    z.e = 0;
    
}/* long32tofp */

// x should be log(number)
// this returns a string with that number, even if it wouldn't fit into an fp type
// note that if the log to base 10 is an integer, the mantissa of 10^theLog is 1
// as the precision goes up the log10 develops a decimal part which then grow as
// (decPrec - length of exponent) as DecPrec increases from length of exponent
// so long as zFrac is zero we don't have to adjust output
// when zFrac>0 we must adjust the length of theMantissa
char* GetNumberFromLog(const fp& x)
{
	char	*outString, *theMantissa, *theExponent;
	fp		z, ten, zExponent, zFrac, zMantissa;
	INT32	expLength, mantissaLength, saveDecPrec, newDecPrec;
	
	if(!x.i.n)
	{
		outString = (char*)malloc(2*sizeof(char));
		outString[0] = '1';
		outString[1] = 0;
		return outString;
	}
	
	ten = 10.;
	
	z = x/log(ten); // log to base 10
	
	floor(zExponent, z);
	
	zFrac = z - zExponent;
	
	zMantissa = pow(ten, zFrac);
	
	saveDecPrec = getDecPrec();
	
	theMantissa = fpToStr(zMantissa, saveDecPrec);
	
	theExponent = fpToStr(zExponent, saveDecPrec);
	
	mantissaLength = strlen(theMantissa);
	expLength = strlen(theExponent);
	
	if(!zFrac.i.n)
	{
		// we have an incomplete exponent so it's 1 x 10^exponent
		outString = (char*)malloc((mantissaLength+expLength+10)*sizeof(char));
		strcpy(outString, theMantissa);
		strcat(outString, "e");
		strcat(outString, theExponent);
		return outString;
	}
	
	// now zFrac>0 and we have a complete exponent, so we need to investigate blockPrec
	
	newDecPrec = 8*getBlockPrec();
	
	newDecPrec = newDecPrec - expLength;
	if(zExponent<0)
		newDecPrec++;
	
	if(newDecPrec<2)
	{
		outString = (char*)malloc((expLength+10)*sizeof(char));
		strcpy(outString, "1");
		strcat(outString, "e");
		strcat(outString, theExponent);
		return outString;
	}
	
	// now newDecPrec>=2 so redo theMantissa
	
	theMantissa = fpToStr(zMantissa, newDecPrec);
	mantissaLength = strlen(theMantissa);
	
	outString = (char*)malloc((mantissaLength+expLength+10)*sizeof(char));
	
	strcpy(outString, theMantissa);
	strcat(outString, "e");
	strcat(outString, theExponent);
	
	return outString;
	
}/* GetNumberFromLog */

