drygdryg/calc
1
0
Fork 0
mirror of https://github.com/lcn2/calc.git synced 2025-08-16 01:03:29 +03:00
Code Issues Packages Projects Releases Wiki Activity
Files
c78a89386249916ac1dfad9bdad63d9790b64f02
calc/comfunc.c
Landon Curt Noll c78a893862 add exterior trigonometric functions 
Added the following new trigonometric functions:

    exsec(x [,eps])		exterior trigonometric secant
    aexsec(x [,eps])		inverse exterior trigonometric secant
    excsc(x [,eps])		exterior trigonometric cosecant
    aexcsc(x [,eps])		inverse exterior trigonometric cosecant

Added to test 95dd and test9500.trigeq.cal to the calc regression test
suite to perform extensive test of trigonometric functions.

Added to test 34dd, some if the missing inverse trigonometric tests.
2023-10-01 23:12:21 -07:00

2620 lines
51 KiB
C
Raw Blame History

/*
* 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 "errtbl.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;
}
/*
* c_tan - complex trigonometric tangent
*
* This uses the formula:
*
* tan(x) = sin(x) / cos(x)
*
* given:
* c argument to the trigonometric function
* epsilon precision of the trigonometric calculation
*
* returns:
* != NULL ==> allocated pointer to COMPLEX result
* NULL ==> invalid trigonometric argument
*
* NOTE: When the trigonometric result is returned as non-NULL result,
* the value may be a real value. The caller may wish to:
*
* COMPLEX *c; (* return result of this function *)
* NUMBER *q; (* COMPLEX result when c is a real number *)
*
* if (c == NULL) {
* math_error("... some error message");
* not_reached();
* }
* if (cisreal(c)) {
* q = c_to_q(c, ok_to_free);
* }
*/
COMPLEX *
c_tan(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *denom; /* trigonometric identity numerator */
COMPLEX *numer; /* trigonometric identity denominator */
COMPLEX *res; /* trigonometric result */
/*
* firewall - check args
*/
if (c == NULL) {
return NULL;
}
if (check_epsilon(epsilon) == false) {
return NULL;
}
/*
* evaluate the cos(x) denominator
*
* Return NULL if cos(x) failed or we otherwise divide by zero.
*/
denom = c_cos(c, epsilon);
if (denom == NULL || ciszero(denom)) {
return NULL;
}
/*
* evaluate the sin(x) numerator
*
* Return NULL if sin(x) failed.
*/
numer = c_sin(c, epsilon);
if (numer == NULL) {
comfree(denom);
return NULL;
}
/*
* catch the special case of numerator of 0
*/
if (ciszero(numer)) {
comfree(denom);
comfree(numer);
return clink(&_czero_);
}
/*
* compute the trigonometric function value
*/
res = c_div(numer, denom);
comfree(denom);
comfree(numer);
/*
* return the trigonometric result
*/
return res;
}
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;
}
/*
* c_cot - complex trigonometric cotangent
*
* This uses the formula:
*
* cot(x) = cos(x) / sin(x)
*
* given:
* c argument to the trigonometric function
* epsilon precision of the trigonometric calculation
*
* returns:
* != NULL ==> allocated pointer to COMPLEX result
* NULL ==> invalid trigonometric argument
*
* NOTE: When the trigonometric result is returned as non-NULL result,
* the value may be a real value. The caller may wish to:
*
* COMPLEX *c; (* return result of this function *)
* NUMBER *q; (* COMPLEX result when c is a real number *)
*
* if (c == NULL) {
* math_error("... some error message");
* not_reached();
* }
* if (cisreal(c)) {
* q = c_to_q(c, ok_to_free);
* }
*/
COMPLEX *
c_cot(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *denom; /* trigonometric identity numerator */
COMPLEX *numer; /* trigonometric identity denominator */
COMPLEX *res; /* trigonometric result */
/*
* firewall - check args
*/
if (c == NULL) {
return NULL;
}
if (check_epsilon(epsilon) == false) {
return NULL;
}
/*
* evaluate the sin(x) denominator
*
* Return NULL if sin(x) failed or we otherwise divide by zero.
*/
denom = c_sin(c, epsilon);
if (denom == NULL || ciszero(denom)) {
return NULL;
}
/*
* evaluate the cos(x) numerator
*
* Return NULL if cos(x) failed.
*/
numer = c_cos(c, epsilon);
if (numer == NULL) {
comfree(denom);
return NULL;
}
/*
* catch the special case of numerator of 0
*/
if (ciszero(numer)) {
comfree(denom);
comfree(numer);
return clink(&_czero_);
}
/*
* compute the trigonometric function value
*/
res = c_div(numer, denom);
comfree(denom);
comfree(numer);
/*
* return the trigonometric result
*/
return res;
}
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;
}
/*
* c_sec - complex trigonometric tangent
*
* This uses the formula:
*
* sec(x) = 1 / cos(x)
*
* given:
* c argument to the trigonometric function
* epsilon precision of the trigonometric calculation
*
* returns:
* != NULL ==> allocated pointer to COMPLEX result
* NULL ==> invalid trigonometric argument
*
* NOTE: When the trigonometric result is returned as non-NULL result,
* the value may be a real value. The caller may wish to:
*
* COMPLEX *c; (* return result of this function *)
* NUMBER *q; (* COMPLEX result when c is a real number *)
*
* if (c == NULL) {
* math_error("... some error message");
* not_reached();
* }
* if (cisreal(c)) {
* q = c_to_q(c, ok_to_free);
* }
*/
COMPLEX *
c_sec(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *denom; /* trigonometric identity numerator */
COMPLEX *res; /* trigonometric result */
/*
* firewall - check args
*/
if (c == NULL) {
return NULL;
}
if (check_epsilon(epsilon) == false) {
return NULL;
}
/*
* evaluate the cos(x) denominator
*
* Return NULL if cos(x) failed or we otherwise divide by zero.
*/
denom = c_cos(c, epsilon);
if (denom == NULL || ciszero(denom)) {
return NULL;
}
/*
* compute the trigonometric function value
*/
res = c_div(&_cone_, denom);
comfree(denom);
/*
* return the trigonometric result
*/
return res;
}
COMPLEX *
c_asec(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *tmp1, *tmp2;
tmp1 = c_inv(c);
tmp2 = c_acos(tmp1, epsilon);
comfree(tmp1);
return tmp2;
}
/*
* c_sec - complex trigonometric cosecant
*
* This uses the formula:
*
* csc(x) = 1 / sin(x)
*
* given:
* c argument to the trigonometric function
* epsilon precision of the trigonometric calculation
*
* returns:
* != NULL ==> allocated pointer to COMPLEX result
* NULL ==> invalid trigonometric argument
*
* NOTE: When the trigonometric result is returned as non-NULL result,
* the value may be a real value. The caller may wish to:
*
* COMPLEX *c; (* return result of this function *)
* NUMBER *q; (* COMPLEX result when c is a real number *)
*
* if (c == NULL) {
* math_error("... some error message");
* not_reached();
* }
* if (cisreal(c)) {
* q = c_to_q(c, ok_to_free);
* }
*/
COMPLEX *
c_csc(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *denom; /* trigonometric identity numerator */
COMPLEX *res; /* trigonometric result */
/*
* firewall - check args
*/
if (c == NULL) {
return NULL;
}
if (check_epsilon(epsilon) == false) {
return NULL;
}
/*
* evaluate the sin(x) denominator
*
* Return NULL if sin(x) failed or we otherwise divide by zero.
*/
denom = c_sin(c, epsilon);
if (denom == NULL || ciszero(denom)) {
return NULL;
}
/*
* compute the trigonometric function value
*/
res = c_div(&_cone_, denom);
comfree(denom);
/*
* return the trigonometric result
*/
return res;
}
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 - COMPLEX valued versed trigonometric sine
*
* 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 cosine for complex versin");
not_reached();
}
r = c_sub(&_cone_, ctmp);
/*
* return trigonometric result
*/
comfree(ctmp);
return r;
}
/*
* c_aversin - COMPLEX valued inverse versed trigonometric sine
*
* 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 *ctmp; /* 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
*/
ctmp = c_sub(&_cone_, c);
r = c_acos(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_coversin - COMPLEX valued coversed trigonometric sine
*
* This uses the formula:
*
* coversin(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_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 sine for complex coversin");
not_reached();
}
r = c_sub(&_cone_, ctmp);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_acoversin - COMPLEX valued inverse coversed trigonometric sine
*
* 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 *ctmp; /* 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
*/
ctmp = c_sub(&_cone_, c);
r = c_asin(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_vercos - COMPLEX valued versed trigonometric cosine
*
* This uses the formula:
*
* vercos(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_vercos(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 cosine for complex vercos");
not_reached();
}
r = c_add(&_cone_, ctmp);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_avercos - COMPLEX valued inverse versed trigonometric cosine
*
* This uses the formula:
*
* avercos(x) = acos(x - 1)
*
* 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_avercos(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* 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
*/
ctmp = c_sub(c, &_cone_);
r = c_acos(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_covercos - COMPLEX valued coversed trigonometric cosine
*
* This uses the formula:
*
* covercos(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_covercos(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 sine for complex covercos");
not_reached();
}
r = c_add(&_cone_, ctmp);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_acovercos - COMPLEX valued inverse coversed trigonometric cosine
*
* This uses the formula:
*
* acovercos(x) = asin(x - 1)
*
* 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_acovercos(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* 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
*/
ctmp = c_sub(c, &_cone_);
r = c_asin(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_haversin - COMPLEX valued half versed trigonometric sine
*
* This uses the formula:
*
* haversin(x) = versin(x) / 2
*
* 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_haversin(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_versin(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex versed sine for complex haversin");
not_reached();
}
r = c_divq(ctmp, &_qtwo_);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_ahaversin - COMPLEX valued inverse half versed trigonometric sine
*
* This uses the formula:
*
* ahaversin(x) = acos(1 - 2*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_ahaversin(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* argument to inverse trig function */
COMPLEX *x2; /* twice x */
/*
* 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
*/
x2 = c_mulq(c, &_qtwo_);
ctmp = c_sub(&_cone_, x2);
comfree(x2);
r = c_acos(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_hacoversin - COMPLEX valued half coversed trigonometric sine
*
* This uses the formula:
*
* hacoversin(x) = coversin(x) / 2
*
* 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_hacoversin(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_coversin(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex coversed sine for complex hacoversin");
not_reached();
}
r = c_divq(ctmp, &_qtwo_);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_ahacoversin - COMPLEX valued inverse half coversed trigonometric sine
*
* This uses the formula:
*
* ahacoversin(x) = asin(1 - 2*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_ahacoversin(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* argument to inverse trig function */
COMPLEX *x2; /* twice x */
/*
* 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
*/
x2 = c_mulq(c, &_qtwo_);
ctmp = c_sub(&_cone_, x2);
comfree(x2);
r = c_asin(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_havercos - COMPLEX valued half coversed trigonometric sine
*
* This uses the formula:
*
* havercos(x) = vercos(x) / 2
*
* 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_havercos(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_vercos(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex versed cosine for complex havercos");
not_reached();
}
r = c_divq(ctmp, &_qtwo_);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_ahavercos - COMPLEX valued inverse half coversed trigonometric sine
*
* This uses the formula:
*
* ahavercos(x) = acos(2*x - 1)
*
* 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_ahavercos(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* argument to inverse trig function */
COMPLEX *x2; /* twice x */
/*
* 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
*/
x2 = c_mulq(c, &_qtwo_);
ctmp = c_sub(&_cone_, x2);
comfree(x2);
r = c_acos(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_hacovercos - COMPLEX valued half coversed trigonometric cosine
*
* This uses the formula:
*
* hacovercos(x) = covercos(x) / 2
*
* 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_hacovercos(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_covercos(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex coversed cosine for complex hacovercos");
not_reached();
}
r = c_divq(ctmp, &_qtwo_);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_ahacovercos - COMPLEX valued inverse half coversed trigonometric cosine
*
* This uses the formula:
*
* ahacovercos(x) = asin(2*x - 1)
*
* 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_ahacovercos(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* argument to inverse trig function */
COMPLEX *x2; /* twice x */
/*
* 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
*/
x2 = c_mulq(c, &_qtwo_);
ctmp = c_sub(&_cone_, x2);
comfree(x2);
r = c_asin(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_exsec - COMPLEX valued exterior trigonometric secant
*
* This uses the formula:
*
* exsec(x) = sec(x) - 1
*
* 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_exsec(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* return COMPLEX value */
COMPLEX *ctmp; /* complex sec(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_sec(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex cosine for complex exsec");
not_reached();
}
r = c_sub(ctmp, &_cone_);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_aexsec - COMPLEX valued inverse exterior trigonometric secant
*
* This uses the formula:
*
* aexsec(x) = asec(x + 1)
*
* 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_aexsec(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* 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
*/
ctmp = c_addq(c, &_qone_);
r = c_asec(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
/*
* c_excsc - COMPLEX valued exterior trigonometric cosecant
*
* This uses the formula:
*
* excsc(x) = csc(x) - 1
*
* 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_excsc(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_csc(c, epsilon);
if (ctmp == NULL) {
math_error("Failed to compute complex sine for complex excsc");
not_reached();
}
r = c_sub(ctmp, &_cone_);
comfree(ctmp);
/*
* return trigonometric result
*/
return r;
}
/*
* c_aexcsc - COMPLEX valued inverse exterior trigonometric cosecant
*
* This uses the formula:
*
* aexcsc(x) = acsc(x + 1)
*
* 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_aexcsc(COMPLEX *c, NUMBER *epsilon)
{
COMPLEX *r; /* inverse trig value result */
COMPLEX *ctmp; /* 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
*/
ctmp = c_addq(c, &_qone_);
r = c_acsc(ctmp, epsilon);
comfree(ctmp);
/*
* return inverse trigonometric result
*/
return r;
}
Reference in New Issue View Git Blame Copy Permalink