/*
 *  fpFuncs.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;


// z = floor(x)
void floor(fp& z, const fp& x)
{
	fp			xt;
	
	equate(xt, x);
	
	if(xt.e>=0) // it already is an integer
	{
		equate(z, xt);
		return;
	}
	
	// now xt.e<0
	mbShiftRight(xt.i, xt.i, abs(xt.e));
	xt.e = 0;
	
	if(xt.i.n<0)
		xt = xt - 1;
	
	if(xt.i.n==0 && x.i.n<0)
		xt = xt - 1;
	
	equate(z, xt);
	
}/* floor */


fp floor(const fp& x)
{
	fp		z;
	
	floor(z, x);
	return z;
	
}/* floor */


// z = ceil(x)
void ceil(fp& z, const fp& x)
{
	fp			xt;
	
	equate(xt, x);
	
	if(xt.e>=0) // it already is an integer
	{
		equate(z, xt);
		return;
	}
	
	// now xt.e<0
	mbShiftRight(xt.i, xt.i, abs(xt.e));
	xt.e = 0;
	
	if(xt.i.n>0)
		xt = xt + 1;
	
	if(xt.i.n==0 && x.i.n>0)
		xt = xt + 1;
	
	equate(z, xt);
	
}/* ceil */


fp ceil(const fp& x)
{
	fp		z;
	
	ceil(z, x);
	return z;
	
}/* ceil */


bool isItInteger(const fp& x)
{
	fp			xt;
	INT32		i;
	bool		isInteger;
	
	equate(xt, x);
	fpNormalize(xt);
	
	if(xt.e>=0)
		return true;
	
	// now xt.e<0 so might or might not have an integer
	
	if(xt.e<=-32)
		return false; // hex point is to left of block zero, which is not zero
		
	// now -31 <= xt.e <= -1
	// so using x.e block zero is of form ...xxx.xxxxx... say, where x is a binary bit
	// the binary point is abs(x.e) from the right
	// so we need to check the bits to the right of the binary point
	// if all are zero, we have an integer; if not, we don't
	
	isInteger = true;
	for(i=0;i<abs(xt.e);i++)
	{
		if(xt.i.b[0] & 1)
		{
			isInteger = false;
			break;
		}
		xt.i.b[0] = (xt.i.b[0]>>1);
	}
	
	return isInteger;
	
}/* isItInteger */


// assumes x is an integer disguised as an fp
bool isItEven(const fp& x)
{
	mb			xt;
		
	if(x.e>0)
		return true;
	
	equate(xt, x.i);
	
	if(xt.b[0] & 1)
		return false;
	
	return true;
	
}/* isItEven */


// converts x to a INT32; it assumes x is in proper range
INT32 to_long(const fp& x)
{
	fp					sx;
	INT32				z;
	
	equate(sx, x);
	
	if(sx.e>=0)
	{
		mbShiftLeft(sx.i, sx.i, sx.e);
	}
	else
	{
		mbShiftRight(sx.i, sx.i, abs(x.e));
	}
	
	if(!sx.i.n)
		return 0;
	
	z = sx.i.b[0];
	if(x.i.n<0)
		z = -z;
	
	return z;
	
}/* to_long */


// z = x^2
void sqr(fp& z, const fp& x)
{
	mul(z, x, x);
	
}/* sqr */


fp sqr(const fp& x)
{
	fp		z;
	
	sqr(z, x);
	return z;
	
}/* sqr */


// z = sqrt(x) must be calculated accurate to blockPrec
void sqrt(fp& z, const fp& x)
{
	fp				s;
	INT32			numBits;
	
	if(x.i.n<0)
		return;
	
	if(!x.i.n)
	{
		z = 0;
		return;
	}
	
	blockPrec*=2;
	
	// initialize s
	equate(s, x);
	
	// bring the number of blocks of s to at least blockPrec
	if(s.i.n<blockPrec)
	{
		numBits = (blockPrec - s.i.n)*blockBits;
		// shift left by numBits
		mbShiftLeft(s.i, s.i, numBits);
		s.e-=numBits;
	}
	
	if(2*(s.e/2)!=s.e)
	{
		// s.e is odd, make it even
		mbShiftLeft(s.i, s.i, 1);
		s.e--;
	}
	
	s.e/=2;
	
	sqrt(s.i, s.i);
		
	blockPrec/=2;
	
	equate(z, s);
	
	return;
	
}/* sqrt */


