comfunc.c

Calc Software Package for Number Calc

/*
 * comfunc - extended precision complex arithmetic non-primitive routines
 *
 * Copyright (C) 1999  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.6 $
 * @(#) $Id: comfunc.c,v 29.6 2006/05/20 08:43:55 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/comfunc.c,v $
 *
 * Under source code control:	1990/02/15 01:48:13
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */


#include "config.h"
#include "cmath.h"

/*
 * cache the natural logarithm of 10
 */
static COMPLEX *cln_10 = NULL;
static NUMBER *cln_10_epsilon = NULL;
static NUMBER _q10_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
static NUMBER _q0_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
COMPLEX _cten_ = { &_q10_, &_q0_, 1 };


/*
 * Compute the result of raising a complex number to an integer power.
 *
 * given:
 *	c		complex number to be raised
 *	q		power to raise it to
 */
COMPLEX *
c_powi(COMPLEX *c, NUMBER *q)
{
	COMPLEX *tmp, *res;	/* temporary values */
	long power;		/* power to raise to */
	FULL bit;		/* current bit value */
	int sign;

	if (qisfrac(q)) {
		math_error("Raising number to non-integral power");
		/*NOTREACHED*/
	}
	if (zge31b(q->num)) {
		math_error("Raising number to very large power");
		/*NOTREACHED*/
	}
	power = ztolong(q->num);
	if (ciszero(c) && (power == 0)) {
		math_error("Raising zero to zeroth power");
		/*NOTREACHED*/
	}
	sign = 1;
	if (qisneg(q))
		sign = -1;
	/*
	 * Handle some low powers specially
	 */
	if (power <= 4) {
		switch ((int) (power * sign)) {
		case 0:
			return clink(&_cone_);
		case 1:
			return clink(c);
		case -1:
			return c_inv(c);
		case 2:
			return c_square(c);
		case -2:
			tmp = c_square(c);
			res = c_inv(tmp);
			comfree(tmp);
			return res;
		case 3:
			tmp = c_square(c);
			res = c_mul(c, tmp);
			comfree(tmp);
			return res;
		case 4:
			tmp = c_square(c);
			res = c_square(tmp);
			comfree(tmp);
			return res;
		}
	}
	/*
	 * Compute the power by squaring and multiplying.
	 * This uses the left to right method of power raising.
	 */
	bit = TOPFULL;
	while ((bit & power) == 0)
		bit >>= 1L;
	bit >>= 1L;
	res = c_square(c);
	if (bit & power) {
		tmp = c_mul(res, c);
		comfree(res);
		res = tmp;
	}
	bit >>= 1L;
	while (bit) {
		tmp = c_square(res);
		comfree(res);
		res = tmp;
		if (bit & power) {
			tmp = c_mul(res, c);
			comfree(res);
			res = tmp;
		}
		bit >>= 1L;
	}
	if (sign < 0) {
		tmp = c_inv(res);
		comfree(res);
		res = tmp;
	}
	return res;
}


/*
 * Calculate the square root of a complex number to specified accuracy.
 * Type of rounding of each component specified by R as for qsqrt().
 */
