commath.c

/*
 * commath - extended precision complex arithmetic primitive routines
 *
 * Copyright (C) 1999-2007,2021-2023  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:10
 * File existed as early as:    before 1990
 *
 * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
 */


#include "cmath.h"


#include "errtbl.h"
#include "banned.h"     /* include after system header <> includes */


COMPLEX _czero_ =               { &_qzero_, &_qzero_, 1 };
COMPLEX _cone_ =                { &_qone_, &_qzero_, 1 };
COMPLEX _conei_ =               { &_qzero_, &_qone_, 1 };

STATIC COMPLEX _cnegone_ =      { &_qnegone_, &_qzero_, 1 };


/*
 * cmappr - complex multiple approximation
 *
 * Approximate a number to nearest multiple of a given real number. Whether
 * rounding is down, up, etc. is determined by rnd.
 *
 * This function is useful to round a result to the nearest epsilon:
 *
 *      COMPLEX *c;             (* complex number to round to nearest epsilon *)
 *      NUMBER *eps;            (* epsilon rounding precision *)
 *      COMPLEX *res;           (* c rounded to nearest epsilon *)
 *      long rnd = 24L;         (* a common rounding mode *)
 *      bool ok_to_free;        (* true ==> free c, false ==> do not free c *)
 *
 *      ...
 *
 *      res = cmappr(c, eps, ok_to_free);
 *
 * given:
 *      c       pointer to COMPLEX value to round
 *      e       pointer to NUMBER multiple
 *      rnd     rounding mode
 *      cfree   true ==> free c, false ==> do not free c
 *
 * returns:
 *      allocated pointer to COMPLEX multiple of e approximation of c
 */
COMPLEX *
cmappr(COMPLEX *c, NUMBER *e, long rnd, bool cfree)
{
        COMPLEX *r;     /* COMPLEX multiple of e approximation of c */

        /*
         * firewall
         */
        if (c == NULL) {
                math_error("%s: c is NULL", __func__);
                not_reached();
        }
        if (e == NULL) {
                math_error("%s: e is NULL", __func__);
                not_reached();
        }

        /*
         * allocate return result
         */
        r = comalloc();

        /*
         * round c to multiple of e
         */
        qfree(r->real);
        r->real = qmappr(c->real, e, rnd);
        qfree(r->imag);
        r->imag = qmappr(c->imag, e, rnd);

        /*
         * free c if requested
         */
        if (cfree == true) {
                comfree(c);
        }

        /*
          * return the allocated multiple of e approximation of c
         */
        return r;
}


/*
 * Add two complex numbers.
 */
COMPLEX *
c_add(COMPLEX *c1, COMPLEX *c2)
{
        COMPLEX *r;

        if (ciszero(c1))
                return clink(c2);
        if (ciszero(c2))
                return clink(c1);
        r = comalloc();
        if (!qiszero(c1->real) || !qiszero(c2->real)) {
                qfree(r->real);
                r->real = qqadd(c1->real, c2->real);
        }
        if (!qiszero(c1->imag) || !qiszero(c2->imag)) {
                qfree(r->imag);
                r->imag = qqadd(c1->imag, c2->imag);
        }
        return r;
}


/*
 * Subtract two complex numbers.
 */
COMPLEX *
c_sub(COMPLEX *c1, COMPLEX *c2)
{
        COMPLEX *r;

        if ((c1->real == c2->real) && (c1->imag == c2->imag))
                return clink(&_czero_);
        if (ciszero(c2))
                return clink(c1);
        r = comalloc();
        if (!qiszero(c1->real) || !qiszero(c2->real)) {
                qfree(r->real);
                r->real = qsub(c1->real, c2->real);
        }
        if (!qiszero(c1->imag) || !qiszero(c2->imag)) {
                qfree(r->imag);
                r->imag = qsub(c1->imag, c2->imag);
        }
        return r;
}


/*
 * Multiply two complex numbers.
 * This saves one multiplication over the obvious algorithm by
 * trading it for several extra additions, as follows.  Let
 *      q1 = (a + b) * (c + d)
 *      q2 = a * c
 *      q3 = b * d
 * Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
 */