fp sqrt(const fp& x)
{
	fp		z;
	
	sqrt(z, x);
	return z;
	
}/* sqrt */


// z = x^y
bool pow(fp& z, const fp& x, const fp& y)
{
	static fp			one;
	fp					xt, yt, zt;
	static bool			initGood=false;
	bool				isNegative;
	INT32				yInt;
	
	
	if(!initGood)
	{
		init(one, 1.);
		//init(lg2fpMax, "2147483647.");
		initGood = true;
	}
	
	if(!y.i.n)
	{
		z = one; // so 0^0 = 1
		return true;
	}
	
	if(!x.i.n)
	{
		z.i.n = 0; // the caller has to check for y.i.n<0 for this case
		return true;
	}
	
	if(!compare(x, one))
	{
		equate(z, one);
		return true;
	}
	
	if(equate(yInt, y))
	{
		// y is an integer
		/*
		if(abs(y)*(1+Lg2(x))>lg2fpMax)
		{
			z = 0;
			return false;
		}
		 */
		pow(z, x, yInt);
		return true;
	}
	
	equate(xt, x);
	
	if(xt.i.n<0)
	{
		isNegative = true;
		xt.i.n = -xt.i.n;
	}
	else
		isNegative = false;
	
	if(!isItInteger(y) && isNegative)
	{
		z = 0;
		return false;
	}
	
	if(!y.i.n)
	{
		equate(z, one);
		return true;
	}
	
	/*
	if(abs(y)*(1+Lg2(x))>lg2fpMax)
	{
		z = 0;
		return false;
	}
	 */
	
		
	if(isItInteger(y))
	{
		// it is up to caller to check that abs(y.i.n)=1
		// and that there will not be an overflow
		power(z, x, to_long(y));
		return true;
	}
	
	blockPrec++;
	
	log(zt, xt);
	mul(zt, zt, y);
	exp(zt, zt);
	
	equate(z, zt);
	
	blockPrec--;
	
	fpNormalize(z);
		
	return true;
	
}/* pow */


fp pow(const fp& x, const fp& y)
{
	fp		z;
	
	pow(z, x, y);
		
	return z;
	
}/* pow */


bool pow(fp& z, const fp& x, double y)
{
	fp			yt;
		
	equate(yt, y);
	
	return pow(z, x, yt);
	
}/* pow */


fp pow(const fp& x, double y)
{
	fp		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


bool pow(fp& z, double x, const fp& y)
{
	fp			xt;
	
	equate(xt, x);
	
	return pow(z, xt, y);
	
}/* pow */


fp pow(double x, const fp& y)
{
	fp		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


bool pow(fp& z, const fp& x, INT32 y)
{
	INT32				i, n, sy;
	static double		lg2fpMax = 2147483647.;
	static UINT32		mask; // 0x80000000
	static fp			one;
	fp					zt;
	static bool			initGood=false;
	bool				isNegative;
	
	if(!initGood)
	{
		one = 1.;
		mask = 1;
		mask = (mask<<31); // 0x80000000
		initGood = true;
	}
	
	if(x==one)
	{
		z = one;
		return true;
	}
	
	if(!x)
	{
		z = 0.;
		return true;
	}
	
	if(!y)
	{
		z = one;
		return true;
	}
	
	if(y>0)
	{
		isNegative = false;
		sy = y;
	}
	else
	{
		isNegative = true;
		sy = -y;
	}
	
	if(abs(y)*Lg2(x)>lg2fpMax)
	{
		z = 0;
		return false;
	}
	
	zt = one;
	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);
	}
	
	z = zt;
	if(isNegative)
		z = 1/z;
	
	return true;
	
}/* pow */


