comfunc.c

Calc Software Package for Number Calc

	qfree(r->real);
	r->real = qmappr(ctmp1->real, epsilon, 24L);
	qfree(r->imag);
	r->imag = qmappr(ctmp1->imag, epsilon, 24L);
	comfree(ctmp1);
	return r;
}


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

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for csin");
		/*NOTREACHED*/
	}
	if (ciszero(c))
		return clink(&_czero_);
	n = qilog2(epsilon);
	ctmp1 = comalloc();
	neg = qisneg(c->imag);
	qfree(ctmp1->real);
	qfree(ctmp1->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_sub(ctmp2, ctmp1);
	comfree(ctmp1);
	comfree(ctmp2);
	ctmp1 = c_scale(ctmp3, -1);
	comfree(ctmp3);
	r = comalloc();
	qtmp = neg ? qlink(ctmp1->imag) : qneg(ctmp1->imag);
	qfree(r->real);
	r->real = qmappr(qtmp, epsilon, 24L);
	qfree(qtmp);
	qtmp = neg ? qneg(ctmp1->real) : qlink(ctmp1->real);
	qfree(r->imag);
	r->imag = qmappr(qtmp, epsilon, 24L);
	qfree(qtmp);
	comfree(ctmp1);
	return r;
}


COMPLEX *
c_cosh(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;

	tmp1 = c_exp(c, epsilon);
	if (tmp1 == NULL)
		return NULL;
	tmp2 = c_neg(c);
	tmp3 = c_exp(tmp2, epsilon);
	comfree(tmp2);
	if (tmp3 == NULL)
		return NULL;
	tmp2 = c_add(tmp1, tmp3);
	comfree(tmp1);
	comfree(tmp3);
	tmp1 = c_scale(tmp2, -1);
	comfree(tmp2);
	return tmp1;
}


COMPLEX *
c_sinh(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;

	tmp1 = c_exp(c, epsilon);
	if (tmp1 == NULL)
		return NULL;
	tmp2 = c_neg(c);
	tmp3 = c_exp(tmp2, epsilon);
	comfree(tmp2);
	if (tmp3 == NULL)
		return NULL;
	tmp2 = c_sub(tmp1, tmp3);
	comfree(tmp1);
	comfree(tmp3);
	tmp1 = c_scale(tmp2, -1);
	comfree(tmp2);
	return tmp1;
}


COMPLEX *
c_asin(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_mul(&_conei_, c);
	tmp2 = c_asinh(tmp1, epsilon);
	comfree(tmp1);
	tmp1 = c_div(tmp2, &_conei_);
	comfree(tmp2);
	return tmp1;
}


COMPLEX *
c_acos(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_square(c);
	tmp2 = c_sub(&_cone_, tmp1);
	comfree(tmp1);
	tmp1 = c_sqrt(tmp2, epsilon, 24);
	comfree(tmp2);
	tmp2 = c_mul(&_conei_, tmp1);
	comfree(tmp1);
	tmp1 = c_add(c, tmp2);
	comfree(tmp2);
	tmp2 = c_ln(tmp1, epsilon);
	comfree(tmp1);
	tmp1 = c_div(tmp2, &_conei_);
	comfree(tmp2);
	return tmp1;
}


COMPLEX *
c_asinh(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;
	BOOL neg;

	neg = qisneg(c->real);
	tmp1 = neg ? c_neg(c) : clink(c);
	tmp2 = c_square(tmp1);
	tmp3 = c_add(&_cone_, tmp2);
	comfree(tmp2);
	tmp2 = c_sqrt(tmp3, epsilon, 24);
	comfree(tmp3);
	tmp3 = c_add(tmp2, tmp1);
	comfree(tmp1);
	comfree(tmp2);
	tmp1 = c_ln(tmp3, epsilon);
	comfree(tmp3);
	if (neg) {
		tmp2 = c_neg(tmp1);
		comfree(tmp1);
		return tmp2;
	}
	return tmp1;
}


