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


#include "fp.h"

INT32	blockPrec=5;

using namespace std;


void setBlockPrec(INT32 prec)
{
	blockPrec = prec;
	
}/* SetBlockPrec */


UINT32 getBlockPrec()
{
	return blockPrec;
	
}/* getBlockPrec */


// don't use x=abs(x) with a static variable, it will cause a bus error.
fp abs(const fp& x)
{
	fp				xt;
	
	equate(xt, x);
	if(xt.i.n<0)
		xt.i.n = -xt.i.n;
	return xt;
	
}/* abs */


// this removes leading and trailing zeroblocks from the integer part and then if has more blocks than blockPrec, it trims down
// we need to make sure that the highest order block is as large as possible by using shiftLeft
void fpNormalize(fp& x)
{
	static fp			one;  // one(1.) gives recursive error
	static UINT64		mask1; // 0x100000000
	static UINT64		mask2; // 0xFFFFFFFF
	static bool			initGood=false;
	INT32				numBlocks, numZeroBlocks, numZeroBits, bitsShift;
	bool				isAllFs, isNegative;
	INT32				i;
	UINT64				cx, carry;
	
	if(!initGood)
	{
		one.i.n = one.i.nn = 1;
		one.e = 0;
		one.i.b = (UINT32*)malloc(one.i.nn*sizeof(UINT32));
		if(!one.i.b)
			exit(1);
		one.i.b[0] = 1;
		mask1 = 1;
		mask1 = (mask1<<32);
		mask2 = mask1 - 1;
		initGood = true;
	}
	
	if(x.i.n==0)
	{
		x.e = 0;
		return;
	}
	
	// remove leading zero blocks
	if(x.i.n<0)
	{
		isNegative = true;
		x.i.n = -x.i.n;
	}
	else
		isNegative = false;
	
	while(x.i.n>0 && x.i.b[x.i.n-1]==0)
		x.i.n--;
	
	if(!x.i.n)
		update(x.i);
	
	if(x.i.n==0)
	{
		x.e = 0;
		return;
	}
	
	if(x.i.n<=blockPrec)
	{
		if(isNegative)
			x.i.n = -x.i.n;
		return;
	}
	
	// shift left to get largest high order block

	numZeroBits = alignLeft(x.i);
	x.e-=numZeroBits;
	
	// remove trailing zero blocks
	numZeroBlocks = 0;
	for(i=0;i<x.i.n;i++)
	{
		if(!x.i.b[i])
			numZeroBlocks++;
		else
			break;
	}
	
	if(numZeroBlocks)
	{
		bitsShift = blockBits*numZeroBlocks;
		mbShiftRight(x.i, x.i, bitsShift);
		x.e+=bitsShift;
	}
	
	if(x.i.n==0)
	{
		x.e = 0;
		return;
	}
	
	if(x.i.n<=blockPrec)
	{
		if(isNegative)
			x.i.n = -x.i.n;
		return;
	}
	
	// now x.i.n>blockPrec so we need to round up and trim
	// but have to worry about the leading blocks some or all
	// being FFFFFFFF.
	// the lowest block just below the blockPrec ones has index x.i.n-blockPrec-1
	// if its value is >0x7FFFFFFF we need to round, else just trim
	cx = x.i.b[x.i.n-blockPrec-1];
	if(cx>0x7FFFFFFF)
	{
		// special case of all blockPrec values being FFFFFFFF
		isAllFs = true;
		for(i=x.i.n-1;i>=x.i.n-blockPrec;i--)
			if(x.i.b[i]!=mask2)
			{
				isAllFs = false;
				break;
			}
		if(isAllFs)
		{
			x.e+=(blockBits*(x.i.n));
			x.i.n = 1;
			x.i.b[0] = 1;
			if(isNegative)
				x.i.n = -x.i.n;
			return;
		}
		// now must add 1 to blockPrec part
		carry = 1;
		for(i=x.i.n-blockPrec;i<x.i.n;i++)
		{
			cx = x.i.b[i];
			cx = cx + carry;
			carry = ((cx & mask1)>>32);
			x.i.b[i] = (cx & mask2);
			if(!carry)
				break;
		}
		// trim trailing zeroes
		numZeroBlocks = 0;
		for(i=0;i<x.i.n;i++)
		{
			if(!x.i.b[i])
				numZeroBlocks++;
			else
				break;
		}
		
		if(numZeroBlocks)
		{
			bitsShift = blockBits*numZeroBlocks;
			mbShiftRight(x.i, x.i, bitsShift);
			x.e+=bitsShift;
		}
		
		if(x.i.n==0)
		{
			x.e = 0;
			return;
		}
	}
	
	// now trim if necessary
	numBlocks = x.i.n - blockPrec;
	if(numBlocks>0)
	{
		bitsShift = blockBits*numBlocks;
		mbShiftRight(x.i, x.i, bitsShift);
		x.e+=bitsShift;
	}
	
	if(isNegative)
		x.i.n = -x.i.n;
	
}/* fpNormalize */


