comfunc.c

/*
 * comfunc - extended precision complex arithmetic non-primitive routines
 *
 * Copyright (C) 1999-2007,2021-2023,2025  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, conf->triground);
    }
    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)).
 *
 * IMPORTANT NOTE: This function will return NULL when
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_EXP_3 error.
 */
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, conf->triground);
    qfree(tmp2);
    tmp2 = qmul(tmp1, sin);
    qfree(tmp1);
    qfree(sin);
    qfree(r->imag);
    r->imag = qmappr(tmp2, epsilon, conf->triground);
    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;
 *
 * IMPORTANT NOTE: This function will return NULL when
 *                 c_exp(COMPLEX *c, NUMBER *epsilon) return NULL because
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_COS_3 error.
 */
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, conf->triground);
    qfree(r->imag);
    r->imag = qmappr(ctmp1->imag, epsilon, conf->triground);
    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).
 *
 * IMPORTANT NOTE: This function will return NULL when
 *                 c_exp(COMPLEX *c, NUMBER *epsilon) return NULL because
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_SIN_3 error.
 */
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, conf->triground);
    qfree(qtmp);
    qtmp = neg ? qneg(ctmp1->real) : qlink(ctmp1->real);
    qfree(r->imag);
    r->imag = qmappr(qtmp, epsilon, conf->triground);
    qfree(qtmp);
    comfree(ctmp1);
    return r;
}

/*
 * IMPORTANT NOTE: This function will return NULL when
 *                 c_exp(COMPLEX *c, NUMBER *epsilon) return NULL because
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_COSH_3 or E_TANH_4 or E_COTH_3 or E_SECH_3 error.
 */
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;
}

/*
 * IMPORTANT NOTE: This function will return NULL when
 *                 c_exp(COMPLEX *c, NUMBER *epsilon) return NULL because
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_SINH_3 or E_TANH_3 or E_COTH_4 or E_CSCH_3 error.
 */
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, conf->triground);
    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, conf->triground);
    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, conf->triground);
    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;
}

/*
 * IMPORTANT NOTE: This function will return NULL when
 *                 c_exp(COMPLEX *c, NUMBER *epsilon) return NULL because
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_GD_3 error.
 */
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, conf->triground);
    qfree(tmp);
    tmp = qmul(q1, sin);
    qfree(sin);
    qfree(r->imag);
    r->imag = qmappr(tmp, epsilon, conf->triground);
    qfree(tmp);
    return r;
}

/*
 * Raise one complex number to the power of another one to within the
 * specified error.
 *
 * IMPORTANT NOTE: This function will return NULL when
 *                 c_exp(COMPLEX *c, NUMBER *epsilon) return NULL because
 *                 qexp(NUMBER *q, NUMBER *epsilon) returns NULL because
 *                 "Too-large re(argument) for exp" resulting
 *                 in an E_POWER_4 error.
 *
 * This function also returns NULL when epsilon is invalid,
 * or if we are given a non-positive real exponent of zero for complex power.
 */