COMPLEX *
c_acosh(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_square(c);
	tmp2 = c_sub(tmp1, &_cone_);
	comfree(tmp1);
	tmp1 = c_sqrt(tmp2, epsilon, 24);
	comfree(tmp2);
	tmp2 = c_add(c, tmp1);
	comfree(tmp1);
	tmp1 = c_ln(tmp2, epsilon);
	comfree(tmp2);
	return tmp1;
}


COMPLEX *
c_atan(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;

	if (qiszero(c->real) && qisunit(c->imag))
		return NULL;
	tmp1 = c_sub(&_conei_, c);
	tmp2 = c_add(&_conei_, c);
	tmp3 = c_div(tmp1, tmp2);
	comfree(tmp1);
	comfree(tmp2);
	tmp1 = c_ln(tmp3, epsilon);
	comfree(tmp3);
	tmp2 = c_scale(tmp1, -1);
	comfree(tmp1);
	tmp1 = c_div(tmp2, &_conei_);
	comfree(tmp2);
	return tmp1;
}


COMPLEX *
c_acot(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;

	if (qiszero(c->real) && qisunit(c->imag))
		return NULL;
	tmp1 = c_add(c, &_conei_);
	tmp2 = c_sub(c, &_conei_);
	tmp3 = c_div(tmp1, tmp2);
	comfree(tmp1);
	comfree(tmp2);
	tmp1 = c_ln(tmp3, epsilon);
	comfree(tmp3);
	tmp2 = c_scale(tmp1, -1);
	comfree(tmp1);
	tmp1 = c_div(tmp2, &_conei_);
	comfree(tmp2);
	return tmp1;
}

COMPLEX *
c_asec(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_inv(c);
	tmp2 = c_acos(tmp1, epsilon);
	comfree(tmp1);
	return tmp2;
}

COMPLEX *
c_acsc(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_inv(c);
	tmp2 = c_asin(tmp1, epsilon);
	comfree(tmp1);
	return tmp2;
}


COMPLEX *
c_atanh(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;

	if (qiszero(c->imag) && qisunit(c->real))
		return NULL;
	tmp1 = c_add(&_cone_, c);
	tmp2 = c_sub(&_cone_, c);
	tmp3 = c_div(tmp1, tmp2);
	comfree(tmp1);
	comfree(tmp2);
	tmp1 = c_ln(tmp3, epsilon);
	comfree(tmp3);
	tmp2 = c_scale(tmp1, -1);
	comfree(tmp1);
	return tmp2;
}


COMPLEX *
c_acoth(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;

	if (qiszero(c->imag) && qisunit(c->real))
		return NULL;
	tmp1 = c_add(c, &_cone_);
	tmp2 = c_sub(c, &_cone_);
	tmp3 = c_div(tmp1, tmp2);
	comfree(tmp1);
	comfree(tmp2);
	tmp1 = c_ln(tmp3, epsilon);
	comfree(tmp3);
	tmp2 = c_scale(tmp1, -1);
	comfree(tmp1);
	return tmp2;
}

COMPLEX *
c_asech(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_inv(c);
	tmp2 = c_acosh(tmp1, epsilon);
	comfree(tmp1);
	return tmp2;
}

COMPLEX *
c_acsch(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_inv(c);
	tmp2 = c_asinh(tmp1, epsilon);
	comfree(tmp1);
	return tmp2;
}