fp pow(const fp& x, INT32 y)
{
	fp			z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


void ComputeE(fp& z)
{
	static fp			one, e;
	fp					s, s1, t;
	static INT32		blockPrecCurr;
	static bool			initGood=false;
	INT32				i;
	
	if(!initGood)
	{
		init(one, 1.);
		blockPrecCurr = 0;
		initGood = true;
	}
	
	if(blockPrec<=blockPrecCurr)
	{
		equate(z, e);
		fpNormalize(z);
		return;
	}
	
	blockPrecCurr = blockPrec;
	
	equate(s, one);
	equate(t, one);
	
	for(i=2;;i++)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		div(t, t, i);
	}
	
	equate(e, s);
	
	equate(z, s);
	
}/* ComputeE */


// uses Bailey-Borwein-Plouffe algorithm
void ComputePi(fp& z)
{
	static fp				one, s47(47.), s15(15.), oneSixteenth, pi;
	fp						s, s1, t, t1, t2;
	static INT32			blockPrecCurr;
	static bool				initGood=false;
	UINT32					i;
	
	if(!initGood)  // the calculation of pi is much slower if we don't init like this
	{
		init(one, 1.);
		init(oneSixteenth, 1.);
		oneSixteenth.e = -4;   // if use oneSixteenth(.0625) get a 5X slowdown!
		blockPrecCurr = 0;
		initGood = true;
	}
	
	if(blockPrec<=blockPrecCurr)
	{
		equate(z, pi);
		fpNormalize(z);
		return;
	}
	
	blockPrecCurr = blockPrec;
	
	mul(s, one, s47);
	div(s, s, s15); // i = 0 term
	
	equate(s1, s);
	
	equate(t, one); // powerFactor
	
	for(i=1; ;i++)
	{
		mul(t, t, oneSixteenth);
		mul(t1, one, 4);
		div(t1, t1, 8*i+1);
		mul(t2, one, 2);
		div(t2, t2, 8*i+4);
		sub(t1, t1, t2);
		div(t2, one, 8*i+5);
		sub(t1, t1, t2);
		div(t2, one, 8*i+6);
		sub(t1,t1,t2);
		mul(t1, t1, t);
		add(s, s, t1);
		if(!compare(s, s1))
			break;
		equate(s1, s);
	}
	
	equate(pi, s);
	
	equate(z, s);
	
}/* ComputePi */


fp Pi()
{
	fp		z;
	
	ComputePi(z);
	return z;
	
}/* Pi */


// exp emulates Victor Shoup's NTL method
void exp(fp& z, const fp& x)
{
	static fp		one;
	fp				f, nn, e, t, t1, t2, s, s1;
	static bool		initGood=false;
	INT32			i, n;
	
	if(!initGood)
	{
		init(one, 1.);
		initGood = true;
	}
	
	if(abs(Lg2(x))>.69*0x7FFFFFFF)
		return; //avoids overflow
		
	floor(nn, x);
	
	sub(f, x, nn);

	// convert nn to INT32
	n = to_long(nn);
	
	// step 1: calculate t1 = e^n by repeated squaring
	
	ComputeE(e);
	
	power(t1, e, n);
	
	// step 2: calculate t2 = e^f using Taylor series expansion
	
	s.i.n = 0;
	s.e = 0;
	equate(t, one);
	
	for(i=1; ;i++)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t, t, f);
		div(t, t, i);
	}
	
	equate(t2, s);
	
	mul(s, t1, t2);
	
	equate(z, s);
	
}/* exp */


fp exp(const fp& x)
{
	fp		z;
	
	exp(z, x);
	return z;
	
}/* exp */


// returns approximate log(abs(x)) to base 2
INT32 Lg2(const fp& x)
{
	static fp		xt;
	static UINT32	mask;
	static bool		initGood=false;
	INT32			i, xPowerOfTwo, numZeroBits;
	UINT32			cs;
	
	if(!initGood)
	{
		init(xt);
		mask = 1;
		mask = (mask<<31); // 0x80000000 (which we can't use directly since it's a negative INT32!
		initGood = true;
	}
	
	if(!x.i.n)
		return 0; // bad input
	
	equate(xt, x);
	xt.i.n = abs(xt.i.n);
	
	xPowerOfTwo = xt.e + blockBits*(xt.i.n-1);
	
	// look at high blocks
	numZeroBits = 0;
	cs = xt.i.b[xt.i.n-1];
	for(i=0;i<blockBits;i++)
	{
		if(!(cs & mask))
		{
			numZeroBits++;
			cs = (cs<<1);
		}
		else
			break;
	}
	
	return xPowerOfTwo + blockBits - numZeroBits - 1;
	
}/* Lg2 */