// for a INT32, counts the number of bits after and including the leading 1 bit
INT32 NumBits(INT32 x)
{
	static UINT32	mask;
	static bool		initGood=false;
	INT32			theBits;
	
	if(!initGood)
	{
		mask = 1;
		mask = (mask<<31); // 0x80000000 (which we can't use directly since it's a negative INT32!
		initGood = true;
	}
	
	theBits = 31;
	
	x = abs(x);
	
	x = (x<<1); // gets rid of sign bit
	
	while(!(x & mask) && theBits)
	{
		x = (x<<1);
		theBits--;
	}
	
	return theBits;
	
}/* NumBits */


// for a UINT32, counts the number of bits after and including the leading 1 bit
INT32 NumBits(UINT32 x)
{
	static UINT32	mask;
	static bool		initGood=false;
	INT32			theBits;
	
	if(!initGood)
	{
		mask = 1;
		mask = (mask<<31); // 0x80000000 (which we can't use directly since it's a negative INT32!
		initGood = true;
	}
	
	theBits = 32;
	
	while(!(x & mask) && theBits)
	{
		x = (x<<1);
		theBits--;
	}
	
	return theBits;
	
}/* NumBits for UINT32 */


// compares abs(x) to abs(y); returns:
// +1 if abs(x)>abs(y)
// 0  if abs(x)=abs(y)
// -1 if abs(x)<abs(y)
static INT32 fpCompareAbs(const fp& x, const fp& y)
{
	static fp		xt, yt;
	INT32			xPowerOfTwo, yPowerOfTwo;
	INT32			i, bitShift;
	
	if(!x.i.n && !y.i.n)
		return 0;
	
	if(!y.i.n)
		return 1;
	
	if(!x.i.n)
		return -1;
	
	equateWithCare(xt, x);
	equateWithCare(yt, y);
	xt.i.n = abs(xt.i.n);
	yt.i.n = abs(yt.i.n);
	
	fpNormalize(xt);
	fpNormalize(yt);
	// most sig blocks are not zero
	
	// coarse compare first
	xPowerOfTwo = xt.e + blockBits*(xt.i.n-1) + NumBits(xt.i.b[xt.i.n-1]);
	yPowerOfTwo = yt.e + blockBits*(yt.i.n-1) + NumBits(yt.i.b[yt.i.n-1]);
	
	if(xPowerOfTwo>yPowerOfTwo)
		return 1;
	
	if(xPowerOfTwo<yPowerOfTwo)
		return -1;
	
	// now xPowerOfTwo = yPowerOfTwo
	// we need the most sig bit of both xt and yt to be highest order bit of highest order block
	// so we can compare block by block from highest to lowest
														
	bitShift = blockBits - NumBits(xt.i.b[xt.i.n-1]);
	mbShiftLeft(xt.i, xt.i, bitShift);
	xt.e-=bitShift; // not needed
	
	bitShift = blockBits - NumBits(yt.i.b[yt.i.n-1]);
	mbShiftLeft(yt.i, yt.i, bitShift);
	yt.e-=bitShift; // not needed
	
	if(xt.i.n>=yt.i.n)
	{
		for(i=xt.i.n-1;i>=xt.i.n-yt.i.n;i--)
		{
			if(xt.i.b[i]>yt.i.b[yt.i.n-xt.i.n+i])
				return 1;
			else
				if(xt.i.b[i]<yt.i.b[yt.i.n-xt.i.n+i])
					return -1;
		}
		
		for(i=xt.i.n-yt.i.n-1;i>=0;i--)
		{
			if(xt.i.b[i]>0)
				return 1;
		}
		
		return 0;
	}
	else
	{
		for(i=yt.i.n-1;i>=yt.i.n-xt.i.n;i--)
		{
			if(yt.i.b[i]>xt.i.b[xt.i.n-yt.i.n+i])
				return -1;
			else
				if(yt.i.b[i]<xt.i.b[xt.i.n-yt.i.n+i])
					return 1;
		}
		
		for(i=yt.i.n-xt.i.n-1;i>=0;i--)
		{
			if(yt.i.b[i]>0)
				return -1;
		}
		
		return 0;
	}
	
	
}/* fpCompareAbs */