COMPLEX *
c_gd(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2, *tmp3;
	NUMBER *q1, *q2;
	NUMBER *sin, *cos;
	NUMBER *eps;
	int n, n1;
	BOOL neg;

	if (cisreal(c)) {
		q1 = qscale(c->real, -1);
		eps = qscale(epsilon, -1);
		q2 = qtanh(q1, eps);
		qfree(q1);
		q1 = qatan(q2, eps);
		qfree(eps);
		qfree(q2);
		tmp1 = comalloc();
		qfree(tmp1->real);
		tmp1->real = qscale(q1, 1);
		qfree(q1);
		return tmp1;
	}
	if (qiszero(c->real)) {
		n = - qilog2(epsilon);
		qsincos(c->imag, n + 8, &sin, &cos);
		if (qiszero(cos) || (n1 = -qilog2(cos)) > n) {
			qfree(sin);
			qfree(cos);
			return NULL;
		}
		neg = qisneg(sin);
		q1 = neg ? qsub(&_qone_, sin) : qqadd(&_qone_, sin);
		qfree(sin);
		if (n1 > 8) {
			qfree(q1);
			qfree(cos);
			qsincos(c->imag, n + n1, &sin, &cos);
			q1 = neg ? qsub(&_qone_, sin) : qqadd(&_qone_, sin);
			qfree(sin);
		}
		q2 = qqdiv(q1, cos);
		qfree(q1);
		q1 = qln(q2, epsilon);
		qfree(q2);
		if (neg) {
			q2 = qneg(q1);
			qfree(q1);
			q1 = q2;
		}
		tmp1 = comalloc();
		qfree(tmp1->imag);
		tmp1->imag = q1;
		if (qisneg(cos)) {
			qfree(tmp1->real);
			q1 = qpi(epsilon);
			if (qisneg(c->imag)) {
				q2 = qneg(q1);
				qfree(q1);
				q1 = q2;
			}
			tmp1->real = q1;
		}
		qfree(cos);
		return tmp1;
	}
	neg = qisneg(c->real);
	tmp1 = neg ? c_neg(c) : clink(c);
	tmp2 = c_exp(tmp1, epsilon);
	comfree(tmp1);
	if (tmp2 == NULL)
		return NULL;
	tmp1 = c_mul(&_conei_, tmp2);
	tmp3 = c_add(&_conei_, tmp2);
	comfree(tmp2);
	tmp2 = c_add(tmp1, &_cone_);
	comfree(tmp1);
	if (ciszero(tmp2) || ciszero(tmp3)) {
		comfree(tmp2);
		comfree(tmp3);
		return NULL;
	}
	tmp1 = c_div(tmp2, tmp3);
	comfree(tmp2);
	comfree(tmp3);
	tmp2 = c_ln(tmp1, epsilon);
	comfree(tmp1);
	tmp1 = c_div(tmp2, &_conei_);
	comfree(tmp2);
	if (neg) {
		tmp2 = c_neg(tmp1);
		comfree(tmp1);
		return tmp2;
	}
	return tmp1;
}


COMPLEX *
c_agd(COMPLEX *c, NUMBER *epsilon)
{
	COMPLEX *tmp1, *tmp2;

	tmp1 = c_mul(&_conei_, c);
	tmp2 = c_gd(tmp1, epsilon);
	comfree(tmp1);
	if (tmp2 == NULL)
		return NULL;
	tmp1 = c_div(tmp2, &_conei_);
	comfree(tmp2);
	return tmp1;
}


/*
 * Convert a number from polar coordinates to normal complex number form
 * within the specified accuracy.  This produces the value:
 *	q1 * cos(q2) + q1 * sin(q2) * i.
 */
COMPLEX *
c_polar(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
{
	COMPLEX *r;
	NUMBER *tmp, *cos, *sin;
	long m, n;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for cpolar");
		/*NOTREACHED*/
	}
	if (qiszero(q1))
		return clink(&_czero_);
	m = qilog2(q1) + 1;
	n = qilog2(epsilon);
	if (m < n)
		return clink(&_czero_);
	r = comalloc();
	if (qiszero(q2)) {
		qfree(r->real);
		r->real = qlink(q1);
		return r;
	}
	qsincos(q2, m - n + 2, &sin, &cos);
	tmp = qmul(q1, cos);
	qfree(cos);
	qfree(r->real);
	r->real = qmappr(tmp, epsilon, 24L);
	qfree(tmp);
	tmp = qmul(q1, sin);
	qfree(sin);
	qfree(r->imag);
	r->imag = qmappr(tmp, epsilon, 24L);
	qfree(tmp);
	return r;
}


/*
 * Raise one complex number to the power of another one to within the
 * specified error.
 */