// z = ln(2)
void ComputeLn2(fp& z)
{
	static fp		pointFive, Ln2;
	fp				s, s1, t, t1;
	static INT32	blockPrecCurr;
	static bool		initGood=false;
	INT32			i;
	
	if(!initGood)
	{
		blockPrecCurr = 0;
		init(pointFive, 1.);
		pointFive.e = -1;
		initGood = true;
	}
	
	if(blockPrec<=blockPrecCurr)
	{
		equate(z, Ln2);
		fpNormalize(z);
		return;
	}
	
	blockPrecCurr = blockPrec;
	
	equate(t, pointFive);
	equate(t1, pointFive);
	s.i.n = 0; // this was initially missing!
	s.e = 0;
	
	for (i=2; ; i++)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t1, t1, pointFive);
		div(t, t1, i);
	}
	
	equate(z, s);
	
	equate(Ln2, s);
	
}/* ComputeLn2 */


// a fast z = 2^n
void power2(fp& z, INT32 n)
{
	static fp		one;
	static bool		initGood=false;
	
	if(!initGood)
	{
		init(one, 1.);
		initGood = true;
	}
	
	one.e = n; // one is now 2^n
	
	equate(z, one);
	
}/* power2 */


fp power2(INT32 n)
{
	fp		z;
	
	power2(z, n);
	return z;
	
}/* power2 */


// log emulates Victor Shoup's NTL method
void log(fp& z, const fp& x)
{
	static fp		pointFive, pointSevenFive, onePointFive, one, two, y, s, s1, t, t1;
	static bool		initGood=false;
	INT32			i, n, result;
	
	if(!initGood)
	{
		init(one, 1.);
		init(two, 1.);
		two.e = 1; // should be faster in use than init(two, 2.)
		init(pointFive, 1.);
		pointFive.e = -1;
		init(onePointFive, 1.5);
		init(pointSevenFive, .75);
		// the above are all exact, so blockPrec doesn't apply
		
		init(y);
		init(s);
		init(s1);
		init(t);
		init(t1);
		initGood = true;
	}
	
	if(x.i.n<=0)
		return;
	
	// re-write x = 2^n * (1 - y), where -1/2 < y < 1/4  (so 3/4 < 1-y < 3/2)
	result = compare(x, pointSevenFive);
	
	if((result==1 || result==0) && compare(x, onePointFive)==-1)
	{
		n = 0;
		sub(y, one, x);
	}
	else
	{
		n = Lg2(x);
		power2(t, -n);
		mul(t, t, x);
		
		while(compare(t, onePointFive)==1)
		{
			mul(t, t, pointFive);
			n++;
		}
		
		sub(y, one, t);
	}
	
	// compute s = - ln(1-y) by power series expansion
	
	s.i.n = 0;
	s.e = 0;
	
	equate(t, y);
	equate(t1, y);
	
	for(i=2; ;i++)
	{	
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t1, t1, y);
		div(t, t1, i);
	}
	
	if(n==0)
	{
		t.i.n = 0;
		t.e = 0;
	}
	else
	{
		ComputeLn2(t);
		mul(t, t, abs(n));
		if(n<0)
			t.i.n = -t.i.n;
	}
	
	sub(z, t, s);
	
}/* log */


fp log(const fp& x)
{
	fp		z;
	
	log(z, x);
	return z;
	
}/* log */


