# 4.4BSD/usr/src/contrib/calc-1.26.4/commath.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 primitive routines
*/

#include <stdio.h>
#include "math.h"

COMPLEX _czero_ =		{ &_qzero_, &_qzero_, 1 };
COMPLEX _cone_ =		{ &_qone_, &_qzero_, 1 };
static COMPLEX _cnegone_ =	{ &_qnegone_, &_qzero_, 1 };

#if 0
COMPLEX _conei_ =	{ &_qzero_, &_qone_, 1 };
#endif

/*
* Free list for complex numbers.
*/
static FREELIST freelist = {
sizeof(COMPLEX),	/* size of an item */
100			/* number of free items to keep */
};

/*
*/
COMPLEX *
COMPLEX *c1, *c2;
{
COMPLEX *r;

if (ciszero(c1))
if (ciszero(c2))
r = comalloc();
if (!qiszero(c1->real) || !qiszero(c2->real))
if (!qiszero(c1->imag) || !qiszero(c2->imag))
return r;
}

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

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

/*
* Multiply two complex numbers.
* This saves one multiplication over the obvious algorithm by
*	q1 = (a + b) * (c + d)
*	q2 = a * c
*	q3 = b * d
* Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
*/
COMPLEX *
cmul(c1, c2)
COMPLEX *c1, *c2;
{
COMPLEX *r;
NUMBER *q1, *q2, *q3, *q4;

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

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

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

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

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

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

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

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

if (ciszero(c))
r = comalloc();
if (!qiszero(c->real))
r->real = qneg(c->real);
if (!qiszero(c->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 *
cint(c)
COMPLEX *c;
{
COMPLEX *r;

if (cisint(c))
r = comalloc();
r->real = qint(c->real);
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 *
cfrac(c)
COMPLEX *c;
{
COMPLEX *r;

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

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

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

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

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

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

if (cisreal(c))
r = comalloc();
return r;
}
#endif

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

if (qiszero(q))
r = comalloc();
return r;
}

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

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

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

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

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

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

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

if (qiszero(q))
if (qisone(q))
if (qisnegone(q))
return cneg(c);
r = comalloc();
r->real = qmul(c->real, q);
r->imag = qmul(c->imag, q);
return r;
}

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

if (qiszero(q))
error("Division by zero");
if (qisone(q))
if (qisnegone(q))
return cneg(c);
r = comalloc();
r->real = qdiv(c->real, q);
r->imag = qdiv(c->imag, q);
return r;
}

/*
* Take the integer quotient of a complex number by a real number.
* This is defined to be the result of doing the quotient for each component.
*/
COMPLEX *
cquoq(c, q)
COMPLEX *c;
NUMBER *q;
{
COMPLEX *r;

if (qiszero(q))
error("Division by zero");
r = comalloc();
r->real = qquo(c->real, q);
r->imag = qquo(c->imag, q);
return r;
}

/*
* Take the modulus of a complex number by a real number.
* This is defined to be the result of doing the modulo for each component.
*/
COMPLEX *
cmodq(c, q)
COMPLEX *c;
NUMBER *q;
{
COMPLEX *r;

if (qiszero(q))
error("Division by zero");
r = comalloc();
r->real = qmod(c->real, q);
r->imag = qmod(c->imag, q);
return r;
}

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

if (qiszero(q1) && qiszero(q2))
r = comalloc();
if (!qiszero(q1))
if (!qiszero(q2))
return r;
}
#endif

/*
* Compare two complex numbers for equality, returning FALSE if they are equal,
* and TRUE if they differ.
*/
BOOL
ccmp(c1, c2)
COMPLEX *c1, *c2;
{
BOOL i;

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

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

r = (COMPLEX *) allocitem(&freelist);
if (r == NULL)
error("Cannot allocate complex number");
return r;
}

/*
* Free a complex number.
*/
void
comfree(c)
COMPLEX *c;
{