COMPLEX *
c_sqrt(COMPLEX *c, NUMBER *epsilon, long R)
{
	COMPLEX *r;
	NUMBER *es, *aes, *bes, *u, *v, qtemp;
	NUMBER *ntmp;
	ZVALUE g, a, b, d, aa, cc;
	ZVALUE tmp1, tmp2, tmp3, mul1, mul2;
	long s1, s2, s3, up1, up2;
	int imsign, sign;

	if (ciszero(c))
		return clink(c);
	if (cisreal(c)) {
		r = comalloc();
		if (!qisneg(c->real)) {
			qfree(r->real);
			r->real = qsqrt(c->real, epsilon, R);
			return r;
		}
		ntmp = qneg(c->real);
		qfree(r->imag);
		r->imag = qsqrt(ntmp, epsilon, R);
		qfree(ntmp);
		return r;
	}

	up1 = up2 = 0;
	sign = (R & 64) != 0;
	imsign = c->imag->num.sign;
	es = qsquare(epsilon);
	aes = qqdiv(c->real, es);
	bes = qqdiv(c->imag, es);
	qfree(es);
	zgcd(aes->den, bes->den, &g);
	zequo(bes->den, g, &tmp1);
	zmul(aes->num, tmp1, &a);
	zmul(aes->den, tmp1, &tmp2);
	zshift(tmp2, 1, &d);
	zfree(tmp1);
	zfree(tmp2);
	zequo(aes->den, g, &tmp1);
	zmul(bes->num, tmp1, &b);
	zfree(tmp1);
	zfree(g);
	qfree(aes);
	qfree(bes);
	zsquare(a, &tmp1);
	zsquare(b, &tmp2);
	zfree(b);
	zadd(tmp1, tmp2, &tmp3);
	zfree(tmp1);
	zfree(tmp2);
	if (R & 16) {
		zshift(tmp3, 4, &tmp1);
		zfree(tmp3);
		zshift(a, 2, &aa);
		zfree(a);
		s1 = zsqrt(tmp1, &cc, 16);
		zfree(tmp1);
		zadd(cc, aa, &tmp1);
		if (s1 == 0 && R & 32) {
			zmul(tmp1, d, &tmp2);
			zfree(tmp1);
			s2 = zsqrt(tmp2, &tmp3, 16);
			zfree(tmp2);
			if (s2 == 0) {
				aes = qalloc();
				zshift(d, 1, &tmp1);
				zreduce(tmp3, tmp1, &aes->num, &aes->den);
				zfree(tmp1);
				zfree(tmp3);
				zfree(aa);
				zfree(cc);
				zfree(d);
				r = comalloc();
				qtemp = *aes;
				qtemp.num.sign = sign;
				qfree(r->real);
				r->real = qmul(&qtemp, epsilon);
				qfree(aes);
				bes = qscale(r->real, 1);
				qtemp = *bes;
				qtemp.num.sign = sign ^ imsign;
				qfree(r->imag);
				r->imag = qqdiv(c->imag, &qtemp);
				qfree(bes);
				return r;
			}
			s3 = zquo(tmp3, d, &tmp1, s2 < 0);
		} else {
			s2 = zquo(tmp1, d, &tmp3, s1 ? (s1 < 0) : 16);
			zfree(tmp1);
			s3 = zsqrt(tmp3,&tmp1,(s1||s2) ? (s1<0 || s2<0) : 16);
		}
		zfree(tmp3);
		zshift(tmp1, -1, &mul1);
		if (*tmp1.v & 1)
			up1 = s1 + s2 + s3;
		else
			up1 = -1;
		zfree(tmp1);
		zsub(cc, aa, &tmp1);
		s2 = zquo(tmp1, d, &tmp2, s1 ? (s1 < 0) : 16);
		zfree(tmp1);
		s3 = zsqrt(tmp2, &tmp1, (s1 || s2) ? (s1 < 0 || s2 < 0) : 16);
		zfree(tmp2);
		zshift(tmp1, -1, &mul2);
		if (*tmp1.v & 1)
			up2 = s1 + s2 + s3;
		else
			up2 = -1;
		zfree(tmp1);
		zfree(aa);
	} else {
		s1 = zsqrt(tmp3, &cc, 0);
		zfree(tmp3);
		zadd(cc, a, &tmp1);
		if (s1 == 0 && R & 32) {
			zmul(tmp1, d, &tmp2);
			zfree(tmp1);
			s2 = zsqrt(tmp2, &tmp3, 0);
			zfree(tmp2);
			if (s2 == 0) {
				aes = qalloc();
				zreduce(tmp3, d, &aes->num, &aes->den);
				zfree(tmp3);
				zfree(a);
				zfree(cc);
				zfree(d);
				r = comalloc();
				qtemp = *aes;
				qtemp.num.sign = sign;
				qfree(r->real);
				r->real = qmul(&qtemp, epsilon);
				qfree(aes);
				bes = qscale(r->real, 1);
				qtemp = *bes;
				qtemp.num.sign = sign ^ imsign;
				qfree(r->imag);
				r->imag = qqdiv(c->imag, &qtemp);
				qfree(bes);
				return r;
			}
			s3 = zquo(tmp3, d, &mul1, 0);
		} else {
			s2 = zquo(tmp1, d, &tmp3, 0);
			zfree(tmp1);
			s3 = zsqrt(tmp3, &mul1, 0);
		}
		up1 = (s1 + s2 + s3) ? 0 : -1;
		zfree(tmp3);
		zsub(cc, a, &tmp1);
		s2 = zquo(tmp1, d, &tmp2, 0);
		zfree(tmp1);
		s3 = zsqrt(tmp2, &mul2, 0);
		up2 = (s1 + s2 + s3) ? 0 : -1;
		zfree(tmp2);
		zfree(a);
	}
	zfree(cc); zfree(d);
	if (up1 == 0) {
		if (R & 8)
			up1 = (long)((R ^ *mul1.v) & 1);
		else
			up1 = (R ^ epsilon->num.sign ^ sign) & 1;
		if (R & 2)
			up1 ^= epsilon->num.sign ^ sign;
		if (R & 4)
			up1 ^= epsilon->num.sign;
	}
	if (up2 == 0) {
		if (R & 8)
			up2 = (long)((R ^ *mul2.v) & 1);
		else
			up2 = (R ^ epsilon->num.sign ^ sign ^ imsign) & 1;
		if (R & 2)
			up2 ^= epsilon->num.sign ^ imsign ^ sign;
		if (R & 4)
			up2 ^= epsilon->num.sign;
	}
	if (up1 > 0) {
		zadd(mul1, _one_, &tmp1);
		zfree(mul1);
		mul1 = tmp1;
	}
	if (up2 > 0) {
		zadd(mul2, _one_, &tmp2);
		zfree(mul2);
		mul2 = tmp2;
	}
	if (ziszero(mul1)) {
		u = qlink(&_qzero_);
	} else {
		mul1.sign = sign ^ epsilon->num.sign;
		u = qalloc();
		zreduce(mul1, epsilon->den, &tmp2, &u->den);
		zmul(tmp2, epsilon->num, &u->num);
		zfree(tmp2);
	}
	zfree(mul1);
	if (ziszero(mul2)) {
		v = qlink(&_qzero_);
	} else {
		mul2.sign = imsign ^ sign ^ epsilon->num.sign;
		v = qalloc();
		zreduce(mul2, epsilon->den, &tmp2, &v->den);
		zmul(tmp2, epsilon->num, &v->num);
		zfree(tmp2);
	}
	zfree(mul2);
	if (qiszero(u) && qiszero(v)) {
		qfree(u);
		qfree(v);
		return clink(&_czero_);
	}
	r = comalloc();
	qfree(r->real);
	qfree(r->imag);
	r->real = u;
	r->imag = v;
	return r;
}