void ComputeLog_10(fp& z)
{
	static fp			ten, LogOfTen;
	static INT32		blockPrecCurr;
	static bool			initGood=false;
	
	if(!initGood)
	{
		init(LogOfTen);
		blockPrecCurr = 0;
		init(ten, 10.);
		initGood = true;
	}
	
	if(blockPrec<=blockPrecCurr)
	{
		equate(z, LogOfTen);
		fpNormalize(z);
		return;
	}
	
	blockPrecCurr = blockPrec;
	
	log(z, ten);
	
	equate(LogOfTen, z);
	
}/* ComputeLog_10 */


// sin emulates Victor Shoup's NTL method with big changes
void sin(fp& z, const fp& x)
{
	static fp		one, three, pointFive;
	fp				f, n, t, t1, t2, s, s1, x2;
	static bool		initGood=false;
	INT32			i;
	
	
	if(!initGood)
	{
		init(one, 1.);
		init(three, 3.);
		init(pointFive, 1.);
		pointFive.e = -1;
		initGood = true;
	}
	
	//if(abs(Lg2(x))>1000)
	//	return; //avoids overflow
		
	if(!x.i.n)
	{
		z = 0;
		return;
	}
	
	mul(x2, x, x); // x*x
	
	if(compare(x2, three)==-1)
	{
		equate(f, x);
	}
	else
	{
		div(t1, x, Pi());
		floor(n, t1);
		sub(f, t1, n);
		if(compare(f, pointFive)==1)
		{
			add(n, n, one);
			sub(f, t1, n);
		}
		
		mul(f, Pi(), f);
		
		if(n.i.n && !isItEven(n))
		{
			f.i.n = -f.i.n;
		}
	}
	
	equate(t, f);
	s.i.n = 0;
	s.e = 0;
	fp		fSqr;
	fSqr = f*f;
	for(i=3; ;i+=2)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t, t, fSqr);
		div(t, t, i*(i-1));
		t.i.n = -t.i.n;
	}
	
	equate(z, s);
	
	return;
	
}/* sin */


fp sin(const fp& x)
{
	fp		z;
	
	sin(z, x);
	return z;
	
}/* sin */


// cos emulates Victor Shoup's NTL method
void cos(fp& z, const fp& x)
{
	static fp		pointFive, one, three;
	fp				f, n, t, t1, s, s1, temp;
	static bool		initGood=false;
	INT32			i, bp;
	
	if(!initGood)
	{
		init(pointFive, 1.);
		init(three, 3.);
		pointFive.e = -1;
		init(one, 1.);
		initGood = true;
	}
	
	//if(abs(Lg2(x))>1000)
	//	return; //avoids overflow
		
	if(!x.i.n)
	{
		equate(z, one);
		return;
	}
	
	bp = blockPrec;
	
	blockPrec++;
	
	div(t1, x, Pi());
	floor(n, t1); 
	add(temp, n, pointFive);
	sub(f, t1, temp);
	mul(f, Pi(), f);
	
	if(!n.i.n || isItEven(n))
	{
		f.i.n = -f.i.n;
	}
	
	equate(t, f);
	s.i.n = 0;
	s.e = 0;
	
	fp		fSqr;
	fSqr = f*f;
	
	for(i=3; ;i+=2)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t, t, fSqr);
		div(t, t, i*(i-1));
		t.i.n = -t.i.n;
	}
	
	blockPrec = bp;
	
	equate(z, s);
	
	fpNormalize(z);
	
	return;
	
}/* cos */


fp cos(const fp& x)
{
	fp		z;
	
	cos(z, x);
	return z;
	
}/* cos */


void tan(fp& z, const fp& x)
{
	fp				myS, myC;
		
	cos(myC, x);
	
	if(!myC.i.n)
		return;
	
	sin(myS, x);
	
	div(z, myS, myC);
	
	return;
	
}/* tan */


fp tan(const fp& x)
{
	fp		z;
	
	tan(z, x);
	return z;
	
}/* tan */


