qtrans.c

Calc Software Package for Number Calc

	if (qisneg(q1)) {
		math_error("Negative base for qpower");
		/*NOTREACHED*/
	}
	if (qisone(q2))
		return qmappr(q1, epsilon, 24);
	if (zrel(q1->num, q1->den) < 0) {
		q1tmp = qinv(q1);
		q2tmp = qneg(q2);
	} else {
		q1tmp = qlink(q1);
		q2tmp = qlink(q2);
	}
	if (qisone(q2tmp)) {
		qfree(q2tmp);
		q2tmp = qmappr(q1tmp, epsilon, 24);
		qfree(q1tmp);
		return q2tmp;
	}
	m = qilog2(q1tmp);
	n = qilog2(epsilon);
	if (qisneg(q2tmp)) {
		if (m > 0) {
			tmp1 = itoq(m);
			tmp2 = qmul(tmp1, q2tmp);
			m = qtoi(tmp2);
		} else {
			tmp1 = qdec(q1tmp);
			tmp2 = qqdiv(tmp1, q1tmp);
			qfree(tmp1);
			tmp1 = qmul(tmp2, q2tmp);
			qfree(tmp2);
			tmp2 = qmul(tmp1, &_qlge_);
			m = qtoi(tmp2);
		}
	} else {
		if (m > 0) {
			tmp1 = itoq(m + 1);
			tmp2 = qmul(tmp1, q2tmp);
			m = qtoi(tmp2);
		} else {
			tmp1 = qdec(q1tmp);
			tmp2 = qmul(tmp1, q2tmp);
			qfree(tmp1);
			tmp1 = qmul(tmp2, &_qlge_);
			m = qtoi(tmp1);
		}
	}
	qfree(tmp1);
	qfree(tmp2);
	if (m > (1 << 30)) {
		qfree(q1tmp);
		qfree(q2tmp);
		return NULL;
	}
	m += 1;
	if (m < n) {
		qfree(q1tmp);
		qfree(q2tmp);
		return qlink(&_qzero_);
	}
	tmp1 = qqdiv(epsilon, q2tmp);
	tmp2 = qscale(tmp1, -m - 4);
	epsilon2 = qqabs(tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qln(q1tmp, epsilon2);
	qfree(epsilon2);
	tmp2 = qmul(tmp1, q2tmp);
	qfree(tmp1);
	qfree(q1tmp);
	qfree(q2tmp);
	if (qisneg(tmp2)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		tmp2 = qexprel(tmp1, m - n + 3);
		if (tmp2 == NULL) {
			qfree(tmp1);
			return NULL;
		}
		qfree(tmp1);
		tmp1 = qinv(tmp2);
	} else {
		tmp1 = qexprel(tmp2, m - n + 3) ;
		if (tmp1 == NULL) {
			qfree(tmp2);
			return NULL;
		}
	}
	qfree(tmp2);
	tmp2 = qmappr(tmp1, epsilon, 24L);
	qfree(tmp1);
	return tmp2;
}


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

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for root");
		/*NOTREACHED*/
	}
	if (qisneg(q2) || qiszero(q2) || qisfrac(q2)) {
		math_error("Taking bad root of number");
		/*NOTREACHED*/
	}
	if (qiszero(q1) || qisone(q1) || qisone(q2))
		return qlink(q1);
	if (qistwo(q2))
		return qsqrt(q1, epsilon, 24L);
	neg = qisneg(q1);
	if (neg) {
		if (ziseven(q2->num)) {
			math_error("Taking even root of negative number");
			/*NOTREACHED*/
		}
		q1 = qqabs(q1);
	}
	tmp2 = qinv(q2);
	tmp1 = qpower(q1, tmp2, epsilon);
	qfree(tmp2);
	if (tmp1 == NULL)
		return NULL;
	if (neg) {
		tmp2 = qneg(tmp1);
		qfree(tmp1);
		tmp1 = tmp2;
	}
	return tmp1;
}


/* Calculate the hyperbolic cosine function to the nearest or next to
 * nearest multiple of epsilon.
 * This is calculated using cosh(x) = (exp(x) + 1/exp(x))/2;
 */
NUMBER *
qcosh(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;

	epsilon1 = qscale(epsilon, -2);
	tmp1 = qqabs(q);
	tmp2 = qexp(tmp1, epsilon1);
	qfree(tmp1);
	qfree(epsilon1);
	if (tmp2 == NULL)
		return NULL;
	tmp1 = qinv(tmp2);
	tmp3 = qqadd(tmp1, tmp2);
	qfree(tmp1)
	qfree(tmp2)
	tmp1 = qscale(tmp3, -1);
	qfree(tmp3);
	tmp2 = qmappr(tmp1, epsilon, 24L);
	qfree(tmp1);
	return tmp2;
}


/*
 * Calculate the hyperbolic sine to the nearest or next to nearest
 * multiple of epsilon.
 * This is calculated using sinh(x) = (exp(x) - 1/exp(x))/2.
 */