COMPLEX *
c_power(COMPLEX *c1, COMPLEX *c2, NUMBER *epsilon)
{
	COMPLEX *ctmp1, *ctmp2;
	long k1, k2, k, m1, m2, m, n;
	NUMBER *a2b2, *qtmp1, *qtmp2, *epsilon1;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for cpower");
		/*NOTREACHED*/
	}
	if (ciszero(c1)) {
		if (qisneg(c2->real) || qiszero(c2->real)) {
			math_error ("Non-positive exponent of zero");
			/*NOTREACHED*/
		}
		return clink(&_czero_);
	}
	n = qilog2(epsilon);
	m1 = m2 = -1000000;
	k1 = k2 = 0;
	if (!qiszero(c2->real)) {
		qtmp1 = qsquare(c1->real);
		qtmp2 = qsquare(c1->imag);
		a2b2 = qqadd(qtmp1, qtmp2);
		qfree(qtmp1);
		qfree(qtmp2);
		m1 = qilog2(c2->real);
		epsilon1 = qbitvalue(-m1 - 1);
		qtmp1 = qln(a2b2, epsilon1);
		qfree(epsilon1);
		qfree(a2b2);
		qtmp2 = qmul(qtmp1, c2->real);
		qfree(qtmp1);
		qtmp1 = qmul(qtmp2, &_qlge_);
		qfree(qtmp2);
		k1 = qtoi(qtmp1);
		qfree(qtmp1);
	}
	if (!qiszero(c2->imag)) {
		m2 = qilog2(c2->imag);
		epsilon1 = qbitvalue(-m2 - 1);
		qtmp1 = qatan2(c1->imag, c1->real, epsilon1);
		qfree(epsilon1);
		qtmp2 = qmul(qtmp1, c2->imag);
		qfree(qtmp1);
		qtmp1 = qscale(qtmp2, -1);
		qfree(qtmp2);
		qtmp2 = qmul(qtmp1, &_qlge_);
		qfree(qtmp1);
		k2 = qtoi(qtmp2);
		qfree(qtmp2);
	}
	m = (m2 > m1) ? m2 : m1;
	k = k1 - k2 + 1;
	if (k < n)
		return clink(&_czero_);
	epsilon1 = qbitvalue(n - k - m - 2);
	ctmp1 = c_ln(c1, epsilon1);
	qfree(epsilon1);
	ctmp2 = c_mul(ctmp1, c2);
	comfree(ctmp1);
	ctmp1 = c_exp(ctmp2, epsilon);
	comfree(ctmp2);
	return ctmp1;
}


/*
 * Print a complex number in the current output mode.
 */
void
comprint(COMPLEX *c)
{
	NUMBER qtmp;

	if (conf->outmode == MODE_FRAC) {
		cprintfr(c);
		return;
	}
	if (!qiszero(c->real) || qiszero(c->imag))
		qprintnum(c->real, MODE_DEFAULT);
	qtmp = c->imag[0];
	if (qiszero(&qtmp))
		return;
	if (!qiszero(c->real) && !qisneg(&qtmp))
		math_chr('+');
	if (qisneg(&qtmp)) {
		math_chr('-');
		qtmp.num.sign = 0;
	}
	qprintnum(&qtmp, MODE_DEFAULT);
	math_chr('i');
}


/*
 * Print a complex number in rational representation.
 * Example:  2/3-4i/5
 */
void
cprintfr(COMPLEX *c)
{
	NUMBER *r;
	NUMBER *i;

	r = c->real;
	i = c->imag;
	if (!qiszero(r) || qiszero(i))
		qprintfr(r, 0L, FALSE);
	if (qiszero(i))
		return;
	if (!qiszero(r) && !qisneg(i))
		math_chr('+');
	zprintval(i->num, 0L, 0L);
	math_chr('i');
	if (qisfrac(i)) {
		math_chr('/');
		zprintval(i->den, 0L, 0L);
	}
}


NUMBER *
c_ilog(COMPLEX *c, ZVALUE base)
{
	NUMBER *qr, *qi;

	qr = qilog(c->real, base);
	qi = qilog(c->imag, base);

	if (qr == NULL) {
		if (qi == NULL)
			return NULL;
		return qi;
	}
	if (qi == NULL)
		return qr;
	if (qrel(qr, qi) >= 0) {
		qfree(qi);
		return qr;
	}
	qfree(qr);
	return qi;
}