usr/src/contrib/calc-2.9.3t6/comfunc.c

/*
 * Copyright (c) 1993 David I. Bell
 * Permission is granted to use, distribute, or modify this source,
 * provided that this copyright notice remains intact.
 *
 * Extended precision complex arithmetic non-primitive routines
 */

#include "cmath.h"


/*
 * Round a complex number to the specified number of decimal places.
 * This simply means to round each of the components of the number.
 * Zero decimal places means round to the nearest complex integer.
 */
COMPLEX *
cround(c, places)
        COMPLEX *c;
        long places;
{
        COMPLEX *res;           /* result */

        res = comalloc();
        res->real = qround(c->real, places);
        res->imag = qround(c->imag, places);
        return res;
}


/*
 * Round a complex number to the specified number of binary decimal places.
 * This simply means to round each of the components of the number.
 * Zero binary places means round to the nearest complex integer.
 */
COMPLEX *
cbround(c, places)
        COMPLEX *c;
        long places;
{
        COMPLEX *res;           /* result */

        res = comalloc();
        res->real = qbround(c->real, places);
        res->imag = qbround(c->imag, places);
        return res;
}


/*
 * Compute the result of raising a complex number to an integer power.
 */
COMPLEX *
cpowi(c, q)
        COMPLEX *c;             /* complex number to be raised */
        NUMBER *q;              /* power to raise it to */
{
        COMPLEX *tmp, *res;     /* temporary values */
        long power;             /* power to raise to */
        unsigned long bit;      /* current bit value */
        int sign;

        if (qisfrac(q))
                math_error("Raising number to non-integral power");
        if (zisbig(q->num))
                math_error("Raising number to very large power");
        power = (zistiny(q->num) ? z1tol(q->num) : z2tol(q->num));
        if (ciszero(c) && (power == 0))
                math_error("Raising zero to zeroth power");
        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 cinv(c);
                        case 2:
                                return csquare(c);
                        case -2:
                                tmp = csquare(c);
                                res = cinv(tmp);
                                comfree(tmp);
                                return res;
                        case 3:
                                tmp = csquare(c);
                                res = cmul(c, tmp);
                                comfree(tmp);
                                return res;
                        case 4:
                                tmp = csquare(c);
                                res = csquare(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 = csquare(c);
        if (bit & power) {
                tmp = cmul(res, c);
                comfree(res);
                res = tmp;
        }
        bit >>= 1L;
        while (bit) {
                tmp = csquare(res);
                comfree(res);
                res = tmp;
                if (bit & power) {
                        tmp = cmul(res, c);
                        comfree(res);
                        res = tmp;
                }
                bit >>= 1L;
        }
        if (sign < 0) {
                tmp = cinv(res);
                comfree(res);
                res = tmp;
        }
        return res;
}


/*
 * Calculate the square root of a complex number, with each component
 * within the specified error.  If the number is a square, then the error
 * is zero.  For sqrt(a + bi), this calculates:
 *      R = sqrt(a^2 + b^2)
 *      U = sqrt((R + abs(a))/2)
 *      V = b/(2 * U)
 *      then sqrt(a + bi) = U + Vi if a >= 0,
 *      or abs(V) + sgn(b) * U  if a < 0
 */
COMPLEX *
csqrt(c, epsilon)
        COMPLEX *c;
        NUMBER *epsilon;
{
        COMPLEX *r;
        NUMBER *A, *B, *R, *U, *V, *tmp1, *tmp2, *epsilon2;
        long m, n;

        if (ciszero(c) || cisone(c))
                return clink(c);
        if (cisreal(c)) {
                r = comalloc();
                if (!qisneg(c->real)) {
                        r->real = qsqrt(c->real, epsilon);
                        return r;
                }
                tmp1 = qneg(c->real);
                r->imag = qsqrt(tmp1, epsilon);
                qfree(tmp1);
                return r;
        }

        A = qlink(c->real);
        B = qlink(c->imag);
        n = zhighbit(B->num) - zhighbit(B->den);
        if (!qiszero(A)) {
                m = zhighbit(A->num) - zhighbit(A->den);
                if (m > n)
                        n = m;
        }
        epsilon2 = qscale(epsilon, n/2);
        R = qhypot(A, B, epsilon2);
        qfree(epsilon2);
        if (qisneg(A))
                tmp1 = qsub(R, A);
        else
                tmp1 = qadd(R, A);
        qfree(A);
        tmp2 = qscale(tmp1, -1L);
        qfree(tmp1);
        U = qsqrt(tmp2, epsilon);
        qfree(tmp2);
        qfree(R);
        if (qiszero(U)) {
                qfree(B);
                qfree(U);
                return clink(&_czero_);
        }
        tmp1 = qdiv(B, U);
        V = qscale(tmp1, -1L);
        qfree(tmp1);
        r = comalloc();
        if (qisneg(c->real)) {
                if (qisneg(B)) {
                        tmp1 = qneg(U);
                        qfree(U);
                        U = tmp1;
                        tmp2 = qabs(V);
                        qfree(V);
                        V = tmp2;
                }
                r->real = V;
                r->imag = U;
        } else {
                r->real = U;
                r->imag = V;
        }
        qfree(B);
        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 *
croot(c, q, epsilon)
        COMPLEX *c;
        NUMBER *q, *epsilon;
{
        COMPLEX *r;
        NUMBER *a2pb2, *root, *tmp1, *tmp2, *epsilon2;

        if (qisneg(q) || qiszero(q) || qisfrac(q))
                math_error("Taking bad root of complex number");
        if (cisone(c) || qisone(q))
                return clink(c);
        if (qistwo(q))
                return csqrt(c, epsilon);
        r = comalloc();
        if (cisreal(c) && !qisneg(c->real)) {
                r->real = qroot(c->real, q, epsilon);
                return r;
        }
        /*
         * Calculate the root using the formula:
         *      croot(a + bi, n) =
         *              cpolar(qroot(a^2 + b^2, 2 * n), qatan2(b, a) / n).
         */
        epsilon2 = qscale(epsilon, -8L);
        tmp1 = qsquare(c->real);
        tmp2 = qsquare(c->imag);
        a2pb2 = qadd(tmp1, tmp2);
        qfree(tmp1);
        qfree(tmp2);
        tmp1 = qscale(q, 1L);
        root = qroot(a2pb2, tmp1, epsilon2);
        qfree(a2pb2);
        qfree(tmp1);
        tmp1 = qatan2(c->imag, c->real, epsilon2);
        qfree(epsilon2);
        tmp2 = qdiv(tmp1, q);
        qfree(tmp1);
        r = cpolar(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 *
cexp(c, epsilon)
        COMPLEX *c;
        NUMBER *epsilon;
{
        COMPLEX *r;
        NUMBER *tmp1, *tmp2, *epsilon2;

        if (ciszero(c))
                return clink(&_cone_);
        r = comalloc();
        if (cisreal(c)) {
                r->real = qexp(c->real, epsilon);
                return r;
        }
        epsilon2 = qscale(epsilon, -2L);
        r->real = qcos(c->imag, epsilon2);
        r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
        if (qiszero(c->real)) {
                qfree(epsilon2);
                return r;
        }
        tmp1 = qexp(c->real, epsilon2);
        qfree(epsilon2);
        tmp2 = qmul(r->real, tmp1);
        qfree(r->real);
        r->real = tmp2;
        tmp2 = qmul(r->imag, tmp1);
        qfree(r->imag);
        qfree(tmp1);
        r->imag = 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 *
cln(c, epsilon)
        COMPLEX *c;
        NUMBER *epsilon;
{
        COMPLEX *r;
        NUMBER *a2b2, *tmp1, *tmp2;

        if (ciszero(c))
                math_error("Logarithm of zero");
        if (cisone(c))
                return clink(&_czero_);
        r = comalloc();
        if (cisreal(c) && !qisneg(c->real)) {
                r->real = qln(c->real, epsilon);
                return r;
        }
        tmp1 = qsquare(c->real);
        tmp2 = qsquare(c->imag);
        a2b2 = qadd(tmp1, tmp2);
        qfree(tmp1);
        qfree(tmp2);
        tmp1 = qln(a2b2, epsilon);
        qfree(a2b2);
        r->real = qscale(tmp1, -1L);
        qfree(tmp1);
        r->imag = qatan2(c->imag, c->real, epsilon);
        return r;
}


/*
 * Calculate the complex cosine within the specified accuracy.
 * This uses the formula:
 *      cos(a + bi) = cos(a) * cosh(b) - sin(a) * sinh(b) * i.
 */
COMPLEX *
ccos(c, epsilon)
        COMPLEX *c;
        NUMBER *epsilon;
{
        COMPLEX *r;
        NUMBER *cosval, *coshval, *tmp1, *tmp2, *tmp3, *epsilon2;
        int negimag;

        if (ciszero(c))
                return clink(&_cone_);
        r = comalloc();
        if (cisreal(c)) {
                r->real = qcos(c->real, epsilon);
                return r;
        }
        if (qiszero(c->real)) {
                r->real = qcosh(c->imag, epsilon);
                return r;
        }
        epsilon2 = qscale(epsilon, -2L);
        coshval = qcosh(c->imag, epsilon2);
        cosval = qcos(c->real, epsilon2);
        negimag = !_sinisneg_;
        if (qisneg(c->imag))
                negimag = !negimag;
        r->real = qmul(cosval, coshval);
        /*
         * Calculate the imaginary part using the formula:
         *      sin(a) * sinh(b) = sqrt((1 - a^2) * (b^2 - 1)).
         */
        tmp1 = qsquare(cosval);
        qfree(cosval);
        tmp2 = qdec(tmp1);
        qfree(tmp1);
        tmp1 = qneg(tmp2);
        qfree(tmp2);
        tmp2 = qsquare(coshval);
        qfree(coshval);
        tmp3 = qdec(tmp2);
        qfree(tmp2);
        tmp2 = qmul(tmp1, tmp3);
        qfree(tmp1);
        qfree(tmp3);
        r->imag = qsqrt(tmp2, epsilon2);
        qfree(tmp2);
        qfree(epsilon2);
        if (negimag) {
                tmp1 = qneg(r->imag);
                qfree(r->imag);
                r->imag = tmp1;
        }
        return r;
}


/*
 * Calculate the complex sine within the specified accuracy.
 * This uses the formula:
 *      sin(a + bi) = sin(a) * cosh(b) + cos(a) * sinh(b) * i.
 */
COMPLEX *
csin(c, epsilon)
        COMPLEX *c;
        NUMBER *epsilon;
{
        COMPLEX *r;

        NUMBER *cosval, *coshval, *tmp1, *tmp2, *epsilon2;

        if (ciszero(c))
                return clink(&_czero_);
        r = comalloc();
        if (cisreal(c)) {
                r->real = qsin(c->real, epsilon);
                return r;
        }
        if (qiszero(c->real)) {
                r->imag = qsinh(c->imag, epsilon);
                return r;
        }
        epsilon2 = qscale(epsilon, -2L);
        coshval = qcosh(c->imag, epsilon2);
        cosval = qcos(c->real, epsilon2);
        tmp1 = qlegtoleg(cosval, epsilon2, _sinisneg_);
        r->real = qmul(tmp1, coshval);
        qfree(tmp1);
        tmp1 = qsquare(coshval);
        qfree(coshval);
        tmp2 = qdec(tmp1);
        qfree(tmp1);
        tmp1 = qsqrt(tmp2, epsilon2);
        qfree(tmp2);
        r->imag = qmul(tmp1, cosval);
        qfree(tmp1);
        qfree(cosval);
        if (qisneg(c->imag)) {
                tmp1 = qneg(r->imag);
                qfree(r->imag);
                r->imag = tmp1;
        }
        return r;
}


/*
 * 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 *
cpolar(q1, q2, epsilon)
        NUMBER *q1, *q2, *epsilon;
{
        COMPLEX *r;
        NUMBER *tmp, *epsilon2;
        long scale;

        r = comalloc();
        if (qiszero(q1) || qiszero(q2)) {
                r->real = qlink(q1);
                return r;
        }
        epsilon2 = epsilon;
        if (!qisunit(q1)) {
                scale = zhighbit(q1->num) - zhighbit(q1->den) + 1;
                if (scale > 0)
                        epsilon2 = qscale(epsilon, -scale);
        }
        r->real = qcos(q2, epsilon2);
        r->imag = qlegtoleg(r->real, epsilon2, _sinisneg_);
        if (epsilon2 != epsilon)
                qfree(epsilon2);
        if (qisone(q1))
                return r;
        tmp = qmul(r->real, q1);
        qfree(r->real);
        r->real = tmp;
        tmp = qmul(r->imag, q1);
        qfree(r->imag);
        r->imag = tmp;
        return r;
}


/*
 * Raise one complex number to the power of another one to within the
 * specified error.
 */
COMPLEX *
cpower(c1, c2, epsilon)
        COMPLEX *c1, *c2;
        NUMBER *epsilon;
{
        COMPLEX *tmp1, *tmp2;
        NUMBER *epsilon2;

        if (cisreal(c2) && qisint(c2->real))
                return cpowi(c1, c2->real);
        if (cisone(c1) || ciszero(c1))
                return clink(c1);
        epsilon2 = qscale(epsilon, -4L);
        tmp1 = cln(c1, epsilon2);
        tmp2 = cmul(tmp1, c2);
        comfree(tmp1);
        tmp1 = cexp(tmp2, epsilon);
        comfree(tmp2);
        qfree(epsilon2);
        return tmp1;
}


/*
 * Return a trivial hash value for a complex number.
 */
HASH
chash(c)
        COMPLEX *c;
{
        HASH hash;

        hash = qhash(c->real);
        if (!cisreal(c))
                hash += qhash(c->imag) * 2000029;
        return hash;
}


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

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


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

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

/* END CODE */