/*
	assumes x an y are normalized
	returns:
	1 if x>y
	0 if x=y
	-1 if x<y
*/
INT32 compare(const fp& x, const fp& y)
{
	INT32	xSign, ySign;
	INT32	result;
	
	if(x.i.n>=0)
		xSign = 1;
	else
		xSign = -1;
	
	if(y.i.n>=0)
		ySign = 1;
	else
		ySign = -1;
	
	result = fpCompareAbs(x, y);
	
	if(!result)
	{
		if(xSign==1 && ySign==-1)
			return 1;
	
		if(xSign==-1 && ySign==1)
			return -1;
		
		return 0;
	}
	
	if(result==1)
	{
		if(xSign==1)
			return 1;
		else
			return -1;
	}
	
	// now result = -1
	if(ySign==1)
		return -1;
	else
		return 1;
	
}/* compare */



// z = x + y
void add(fp& z, const fp& x, const fp& y)
{
	static fp		xt, yt, zt;
	INT32			eDiff, blocksShift, bitsLeft, bitsRight, restBits;
	INT32			s;
	INT64			lg2x, lg2y;
	
	if(!x.i.n)
	{
		equate(z, y);
		return;
	}
	
	if(!y.i.n)
	{
		equate(z, x);
		return;
	}
	
	equateWithCare(xt.i, x.i);
	xt.e = x.e;
	
	equateWithCare(yt.i, y.i);
	yt.e = y.e;
	
	
	// obvious cases
	if(xt==yt)
	{
		xt.e++; // multiply by 2
		equate(z, xt);
		fpNormalize(z);
		return;
	}
	
	if(xt==-yt)
	{
		z = 0;
		return;
	}
	
	lg2x = Lg2(xt);
	lg2y = Lg2(yt);
	
	if(abs(lg2x-lg2y)>((blockPrec+1)*blockBits))
	{
		if(lg2x>lg2y)
		{
			z = xt;
			fpNormalize(z);
			return;
		}
		else
		{
			z = yt;
			fpNormalize(z);
			return;
		}
	}
	
	
	if(xt.e==yt.e)
	{
		add(zt.i, xt.i, yt.i);
		zt.e = xt.e;
		equate(z, zt);
		fpNormalize(z);
		return;
	}
	
	if(xt.e>yt.e)
	{
		eDiff = xt.e - yt.e;
		blocksShift = eDiff/blockBits;
		restBits = eDiff - blockBits*blocksShift;
		mbShiftLeft(xt.i, xt.i, restBits);
			xt.e-=restBits;
	}
	else
	{
		eDiff = yt.e - xt.e;
		blocksShift = eDiff/blockBits;
		restBits = eDiff - blockBits*blocksShift;
		mbShiftLeft(yt.i, yt.i, restBits);
			yt.e-=restBits;
	}
	
	if(xt.e==yt.e)
	{
		add(zt.i, xt.i, yt.i);
		zt.e = xt.e;
		equate(z, zt);
		fpNormalize(z);
		return;
	}
	
	// the exponents now differ by an integral multiple of blockBits
	// For unequal exponents we must always decrease the higher exponent down to match the lower
	// But worry about exceeding blockPrec
	if(xt.e>yt.e)
	{
		eDiff = xt.e - yt.e;
		
		blocksShift = eDiff/blockBits;
		
		if(abs(xt.i.n)>=blockPrec)
		{
			bitsRight = blockBits*blocksShift;
			mbShiftRight(yt.i, yt.i, bitsRight);
			yt.e+=bitsRight;
		}
		else
		{
			s = blockPrec - abs(xt.i.n);
			if(s>=blocksShift)
			{
				bitsLeft = blockBits*blocksShift;
				mbShiftLeft(xt.i, xt.i, bitsLeft);
				xt.e-=bitsLeft;
			}
			else
			{
				bitsLeft = blockBits*s;
				bitsRight = blockBits*(blocksShift - s);
				mbShiftLeft(xt.i, xt.i, bitsLeft);
				xt.e-=bitsLeft;
				mbShiftRight(yt.i, yt.i, bitsRight);
				yt.e+=bitsRight;
			}
		}
	}
	else
	{
		eDiff = yt.e - xt.e;
		
		blocksShift = eDiff/blockBits;
		
		if(abs(yt.i.n)>=blockPrec)
		{
			bitsRight = blockBits*blocksShift;
			mbShiftRight(xt.i, xt.i, bitsRight);
			xt.e+=bitsRight;
		}
		else
		{
			s = blockPrec - abs(yt.i.n);
			if(s>=blocksShift)
			{
				bitsLeft = blockBits*blocksShift;
				mbShiftLeft(yt.i, yt.i, bitsLeft);
				yt.e-=bitsLeft;
			}
			else
			{
				bitsLeft = blockBits*s;
				bitsRight = blockBits*(blocksShift - s);
				mbShiftLeft(yt.i, yt.i, bitsLeft);
				yt.e-=bitsLeft;
				mbShiftRight(xt.i, xt.i, bitsRight);
				xt.e+=bitsRight;
			}
		}
	}
	
	// now can add integer parts
	add(zt.i, xt.i, yt.i);
	zt.e = xt.e;
	
	equate(z, zt);
	
	fpNormalize(z);
	
}/* add */


