mirror of
https://github.com/lcn2/calc.git
synced 2025-08-16 01:03:29 +03:00
86f1d9e029
Added new aversin(x, [,eps]) for inverse versed sine and acoversin(x, [,eps]) for inverse coversed sine. Improved trig function help files to reference use of complex arguments that while supported were not documented. Removed old Makefile testing rules for make dbx and make gdb. Improved "make run" to execute calccalc using shared libraries from the local directory, and with reading of the startup scripts disabled. Changed "make prep" to perform various tests that are used to help verify that calc is ready for a release. Added Makefile testing rule testfuncsort to check for the sort of the builtin function list. Changed the order that builtin functions are listed by "show builtin" and the help/builtin to match the sorting of "LANG=C LC_ALL=C sort -d -u".
1520 lines
29 KiB
C
1520 lines
29 KiB
C
/*
|
|
* comfunc - extended precision complex arithmetic non-primitive routines
|
|
*
|
|
* Copyright (C) 1999-2007,2021-2023 David I. Bell, Landon Curt Noll 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.
|
|
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
*
|
|
* 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"
|
|
|
|
|
|
#include "attribute.h"
|
|
#include "banned.h" /* include after system header <> includes */
|
|
|
|
|
|
/*
|
|
* cache the natural logarithm of 10
|
|
*/
|
|
STATIC COMPLEX *cln_10 = NULL;
|
|
STATIC NUMBER *cln_10_epsilon = NULL;
|
|
STATIC COMPLEX *cln_2 = NULL;
|
|
STATIC NUMBER *cln_2_epsilon = NULL;
|
|
STATIC NUMBER _q10_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
STATIC NUMBER _q2_ = { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
STATIC NUMBER _q0_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
|
|
COMPLEX _cten_ = { &_q10_, &_q0_, 1 };
|
|
COMPLEX _ctwo_ = { &_q2_, &_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");
|
|
not_reached();
|
|
}
|
|
if (zge31b(q->num)) {
|
|
math_error("Raising number to very large power");
|
|
not_reached();
|
|
}
|
|
power = ztolong(q->num);
|
|
if (ciszero(c) && (power == 0)) {
|
|
math_error("Raising zero to zeroth power");
|
|
not_reached();
|
|
}
|
|
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");
|
|
not_reached();
|
|
}
|
|
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 (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon value for complex exp");
|
|
not_reached();
|
|
}
|
|
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 (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon value for complex ln");
|
|
not_reached();
|
|
}
|
|
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);
|
|
/* quick return for 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) == true) {
|
|
/* replaced cached 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 base 2 logarithm by:
|
|
*
|
|
* log(c) = ln(c) / ln(2)
|
|
*/
|
|
COMPLEX *
|
|
c_log2(COMPLEX *c, NUMBER *epsilon)
|
|
{
|
|
int need_new_cln_2 = true; /* false => use cached cln_2 value */
|
|
COMPLEX *ln_c; /* ln(x) */
|
|
COMPLEX *ret; /* base 2 of x */
|
|
|
|
/*
|
|
* compute ln(c) first
|
|
*/
|
|
ln_c = c_ln(c, epsilon);
|
|
/* quick return for log(1) == 0 */
|
|
if (ciszero(ln_c)) {
|
|
return ln_c;
|
|
}
|
|
|
|
/*
|
|
* save epsilon for ln(2) if needed
|
|
*/
|
|
if (cln_2_epsilon == NULL) {
|
|
/* first time call */
|
|
cln_2_epsilon = qcopy(epsilon);
|
|
} else if (qcmp(cln_2_epsilon, epsilon) == true) {
|
|
/* replaced cached value with epsilon arg */
|
|
qfree(cln_2_epsilon);
|
|
cln_2_epsilon = qcopy(epsilon);
|
|
} else if (cln_2 != NULL) {
|
|
/* the previously computed ln(2) is OK to use */
|
|
need_new_cln_2 = false;
|
|
}
|
|
|
|
/*
|
|
* compute ln(2) if needed
|
|
*/
|
|
if (need_new_cln_2 == true) {
|
|
if (cln_2 != NULL) {
|
|
comfree(cln_2);
|
|
}
|
|
cln_2 = c_ln(&_ctwo_, cln_2_epsilon);
|
|
}
|
|
|
|
/*
|
|
* return ln(c) / ln(2)
|
|
*/
|
|
ret = c_div(ln_c, cln_2);
|
|
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 (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon value for complex cos");
|
|
not_reached();
|
|
}
|
|
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();
|
|
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 (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon value for complex sin");
|
|
not_reached();
|
|
}
|
|
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 (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon value for complex polar");
|
|
not_reached();
|
|
}
|
|
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 (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon value for complex power");
|
|
not_reached();
|
|
}
|
|
if (ciszero(c1)) {
|
|
if (cisreal(c2) && qisneg(c2->real)) {
|
|
math_error ("Non-positive real exponent of zero for complex power");
|
|
not_reached();
|
|
}
|
|
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, conf->outdigits);
|
|
qtmp = c->imag[0];
|
|
if (qiszero(&qtmp))
|
|
return;
|
|
if (conf->complex_space) {
|
|
math_chr(' ');
|
|
}
|
|
if (!qiszero(c->real) && !qisneg(&qtmp))
|
|
math_chr('+');
|
|
if (qisneg(&qtmp)) {
|
|
math_chr('-');
|
|
qtmp.num.sign = 0;
|
|
}
|
|
if (conf->complex_space) {
|
|
math_chr(' ');
|
|
}
|
|
qprintnum(&qtmp, MODE_DEFAULT, conf->outdigits);
|
|
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)) {
|
|
if (conf->fraction_space) {
|
|
math_chr(' ');
|
|
}
|
|
math_chr('/');
|
|
if (conf->fraction_space) {
|
|
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;
|
|
}
|
|
|
|
|
|
/*
|
|
* c_versin - versed sine for COMPLEX values
|
|
*
|
|
* This uses the formula:
|
|
*
|
|
* versin(x) = 1 - cos(x)
|
|
*
|
|
* given:
|
|
* q complex value to pass to the trig function
|
|
* epsilon error tolerance / precision for trig calculation
|
|
*
|
|
* returns:
|
|
* complex value result of trig function on q with error epsilon
|
|
*/
|
|
COMPLEX *
|
|
c_versin(COMPLEX *c, NUMBER *epsilon)
|
|
{
|
|
COMPLEX *r; /* return COMPLEX value */
|
|
COMPLEX *ctmp; /* complex cos(c) */
|
|
|
|
/*
|
|
* firewall
|
|
*/
|
|
if (c == NULL) {
|
|
math_error("%s: c is NULL", __func__);
|
|
not_reached();
|
|
}
|
|
if (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon arg for %s", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/*
|
|
* calculate complex trig function value
|
|
*/
|
|
ctmp = c_cos(c, epsilon);
|
|
if (ctmp == NULL) {
|
|
math_error("Failed to compute complex cos for complex versin");
|
|
not_reached();
|
|
}
|
|
r = c_sub(&_cone_, ctmp);
|
|
|
|
/*
|
|
* return complex 1 - cos(x)
|
|
*/
|
|
comfree(ctmp);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* c_aversin - inverse versed sine for COMPLEX values
|
|
*
|
|
* This uses the formula:
|
|
*
|
|
* aversin(x) = acos(1 - x)
|
|
*
|
|
* given:
|
|
* q complex value to pass to the trig function
|
|
* epsilon error tolerance / precision for trig calculation
|
|
*
|
|
* returns:
|
|
* complex value result of trig function on q with error epsilon
|
|
*/
|
|
COMPLEX *
|
|
c_aversin(COMPLEX *c, NUMBER *epsilon)
|
|
{
|
|
COMPLEX *r; /* inverse trig value result */
|
|
COMPLEX *x; /* argument to inverse trig function */
|
|
|
|
/*
|
|
* firewall
|
|
*/
|
|
if (c == NULL) {
|
|
math_error("%s: c is NULL", __func__);
|
|
not_reached();
|
|
}
|
|
if (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon arg for %s", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/*
|
|
* calculate complex inverse trig function value
|
|
*/
|
|
x = c_sub(&_cone_, c);
|
|
r = c_acos(x, epsilon);
|
|
comfree(x);
|
|
|
|
/*
|
|
* return complex acos(1 - x)
|
|
*/
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* c_coversin - coversed sine for COMPLEX values
|
|
*
|
|
* This uses the formula:
|
|
*
|
|
* coversin(x) = 1 - sin(x)
|
|
*
|
|
* given:
|
|
* q complex value to pass to the trig function
|
|
* epsilon error tolerance / precision for trig calculation
|
|
*
|
|
* returns:
|
|
* complex value result of trig function on q with error epsilon
|
|
*/
|
|
COMPLEX *
|
|
c_coversin(COMPLEX *c, NUMBER *epsilon)
|
|
{
|
|
COMPLEX *r; /* return COMPLEX value */
|
|
COMPLEX *ctmp; /* complex sin(c) */
|
|
|
|
/*
|
|
* firewall
|
|
*/
|
|
if (c == NULL) {
|
|
math_error("%s: c is NULL", __func__);
|
|
not_reached();
|
|
}
|
|
if (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon arg for %s", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/*
|
|
* calculate complex trig function value
|
|
*/
|
|
ctmp = c_sin(c, epsilon);
|
|
if (ctmp == NULL) {
|
|
math_error("Failed to compute complex sin for complex coversin");
|
|
not_reached();
|
|
}
|
|
r = c_sub(&_cone_, ctmp);
|
|
|
|
/*
|
|
* return complex 1 - sin(x)
|
|
*/
|
|
comfree(ctmp);
|
|
return r;
|
|
}
|
|
|
|
|
|
/*
|
|
* c_acoversin - inverse versed sine for COMPLEX values
|
|
*
|
|
* This uses the formula:
|
|
*
|
|
* acoversin(x) = asin(1 - x)
|
|
*
|
|
* given:
|
|
* q complex value to pass to the trig function
|
|
* epsilon error tolerance / precision for trig calculation
|
|
*
|
|
* returns:
|
|
* complex value result of trig function on q with error epsilon
|
|
*/
|
|
COMPLEX *
|
|
c_acoversin(COMPLEX *c, NUMBER *epsilon)
|
|
{
|
|
COMPLEX *r; /* inverse trig value result */
|
|
COMPLEX *x; /* argument to inverse trig function */
|
|
|
|
/*
|
|
* firewall
|
|
*/
|
|
if (c == NULL) {
|
|
math_error("%s: c is NULL", __func__);
|
|
not_reached();
|
|
}
|
|
if (check_epsilon(epsilon) == false) {
|
|
math_error("Invalid epsilon arg for %s", __func__);
|
|
not_reached();
|
|
}
|
|
|
|
/*
|
|
* calculate complex inverse trig function value
|
|
*/
|
|
x = c_sub(&_cone_, c);
|
|
r = c_asin(x, epsilon);
|
|
comfree(x);
|
|
|
|
/*
|
|
* return complex asin(1 - x)
|
|
*/
|
|
return r;
|
|
}
|