drygdryg/calc
1
0
Fork 0
mirror of https://github.com/lcn2/calc.git synced 2025-08-16 01:03:29 +03:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
calc/commath.c
Landon Curt Noll db77e29a23 convert ASCII TABs to ASCII SPACEs 
Converted all ASCII tabs to ASCII spaces using a 8 character
tab stop, for all files, except for all Makefiles (plus rpm.mk).
The `git diff -w` reports no changes.
2024-07-11 22:03:52 -07:00

792 lines
17 KiB
C
Raw Permalink Blame History

/*
* 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);
}
Reference in New Issue View Git Blame Copy Permalink