void atan(fp& z, const fp& x)
{
	static fp		one, pointFive, pointZeroOne;
	fp				theTan, yPrev, y, y2, currentTerm, pi, temp;
	static bool		initGood=false;
	INT32			numReductions;
	bool			isReciprocal, isNegative;
	INT32			i;
	
	if(!initGood)
	{
		init(pointFive, 1.);
		pointFive.e = -1;
		init(pointZeroOne, .01);
		init(one, 1.);
		initGood = true;
	}
	
	numReductions = 0;
	
	if(!x.i.n)
	{
		z.i.n = 0;
		return;
	}
	
	blockPrec+=2;
	
	equate(theTan, x);
	
	if(theTan.i.n<0)
	{
		theTan.i.n = -theTan.i.n;
		isNegative = true;
	}
	else
		isNegative = false;
	
	if(!compare(theTan, one))
	{
		ComputePi(pi);
		blockPrec-=2;
		div(z, pi, 4);
		if(isNegative)
			z.i.n = -z.i.n;
		return;
	}
	
	if(compare(theTan, one)==1)
	{
		div(theTan, one, theTan);
		isReciprocal = true;
	}
	else
		isReciprocal = false;
	
	// divides arctan by a factor of two for each iteration
	while(compare(theTan, pointZeroOne)==1)
	{
		div(y, one, theTan);
		mul(y2, y, y);
		add(y2, y2, one);
		sqrt(y2, y2);
		add(y2, y2, y);
		div(theTan, one, y2);
		numReductions++;
	}
	
	// series method
	// atan z = z - z^3/3 + z^5/5 - z^7/7 + ...
	
	equate(y, theTan);
	equate(yPrev, theTan);
	equate(currentTerm, theTan);
	
	mul(y2, y, y);
	y2.i.n = -y2.i.n;
	
	for(i=0; ;i+=2)
	{
		mul(currentTerm, currentTerm, y2);
		div(temp, currentTerm, i+3);
		add(y, y, temp);
		if(!compare(y, yPrev))
			break;
		equate(yPrev, y);
	}
	
	equate(temp, one);
	temp.e+=numReductions;
	
	mul(y, y, temp);
	
	if(isReciprocal)
	{
		ComputePi(pi);
		pi.e--; // pi/2
		sub(y, pi, y);
	}
	
	if(isNegative)
		y.i.n = -y.i.n;
	
	equate(z, y);
	
	blockPrec-=2;
	
	fpNormalize(z);
	
}/* atan */


fp atan(const fp& x)
{
	fp		z;
	
	atan(z, x);
	return z;
	
}/* atan */

// z = atan(y/x) in correct quadrant
void atan2(fp& z, const fp& y, const fp& x)
{
	fp				myPi, myTan;

	if(!x.i.n && !y.i.n)
	{
		z.i.n = 0; // actually undetermined
		z.e = 0;
		return; 
	}
	
	ComputePi(myPi);
	
	if(!x.i.n)
	{
		equate(z, myPi);
		div(z, z, 2.);
		if(y.i.n<0)
			z.i.n = -z.i.n;
		return;
	}
	
	if(!y.i.n)
	{
		if(x.i.n>0)
		{
			z.i.n = 0;
			z.e = 0;
		}
		else
			equate(z, myPi);
		return;
	}
	
	div(myTan, y, x);
	//myTan = abs(myTan); // causes bus error, I suppose because myTan is a static variable
						// and abs(myTan) is a different variable
	if(myTan.i.n<0)
		myTan.i.n = -myTan.i.n;
	atan(z, myTan);
	
	if(x.i.n<0 && y.i.n>=0)
	{
		sub(z, myPi, z);
		return;
	}
	
	if(x.i.n<0 && y.i.n<0)
	{
		sub(z, z, myPi);
		return;
	}
	
	if(x.i.n>=0 && y.i.n<0)
		z.i.n = -z.i.n;
	
}/* atan2 */


fp atan2(const fp& y, const fp& x)
{
	fp		z;
	
	atan2(z, y, x);
	return z;
	
}/* atan2 */


