/* * 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();