/*
 * Copyright (c) 1994 David I. Bell
 * Permission is granted to use, distribute, or modify this source,
 * provided that this copyright notice remains intact.
 *
 * Transcendental functions for real numbers.
 * These are sin, cos, exp, ln, power, cosh, sinh.
 */

#include "qmath.h"

static NUMBER *qcosimpl MATH_PROTO((NUMBER *q, NUMBER *epsilon, BOOL *sinisnegPtr));

/*
 * Calculate the cosine of a number with an accuracy within epsilon.
 * This also saves the sign of the corresponding sin function.
 */
NUMBER *
qcos(q, epsilon)
	NUMBER *q, *epsilon;
{
    BOOL sinisneg = qisneg(q);

    return qcosimpl(q, epsilon, &sinisneg);
}

static NUMBER *
qcosimpl(q, epsilon, sinisnegPtr)
	NUMBER *q, *epsilon;
	BOOL *sinisnegPtr;
{
	NUMBER *term, *sum, *qsq, *epsilon2, *tmp;
	FULL n, i;
	long scale, bits, bits2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for cosine");
	if (qiszero(q))
		return qlink(&_qone_);
	bits = qprecision(epsilon) + 1;
	epsilon = qscale(epsilon, -4L);
	/*
	 * If the argument is larger than one, then divide it by a power of two
	 * so that it is one or less.  This will make the series converge quickly.
	 * We will extrapolate the result for the original argument afterwards.
	 */
	scale = zhighbit(q->num) - zhighbit(q->den) + 1;
	if (scale < 0)
		scale = 0;
	if (scale > 0) {
		q = qscale(q, -scale);
		tmp = qscale(epsilon, -scale);
		qfree(epsilon);
		epsilon = tmp;
	}
	epsilon2 = qscale(epsilon, -4L);
	qfree(epsilon);
	bits2 = qprecision(epsilon2) + 10;
	/*
	 * Now use the Taylor series expansion to calculate the cosine.
	 * Keep using approximations so that the fractions don't get too large.
	 */
	qsq = qsquare(q);
	if (scale > 0)
		qfree(q);
	term = qlink(&_qone_);
	sum = qlink(&_qone_);
	n = 0;
	while (qrel(term, epsilon2) > 0) {
		i = ++n;
		i *= ++n;
		tmp = qmul(term, qsq);
		qfree(term);
		term = qdivi(tmp, (long) i);
		qfree(tmp);
		tmp = qbround(term, bits2);
		qfree(term);
		term = tmp;
		if (n & 2)
			tmp = qsub(sum, term);
		else
			tmp = qadd(sum, term);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	qfree(term);
	qfree(qsq);
	qfree(epsilon2);
	/*
	 * Now scale back up to the original value of x by using the formula:
	 *	cos(2 * x) = 2 * (cos(x) ^ 2) - 1.
	 */
	while (--scale >= 0) {
		if (qisneg(sum))
			*sinisnegPtr = !(*sinisnegPtr);
		tmp = qsquare(sum);
		qfree(sum);
		sum = qscale(tmp, 1L);
		qfree(tmp);
		tmp = qdec(sum);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	tmp = qbround(sum, bits);
	qfree(sum);
	return tmp;
}


/*
 * Calculate the sine of a number with an accuracy within epsilon.
 * This is calculated using the formula:
 *	sin(x)^2 + cos(x)^2 = 1.
 * The only tricky bit is resolving the sign of the result.
 * Future: Use sin(3*x) = 3*sin(x) - 4*sin(x)^3.
 */
NUMBER *
qsin(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;
	BOOL sinisneg = qisneg(q);

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for sine");
	if (qiszero(q))
		return qlink(q);
	epsilon2 = qsquare(epsilon);
	tmp1 = qcosimpl(q, epsilon2, &sinisneg);
	qfree(epsilon2);
	tmp2 = qlegtoleg(tmp1, epsilon, sinisneg);
	qfree(tmp1);
	return tmp2;
}


/*
 * Calculate the tangent function.
 */
NUMBER *
qtan(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *cosval, *sinval, *epsilon2, *tmp, *res;
	BOOL sinisneg = qisneg(q);

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for tangent");
	if (qiszero(q))
		return qlink(q);
	epsilon2 = qsquare(epsilon);
	cosval = qcosimpl(q, epsilon2, &sinisneg);
	sinval = qlegtoleg(cosval, epsilon2, sinisneg);
	qfree(epsilon2);
	tmp = qdiv(sinval, cosval);
	qfree(cosval);
	qfree(sinval);
	res = qbround(tmp, qprecision(epsilon) + 1);
	qfree(tmp);
	return res;
}


/*
 * Calculate the arcsine function.
 * The result is in the range -pi/2 to pi/2.
 */
NUMBER *
qasin(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *sum, *term, *epsilon2, *qsq, *tmp;
	FULL n, i;
	long bits, bits2;
	int neg;
	NUMBER mulnum;
	HALF numval[2];
	HALF denval[2];

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for arcsine");
	if (qiszero(q))
		return qlink(&_qzero_);
	if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
		math_error("Argument too large for asin");
	neg = qisneg(q);
	q = qabs(q);
	epsilon = qscale(epsilon, -4L);
	epsilon2 = qscale(epsilon, -4L);
	mulnum.num.sign = 0;
	mulnum.num.len = 1;
	mulnum.num.v = numval;
	mulnum.den.sign = 0;
	mulnum.den.len = 1;
	mulnum.den.v = denval;
	/*
	 * If the argument is too near one (we use .5) then reduce the
	 * argument to a more accurate range using the formula:
	 *	asin(x) = 2 * asin(sqrt((1 - sqrt(1 - x^2)) / 2)).
	 */
	if (qrel(q, &_qonehalf_) > 0) {
		sum = qlegtoleg(q, epsilon2, FALSE);
		qfree(q);
		tmp = qsub(&_qone_, sum);
		qfree(sum);
		sum = qscale(tmp, -1L);
		qfree(tmp);
		tmp = qsqrt(sum, epsilon2);
		qfree(sum);
		qfree(epsilon2);
		sum = qasin(tmp, epsilon);
		qfree(tmp);
		qfree(epsilon);
		tmp = qscale(sum, 1L);
		qfree(sum);
		if (neg) {
			sum = qneg(tmp);
			qfree(tmp);
			tmp = sum;
		}
		return tmp;
	}
	/*
	 * Argument is between zero and .5, so use the series.
	 */
	tmp = qscale(epsilon, -4L);
	qfree(epsilon);
	epsilon = tmp;
	qfree(epsilon2);
	epsilon2 = qscale(epsilon, -4L);
	bits = qprecision(epsilon) + 1;
	bits2 = bits + 10;
	sum = qlink(q);
	term = qlink(q);
	qsq = qsquare(q);
	qfree(q);
	n = 1;
	while (qrel(term, epsilon2) > 0) {
		/*
		 * Maximum value for n is 33207, which occurs at q == 0.5 and
		 * mp_precision at the maximum value of 10000.  This value of
		 * n is small enough that (n+1)*(n+2) will fit into two 16-bit
		 * unsigned values, hence numval and denval don't need length
		 * longer than 2.
		 */
		i = n * n;
		numval[0] = i & BASE1;
		if (i >= BASE) {
			numval[1] = i / BASE;
			mulnum.num.len = 2;
		}
		i = (n + 1) * (n + 2);
		denval[0] = i & BASE1;
		if (i >= BASE) {
			denval[1] = i / BASE;
			mulnum.den.len = 2;
		}
		tmp = qmul(term, qsq);
		qfree(term);
		term = qmul(tmp, &mulnum);
		qfree(tmp);
		tmp = qbround(term, bits2);
		qfree(term);
		term = tmp;
		tmp = qadd(sum, term);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
		n += 2;
	}
	qfree(epsilon);
	qfree(epsilon2);
	qfree(term);
	qfree(qsq);
	tmp = qbround(sum, bits);
	qfree(sum);
	if (neg) {
		term = qneg(tmp);
		qfree(tmp);
		tmp = term;
	}
	return tmp;
}


/*
 * Calculate the acos function.
 * The result is in the range 0 to pi.
 */
NUMBER *
qacos(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for arccosine");
	if (qisone(q))
		return qlink(&_qzero_);
	if ((qrel(q, &_qone_) > 0) || (qrel(q, &_qnegone_) < 0))
		math_error("Argument too large for acos");
	/*
	 * Calculate the result using the formula:
	 *	acos(x) = asin(sqrt(1 - x^2)).
	 * The formula is only good for positive x, so we must fix up the
	 * result for negative values.
	 */
	epsilon2 = qscale(epsilon, -8L);
	tmp1 = qlegtoleg(q, epsilon2, FALSE);
	qfree(epsilon2);
	tmp2 = qasin(tmp1, epsilon);
	qfree(tmp1);
	if (!qisneg(q))
		return tmp2;
	/*
	 * For negative values, we need to subtract the asin from pi.
	 */
	tmp1 = qpi(epsilon);
	tmp3 = qsub(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qbround(tmp3, qprecision(epsilon) + 1);
	qfree(tmp3);
	return tmp1;
}


/*
 * Calculate the arctangent function with a accuracy less than epsilon.
 * This uses the formula:
 *	atan(x) = asin(sqrt(x^2 / (x^2 + 1))).
 */
NUMBER *
qatan(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for arctangent");
	if (qiszero(q))
		return qlink(&_qzero_);
	tmp1 = qsquare(q);
	tmp2 = qinc(tmp1);
	tmp3 = qdiv(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	epsilon2 = qscale(epsilon, -8L);
	tmp1 = qsqrt(tmp3, epsilon2);
	qfree(epsilon2);
	qfree(tmp3);
	tmp2 = qasin(tmp1, epsilon);
	qfree(tmp1);
	if (qisneg(q)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		tmp2 = tmp1;
	}
	return tmp2;
}


/*
 * Calculate the angle which is determined by the point (x,y).
 * This is the same as arctan for non-negative x, but gives the correct
 * value for negative x.  By convention, y is the first argument.
 * For example, qatan2(1, -1) = 3/4 * pi.
 */
NUMBER *
qatan2(qy, qx, epsilon)
	NUMBER *qy, *qx, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for atan2");
	if (qiszero(qy) && qiszero(qx)) {
		/* conform to 4.3BSD ANSI/IEEE 754-1985 math lib */
		return qlink(&_qzero_);
	}
	/*
	 * If the point is on the negative real axis, then the answer is pi.
	 */
	if (qiszero(qy) && qisneg(qx))
		return qpi(epsilon);
	/*
	 * If the point is in the right half plane, then use the normal atan.
	 */
	if (!qisneg(qx) && !qiszero(qx)) {
		if (qiszero(qy))
			return qlink(&_qzero_);
		tmp1 = qdiv(qy, qx);
		tmp2 = qatan(tmp1, epsilon);
		qfree(tmp1);
		return tmp2;
	}
	/*
	 * The point is in the left half plane.  Calculate the angle by finding
	 * the atan of half the angle using the formula:
	 *	atan2(y,x) = 2 * atan((sqrt(x^2 + y^2) - x) / y).
	 */
	epsilon2 = qscale(epsilon, -4L);
	tmp1 = qhypot(qx, qy, epsilon2);
	tmp2 = qsub(tmp1, qx);
	qfree(tmp1);
	tmp1 = qdiv(tmp2, qy);
	qfree(tmp2);
	tmp2 = qatan(tmp1, epsilon2);
	qfree(tmp1);
	qfree(epsilon2);
	tmp1 = qscale(tmp2, 1L);
	qfree(tmp2);
	return tmp1;
}


/*
 * Calculate the value of pi to within the required epsilon.
 * This uses the following formula which only needs integer calculations
 * except for the final operation:
 *	pi = 1 / SUMOF(comb(2 * N, N) ^ 3 * (42 * N + 5) / 2 ^ (12 * N + 4)),
 * where the summation runs from N=0.  This formula gives about 6 bits of
 * accuracy per term.  Since the denominator for each term is a power of two,
 * we can simply use shifts to sum the terms.  The combinatorial numbers
 * in the formula are calculated recursively using the formula:
 *	comb(2*(N+1), N+1) = 2 * comb(2 * N, N) * (2 * N + 1) / N.
 */
NUMBER *
qpi(epsilon)
	NUMBER *epsilon;
{
	ZVALUE comb;			/* current combinatorial value */
	ZVALUE sum;			/* current sum */
	ZVALUE tmp1, tmp2;
	NUMBER *r, *t1, qtmp;
	long shift;			/* current shift of result */
	long N;				/* current term number */
	long bits;			/* needed number of bits of precision */
	long t;

	if (qiszero(epsilon) || qisneg(epsilon))
		math_error("Bad epsilon value for pi");
	bits = qprecision(epsilon) + 4;
	comb = _one_;
	itoz(5L, &sum);
	N = 0;
	shift = 4;
	do {
		t = 1 + (++N & 0x1);
		(void) zdivi(comb, N / (3 - t), &tmp1);
		zfree(comb);
		zmuli(tmp1, t * (2 * N - 1), &comb);
		zfree(tmp1);
		zsquare(comb, &tmp1);
		zmul(comb, tmp1, &tmp2);
		zfree(tmp1);
		zmuli(tmp2, 42 * N + 5, &tmp1);
		zfree(tmp2);
		zshift(sum, 12L, &tmp2);
		zfree(sum);
		zadd(tmp1, tmp2, &sum);
		t = zhighbit(tmp1);
		zfree(tmp1);
		zfree(tmp2);
		shift += 12;
	} while ((shift - t) < bits);
	zfree(comb);	/* tp: fix leak */
	qtmp.num = _one_;
	qtmp.den = sum;
	t1 = qscale(&qtmp, shift);
	zfree(sum);
	r = qbround(t1, bits);
	qfree(t1);
	return r;
}


/*
 * Calculate the exponential function with a relative accuracy less than
 * epsilon.
 */
NUMBER *
qexp(q, epsilon)
	NUMBER *q, *epsilon;
{
	long scale;
	FULL n;
	long bits, bits2;
	NUMBER *sum, *term, *qs, *epsilon2, *tmp;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for exp");
	if (qiszero(q))
		return qlink(&_qone_);
	epsilon = qscale(epsilon, -4L);
	/*
	 * If the argument is larger than one, then divide it by a power of two
	 * so that it is one or less.  This will make the series converge quickly.
	 * We will extrapolate the result for the original argument afterwards.
	 * Also make the argument non-negative.
	 */
	qs = qabs(q);
	scale = zhighbit(q->num) - zhighbit(q->den) + 1;
	if (scale < 0)
		scale = 0;
	if (scale > 0) {
		if (scale > 100000)
			math_error("Very large argument for exp");
		tmp = qscale(qs, -scale);
		qfree(qs);
		qs = tmp;
		tmp = qscale(epsilon, -scale);
		qfree(epsilon);
		epsilon = tmp;
	}
	epsilon2 = qscale(epsilon, -4L);
	bits = qprecision(epsilon) + 1;
	bits2 = bits + 10;
	qfree(epsilon);
	/*
	 * Now use the Taylor series expansion to calculate the exponential.
	 * Keep using approximations so that the fractions don't get too large.
	 */
	sum = qlink(&_qone_);
	term = qlink(&_qone_);
	n = 0;
	while (qrel(term, epsilon2) > 0) {
		n++;
		tmp = qmul(term, qs);
		qfree(term);
		term = qdivi(tmp, (long) n);
		qfree(tmp);
		tmp = qbround(term, bits2);
		qfree(term);
		term = tmp;
		tmp = qadd(sum, term);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	qfree(term);
	qfree(qs);
	qfree(epsilon2);
	/*
	 * Now repeatedly square the answer to get the final result.
	 * Then invert it if the original argument was negative.
	 */
	while (--scale >= 0) {
		tmp = qsquare(sum);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	tmp = qbround(sum, bits);
	qfree(sum);
	if (qisneg(q)) {
		sum = qinv(tmp);
		qfree(tmp);
		tmp = sum;
	}
	return tmp;
}


/*
 * Calculate the natural logarithm of a number accurate to the specified
 * epsilon.
 */
NUMBER *
qln(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *term, *term2, *sum, *epsilon2, *tmp1, *tmp2, *maxr;
	long shift, bits, bits2;
	int j, k;
	FULL n;
	BOOL neg;

	if (qisneg(q) || qiszero(q))
		math_error("log of non-positive number");
	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon for ln");
	if (qisone(q))
		return qlink(&_qzero_);
	/*
	 * If the number is less than one, invert it and remember that
	 * the result is to be negative.
	 */
	neg = FALSE;
	if (zrel(q->num, q->den) < 0) {
		neg = TRUE;
		q = qinv(q);
	} else
		q = qlink(q);
	j = BASEB;
	k = zhighbit(q->num) - zhighbit(q->den) + 1;
	while (k >>= 1)
		j++;
	epsilon2 = qscale(epsilon, -j);
	bits = qprecision(epsilon) + 1;
	bits2 = qprecision(epsilon2) + 5;
	/*
	 * By repeated square-roots scale number down to a value close
	 * to 1 so that Taylor series to be used will converge rapidly.
	 * The effect of scaling will be reversed by a later shift.
	 */
	maxr = iitoq(BASE + 1, BASE);
	shift = 1;
	while (qrel(q, maxr) > 0) {
		tmp1 = qsqrt(q, epsilon2);
		qfree(q);
		q = tmp1;
		shift++;
	}
	qfree(maxr);
	/*
	 * Calculate a value which will always converge using the formula:
	 *	ln((1+x)/(1-x)) = ln(1+x) - ln(1-x).
	 */
	tmp1 = qdec(q);
	tmp2 = qinc(q);
	term = qdiv(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	qfree(q);
	/*
	 * Now use the Taylor series expansion to calculate the result.
	 */
	n = 1;
	term2 = qsquare(term);
	sum = qlink(term);
	while (qrel(term, epsilon2) > 0) {
		n += 2;
		tmp1 = qmul(term, term2);
		qfree(term);
		term = qbround(tmp1, bits2);
		qfree(tmp1);
		tmp1 = qdivi(term, (long) n);
		tmp2 = qadd(sum, tmp1);
		qfree(tmp1);
		qfree(sum);
		sum = qbround(tmp2, bits2);
		qfree(tmp2);
	}
	qfree(epsilon2);
	qfree(term);
	qfree(term2);
	/*
	 * Calculate the final result by multiplying by the proper power
	 * of two to undo the square roots done at the top, and possibly
	 * negating the result.
	 */
	tmp1 = qscale(sum, shift);
	qfree(sum);
	sum = qbround(tmp1, bits);
	qfree(tmp1);
	if (neg) {
		tmp1 = qneg(sum);
		qfree(sum);
		sum = tmp1;
	}
	return sum;
}


/*
 * Calculate the result of raising one number to the power of another.
 * The result is calculated to within the specified relative error.
 */
NUMBER *
qpower(q1, q2, epsilon)
	NUMBER *q1, *q2, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;

	if (qisint(q2))
		return qpowi(q1, q2);
	epsilon2 = qscale(epsilon, -4L);
	tmp1 = qln(q1, epsilon2);
	tmp2 = qmul(tmp1, q2);
	qfree(tmp1);
	tmp1 = qexp(tmp2, epsilon);
	qfree(tmp2);
	qfree(epsilon2);
	return tmp1;
}


/*
 * Calculate the Kth root of a number to within the specified accuracy.
 */
NUMBER *
qroot(q1, q2, epsilon)
	NUMBER *q1, *q2, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;
	int neg;

	if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
		math_error("Taking bad root of number");
	if (qiszero(q1) || qisone(q1) || qisone(q2))
		return qlink(q1);
	if (qistwo(q2))
		return qsqrt(q1, epsilon);
	neg = qisneg(q1);
	if (neg) {
		if (ziseven(q2->num))
			math_error("Taking even root of negative number");
		q1 = qabs(q1);
	}
	epsilon2 = qscale(epsilon, -4L);
	tmp1 = qln(q1, epsilon2);
	tmp2 = qdiv(tmp1, q2);
	qfree(tmp1);
	tmp1 = qexp(tmp2, epsilon);
	qfree(tmp2);
	qfree(epsilon2);
	if (neg) {
		tmp2 = qneg(tmp1);
		qfree(tmp1);
		tmp1 = tmp2;
	}
	return tmp1;
}


/*
 * Calculate the hyperbolic cosine function with a relative accuracy less
 * than epsilon.  This is defined by:
 *	cosh(x) = (exp(x) + exp(-x)) / 2.
 */
NUMBER *
qcosh(q, epsilon)
	NUMBER *q, *epsilon;
{
	long scale;
	FULL n;
	FULL m;
	long bits, bits2;
	NUMBER *sum, *term, *qs, *epsilon2, *tmp;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for exp");
	if (qiszero(q))
		return qlink(&_qone_);
	epsilon = qscale(epsilon, -4L);
	/*
	 * If the argument is larger than one, then divide it by a power of two
	 * so that it is one or less.  This will make the series converge quickly.
	 * We will extrapolate the result for the original argument afterwards.
	 */
	qs = qabs(q);
	scale = zhighbit(q->num) - zhighbit(q->den) + 1;
	if (scale < 0)
		scale = 0;
	if (scale > 0) {
		if (scale > 100000)
			math_error("Very large argument for exp");
		tmp = qscale(qs, -scale);
		qfree(qs);
		qs = tmp;
		tmp = qscale(epsilon, -scale);
		qfree(epsilon);
		epsilon = tmp;
	}
	epsilon2 = qscale(epsilon, -4L);
	bits = qprecision(epsilon) + 1;
	bits2 = bits + 10;
	qfree(epsilon);
	tmp = qsquare(qs);
	qfree(qs);
	qs = tmp;
	/*
	 * Now use the Taylor series expansion to calculate the exponential.
	 * Keep using approximations so that the fractions don't get too large.
	 */
	sum = qlink(&_qone_);
	term = qlink(&_qone_);
	n = 0;
	while (qrel(term, epsilon2) > 0) {
		m = ++n;
		m *= ++n;
		tmp = qmul(term, qs);
		qfree(term);
		term = qdivi(tmp, (long) m);
		qfree(tmp);
		tmp = qbround(term, bits2);
		qfree(term);
		term = tmp;
		tmp = qadd(sum, term);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	qfree(term);
	qfree(qs);
	qfree(epsilon2);
	/*
	 * Now bring the number back up into range to get the final result.
	 * This uses the formula:
	 *	cosh(2 * x) = 2 * cosh(x)^2 - 1.
	 */
	while (--scale >= 0) {
		tmp = qsquare(sum);
		qfree(sum);
		sum = qscale(tmp, 1L);
		qfree(tmp);
		tmp = qdec(sum);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	tmp = qbround(sum, bits);
	qfree(sum);
	return tmp;
}


/*
 * Calculate the hyperbolic sine with an accurary less than epsilon.
 * This is calculated using the formula:
 *	cosh(x)^2 - sinh(x)^2 = 1.
 */
NUMBER *
qsinh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for sinh");
	if (qiszero(q))
		return qlink(q);
	epsilon = qscale(epsilon, -4L);
	tmp1 = qcosh(q, epsilon);
	tmp2 = qsquare(tmp1);
	qfree(tmp1);
	tmp1 = qdec(tmp2);
	qfree(tmp2);
	tmp2 = qsqrt(tmp1, epsilon);
	qfree(tmp1);
	if (qisneg(q)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		tmp2 = tmp1;
	}
	qfree(epsilon);
	return tmp2;
}


/*
 * Calculate the hyperbolic tangent with an accurary less than epsilon.
 * This is calculated using the formula:
 *	tanh(x) = sinh(x) / cosh(x).
 */
NUMBER *
qtanh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *coshval;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for tanh");
	if (qiszero(q))
		return qlink(q);
	epsilon = qscale(epsilon, -4L);
	coshval = qcosh(q, epsilon);
	tmp2 = qsquare(coshval);
	tmp1 = qdec(tmp2);
	qfree(tmp2);
	tmp2 = qsqrt(tmp1, epsilon);
	qfree(tmp1);
	if (qisneg(q)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		tmp2 = tmp1;
	}
	qfree(epsilon);
	tmp1 = qdiv(tmp2, coshval);
	qfree(tmp2);
	qfree(coshval);
	return tmp1;
}

/* END CODE */