void asin(fp& z, const fp& x)
{
	static fp		one;
	fp				theTan, y, temp;
	static bool		initGood=false;
	bool			isNegative;
	
	if(!initGood)
	{
		init(one, 1.);
		initGood = true;
	}
	
	if(!x.i.n)
	{
		z.i.n = 0;
		z.e = 0;
		return;
	}
	
	blockPrec+=2;
	
	equate(temp, x);
	
	if(temp.i.n<0)
	{
		temp.i.n = -temp.i.n;
		isNegative = true;
	}
	else
		isNegative = false;
	
	if(!compare(temp, one))
	{
		ComputePi(y);
		y.e--; // pi/2
		if(isNegative)
			y.i.n = -y.i.n;
		equate(z, y);
		blockPrec-=2;
		fpNormalize(z);
		return;
	}
	
	mul(theTan, temp, temp);
	sub(theTan, one, theTan);
	sqrt(theTan, theTan);
	div(theTan, temp, theTan);
	
	atan(y, theTan);
	
	if(isNegative)
		y.i.n = -y.i.n;
	
	equate(z, y);
	
	blockPrec-=2;
	
	fpNormalize(z);
	
}/* asin */


fp asin(const fp& x)
{
	fp		z;
	
	asin(z, x);
	return z;
	
}/* asin */


void acos(fp& z, const fp& x)
{
	static fp		one;
	fp				theTan, y, pi, temp;
	static bool		initGood=false;
	bool			isNegative;
	
	if(!initGood)
	{
		init(one, 1.);
		initGood = true;
	}
	
	if(!x.i.n)
	{
		ComputePi(z);
		z.e--; // pi/2
		return;
	}
	
	blockPrec+=2;  // was ++
	
	equate(temp, x);
	
	if(temp.i.n<0)
	{
		temp.i.n = -temp.i.n;
		isNegative = true;
	}
	else
		isNegative = false;
	
	mul(theTan, temp, temp);
	sub(theTan, one, theTan);
	sqrt(theTan, theTan);
	div(theTan, theTan, temp);
	
	atan(y, theTan);
	
	if(isNegative)
	{
		ComputePi(pi);
		sub(y, pi, y);
	}
	
	equate(z, y);
	
	blockPrec-=2;
	
	fpNormalize(z);
	
}/* acos */


fp acos(const fp& x)
{
	fp		z;
	
	acos(z, x);
	return z;
	
}/* acos */


void sinh(fp& z, const fp& x)
{
	static fp		pointOne;
	fp				xt, temp, temp1, t2, s, s1, t;
	static bool		initGood=false;
	INT32			i;
	
	if(!initGood)
	{
		init(pointOne, .1);
		initGood = true;
	}
	
	equate(xt, x);
	if(xt.i.n<0)
		xt.i.n = -xt.i.n;
	if(compare(xt, pointOne)==1)
	{
		exp(temp, x);
		div(temp1, 1., temp);
		sub(temp1, temp, temp1);
		div(z, temp1, 2.);
		return;
	}
	
	// now abs(x)<=.1 so use power series
	s.i.n = 0;
	s.e = 0;
	equate(t, 1.);
	mul(t2, x, x);
	
	for (i=2; ;i+=2)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t, t, t2);
		div(t, t, i*(i+1));
	}
	mul(z, x, s);

}/* sinh */


fp sinh(const fp& x)
{
	fp		z;
	
	sinh(z, x);
	return z;
	
}/* sinh */


void cosh(fp& z, const fp& x)
{
	fp			temp;
	
	equate(temp, x);
	exp(temp, temp);
	equate(z, temp);
	div(temp, 1., temp);
	add(z, z, temp);
	div(z, z, 2.);
	
}/* cosh */


fp cosh(const fp& x)
{
	fp		z;
	
	cosh(z, x);
	return z;
	
}/* cosh */


void tanh(fp& z, const fp& x)
{
	fp				temp, temp1;
	
	sinh(temp, x);
	cosh(temp1, x);
	div(z, temp, temp1);

}/* tanh */


fp tanh(const fp& x)
{
	fp		z;
	
	tanh(z, x);
	return z;
	
}/* tanh */