/*
 * Take the Nth root of a complex number, where N is a positive integer.
 * Each component of the result is within the specified error.
 */
COMPLEX *
c_root(COMPLEX *c, NUMBER *q, NUMBER *epsilon)
{
	COMPLEX *r;
	NUMBER *a2pb2, *root, *tmp1, *tmp2, *epsilon2;
	long n, m;

	if (qisneg(q) || qiszero(q) || qisfrac(q)) {
		math_error("Taking bad root of complex number");
		/*NOTREACHED*/
	}
	if (cisone(c) || qisone(q))
		return clink(c);
	if (qistwo(q))
		return c_sqrt(c, epsilon, 24L);
	if (cisreal(c) && !qisneg(c->real)) {
		tmp1 = qroot(c->real, q, epsilon);
		if (tmp1 == NULL)
			return NULL;
		r = comalloc();
		qfree(r->real);
		r->real = tmp1;
		return r;
	}
	/*
	 * Calculate the root using the formula:
	 *	c_root(a + bi, n) =
	 *		c_polar(qroot(a^2 + b^2, 2 * n), qatan2(b, a) / n).
	 */
	n = qilog2(epsilon);
	epsilon2 = qbitvalue(n - 4);
	tmp1 = qsquare(c->real);
	tmp2 = qsquare(c->imag);
	a2pb2 = qqadd(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qscale(q, 1L);
	root = qroot(a2pb2, tmp1, epsilon2);
	qfree(a2pb2);
	qfree(tmp1);
	qfree(epsilon2);
	if (root == NULL)
		return NULL;
	m = qilog2(root);
	if (m < n) {
		qfree(root);
		return clink(&_czero_);
	}
	epsilon2 = qbitvalue(n - m - 4);
	tmp1 = qatan2(c->imag, c->real, epsilon2);
	qfree(epsilon2);
	tmp2 = qqdiv(tmp1, q);
	qfree(tmp1);
	r = c_polar(root, tmp2, epsilon);
	qfree(root);
	qfree(tmp2);
	return r;
}


/*
 * Calculate the complex exponential function to the desired accuracy.
 * We use the formula:
 *	exp(a + bi) = exp(a) * (cos(b) + i * sin(b)).
 */
COMPLEX *
c_exp(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *r;
	NUMBER *sin, *cos, *tmp1, *tmp2, *epsilon1;
	long k, n;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for cexp");
		/*NOTREACHED*/
	}
	if (cisreal(c)) {
		tmp1 = qexp(c->real, epsilon);
		if (tmp1 == NULL)
			return NULL;
		r = comalloc();
		qfree(r->real);
		r->real = qexp(c->real, epsilon);
		return r;
	}
	n = qilog2(epsilon);
	epsilon1 = qbitvalue(n - 2);
	tmp1 = qexp(c->real, epsilon1);
	qfree(epsilon1);
	if (tmp1 == NULL)
		return NULL;
	if (qiszero(tmp1)) {
		qfree(tmp1);
		return clink(&_czero_);
	}
	k = qilog2(tmp1) + 1;
	if (k < n) {
		qfree(tmp1);
		return clink(&_czero_);
	}
	qsincos(c->imag, k - n + 2, &sin, &cos);
	tmp2 = qmul(tmp1, cos);
	qfree(cos);
	r = comalloc();
	qfree(r->real);
	r->real = qmappr(tmp2, epsilon, 24L);
	qfree(tmp2);
	tmp2 = qmul(tmp1, sin);
	qfree(tmp1);
	qfree(sin);
	qfree(r->imag);
	r->imag = qmappr(tmp2, epsilon, 24L);
	qfree(tmp2);
	return r;
}


