c_cpp/etc/calc/commath.c

   1 /*
   2  * commath - extended precision complex arithmetic primitive routines
   3  *
   4  * Copyright (C) 1999-2007  David I. Bell
   5  *
   6  * Calc is open software; you can redistribute it and/or modify it under
   7  * the terms of the version 2.1 of the GNU Lesser General Public License
   8  * as published by the Free Software Foundation.
   9  *
  10  * Calc is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13  * Public License for more details.
  14  *
  15  * A copy of version 2.1 of the GNU Lesser General Public License is
  16  * distributed with calc under the filename COPYING-LGPL.  You should have
  17  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  * @(#) $Revision: 30.1 $
  21  * @(#) $Id: commath.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/RCS/commath.c,v $
  23  *
  24  * Under source code control:   1990/02/15 01:48:10
  25  * File existed as early as:    before 1990
  26  *
  27  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  28  */
  29
  30
  31 #include "cmath.h"
  32
  33
  34 COMPLEX _czero_ =               { &_qzero_, &_qzero_, 1 };
  35 COMPLEX _cone_ =                { &_qone_, &_qzero_, 1 };
  36 COMPLEX _conei_ =               { &_qzero_, &_qone_, 1 };
  37
  38 STATIC COMPLEX _cnegone_ =      { &_qnegone_, &_qzero_, 1 };
  39
  40
  41 /*
  42  * Add two complex numbers.
  43  */
  44 COMPLEX *
  45 c_add(COMPLEX *c1, COMPLEX *c2)
  46 {
  47         COMPLEX *r;
  48
  49         if (ciszero(c1))
  50                 return clink(c2);
  51         if (ciszero(c2))
  52                 return clink(c1);
  53         r = comalloc();
  54         if (!qiszero(c1->real) || !qiszero(c2->real)) {
  55                 qfree(r->real);
  56                 r->real = qqadd(c1->real, c2->real);
  57         }
  58         if (!qiszero(c1->imag) || !qiszero(c2->imag)) {
  59                 qfree(r->imag);
  60                 r->imag = qqadd(c1->imag, c2->imag);
  61         }
  62         return r;
  63 }
  64
  65
  66 /*
  67  * Subtract two complex numbers.
  68  */
  69 COMPLEX *
  70 c_sub(COMPLEX *c1, COMPLEX *c2)
  71 {
  72         COMPLEX *r;
  73
  74         if ((c1->real == c2->real) && (c1->imag == c2->imag))
  75                 return clink(&_czero_);
  76         if (ciszero(c2))
  77                 return clink(c1);
  78         r = comalloc();
  79         if (!qiszero(c1->real) || !qiszero(c2->real)) {
  80                 qfree(r->real);
  81                 r->real = qsub(c1->real, c2->real);
  82         }
  83         if (!qiszero(c1->imag) || !qiszero(c2->imag)) {
  84                 qfree(r->imag);
  85                 r->imag = qsub(c1->imag, c2->imag);
  86         }
  87         return r;
  88 }
  89
  90
  91 /*
  92  * Multiply two complex numbers.
  93  * This saves one multiplication over the obvious algorithm by
  94  * trading it for several extra additions, as follows.  Let
  95  *      q1 = (a + b) * (c + d)
  96  *      q2 = a * c
  97  *      q3 = b * d
  98  * Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
  99  */
 100 COMPLEX *
 101 c_mul(COMPLEX *c1, COMPLEX *c2)
 102 {
 103         COMPLEX *r;
 104         NUMBER *q1, *q2, *q3, *q4;
 105
 106         if (ciszero(c1) || ciszero(c2))
 107                 return clink(&_czero_);
 108         if (cisone(c1))
 109                 return clink(c2);
 110         if (cisone(c2))
 111                 return clink(c1);
 112         if (cisreal(c2))
 113                 return c_mulq(c1, c2->real);
 114         if (cisreal(c1))
 115                 return c_mulq(c2, c1->real);
 116         /*
 117          * Need to do the full calculation.
 118          */
 119         r = comalloc();
 120         q2 = qqadd(c1->real, c1->imag);
 121         q3 = qqadd(c2->real, c2->imag);
 122         q1 = qmul(q2, q3);
 123         qfree(q2);
 124         qfree(q3);
 125         q2 = qmul(c1->real, c2->real);
 126         q3 = qmul(c1->imag, c2->imag);
 127         q4 = qqadd(q2, q3);
 128         qfree(r->real);
 129         r->real = qsub(q2, q3);
 130         qfree(r->imag);
 131         r->imag = qsub(q1, q4);
 132         qfree(q1);
 133         qfree(q2);
 134         qfree(q3);
 135         qfree(q4);
 136         return r;
 137 }
 138
 139
 140 /*
 141  * Square a complex number.
 142  */
 143 COMPLEX *
 144 c_square(COMPLEX *c)
 145 {
 146         COMPLEX *r;
 147         NUMBER *q1, *q2;
 148
 149         if (ciszero(c))
 150                 return clink(&_czero_);
 151         if (cisrunit(c))
 152                 return clink(&_cone_);
 153         if (cisiunit(c))
 154                 return clink(&_cnegone_);
 155         r = comalloc();
 156         if (cisreal(c)) {
 157                 qfree(r->real);
 158                 r->real = qsquare(c->real);
 159                 return r;
 160         }
 161         if (cisimag(c)) {
 162                 qfree(r->real);
 163                 q1 = qsquare(c->imag);
 164                 r->real = qneg(q1);
 165                 qfree(q1);
 166                 return r;
 167         }
 168         q1 = qsquare(c->real);
 169         q2 = qsquare(c->imag);
 170         qfree(r->real);
 171         r->real = qsub(q1, q2);
 172         qfree(q1);
 173         qfree(q2);
 174         qfree(r->imag);
 175         q1 = qmul(c->real, c->imag);
 176         r->imag = qscale(q1, 1L);
 177         qfree(q1);
 178         return r;
 179 }
 180
 181
 182 /*
 183  * Divide two complex numbers.
 184  */
 185 COMPLEX *
 186 c_div(COMPLEX *c1, COMPLEX *c2)
 187 {
 188         COMPLEX *r;
 189         NUMBER *q1, *q2, *q3, *den;
 190
 191         if (ciszero(c2)) {
 192                 math_error("Division by zero");
 193                 /*NOTREACHED*/
 194         }
 195         if ((c1->real == c2->real) && (c1->imag == c2->imag))
 196                 return clink(&_cone_);
 197         r = comalloc();
 198         if (cisreal(c1) && cisreal(c2)) {
 199                 qfree(r->real);
 200                 r->real = qqdiv(c1->real, c2->real);
 201                 return r;
 202         }
 203         if (cisimag(c1) && cisimag(c2)) {
 204                 qfree(r->real);
 205                 r->real = qqdiv(c1->imag, c2->imag);
 206                 return r;
 207         }
 208         if (cisimag(c1) && cisreal(c2)) {
 209                 qfree(r->imag);
 210                 r->imag = qqdiv(c1->imag, c2->real);
 211                 return r;
 212         }
 213         if (cisreal(c1) && cisimag(c2)) {
 214                 qfree(r->imag);
 215                 q1 = qqdiv(c1->real, c2->imag);
 216                 r->imag = qneg(q1);
 217                 qfree(q1);
 218                 return r;
 219         }
 220         if (cisreal(c2)) {
 221                 qfree(r->real);
 222                 qfree(r->imag);
 223                 r->real = qqdiv(c1->real, c2->real);
 224                 r->imag = qqdiv(c1->imag, c2->real);
 225                 return r;
 226         }
 227         q1 = qsquare(c2->real);
 228         q2 = qsquare(c2->imag);
 229         den = qqadd(q1, q2);
 230         qfree(q1);
 231         qfree(q2);
 232         q1 = qmul(c1->real, c2->real);
 233         q2 = qmul(c1->imag, c2->imag);
 234         q3 = qqadd(q1, q2);
 235         qfree(q1);
 236         qfree(q2);
 237         qfree(r->real);
 238         r->real = qqdiv(q3, den);
 239         qfree(q3);
 240         q1 = qmul(c1->real, c2->imag);
 241         q2 = qmul(c1->imag, c2->real);
 242         q3 = qsub(q2, q1);
 243         qfree(q1);
 244         qfree(q2);
 245         qfree(r->imag);
 246         r->imag = qqdiv(q3, den);
 247         qfree(q3);
 248         qfree(den);
 249         return r;
 250 }
 251
 252
 253 /*
 254  * Invert a complex number.
 255  */
 256 COMPLEX *
 257 c_inv(COMPLEX *c)
 258 {
 259         COMPLEX *r;
 260         NUMBER *q1, *q2, *den;
 261
 262         if (ciszero(c)) {
 263                 math_error("Inverting zero");
 264                 /*NOTREACHED*/
 265         }
 266         r = comalloc();
 267         if (cisreal(c)) {
 268                 qfree(r->real);
 269                 r->real = qinv(c->real);
 270                 return r;
 271         }
 272         if (cisimag(c)) {
 273                 q1 = qinv(c->imag);
 274                 qfree(r->imag);
 275                 r->imag = qneg(q1);
 276                 qfree(q1);
 277                 return r;
 278         }
 279         q1 = qsquare(c->real);
 280         q2 = qsquare(c->imag);
 281         den = qqadd(q1, q2);
 282         qfree(q1);
 283         qfree(q2);
 284         qfree(r->real);
 285         r->real = qqdiv(c->real, den);
 286         q1 = qqdiv(c->imag, den);
 287         qfree(r->imag);
 288         r->imag = qneg(q1);
 289         qfree(q1);
 290         qfree(den);
 291         return r;
 292 }
 293
 294
 295 /*
 296  * Negate a complex number.
 297  */
 298 COMPLEX *
 299 c_neg(COMPLEX *c)
 300 {
 301         COMPLEX *r;
 302
 303         if (ciszero(c))
 304                 return clink(&_czero_);
 305         r = comalloc();
 306         if (!qiszero(c->real)) {
 307                 qfree(r->real);
 308                 r->real = qneg(c->real);
 309         }
 310         if (!qiszero(c->imag)) {
 311                 qfree(r->imag);
 312                 r->imag = qneg(c->imag);
 313         }
 314         return r;
 315 }
 316
 317
 318 /*
 319  * Take the integer part of a complex number.
 320  * This means take the integer part of both components.
 321  */
 322 COMPLEX *
 323 c_int(COMPLEX *c)
 324 {
 325         COMPLEX *r;
 326
 327         if (cisint(c))
 328                 return clink(c);
 329         r = comalloc();
 330         qfree(r->real);
 331         r->real = qint(c->real);
 332         qfree(r->imag);
 333         r->imag = qint(c->imag);
 334         return r;
 335 }
 336
 337
 338 /*
 339  * Take the fractional part of a complex number.
 340  * This means take the fractional part of both components.
 341  */
 342 COMPLEX *
 343 c_frac(COMPLEX *c)
 344 {
 345         COMPLEX *r;
 346
 347         if (cisint(c))
 348                 return clink(&_czero_);
 349         r = comalloc();
 350         qfree(r->real);
 351         r->real = qfrac(c->real);
 352         qfree(r->imag);
 353         r->imag = qfrac(c->imag);
 354         return r;
 355 }
 356
 357
 358 /*
 359  * Take the conjugate of a complex number.
 360  * This negates the complex part.
 361  */
 362 COMPLEX *
 363 c_conj(COMPLEX *c)
 364 {
 365         COMPLEX *r;
 366
 367         if (cisreal(c))
 368                 return clink(c);
 369         r = comalloc();
 370         if (!qiszero(c->real)) {
 371                 qfree(r->real);
 372                 r->real = qlink(c->real);
 373         }
 374         qfree(r->imag);
 375         r->imag = qneg(c->imag);
 376         return r;
 377 }
 378
 379
 380 /*
 381  * Return the real part of a complex number.
 382  */
 383 COMPLEX *
 384 c_real(COMPLEX *c)
 385 {
 386         COMPLEX *r;
 387
 388         if (cisreal(c))
 389                 return clink(c);
 390         r = comalloc();
 391         if (!qiszero(c->real)) {
 392                 qfree(r->real);
 393                 r->real = qlink(c->real);
 394         }
 395         return r;
 396 }
 397
 398
 399 /*
 400  * Return the imaginary part of a complex number as a real.
 401  */
 402 COMPLEX *
 403 c_imag(COMPLEX *c)
 404 {
 405         COMPLEX *r;
 406
 407         if (cisreal(c))
 408                 return clink(&_czero_);
 409         r = comalloc();
 410         qfree(r->real);
 411         r->real = qlink(c->imag);
 412         return r;
 413 }
 414
 415
 416 /*
 417  * Add a real number to a complex number.
 418  */
 419 COMPLEX *
 420 c_addq(COMPLEX *c, NUMBER *q)
 421 {
 422         COMPLEX *r;
 423
 424         if (qiszero(q))
 425                 return clink(c);
 426         r = comalloc();
 427         qfree(r->real);
 428         qfree(r->imag);
 429         r->real = qqadd(c->real, q);
 430         r->imag = qlink(c->imag);
 431         return r;
 432 }
 433
 434
 435 /*
 436  * Subtract a real number from a complex number.
 437  */
 438 COMPLEX *
 439 c_subq(COMPLEX *c, NUMBER *q)
 440 {
 441         COMPLEX *r;
 442
 443         if (qiszero(q))
 444                 return clink(c);
 445         r = comalloc();
 446         qfree(r->real);
 447         qfree(r->imag);
 448         r->real = qsub(c->real, q);
 449         r->imag = qlink(c->imag);
 450         return r;
 451 }
 452
 453
 454 /*
 455  * Shift the components of a complex number left by the specified
 456  * number of bits.  Negative values shift to the right.
 457  */
 458 COMPLEX *
 459 c_shift(COMPLEX *c, long n)
 460 {
 461         COMPLEX *r;
 462
 463         if (ciszero(c) || (n == 0))
 464                 return clink(c);
 465         r = comalloc();
 466         qfree(r->real);
 467         qfree(r->imag);
 468         r->real = qshift(c->real, n);
 469         r->imag = qshift(c->imag, n);
 470         return r;
 471 }
 472
 473
 474 /*
 475  * Scale a complex number by a power of two.
 476  */
 477 COMPLEX *
 478 c_scale(COMPLEX *c, long n)
 479 {
 480         COMPLEX *r;
 481
 482         if (ciszero(c) || (n == 0))
 483                 return clink(c);
 484         r = comalloc();
 485         qfree(r->real);
 486         qfree(r->imag);
 487         r->real = qscale(c->real, n);
 488         r->imag = qscale(c->imag, n);
 489         return r;
 490 }
 491
 492
 493 /*
 494  * Multiply a complex number by a real number.
 495  */
 496 COMPLEX *
 497 c_mulq(COMPLEX *c, NUMBER *q)
 498 {
 499         COMPLEX *r;
 500
 501         if (qiszero(q))
 502                 return clink(&_czero_);
 503         if (qisone(q))
 504                 return clink(c);
 505         if (qisnegone(q))
 506                 return c_neg(c);
 507         r = comalloc();
 508         qfree(r->real);
 509         qfree(r->imag);
 510         r->real = qmul(c->real, q);
 511         r->imag = qmul(c->imag, q);
 512         return r;
 513 }
 514
 515
 516 /*
 517  * Divide a complex number by a real number.
 518  */
 519 COMPLEX *
 520 c_divq(COMPLEX *c, NUMBER *q)
 521 {
 522         COMPLEX *r;
 523
 524         if (qiszero(q)) {
 525                 math_error("Division by zero");
 526                 /*NOTREACHED*/
 527         }
 528         if (qisone(q))
 529                 return clink(c);
 530         if (qisnegone(q))
 531                 return c_neg(c);
 532         r = comalloc();
 533         qfree(r->real);
 534         qfree(r->imag);
 535         r->real = qqdiv(c->real, q);
 536         r->imag = qqdiv(c->imag, q);
 537         return r;
 538 }
 539
 540
 541
 542
 543 /*
 544  * Construct a complex number given the real and imaginary components.
 545  */
 546 COMPLEX *
 547 qqtoc(NUMBER *q1, NUMBER *q2)
 548 {
 549         COMPLEX *r;
 550
 551         if (qiszero(q1) && qiszero(q2))
 552                 return clink(&_czero_);
 553         r = comalloc();
 554         qfree(r->real);
 555         qfree(r->imag);
 556         r->real = qlink(q1);
 557         r->imag = qlink(q2);
 558         return r;
 559 }
 560
 561
 562 /*
 563  * Compare two complex numbers for equality, returning FALSE if they are equal,
 564  * and TRUE if they differ.
 565  */
 566 BOOL
 567 c_cmp(COMPLEX *c1, COMPLEX *c2)
 568 {
 569         BOOL i;
 570
 571         i = qcmp(c1->real, c2->real);
 572         if (!i)
 573                 i = qcmp(c1->imag, c2->imag);
 574         return i;
 575 }
 576
 577
 578 /*
 579  * Compare two complex numbers and return a complex number with real and
 580  * imaginary parts -1, 0 or 1 indicating relative values of the real and
 581  * imaginary parts of the two numbers.
 582  */
 583 COMPLEX *
 584 c_rel(COMPLEX *c1, COMPLEX *c2)
 585 {
 586         COMPLEX *c;
 587
 588         c = comalloc();
 589         qfree(c->real);
 590         qfree(c->imag);
 591         c->real = itoq((long) qrel(c1->real, c2->real));
 592         c->imag = itoq((long) qrel(c1->imag, c2->imag));
 593
 594         return c;
 595 }
 596
 597
 598 /*
 599  * Allocate a new complex number.
 600  */
 601 COMPLEX *
 602 comalloc(void)
 603 {
 604         COMPLEX *r;
 605
 606         r = (COMPLEX *) malloc(sizeof(COMPLEX));
 607         if (r == NULL) {
 608                 math_error("Cannot allocate complex number");
 609                 /*NOTREACHED*/
 610         }
 611         r->links = 1;
 612         r->real = qlink(&_qzero_);
 613         r->imag = qlink(&_qzero_);
 614         return r;
 615 }
 616
 617
 618 /*
 619  * Free a complex number.
 620  */
 621 void
 622 comfree(COMPLEX *c)
 623 {
 624         if (--(c->links) > 0)
 625                 return;
 626         qfree(c->real);
 627         qfree(c->imag);
 628         free(c);
 629 }