void asinh(fp& z, const fp& x)
{
	static fp		pointOne, max;
	fp				temp, xt, t2, s, s1, t;
	static bool		initGood=false;
	INT32			i;
	
	if(!initGood)
	{
		init(pointOne, .1);
		init(max, "1e30000");
		initGood = true;
	}
	
	equate(xt, x);
	if(xt.i.n<0)
		xt.i.n = -xt.i.n;
	if(compare(xt, max)==1)
	{
		blockPrec++;
		log(z, 2*xt);
		blockPrec--;
		//add(z, temp, log(2.));
		if(x.i.n<0)
			z.i.n = -z.i.n;
		return;
	}
	
	if(compare(xt, pointOne)==1)
	{
		blockPrec++;
		mul(temp, x, x);
		add(temp, temp, 1.);
		sqrt(temp, temp);
		add(temp, xt, temp);
		log(z, temp);
		if(x.i.n<0)
			z.i.n = -z.i.n;
		blockPrec--;
		return;
	}
	
	// abs(x)<=.1 so use power series
	blockPrec++;
	s.i.n = 0;
	s.e = 0;
	equate(t, 1.);
	mul(t2, x, x);
	t2.i.n = -t2.i.n;
	
	for(i=2; ;i+=2)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t, t, t2);
		mul(t, t, (i-1)*(i-1));
		div(t, t, i*(i+1));
	}
	mul(z, x, s);
	
	blockPrec--;
	
}/* asinh */


fp asinh(const fp& x)
{
	fp		z;
	
	asinh(z, x);
	return z;
	
}/* asinh */


void acosh(fp& z, const fp& x)
{
	static fp		max, one;
	fp				temp;
	static bool		initGood=false;
	
	if(!initGood)
	{
		init(max, "1e30000");
		init(one, 1.);
		initGood = true;
	}
	
	if(compare(x, one)==-1) // cosh >= 1 !
		return;
	
	if(!compare(x, one))
	{
		z.i.n = 0;
		z.e = 0;
		return;
	}
	
	if(compare(x, max)==1)
	{
		blockPrec++;
		log(z, 2*x);
		//add(z, z, log(2));
		blockPrec--;
		return;
	}
	
	blockPrec++;
	mul(temp, x, x);
	sub(temp, temp, one);
	sqrt(temp, temp);
	add(temp, x, temp);
	log(z, temp);
	blockPrec--;
	
}/* acosh */


fp acosh(const fp& x)
{
	fp		z;
	
	acosh(z, x);
	return z;
	
}/* acosh */


void atanh(fp& z, const fp& x)
{
	static fp		one, pointOne;
	fp				xt, temp, temp1, t2, t3, s, s1, t;
	static bool		initGood=false;
	INT32			i, result;
	
	if(!initGood)
	{
		init(one, 1.);
		init(pointOne, .1);
		initGood = true;
	}
	
	xt = x;
	if(xt.i.n<0)
		xt.i.n = -xt.i.n;
	result = compare(xt, one);
	if(result>=0) // abs(tanh)<1
		return;
	
	blockPrec+=2;
	
	result = compare(xt, pointOne);
	if(result>=0)  // xt>=.1
	{
		add(temp, one, xt);
		sub(temp1, one, xt);
		div(temp, temp, temp1);
		log(temp, temp);
		div(z, temp, 2.);
		if(x.i.n<0)
			z.i.n = -z.i.n;
		blockPrec-=2;
		return;
	}
	
	// now xt<.1 so use power series
	
	s.i.n = 0;
	s.e = 0;
	equate(t, one);
	mul(t2, xt, xt);
	equate(t3, one);
	
	for(i=3; ;i+=2)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t3, t3, t2);
		div(t, t3, i);
	}
	mul(z, xt, s);
	if(x.i.n<0)
		z.i.n = -z.i.n;
	
	blockPrec-=2;
	
}/* atanh */


fp atanh(const fp& x)
{
	fp		z;
	
	atanh(z, x);
	return z;
	
}/* atanh */


void factorial(fp& z, const fp& x)
{
	Gamma(z, x+1);

}/* factorial */


fp factorial(const fp& x)
{
	fp		z;
	
	factorial(z, x);
	return z;
	
}/* factorial */


void factorial(fp& z, double x)
{
	Gamma(z, x+1);
	
}/* factorial */


fp factorial(double x)
{
	fp		z;
	
	factorial(z, x);
	
	return z;
	
}/* factorial */


