qtrans.c

Calc Software Package for Number Calc

/*
 * qtrans - transcendental functions for real numbers
 *
 * Copyright (C) 1999-2004  David I. Bell and Ernest Bowen
 *
 * Primary author:  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.7 $
 * @(#) $Id: qtrans.c,v 29.7 2006/05/07 13:04:18 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/qtrans.c,v $
 *
 * Under source code control:	1990/02/15 01:48:22
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * Transcendental functions for real numbers.
 * These are sin, cos, exp, ln, power, cosh, sinh.
 */


#include "qmath.h"

HALF _qlgenum_[] = { 36744 };
HALF _qlgeden_[] = { 25469 };
NUMBER _qlge_ = { { _qlgenum_, 1, 0 }, { _qlgeden_, 1, 0 }, 1, NULL };

/*
 * cache the natural logarithm of 10
 */
static NUMBER *ln_10 = NULL;
static NUMBER *ln_10_epsilon = NULL;

static NUMBER *pivalue[2];
static NUMBER *qexprel(NUMBER *q, long bitnum);

/*
 * Evaluate and store in specified locations the sin and cos of a given
 * number to accuracy corresponding to a specified number of binary digits.
 */
void
qsincos(NUMBER *q, long bitnum, NUMBER **vs, NUMBER **vc)
{
	long n, m, k, h, s, t, d;
	NUMBER *qtmp1, *qtmp2;
	ZVALUE X, cossum, sinsum, mul, ztmp1, ztmp2, ztmp3;

	qtmp1 = qqabs(q);
	h = qilog2(qtmp1);
	qfree(qtmp1);
	k = bitnum + h + 1;
	if (k < 0) {
		*vs = qlink(&_qzero_);
		*vc = qlink(&_qone_);
		return;
	}
	s = k;
	if (k) {
		do {
			t = s;
			s = (s + k/s)/2;
		}
		while (t > s);
	}		/* s is int(sqrt(k)) */
	s++;
	if (s < -h)
		s = -h;
	n = h + s;	/* n is number of squarings that will be required */
	m = bitnum + n;
	while (s > 0) {		/* increasing m by ilog2(s) */
		s >>= 1;
		m++;
	}			/* m is working number of bits */
	qtmp1 = qscale(q, m - n);
	zquo(qtmp1->num, qtmp1->den, &X, 24);
	qfree(qtmp1);
	if (ziszero(X)) {
		zfree(X);
		*vs = qlink(&_qzero_);
		*vc = qlink(&_qone_);
		return;
	}
	zbitvalue(m, &cossum);
	zcopy(X, &sinsum);
	zcopy(X, &mul);
	d = 1;
	for (;;) {
		X.sign = !X.sign;
		zmul(X, mul, &ztmp1);
		zfree(X);
		zshift(ztmp1, -m, &ztmp2);
		zfree(ztmp1);
		zdivi(ztmp2, ++d, &X);
		zfree(ztmp2);
		if (ziszero(X))
			break;
		zadd(cossum, X, &ztmp1);
		zfree(cossum);
		cossum = ztmp1;
		zmul(X, mul, &ztmp1);
		zfree(X);
		zshift(ztmp1, -m, &ztmp2);
		zfree(ztmp1);
		zdivi(ztmp2, ++d, &X);
		zfree(ztmp2);
		if (ziszero(X))
			break;
		zadd(sinsum, X, &ztmp1);
		zfree(sinsum);
		sinsum = ztmp1;
	}
	zfree(X);
	zfree(mul);
	while (n-- > 0) {
		zsquare(cossum, &ztmp1);
		zsquare(sinsum, &ztmp2);
		zsub(ztmp1, ztmp2, &ztmp3);
		zfree(ztmp1);
		zfree(ztmp2);
		zmul(cossum, sinsum, &ztmp1);
		zfree(cossum);
		zfree(sinsum);
		zshift(ztmp3, -m, &cossum);
		zfree(ztmp3);
		zshift(ztmp1, 1 - m, &sinsum);
		zfree(ztmp1);
	}
	h = zlowbit(cossum);
	qtmp1 = qalloc();
	if (m > h) {
		zshift(cossum, -h, &qtmp1->num);
		zbitvalue(m - h, &qtmp1->den);
	} else {
		zshift(cossum, - m, &qtmp1->num);
	}
	zfree(cossum);
	*vc = qtmp1;
	h = zlowbit(sinsum);
	qtmp2 = qalloc();
	if (m > h) {
		zshift(sinsum, -h, &qtmp2->num);
		zbitvalue(m - h, &qtmp2->den);
	} else {
		zshift(sinsum, -m, &qtmp2->num);
	}
	zfree(sinsum);
	*vs = qtmp2;
	return;
}