void add(fp& z, const fp& x, const mb& y)
{
	static fp		yt;
	
	equate(yt, y);
	
	add(z, x, yt);
	
}/* add */


void add(fp& z, const mb& x, const fp& y)
{
	static fp		xt;
	
	equate(xt, x);
	
	add(z, xt, y);
	
}/* add */


void add(fp& z, const fp& x, double y)
{
	static fp		yt;
	
	equateWithCare(yt, y);
	
	add(z, x, yt);
	
}/* add */


void add(fp& z, double x, const fp& y)
{
	static fp		xt;
	
	equateWithCare(xt, x);
	
	add(z, xt, y);
	
}/* add */


// z = x - y
void sub(fp& z, const fp& x, const fp& y)
{
	static fp		xt, yt, zt;
	INT32			eDiff, blocksShift, bitsLeft, bitsRight, restBits;
	INT32			s;
	INT64			lg2x, lg2y;
	
	if(!x.i.n)
	{
		equate(z, y);
		z.i.n = -z.i.n;
		return;
	}
	
	if(!y.i.n)
	{
		equate(z, x);
		return;
	}
	
	equateWithCare(xt.i, x.i);
	xt.e = x.e;
	
	equateWithCare(yt.i, y.i);
	yt.e = y.e;
	
	// obvious cases
	if(yt==-xt)
	{
		xt.e++; // multiply by 2
		equate(z, xt);
		fpNormalize(z);
		return;
	}
	
	if(xt==yt)
	{
		z = 0;
		return;
	}
	
	lg2x = Lg2(xt);
	lg2y = Lg2(yt);
	
	if(abs(lg2x-lg2y)>((blockPrec+1)*blockBits))
	{
		if(lg2x>lg2y)
		{
			z = xt;
			fpNormalize(z);
			return;
		}
		else
		{
			z = -yt;
			fpNormalize(z);
			return;
		}
	}
	
	
	if(xt.e==yt.e)
	{
		sub(zt.i, xt.i, yt.i);
		zt.e = xt.e;
		equate(z, zt);
		fpNormalize(z);
		return;
	}
	
	if(xt.e>yt.e)
	{
		eDiff = xt.e - yt.e;
		blocksShift = eDiff/blockBits;
		restBits = eDiff - blockBits*blocksShift;
		mbShiftLeft(xt.i, xt.i, restBits);
			xt.e-=restBits;
	}
	else
	{
		eDiff = yt.e - xt.e;
		blocksShift = eDiff/blockBits;
		restBits = eDiff - blockBits*blocksShift;
		mbShiftLeft(yt.i, yt.i, restBits);
			yt.e-=restBits;
	}
	
	if(xt.e==yt.e)
	{
		sub(zt.i, xt.i, yt.i);
		zt.e = xt.e;
		equate(z, zt);
		fpNormalize(z);
		return;
	}
	
	// For unequal exponents we must always decrease the higher exponent down to match the lower
	// But worry about exceeding blockPrec
	// If shifting eDiff left causes > blockPrec, we stop at blockPrec and shift the other right (we may hit zero)
	if(xt.e>yt.e)
	{
		eDiff = xt.e - yt.e;
		
		blocksShift = eDiff/blockBits;
		
		if(abs(xt.i.n)>=blockPrec)
		{
			bitsRight = blockBits*blocksShift;
			mbShiftRight(yt.i, yt.i, bitsRight);
			yt.e+=bitsRight;
		}
		else
		{
			s = blockPrec - abs(xt.i.n);
			if(s>=blocksShift)
			{
				bitsLeft = blockBits*blocksShift;
				mbShiftLeft(xt.i, xt.i, bitsLeft);
				xt.e-=bitsLeft;
			}
			else
			{
				bitsLeft = blockBits*s;
				bitsRight = blockBits*(blocksShift - s);
				mbShiftLeft(xt.i, xt.i, bitsLeft);
				xt.e-=bitsLeft;
				mbShiftRight(yt.i, yt.i, bitsRight);
				yt.e+=bitsRight;
			}
		}
	}
	else
	{
		eDiff = yt.e - xt.e;
		
		blocksShift = eDiff/blockBits;
		
		if(abs(yt.i.n)>=blockPrec)
		{
			bitsRight = blockBits*blocksShift;
			mbShiftRight(xt.i, xt.i, bitsRight);
			xt.e+=bitsRight;
		}
		else
		{
			s = blockPrec - abs(yt.i.n);
			if(s>=blocksShift)
			{
				bitsLeft = blockBits*blocksShift;
				mbShiftLeft(yt.i, yt.i, bitsLeft);
				yt.e-=bitsLeft;
			}
			else
			{
				bitsLeft = blockBits*s;
				bitsRight = blockBits*(blocksShift - s);
				mbShiftLeft(yt.i, yt.i, bitsLeft);
				yt.e-=bitsLeft;
				mbShiftRight(xt.i, xt.i, bitsRight);
				xt.e+=bitsRight;
			}
		}
	}
	
	// now can subtract integer parts
	sub(zt.i, xt.i, yt.i);
	zt.e = xt.e;
	
	equate(z, zt);
	
	fpNormalize(z);
	
}/* sub */