NUMBER *
qsinh(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;

	if (qiszero(q))
		return qlink(&_qzero_);
	epsilon1 = qscale(epsilon, -3);
	tmp1 = qqabs(q);
	tmp2 = qexp(tmp1, epsilon1);
	qfree(tmp1);
	qfree(epsilon1);
	if (tmp2 == NULL)
		return NULL;
	tmp1 = qinv(tmp2);
	tmp3 = qispos(q) ? qsub(tmp2, tmp1) : qsub(tmp1, tmp2);
	qfree(tmp1)
	qfree(tmp2)
	tmp1 = qscale(tmp3, -1);
	qfree(tmp3);
	tmp2 = qmappr(tmp1, epsilon, 24L);
	qfree(tmp1);
	return tmp2;
}


/*
 * Calculate the hyperbolic tangent to the nearest or next to nearest
 * multiple of epsilon.
 * This is calculated using the formulae:
 *	tanh(x) = 1 or -1 for very large abs(x)
 *	tanh(x) = (+ or -)(1 - 2 * exp(2 * x))	for large abx(x)
 *	tanh(x) = (exp(2*x) - 1)/(exp(2*x) + 1) otherwise
 */
NUMBER *
qtanh(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3;
	long n, m;

	n = qilog2(epsilon);
	if (n > 0 || qiszero(q))
		return qlink(&_qzero_);
	n = -n;
	tmp1 = qqabs(q);
	tmp2 = qmul(tmp1, &_qlge_);
	m = qtoi(tmp2);		/* exp(|q|) < 2^(m+1) or m == MAXLONG */
	qfree(tmp2);
	if (m > 1 + n/2) {
		qfree(tmp1);
		return q->num.sign ? qlink(&_qnegone_) : qlink(&_qone_);
	}
	tmp2 = qscale(tmp1, 1);
	qfree(tmp1);
	tmp1 = qexprel(tmp2, 2 + n);
	qfree(tmp2);
	if (m > 1 + n/4) {
		tmp2 = qqdiv(&_qtwo_, tmp1);
		qfree(tmp1);
		tmp1 = qsub(&_qone_, tmp2);
		qfree(tmp2);
	} else {
		tmp2 = qdec(tmp1);
		tmp3 = qinc(tmp1);
		qfree(tmp1);
		tmp1 = qqdiv(tmp2, tmp3);
		qfree(tmp2);
		qfree(tmp3);
	}
	tmp2 = qmappr(tmp1, epsilon, 24L);
	qfree(tmp1);
	if (qisneg(q)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		return tmp1;
	}
	return tmp2;
}


/*
 * Hyperbolic cotangent.
 * Calculated using coth(x) = 1 + 2/(exp(2*x) - 1)
 */
NUMBER *
qcoth(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *res;
	long n, k;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for coth");
		/*NOTREACHED*/
	}
	if (qiszero(q)) {
		math_error("Zero argument for coth");
		/*NOTREACHED*/
	}
	tmp1 = qscale(q, 1);
	tmp2 = qqabs(tmp1);
	qfree(tmp1);
	k = qilog2(tmp2);
	n = qilog2(epsilon);
	if (k > 0) {
		tmp1 = qmul(&_qlge_, tmp2);
		k = qtoi(tmp1);
		qfree(tmp1);
	} else {
		k = 2 * k;
	}
	k = 4 - k - n;
	if (k < 4)
		k = 4;
	tmp1 = qexprel(tmp2,  k);
	qfree(tmp2);
	tmp2 = qdec(tmp1);
	qfree(tmp1);
	if (qiszero(tmp2)) {
		math_error("This should not happen ????");
		/*NOTREACHED*/
	}
	tmp1 = qinv(tmp2);
	qfree(tmp2);
	tmp2 = qscale(tmp1, 1);
	qfree(tmp1);
	tmp1 = qinc(tmp2);
	qfree(tmp2);
	if (qisneg(q)) {
		tmp2 = qneg(tmp1);
		qfree(tmp1);
		tmp1 = tmp2;
	}
	res = qmappr(tmp1, epsilon, 24L);
	qfree(tmp1);
	return res;
}