COMPLEX *
c_mul(COMPLEX *c1, COMPLEX *c2)
{
        COMPLEX *r;
        NUMBER *q1, *q2, *q3, *q4;

        if (ciszero(c1) || ciszero(c2))
                return clink(&_czero_);
        if (cisone(c1))
                return clink(c2);
        if (cisone(c2))
                return clink(c1);
        if (cisreal(c2))
                return c_mulq(c1, c2->real);
        if (cisreal(c1))
                return c_mulq(c2, c1->real);
        /*
         * Need to do the full calculation.
         */
        r = comalloc();
        q2 = qqadd(c1->real, c1->imag);
        q3 = qqadd(c2->real, c2->imag);
        q1 = qmul(q2, q3);
        qfree(q2);
        qfree(q3);
        q2 = qmul(c1->real, c2->real);
        q3 = qmul(c1->imag, c2->imag);
        q4 = qqadd(q2, q3);
        qfree(r->real);
        r->real = qsub(q2, q3);
        qfree(r->imag);
        r->imag = qsub(q1, q4);
        qfree(q1);
        qfree(q2);
        qfree(q3);
        qfree(q4);
        return r;
}


/*
 * Square a complex number.
 */
COMPLEX *
c_square(COMPLEX *c)
{
        COMPLEX *r;
        NUMBER *q1, *q2;

        if (ciszero(c))
                return clink(&_czero_);
        if (cisrunit(c))
                return clink(&_cone_);
        if (cisiunit(c))
                return clink(&_cnegone_);
        r = comalloc();
        if (cisreal(c)) {
                qfree(r->real);
                r->real = qsquare(c->real);
                return r;
        }
        if (cisimag(c)) {
                qfree(r->real);
                q1 = qsquare(c->imag);
                r->real = qneg(q1);
                qfree(q1);
                return r;
        }
        q1 = qsquare(c->real);
        q2 = qsquare(c->imag);
        qfree(r->real);
        r->real = qsub(q1, q2);
        qfree(q1);
        qfree(q2);
        qfree(r->imag);
        q1 = qmul(c->real, c->imag);
        r->imag = qscale(q1, 1L);
        qfree(q1);
        return r;
}


/*
 * Divide two complex numbers.
 */
COMPLEX *
c_div(COMPLEX *c1, COMPLEX *c2)
{
        COMPLEX *r;
        NUMBER *q1, *q2, *q3, *den;

        if (ciszero(c2)) {
                math_error("Division by zero");
                not_reached();
        }
        if ((c1->real == c2->real) && (c1->imag == c2->imag))
                return clink(&_cone_);
        r = comalloc();
        if (cisreal(c1) && cisreal(c2)) {
                qfree(r->real);
                r->real = qqdiv(c1->real, c2->real);
                return r;
        }
        if (cisimag(c1) && cisimag(c2)) {
                qfree(r->real);
                r->real = qqdiv(c1->imag, c2->imag);
                return r;
        }
        if (cisimag(c1) && cisreal(c2)) {
                qfree(r->imag);
                r->imag = qqdiv(c1->imag, c2->real);
                return r;
        }
        if (cisreal(c1) && cisimag(c2)) {
                qfree(r->imag);
                q1 = qqdiv(c1->real, c2->imag);
                r->imag = qneg(q1);
                qfree(q1);
                return r;
        }
        if (cisreal(c2)) {
                qfree(r->real);
                qfree(r->imag);
                r->real = qqdiv(c1->real, c2->real);
                r->imag = qqdiv(c1->imag, c2->real);
                return r;
        }
        q1 = qsquare(c2->real);
        q2 = qsquare(c2->imag);
        den = qqadd(q1, q2);
        qfree(q1);
        qfree(q2);
        q1 = qmul(c1->real, c2->real);
        q2 = qmul(c1->imag, c2->imag);
        q3 = qqadd(q1, q2);
        qfree(q1);
        qfree(q2);
        qfree(r->real);
        r->real = qqdiv(q3, den);
        qfree(q3);
        q1 = qmul(c1->real, c2->imag);
        q2 = qmul(c1->imag, c2->real);
        q3 = qsub(q2, q1);
        qfree(q1);
        qfree(q2);
        qfree(r->imag);
        r->imag = qqdiv(q3, den);
        qfree(q3);
        qfree(den);
        return r;
}


/*
 * Invert a complex number.
 */