COMPLEX *
c_power(COMPLEX *c1, COMPLEX *c2, NUMBER *epsilon)
{
    COMPLEX *ctmp1, *ctmp2;
    long k1, k2, m1, m2, m, n;
    NUMBER *a2b2, *qtmp1, *qtmp2, *epsilon1;
    NUMBER *k1_num;
    NUMBER *k2_num;
    NUMBER *tmp;
    NUMBER *tmp2;
    NUMBER *k_num;
    NUMBER *n_num;
    NUMBER *m_num;
    NUMBER *ep1_power_of_2;
    NUMBER *two_num;

    if (check_epsilon(epsilon) == false) {
        return NULL;
    }
    if (ciszero(c1)) {
        if (cisreal(c2) && qisneg(c2->real)) {
            return NULL;
        }
        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;

    /*
     * We want to, in effect, compute:
     *
     *      k = k1 - k2 + 1
     *
     * however we must be careful about "Signed integer overflow" during k1 - k2 or due to + 1.
     */
    k1_num = itoq(k1);
    k2_num = itoq(k2);
    tmp = qsub(k1_num, k2_num);
    qfree(k1_num);
    qfree(k2_num);
    k_num = qinc(tmp);
    qfree(tmp);

    /*
     * We want to, in effect, compute:
     *
     *      if (k < n)
     *              return clink(&_czero_);
     *
     * If that is true, then value is < epsilon, so return complex 0.
     */
    n_num = itoq(n);
    if (qrel(k_num, n_num) < 0) {

        /* value is < epsilon, so return complex 0 */
        qfree(k_num) qfree(n_num) return clink(&_czero_);
    }

    /*
     * We want to, in effect, compute:
     *
     *      epsilon1 = "2 to the power of (n - k - m - 2)"
     *
     * I.e.:
     *
     *      epsilon1 = qbitvalue(-(n - k - m - 2));
     *
     * where "n - k - m - 2" does not lead to a "Signed integer overflow"
     * and where "n - k - m - 2" may not fit into a long.
     */
    tmp = qsub(n_num, k_num);
    qfree(n_num);
    qfree(k_num);
    m_num = itoq(m);
    tmp2 = qsub(tmp, m_num);
    qfree(tmp);
    qfree(m_num);
    two_num = itoq(2);
    ep1_power_of_2 = qsub(tmp2, two_num);
    qfree(tmp2);
    qfree(two_num);
    epsilon1 = qqbitvalue(*ep1_power_of_2);
    qfree(ep1_power_of_2);
    if (epsilon1 == NULL) {
        return NULL;
    }

    /* ctmp1 = natural log of epsilon1 */
    ctmp1 = c_ln(c1, epsilon1);
    qfree(epsilon1);

    /* return complex result */
    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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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:
 *      c           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;
}

/*
 * c_crd - COMPLEX valued trigonometric chord of a unit circle
 *
 * This uses the formula:
 *
 *      crd(x) = 2 * sin(x / 2)
 *
 * given:
 *      c           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_crd(COMPLEX *c, NUMBER *epsilon)
{
    COMPLEX *r;     /* return COMPLEX value */
    COMPLEX *cdiv2; /* complex c/2 */
    COMPLEX *ctmp;  /* complex sin(c/2) */

    /*
     * 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
     */
    cdiv2 = c_divq(c, &_qtwo_);
    ctmp = c_sin(cdiv2, epsilon);
    comfree(cdiv2);
    if (ctmp == NULL) {
        math_error("Failed to compute complex sine for complex crd");
        not_reached();
    }
    r = c_mulq(ctmp, &_qtwo_);
    comfree(ctmp);

    /*
     * return trigonometric result
     */
    return r;
}

/*
 * c_acrd - COMPLEX valued inverse trigonometric chord of a unit circle
 *
 * This uses the formula:
 *
 *      acrd(x) = 2 * asin(x / 2)
 *
 * given:
 *      c           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_acrd(COMPLEX *c, NUMBER *epsilon)
{
    COMPLEX *r;     /* inverse trig value result */
    COMPLEX *cdiv2; /* complex c/2 */
    COMPLEX *ctmp;  /* complex asin(c/2) */

    /*
     * 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
     */
    cdiv2 = c_divq(c, &_qtwo_);
    ctmp = c_asin(cdiv2, epsilon);
    comfree(cdiv2);
    r = c_mulq(ctmp, &_qtwo_);
    comfree(ctmp);

    /*
     * return inverse trigonometric result
     */
    return r;
}

/*
 * c_cas - COMPLEX valued cosine plus sine
 *
 * This uses the formula:
 *
 *      cas(x) = cos(x) + sin(x)
 *
 * given:
 *      c           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_cas(COMPLEX *c, NUMBER *epsilon)
{
    COMPLEX *r;    /* return COMPLEX value */
    COMPLEX *csin; /* complex sin(c) */
    COMPLEX *ccos; /* 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
     */
    csin = c_sin(c, epsilon);
    if (csin == NULL) {
        math_error("Failed to compute complex sine for complex cas");
        not_reached();
    }
    ccos = c_cos(c, epsilon);
    if (ccos == NULL) {
        comfree(csin);
        math_error("Failed to compute complex cosine for complex cas");
        not_reached();
    }
    r = c_add(csin, ccos);
    comfree(csin);
    comfree(ccos);

    /*
     * return trigonometric result
     */
    return r;
}

/*
 * c_cis - COMPLEX valued Euler's formula
 *
 * This uses the formula:
 *
 *      cis(x) = cos(x) + 1i*sin(x)
 *      cis(x) = exp(1i * x)
 *
 * given:
 *      c           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_cis(COMPLEX *c, NUMBER *epsilon)
{
    COMPLEX *r;    /* return COMPLEX value */
    COMPLEX *ctmp; /* 1i * 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_mul(c, &_conei_);
    r = c_exp(ctmp, epsilon);
    comfree(ctmp);
    if (r == NULL) {
        math_error("Failed to compute complex exp for complex cis");
        not_reached();
    }

    /*
     * return trigonometric result
     */
    return r;
}