/*
 * Calculate the cosine of a number to a near multiple of epsilon.
 * This calls qsincos() and discards the value of sin.
 */
NUMBER *
qcos(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *sin, *cos, *res;
	long n;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for cosine");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qone_);
	n = -qilog2(epsilon);
	if (n < 0)
		return qlink(&_qzero_);
	qsincos(q, n + 2, &sin, &cos);
	qfree(sin);
	res = qmappr(cos, epsilon, 24);
	qfree(cos);
	return res;
}



/*
 * This calls qsincos() and discards the value of cos.
 */
NUMBER *
qsin(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *sin, *cos, *res;
	long n;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for sine");
		/*NOTREACHED*/
	}
	n = -qilog2(epsilon);
	if (qiszero(q) || n < 0)
		return qlink(&_qzero_);
	qsincos(q, n + 2, &sin, &cos);
	qfree(cos);
	res = qmappr(sin, epsilon, 24);
	qfree(sin);
	return res;
}


/*
 * Calculate the tangent function.
 */
NUMBER *
qtan(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *sin, *cos, *tan, *res;
	long n, k, m;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for tangent");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(q);
	n = qilog2(epsilon);
	k = (n > 0) ? 4 + n/2 : 4;
	for (;;) {
		qsincos(q, 2 * k - n, &sin, &cos);
		if (qiszero(cos)) {
			qfree(sin);
			qfree(cos);
			k = 2 * k - n + 4;
			continue;
		}
		m = -qilog2(cos);
		if (m < k)
			 break;
		qfree(sin);
		qfree(cos);
		k = m + 1;
	}
	tan = qqdiv(sin, cos);
	qfree(sin);
	qfree(cos);
	res = qmappr(tan, epsilon, 24);
	qfree(tan);
	return res;
}


/*
 * Calculate the cotangent function.
 */
NUMBER *
qcot(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *sin, *cos, *cot, *res;
	long n, k, m;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for cotangent");
		/*NOTREACHED*/
	}
	if (qiszero(q)) {
		math_error("Zero argument for cotangent");
		/*NOTREACHED*/
	}
	k = -qilog2(q);
	n = qilog2(epsilon);
	if (k < 0)
		k = (n > 0) ? n/2 : 0;
	k += 4;
	for (;;) {
		qsincos(q, 2 * k - n, &sin, &cos);
		if (qiszero(sin)) {
			qfree(sin);
			qfree(cos);
			k = 2 * k - n + 4;
			continue;
		}
		m = -qilog2(sin);
		if (m < k)
			 break;
		qfree(sin);
		qfree(cos);
		k = m + 1;
	}
	cot = qqdiv(cos, sin);
	qfree(sin);
	qfree(cos);
	res = qmappr(cot, epsilon, 24);
	qfree(cot);
	return res;
}


/*
 * Calculate the secant function.
 */
NUMBER *
qsec(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *sin, *cos, *sec, *res;
	long n, k, m;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for secant");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qone_);
	n = qilog2(epsilon);
	k = (n > 0) ? 4 + n/2 : 4;
	for (;;) {
		qsincos(q, 2 * k - n, &sin, &cos);
		qfree(sin);
		if (qiszero(cos)) {
			qfree(cos);
			k = 2 * k - n + 4;
			continue;
		}
		m = -qilog2(cos);
		if (m < k)
			 break;
		qfree(cos);
		k = m + 1;
	}
	sec = qinv(cos);
	qfree(cos);
	res = qmappr(sec, epsilon, 24);
	qfree(sec);
	return res;
}


/*
 * Calculate the cosecant function.
 */
NUMBER *
qcsc(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *sin, *cos, *csc, *res;
	long n, k, m;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for cosecant");
		/*NOTREACHED*/
	}
	if (qiszero(q)) {
		math_error("Zero argument for cosecant");
		/*NOTREACHED*/
	}
	k = -qilog2(q);
	n = qilog2(epsilon);
	if (k < 0)
		k = (n > 0) ? n/2 : 0;
	k += 4;
	for (;;) {
		qsincos(q, 2 * k - n, &sin, &cos);
		qfree(cos);
		if (qiszero(sin)) {
			qfree(sin);
			k = 2 * k - n + 4;
			continue;
		}
		m = -qilog2(sin);
		if (m < k)
			 break;
		qfree(sin);
		k = m + 1;
	}
	csc = qinv(sin);
	qfree(sin);
	res = qmappr(csc, epsilon, 24);
	qfree(csc);
	return res;
}
/*
 * Calculate the arcsine function.
 * The result is in the range -pi/2 to pi/2.
 */