COMPLEX *
c_inv(COMPLEX *c)
{
        COMPLEX *r;
        NUMBER *q1, *q2, *den;

        if (ciszero(c)) {
                math_error("Inverting zero");
                not_reached();
        }
        r = comalloc();
        if (cisreal(c)) {
                qfree(r->real);
                r->real = qinv(c->real);
                return r;
        }
        if (cisimag(c)) {
                q1 = qinv(c->imag);
                qfree(r->imag);
                r->imag = qneg(q1);
                qfree(q1);
                return r;
        }
        q1 = qsquare(c->real);
        q2 = qsquare(c->imag);
        den = qqadd(q1, q2);
        qfree(q1);
        qfree(q2);
        qfree(r->real);
        r->real = qqdiv(c->real, den);
        q1 = qqdiv(c->imag, den);
        qfree(r->imag);
        r->imag = qneg(q1);
        qfree(q1);
        qfree(den);
        return r;
}


/*
 * Negate a complex number.
 */
COMPLEX *
c_neg(COMPLEX *c)
{
        COMPLEX *r;

        if (ciszero(c))
                return clink(&_czero_);
        r = comalloc();
        if (!qiszero(c->real)) {
                qfree(r->real);
                r->real = qneg(c->real);
        }
        if (!qiszero(c->imag)) {
                qfree(r->imag);
                r->imag = qneg(c->imag);
        }
        return r;
}


/*
 * Take the integer part of a complex number.
 * This means take the integer part of both components.
 */
COMPLEX *
c_int(COMPLEX *c)
{
        COMPLEX *r;

        if (cisint(c))
                return clink(c);
        r = comalloc();
        qfree(r->real);
        r->real = qint(c->real);
        qfree(r->imag);
        r->imag = qint(c->imag);
        return r;
}


/*
 * Take the fractional part of a complex number.
 * This means take the fractional part of both components.
 */
COMPLEX *
c_frac(COMPLEX *c)
{
        COMPLEX *r;

        if (cisint(c))
                return clink(&_czero_);
        r = comalloc();
        qfree(r->real);
        r->real = qfrac(c->real);
        qfree(r->imag);
        r->imag = qfrac(c->imag);
        return r;
}


/*
 * Take the conjugate of a complex number.
 * This negates the complex part.
 */
COMPLEX *
c_conj(COMPLEX *c)
{
        COMPLEX *r;

        if (cisreal(c))
                return clink(c);
        r = comalloc();
        if (!qiszero(c->real)) {
                qfree(r->real);
                r->real = qlink(c->real);
        }
        qfree(r->imag);
        r->imag = qneg(c->imag);
        return r;
}


/*
 * Return the real part of a complex number.
 */
COMPLEX *
c_real(COMPLEX *c)
{
        COMPLEX *r;

        if (cisreal(c))
                return clink(c);
        r = comalloc();
        if (!qiszero(c->real)) {
                qfree(r->real);
                r->real = qlink(c->real);
        }
        return r;
}


/*
 * c_to_q - convert a real part of a COMPLEX to a NUMBER
 *
 * given:
 *      c       complex number for which the real part will be used
 *      cfree   true ==> free c, false ==> do not free c
 *
 * returns:
 *      allocated NUMBER that the equivalent of the real part of a complex number
 *
 * NOTE: Any imaginary part of the COMPLEX value is ignored.
 *
 * NOTE: To avoid a loss of value, test with cisreal(c) first:
 *
 *      COMPLEX *c;
 *      NUMBER *q;
 *      bool ok_to_free;
 *
 *      if (cisreal(c)) {
 *              q = c_to_q(c, ok_to_free);
 *      }
 */
NUMBER *
c_to_q(COMPLEX *c, bool cfree)
{
        NUMBER *r;      /* allocated NUMBER equivalent to return */

        /*
         * firewall
         */
        if (c == NULL) {
                math_error("%s: c is NULL", __func__);
                not_reached();
        }

        /*
         * allocate a new NUMBER
         */
        r = qalloc();

        /*
         * link in the real part of the COMPLEX value
         */
        r = qlink(c->real);

        /*
         * free c if requested
         */
        if (cfree == true) {
                comfree(c);
        }

        /*
         * return the allocated equivalent NUMBER
         */
        return r;
}


/*
 * q_to_c - convert a NUMBER into an allocated COMPLEX
 *
 * given:
 *      q       NUMBER to be converted
 *
 * returns:
 *      allocated COMPLEX number whose real part is NUMBER and imag part is 0
 */
COMPLEX *
q_to_c(NUMBER *q)
{
        COMPLEX *res;   /* COMPLEX number to return */

        /*
         * allocate complex number
         */
        res = comalloc();

        /*
         * assign NUMBER to real part
         */
        qfree(res->real);
        res->real = qlink(q);

        /*
         * return the allocated equivalent COMPLEX
         */
        return res;
}