/*
 * Calculate the natural logarithm of a complex number within the specified
 * error.  We use the formula:
 *	ln(a + bi) = ln(a^2 + b^2) / 2 + i * atan2(b, a).
 */
COMPLEX *
c_ln(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *r;
	NUMBER *a2b2, *tmp1, *tmp2, *epsilon1;

	if (ciszero(c)) {
		math_error("logarithm of zero");
		/*NOTREACHED*/
	}
	if (cisone(c))
		return clink(&_czero_);
	r = comalloc();
	if (cisreal(c) && !qisneg(c->real)) {
		qfree(r->real);
		r->real = qln(c->real, epsilon);
		return r;
	}
	tmp1 = qsquare(c->real);
	tmp2 = qsquare(c->imag);
	a2b2 = qqadd(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	epsilon1 = qscale(epsilon, 1L);
	tmp1 = qln(a2b2, epsilon1);
	qfree(a2b2);
	qfree(epsilon1);
	qfree(r->real);
	r->real = qscale(tmp1, -1L);
	qfree(tmp1);
	qfree(r->imag);
	r->imag = qatan2(c->imag, c->real, epsilon);
	return r;
}

/*
 * Calculate base 10 logarithm by:
 *
 *	log(c) = ln(c) / ln(10)
 */
COMPLEX *
c_log(COMPLEX *c, NUMBER *epsilon)
{
	int need_new_cln_10 = TRUE;	/* FALSE => use cached cln_10 value */
	COMPLEX *ln_c;			/* ln(x) */
	COMPLEX *ret;			/* base 10 logarithm of x */

	/*
	 * compute ln(c) first
	 */
	ln_c = c_ln(c, epsilon);
	/* log(1) == 0 */
	if (ciszero(ln_c)) {
		return ln_c;
	}

	/*
	 * save epsilon for ln(10) if needed
	 */
	if (cln_10_epsilon == NULL) {
		/* first time call */
		cln_10_epsilon = qcopy(epsilon);
	} else if (qcmp(cln_10_epsilon, epsilon) == FALSE) {
		/* replaced cacheed value with epsilon arg */
		qfree(cln_10_epsilon);
		cln_10_epsilon = qcopy(epsilon);
	} else if (cln_10 != NULL) {
		/* the previously computed ln(2) is OK to use */
		need_new_cln_10 = FALSE;
	}

	/*
	 * compute ln(10) if needed
	 */
	if (need_new_cln_10 == TRUE) {
		if (cln_10 != NULL) {
			comfree(cln_10);
		}
		cln_10 = c_ln(&_cten_, cln_10_epsilon);
	}

	/*
	 * return ln(c) / ln(10)
	 */
	ret = c_div(ln_c, cln_10);
	comfree(ln_c);
	return ret;
}


/*
 * Calculate the complex cosine within the specified accuracy.
 * This uses the formula:
 *	cos(x) = (exp(1i * x) + exp(-1i * x))/2;
 */
COMPLEX *
c_cos(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *r, *ctmp1, *ctmp2, *ctmp3;
	NUMBER *epsilon1;
	long n;
	BOOL neg;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for ccos");
		/*NOTREACHED*/
	}
	n = qilog2(epsilon);
	ctmp1 = comalloc();
	qfree(ctmp1->real);
	qfree(ctmp1->imag);
	neg = qisneg(c->imag);
	ctmp1->real = neg ? qneg(c->imag) : qlink(c->imag);
	ctmp1->imag = neg ? qlink(c->real) : qneg(c->real);
	epsilon1 = qbitvalue(n - 2);
	ctmp2 = c_exp(ctmp1, epsilon1);
	comfree(ctmp1);
	qfree(epsilon1);
	if (ctmp2 == NULL)
		return NULL;
	if (ciszero(ctmp2)) {
		comfree(ctmp2);
		return clink(&_czero_);
	}
	ctmp1 = c_inv(ctmp2);
	ctmp3 = c_add(ctmp2, ctmp1);
	comfree(ctmp1);
	comfree(ctmp2);
	ctmp1 = c_scale(ctmp3, -1);
	comfree(ctmp3);
	r = comalloc();