void sub(fp& z, const fp& x, const mb& y)
{
	static fp		yt;
	
	equate(yt, y);
	
	sub(z, x, yt);
	
}/* sub */


void sub(fp& z, const mb& x, const fp& y)
{
	static fp		xt;
	
	equate(xt, x);
	
	sub(z, xt, y);
	
}/* sub */


void sub(fp& z, const fp& x, double y)
{
	static fp		yt;
	
	equateWithCare(yt, y);
	
	sub(z, x, yt);
	
}/* sub */


void sub(fp& z, double x, const fp& y)
{
	static fp		xt;
	
	equateWithCare(xt, x);
	
	sub(z, xt, y);
	
}/* sub */


// z = x*y
void mul(fp& z, const fp& x, const fp& y)
{
	if(!x.i.n || !y.i.n)
	{
		z = 0;
		return;
	}
	
	mul(z.i, x.i, y.i);
	z.e = x.e + y.e;
	
	fpNormalize(z);
	
	return;
	
}/* mul */


void mul(fp& z, const fp& x, const mb& y)
{
	static fp		yt;
	
	equate(yt, y);
	
	mul(z, x, yt);
	
}/* mul */


void mul(fp& z, const mb& x, const fp& y)
{
	static fp		xt;
	
	equate(xt, x);
	
	mul(z, xt, y);
	
}/* mul */


void mul(fp& z, const fp& x, double y)
{
	static fp		yt;
	
	equateWithCare(yt, y);
	
	mul(z, x, yt);
	
}/* mul */


void mul(fp& z, double x, const fp& y)
{
	static fp		xt;
	
	equateWithCare(xt, x);
	
	mul(z, xt, y);
	
}/* mul */