/*
 * Return the imaginary part of a complex number as a real.
 */
COMPLEX *
c_imag(COMPLEX *c)
{
        COMPLEX *r;

        if (cisreal(c))
                return clink(&_czero_);
        r = comalloc();
        qfree(r->real);
        r->real = qlink(c->imag);
        return r;
}


/*
 * Add a real number to a complex number.
 */
COMPLEX *
c_addq(COMPLEX *c, NUMBER *q)
{
        COMPLEX *r;

        if (qiszero(q))
                return clink(c);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qqadd(c->real, q);
        r->imag = qlink(c->imag);
        return r;
}


/*
 * Subtract a real number from a complex number.
 */
COMPLEX *
c_subq(COMPLEX *c, NUMBER *q)
{
        COMPLEX *r;

        if (qiszero(q))
                return clink(c);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qsub(c->real, q);
        r->imag = qlink(c->imag);
        return r;
}


/*
 * Shift the components of a complex number left by the specified
 * number of bits.  Negative values shift to the right.
 */
COMPLEX *
c_shift(COMPLEX *c, long n)
{
        COMPLEX *r;

        if (ciszero(c) || (n == 0))
                return clink(c);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qshift(c->real, n);
        r->imag = qshift(c->imag, n);
        return r;
}


/*
 * Scale a complex number by a power of two.
 */
COMPLEX *
c_scale(COMPLEX *c, long n)
{
        COMPLEX *r;

        if (ciszero(c) || (n == 0))
                return clink(c);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qscale(c->real, n);
        r->imag = qscale(c->imag, n);
        return r;
}


/*
 * Multiply a complex number by a real number.
 */
COMPLEX *
c_mulq(COMPLEX *c, NUMBER *q)
{
        COMPLEX *r;

        if (qiszero(q))
                return clink(&_czero_);
        if (qisone(q))
                return clink(c);
        if (qisnegone(q))
                return c_neg(c);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qmul(c->real, q);
        r->imag = qmul(c->imag, q);
        return r;
}


/*
 * Divide a complex number by a real number.
 */
COMPLEX *
c_divq(COMPLEX *c, NUMBER *q)
{
        COMPLEX *r;

        if (qiszero(q)) {
                math_error("Division by zero");
                not_reached();
        }
        if (qisone(q))
                return clink(c);
        if (qisnegone(q))
                return c_neg(c);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qqdiv(c->real, q);
        r->imag = qqdiv(c->imag, q);
        return r;
}




/*
 * Construct a complex number given the real and imaginary components.
 */
COMPLEX *
qqtoc(NUMBER *q1, NUMBER *q2)
{
        COMPLEX *r;

        if (qiszero(q1) && qiszero(q2))
                return clink(&_czero_);
        r = comalloc();
        qfree(r->real);
        qfree(r->imag);
        r->real = qlink(q1);
        r->imag = qlink(q2);
        return r;
}


/*
 * Compare two complex numbers for equality, returning false if they are equal,
 * and true if they differ.
 */
bool
c_cmp(COMPLEX *c1, COMPLEX *c2)
{
        bool i;

        i = qcmp(c1->real, c2->real);
        if (!i)
                i = qcmp(c1->imag, c2->imag);
        return i;
}


/*
 * Compare two complex numbers and return a complex number with real and
 * imaginary parts -1, 0 or 1 indicating relative values of the real and
 * imaginary parts of the two numbers.
 */
COMPLEX *
c_rel(COMPLEX *c1, COMPLEX *c2)
{
        COMPLEX *c;

        c = comalloc();
        qfree(c->real);
        qfree(c->imag);
        c->real = itoq((long) qrel(c1->real, c2->real));
        c->imag = itoq((long) qrel(c1->imag, c2->imag));

        return c;
}


/*
 * Allocate a new complex number.
 */
COMPLEX *
comalloc(void)
{
        COMPLEX *r;

        r = (COMPLEX *) malloc(sizeof(COMPLEX));
        if (r == NULL) {
                math_error("Cannot allocate complex number");
                not_reached();
        }
        r->links = 1;
        r->real = qlink(&_qzero_);
        r->imag = qlink(&_qzero_);
        return r;
}


/*
 * Free a complex number.
 */
void
comfree(COMPLEX *c)
{
        if (--(c->links) > 0)
                return;
        qfree(c->real);
        qfree(c->imag);
        free(c);
}