NUMBER *
qasin(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *qtmp1, *qtmp2, *epsilon1;
	ZVALUE ztmp;
	BOOL neg;
	FLAG r;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for asin");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qzero_);
	ztmp = q->num;
	neg = ztmp.sign;
	ztmp.sign = 0;
	r = zrel(ztmp, q->den);
	if (r > 0)
		return NULL;
	if (r == 0) {
		epsilon1 = qscale(epsilon, 1L);
		qtmp2 = qpi(epsilon1);
		qtmp1 = qscale(qtmp2, -1L);
	} else {
		epsilon1 = qscale(epsilon, -2L);
		qtmp1 = qalloc();
		zsquare(q->num, &qtmp1->num);
		zsquare(q->den, &ztmp);
		zsub(ztmp, qtmp1->num, &qtmp1->den);
		zfree(ztmp);
		qtmp2 = qsqrt(qtmp1, epsilon1, 24);
		qfree(qtmp1);
		qtmp1 = qatan(qtmp2, epsilon);
	}
	qfree(qtmp2);
	qfree(epsilon1);
	if (neg) {
		qtmp2 = qneg(qtmp1);
		qfree(qtmp1);
		return(qtmp2);
	}
	return qtmp1;
}



/*
 * Calculate the acos function.
 * The result is in the range 0 to pi.
 */
NUMBER *
qacos(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *q1, *q2, *epsilon1;
	ZVALUE z;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for acos");
		/*NOTREACHED*/
	}
	if (qisone(q))
		return qlink(&_qzero_);
	if (qisnegone(q))
		return qpi(epsilon);

	z = q->num;
	z.sign = 0;
	if (zrel(z, q->den) > 0)
		return NULL;
	epsilon1 = qscale(epsilon, -3L);		/* ??? */
	q1 = qalloc();
	zsub(q->den, q->num, &q1->num);
	zadd(q->den, q->num, &q1->den);
	q2 = qsqrt(q1, epsilon1, 24L);
	qfree(q1);
	qfree(epsilon1);
	epsilon1 = qscale(epsilon, -1L);
	q1 = qatan(q2, epsilon1);
	qfree(epsilon1);
	qfree(q2);
	q2 = qscale(q1, 1L);
	qfree(q1)
	return q2;
}


/*
 * Calculate the arctangent function to the nearest or next to nearest
 * multiple of epsilon.	 Algorithm uses
 *	atan(x) = 2 * atan(x/(1 + sqrt(1+x^2)))
 * to reduce x to a small value and then
 *	atan(x) = x - x^3/3 + ...
 */
NUMBER *
qatan(NUMBER *q, NUMBER *epsilon)
{
	long m, k, i;
	FULL d;
	ZVALUE X, D, DD, sum, mul, term, ztmp1, ztmp2;
	NUMBER *qtmp, *res;
	BOOL sign;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for arctangent");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qzero_);
	m = 12 - qilog2(epsilon);
		/* 4 bits for 4 doublings; 8 for rounding */
	if (m < 8)
		m = 8;		/* m is number of working binary digits */
	qtmp = qscale(q, m);
	zquo(qtmp->num, qtmp->den, &X, 24);
	qfree(qtmp);
	zbitvalue(m, &D);	/* q has become X/D */
	zsquare(D, &DD);
	i = 4;			/* maybe this should be larger */
	while (i-- > 0 && !ziszero(X)) {
		zsquare(X, &ztmp1);
		zadd(ztmp1, DD, &ztmp2);
		zfree(ztmp1);
		zsqrt(ztmp2, &ztmp1, 24L);
		zfree(ztmp2);
		zadd(ztmp1, D, &ztmp2);
		zfree(ztmp1);
		zshift(X, m, &ztmp1);
		zfree(X);
		zquo(ztmp1, ztmp2, &X, 24L);
		zfree(ztmp1);
		zfree(ztmp2);
	}
	zfree(DD);
	zfree(D);
	if (ziszero(X)) {
		zfree(X);
		return qlink(&_qzero_);
	}
	zcopy(X, &sum);
	zsquare(X, &ztmp1);
	zshift(ztmp1, -m, &mul);
	zfree(ztmp1);
	d = 3;
	sign = !X.sign;
	for (;;) {
		if (d > BASE) {
			math_error("Too many terms required for atan");
			/*NOTREACHED*/
		}
		zmul(X, mul, &ztmp1);
		zfree(X);
		zshift(ztmp1, -m, &X);	/* X now (original X)^d */
		zfree(ztmp1);
		zdivi(X, d, &term);
		if (ziszero(term)) {
			zfree(term);
			break;
		}
		term.sign = sign;
		zadd(sum, term, &ztmp1);
		zfree(sum);
		zfree(term);
		sum = ztmp1;
		sign = !sign;
		d += 2;
	}
	zfree(mul);