// z = x/y
// we have to worry about going to blockPrec for unending divisions
bool div(fp& z, const fp& x, const fp& y)
{
	static fp		s;
	INT32			blocksShift, bitsShift;
	
	if(!y.i.n)
		return false;
	
	if(!x.i.n)
	{
		z = 0;
		return true;
	}
	
	equateWithCare(s, x);
	
	// shift s.i left so new abs(s.i.n) = blockPrec + abs(y.i.n) + 1
	// may have to increase abs(s.i.n), but mbShiftLeft takes care of that
	blocksShift = blockPrec - abs(x.i.n) + abs(y.i.n) + 1;
	
	if(blocksShift>0)
	{
		bitsShift = blockBits*blocksShift;
		mbShiftLeft(s.i, s.i, bitsShift);
		s.e-=bitsShift;
	}
	
	div(z.i, s.i, y.i);
	
	z.e = s.e - y.e;
	
	fpNormalize(z);
	
	return true;
	
}/* div */


void div(fp& z, const fp& x, const mb& y)
{
	static fp		yt;
	
	equate(yt, y);
	
	div(z, x, yt);
	
}/* div */


void div(fp& z, const mb& x, const fp& y)
{
	static fp		xt;
	
	equate(xt, x);
	
	div(z, xt, y);
	
}/* div */


bool div(fp& z, const fp& x, double y)
{
	static fp		yt;
	
	equateWithCare(yt, y);
	
	return div(z, x, yt);
	
}/* div */


bool div(fp& z, double x, const fp& y)
{
	static fp		xt;
	
	equateWithCare(xt, x);
	
	return div(z, xt, y);
	
}/* div */


// z = x^y ; y is a INT32
/*
	Uses Victor Shoup's NTL method.
	If bit(i) is ith bit of y, i = 0, 1, ... , 30 (bit(31) is sign bit)
	then abs(y) = Sum(i=0 to 30) bit(i)*2^i
	so x^y = x^(bit(0)) * x^(bit(1)*2) * x^(bit(2)*4) * x^(bit(3)*8) * x^(bit(4)*16) * ...
	Note that if all bit(i)=1 each term in the product is the square of the previous term
	So we remember that and multiply in if bit(i)=1, but don't multiply in if bit(i)=0
	Then Shoup does that in an ingenious way
	NumBits(y) gives the number of bits after and including the high order 1 bit, not counting the sign bit
*/
bool power(fp& z, const fp& x, INT32 y)
{
	static fp			one;
	fp					zt;
	//static double		lg2fpMax = 2147483647.;
	static UINT32		mask; // 0x80000000
	static bool			initGood=false;
	INT32				i, n, sy;
	
	if(!initGood)
	{
		init(zt);
		one.i.n = one.i.nn = 1;
		one.i.b = (UINT32*)malloc(one.i.n*sizeof(UINT32));
		if(!one.i.b)
			exit(1);
		one.i.b[0] = 1;
		one.e = 0;
		mask = 1;
		mask = (mask<<31); // 0x80000000 (which we can't use directly since it's a negative INT32!
		initGood = true;
	}
	
	if(!x.i.n)
	{
		z = 0;
		if(y<=0)
			return false;
		
		return true;
	}
	
	if(y==0)
	{
		equate(z, one);
		return true;
	}
	
	//blockPrec++;
	//cout << "abs(y) = " << abs(y) << endl;
	//cout << "in power Lg2(x) = " << Lg2(x) << endl;  // 3 but should be 25.13
	//cout << "(UINT64)abs(y)*(1+Lg2(x)) = " << (UINT64)abs(y)*(1+Lg2(x)) << endl;  // -1609030302
	
	/*
	if((UINT64)abs(y)*(1+Lg2(x))>lg2fpMax)
	{
		z = 0;
		return false;
	}
	 */
	
	equate(zt, one);
	sy = y;
	if(y<0)
		sy = -y;
	n = NumBits(sy);
	
	sy = (sy<<(32-n));  // go to leading 1 bit
	
	for(i=0;i<n;i++)
	{
		mul(zt, zt, zt); // zt^2
		if(sy & mask)
			mul(zt, zt, x);
		sy = (sy<<1);
	}
	
	if(y<0)
		div(zt, one, zt);
	
	equate(z, zt);
	
	//blockPrec--;
	
	//fpNormalize(z);
	
	return true;
	
}/* power */