NUMBER *
qsech(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *res;
	long n, k;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for sech");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qone_);

	tmp1 = qqabs(q);
	k = 0;
	if (zrel(tmp1->num, tmp1->den) >= 0) {
		tmp2 = qmul(&_qlge_, tmp1);
		k = qtoi(tmp2);
		qfree(tmp2);
	}
	n = qilog2(epsilon);
	if (k + n > 1) {
		qfree(tmp1);
		return qlink(&_qzero_);
	}
	tmp2 = qexprel(tmp1, 4 - k - n);
	qfree(tmp1);
	tmp1 = qinv(tmp2);
	tmp3 = qqadd(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qinv(tmp3);
	qfree(tmp3);
	tmp2 = qscale(tmp1, 1);
	qfree(tmp1);
	res = qmappr(tmp2, epsilon, 24L);
	qfree(tmp2);
	return res;
}


NUMBER *
qcsch(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *res;
	long n, k;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for csch");
		/*NOTREACHED*/
	}
	if (qiszero(q)) {
		math_error("Zero argument for csch");
		/*NOTREACHED*/
	}

	n = qilog2(epsilon);
	tmp1 = qqabs(q);
	if (zrel(tmp1->num, tmp1->den) >= 0) {
		tmp2 = qmul(&_qlge_, tmp1);
		k = qtoi(tmp2);
		qfree(tmp2);
	} else {
		k = 2 * qilog2(tmp1);
	}
	if (k + n >= 1) {
		qfree(tmp1);
		return qlink(&_qzero_);
	}
	tmp2 = qexprel(tmp1, 4 - k - n);
	qfree(tmp1);
	tmp1 = qinv(tmp2);
	if (qisneg(q))
		tmp3 = qsub(tmp1, tmp2);
	else
		tmp3 = qsub(tmp2, tmp1);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qinv(tmp3);
	qfree(tmp3)
	tmp2 = qscale(tmp1, 1);
	qfree(tmp1);
	res = qmappr(tmp2, epsilon, 24L);
	qfree(tmp2);
	return res;
}


/*
 * Compute the hyperbolic arccosine within the specified accuracy.
 * This is calculated using the formula:
 *	acosh(x) = ln(x + sqrt(x^2 - 1)).
 */
NUMBER *
qacosh(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *epsilon1;
	long n;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for acosh");
		/*NOTREACHED*/
	}
	if (qisone(q))
		return qlink(&_qzero_);
	if (zrel(q->num, q->den) < 0)
		return NULL;
	n = qilog2(epsilon);
	epsilon1 = qbitvalue(n - 3);
	tmp1 = qsquare(q);
	tmp2 = qdec(tmp1);
	qfree(tmp1);
	tmp1 = qsqrt(tmp2, epsilon1, 24L);
	qfree(tmp2);
	tmp2 = qqadd(tmp1, q);
	qfree(tmp1);
	tmp1 = qln(tmp2, epsilon1);
	qfree(tmp2);
	qfree(epsilon1);
	tmp2 = qmappr(tmp1, epsilon, 24L);
	qfree(tmp1);
	return tmp2;
}


/*
 * Compute the hyperbolic arcsine within the specified accuracy.
 * This is calculated using the formula:
 *	asinh(x) = ln(x + sqrt(x^2 + 1)).
 */
NUMBER *
qasinh(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *epsilon1;
	long n;
	BOOL neg;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for asinh");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qzero_);
	neg = qisneg(q);
	q = qqabs(q);
	n = qilog2(epsilon);
	epsilon1 = qbitvalue(n - 3);
	tmp1 = qsquare(q);
	tmp2 = qinc(tmp1);
	qfree(tmp1);
	tmp1 = qsqrt(tmp2, epsilon1, 24L);
	qfree(tmp2);
	tmp2 = qqadd(tmp1, q);
	qfree(tmp1);
	tmp1 = qln(tmp2, epsilon1);
	qfree(tmp2);
	qfree(q);
	qfree(epsilon1);
	tmp2 = qmappr(tmp1, epsilon, 24L);
	if (neg) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		return tmp1;
	}
	return tmp2;
}


/*
 * Compute the hyperbolic arctangent within the specified accuracy.
 * This is calculated using the formula:
 *	atanh(x) = ln((1 + x) / (1 - x)) / 2.
 */
NUMBER *
qatanh(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
	ZVALUE z;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for atanh");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qzero_);
	z = q->num;
	z.sign = 0;
	if (zrel(z, q->den) >= 0)
		return NULL;
	tmp1 = qinc(q);
	tmp2 = qsub(&_qone_, q);
	tmp3 = qqdiv(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	epsilon1 = qscale(epsilon, 1L);
	tmp1 = qln(tmp3, epsilon1);
	qfree(tmp3);
	tmp2 = qscale(tmp1, -1L);
	qfree(tmp1);
	qfree(epsilon1);
	return tmp2;
}


/*
 * Inverse hyperbolic secant function
 */
NUMBER *
qasech(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp, *res;

	tmp = qinv(q);
	res = qacosh(tmp, epsilon);
	qfree(tmp);
	return res;
}


/*
 * Inverse hyperbolic cosecant function
 */
NUMBER *
qacsch(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp, *res;

	tmp = qinv(q);
	res = qasinh(tmp, epsilon);
	qfree(tmp);
	return res;
}


/*
 * Inverse hyperbolic cotangent function
 */
NUMBER *
qacoth(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp, *res;

	tmp = qinv(q);
	res = qatanh(tmp, epsilon);
	qfree(tmp);
	return res;
}