home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume26 / calc / part05 / commath.c < prev    next >
C/C++ Source or Header  |  1992-05-09  |  10KB  |  611 lines
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision complex arithmetic primitive routines
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include "math.h"
  11.  
  12.  
  13. COMPLEX _czero_ =        { &_qzero_, &_qzero_, 1 };
  14. COMPLEX _cone_ =        { &_qone_, &_qzero_, 1 };
  15. static COMPLEX _cnegone_ =    { &_qnegone_, &_qzero_, 1 };
  16.  
  17. #if 0
  18. COMPLEX _conei_ =    { &_qzero_, &_qone_, 1 };
  19. #endif
  20.  
  21.  
  22. /*
  23.  * Free list for complex numbers.
  24.  */
  25. static FREELIST freelist = {
  26.     sizeof(COMPLEX),    /* size of an item */
  27.     100            /* number of free items to keep */
  28. };
  29.  
  30.  
  31. /*
  32.  * Add two complex numbers.
  33.  */
  34. COMPLEX *
  35. cadd(c1, c2)
  36.     COMPLEX *c1, *c2;
  37. {
  38.     COMPLEX *r;
  39.  
  40.     if (ciszero(c1))
  41.         return clink(c2);
  42.     if (ciszero(c2))
  43.         return clink(c1);
  44.     r = comalloc();
  45.     if (!qiszero(c1->real) || !qiszero(c2->real))
  46.         r->real = qadd(c1->real, c2->real);
  47.     if (!qiszero(c1->imag) || !qiszero(c2->imag))
  48.         r->imag = qadd(c1->imag, c2->imag);
  49.     return r;
  50. }
  51.  
  52.  
  53. /*
  54.  * Subtract two complex numbers.
  55.  */
  56. COMPLEX *
  57. csub(c1, c2)
  58.     COMPLEX *c1, *c2;
  59. {
  60.     COMPLEX *r;
  61.  
  62.     if ((c1->real == c2->real) && (c1->imag == c2->imag))
  63.         return clink(&_czero_);
  64.     if (ciszero(c2))
  65.         return clink(c1);
  66.     r = comalloc();
  67.     if (!qiszero(c1->real) || !qiszero(c2->real))
  68.         r->real = qsub(c1->real, c2->real);
  69.     if (!qiszero(c1->imag) || !qiszero(c2->imag))
  70.         r->imag = qsub(c1->imag, c2->imag);
  71.     return r;
  72. }
  73.  
  74.  
  75. /*
  76.  * Multiply two complex numbers.
  77.  * This saves one multiplication over the obvious algorithm by
  78.  * trading it for several extra additions, as follows.  Let
  79.  *    q1 = (a + b) * (c + d)
  80.  *    q2 = a * c
  81.  *    q3 = b * d
  82.  * Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
  83.  */
  84. COMPLEX *
  85. cmul(c1, c2)
  86.     COMPLEX *c1, *c2;
  87. {
  88.     COMPLEX *r;
  89.     NUMBER *q1, *q2, *q3, *q4;
  90.  
  91.     if (ciszero(c1) || ciszero(c2))
  92.         return clink(&_czero_);
  93.     if (cisone(c1))
  94.         return clink(c2);
  95.     if (cisone(c2))
  96.         return clink(c1);
  97.     if (cisreal(c2))
  98.         return cmulq(c1, c2->real);
  99.     if (cisreal(c1))
  100.         return cmulq(c2, c1->real);
  101.     /*
  102.      * Need to do the full calculation.
  103.      */
  104.     r = comalloc();
  105.     q2 = qadd(c1->real, c1->imag);
  106.     q3 = qadd(c2->real, c2->imag);
  107.     q1 = qmul(q2, q3);
  108.     qfree(q2);
  109.     qfree(q3);
  110.     q2 = qmul(c1->real, c2->real);
  111.     q3 = qmul(c1->imag, c2->imag);
  112.     q4 = qadd(q2, q3);
  113.     r->real = qsub(q2, q3);
  114.     r->imag = qsub(q1, q4);
  115.     qfree(q1);
  116.     qfree(q2);
  117.     qfree(q3);
  118.     qfree(q4);
  119.     return r;
  120. }
  121.  
  122.  
  123. /*
  124.  * Square a complex number.
  125.  */
  126. COMPLEX *
  127. csquare(c)
  128.     COMPLEX *c;
  129. {
  130.     COMPLEX *r;
  131.     NUMBER *q1, *q2;
  132.  
  133.     if (ciszero(c))
  134.         return clink(&_czero_);
  135.     if (cisrunit(c))
  136.         return clink(&_cone_);
  137.     if (cisiunit(c))
  138.         return clink(&_cnegone_);
  139.     r = comalloc();
  140.     if (cisreal(c)) {
  141.         r->real = qsquare(c->real);
  142.         return r;
  143.     }
  144.     if (cisimag(c)) {
  145.         q1 = qsquare(c->imag);
  146.         r->real = qneg(q1);
  147.         qfree(q1);
  148.         return r;
  149.     }
  150.     q1 = qsquare(c->real);
  151.     q2 = qsquare(c->imag);
  152.     r->real = qsub(q1, q2);
  153.     qfree(q1);
  154.     qfree(q2);
  155.     q1 = qmul(c->real, c->imag);
  156.     r->imag = qscale(q1, 1L);
  157.     qfree(q1);
  158.     return r;
  159. }
  160.  
  161.  
  162. /*
  163.  * Divide two complex numbers.
  164.  */
  165. COMPLEX *
  166. cdiv(c1, c2)
  167.     COMPLEX *c1, *c2;
  168. {
  169.     COMPLEX *r;
  170.     NUMBER *q1, *q2, *q3, *den;
  171.  
  172.     if (ciszero(c2))
  173.         error("Division by zero");
  174.     if ((c1->real == c2->real) && (c1->imag == c2->imag))
  175.         return clink(&_cone_);
  176.     r = comalloc();
  177.     if (cisreal(c1) && cisreal(c2)) {
  178.         r->real = qdiv(c1->real, c2->real);
  179.         return r;
  180.     }
  181.     if (cisimag(c1) && cisimag(c2)) {
  182.         r->real = qdiv(c1->imag, c2->imag);
  183.         return r;
  184.     }
  185.     if (cisimag(c1) && cisreal(c2)) {
  186.         r->imag = qdiv(c1->imag, c2->real);
  187.         return r;
  188.     }
  189.     if (cisreal(c1) && cisimag(c2)) {
  190.         q1 = qdiv(c1->real, c2->imag);
  191.         r->imag = qneg(q1);
  192.         qfree(q1);
  193.         return r;
  194.     }
  195.     if (cisreal(c2)) {
  196.         r->real = qdiv(c1->real, c2->real);
  197.         r->imag = qdiv(c1->imag, c2->real);
  198.         return r;
  199.     }
  200.     q1 = qsquare(c2->real);
  201.     q2 = qsquare(c2->imag);
  202.     den = qadd(q1, q2);
  203.     qfree(q1);
  204.     qfree(q2);
  205.     q1 = qmul(c1->real, c2->real);
  206.     q2 = qmul(c1->imag, c2->imag);
  207.     q3 = qadd(q1, q2);
  208.     qfree(q1);
  209.     qfree(q2);
  210.     r->real = qdiv(q3, den);
  211.     qfree(q3);
  212.     q1 = qmul(c1->real, c2->imag);
  213.     q2 = qmul(c1->imag, c2->real);
  214.     q3 = qsub(q2, q1);
  215.     qfree(q1);
  216.     qfree(q2);
  217.     r->imag = qdiv(q3, den);
  218.     qfree(q3);
  219.     qfree(den);
  220.     return r;
  221. }
  222.  
  223.  
  224. /*
  225.  * Invert a complex number.
  226.  */
  227. COMPLEX *
  228. cinv(c)
  229.     COMPLEX *c;
  230. {
  231.     COMPLEX *r;
  232.     NUMBER *q1, *q2, *den;
  233.  
  234.     if (ciszero(c))
  235.         error("Inverting zero");
  236.     r = comalloc();
  237.     if (cisreal(c)) {
  238.         r->real = qinv(c->real);
  239.         return r;
  240.     }
  241.     if (cisimag(c)) {
  242.         q1 = qinv(c->imag);
  243.         r->imag = qneg(q1);
  244.         qfree(q1);
  245.         return r;
  246.     }
  247.     q1 = qsquare(c->real);
  248.     q2 = qsquare(c->imag);
  249.     den = qadd(q1, q2);
  250.     qfree(q1);
  251.     qfree(q2);
  252.     r->real = qdiv(c->real, den);
  253.     q1 = qdiv(c->imag, den);
  254.     r->imag = qneg(q1);
  255.     qfree(q1);
  256.     qfree(den);
  257.     return r;
  258. }
  259.  
  260.  
  261. /*
  262.  * Negate a complex number.
  263.  */
  264. COMPLEX *
  265. cneg(c)
  266.     COMPLEX *c;
  267. {
  268.     COMPLEX *r;
  269.  
  270.     if (ciszero(c))
  271.         return clink(&_czero_);
  272.     r = comalloc();
  273.     if (!qiszero(c->real))
  274.         r->real = qneg(c->real);
  275.     if (!qiszero(c->imag))
  276.         r->imag = qneg(c->imag);
  277.     return r;
  278. }
  279.  
  280.  
  281. /*
  282.  * Take the integer part of a complex number.
  283.  * This means take the integer part of both components.
  284.  */
  285. COMPLEX *
  286. cint(c)
  287.     COMPLEX *c;
  288. {
  289.     COMPLEX *r;
  290.  
  291.     if (cisint(c))
  292.         return clink(c);
  293.     r = comalloc();
  294.     r->real = qint(c->real);
  295.     r->imag = qint(c->imag);
  296.     return r;
  297. }
  298.  
  299.  
  300. /*
  301.  * Take the fractional part of a complex number.
  302.  * This means take the fractional part of both components.
  303.  */
  304. COMPLEX *
  305. cfrac(c)
  306.     COMPLEX *c;
  307. {
  308.     COMPLEX *r;
  309.  
  310.     if (cisint(c))
  311.         return clink(&_czero_);
  312.     r = comalloc();
  313.     r->real = qfrac(c->real);
  314.     r->imag = qfrac(c->imag);
  315.     return r;
  316. }
  317.  
  318.  
  319. #if 0
  320. /*
  321.  * Take the conjugate of a complex number.
  322.  * This negates the complex part.
  323.  */
  324. COMPLEX *
  325. cconj(c)
  326.     COMPLEX *c;
  327. {
  328.     COMPLEX *r;
  329.  
  330.     if (cisreal(c))
  331.         return clink(c);
  332.     r = comalloc();
  333.     if (!qiszero(c->real))
  334.         r->real = qlink(c->real);
  335.     r->imag = qneg(c->imag);
  336.     return r;
  337. }
  338.  
  339.  
  340. /*
  341.  * Return the real part of a complex number.
  342.  */
  343. COMPLEX *
  344. creal(c)
  345.     COMPLEX *c;
  346. {
  347.     COMPLEX *r;
  348.  
  349.     if (cisreal(c))
  350.         return clink(c);
  351.     r = comalloc();
  352.     if (!qiszero(c->real))
  353.         r->real = qlink(c->real);
  354.     return r;
  355. }
  356.  
  357.  
  358. /*
  359.  * Return the imaginary part of a complex number as a real.
  360.  */
  361. COMPLEX *
  362. cimag(c)
  363.     COMPLEX *c;
  364. {
  365.     COMPLEX *r;
  366.  
  367.     if (cisreal(c))
  368.         return clink(&_czero_);
  369.     r = comalloc();
  370.     r->real = qlink(c->imag);
  371.     return r;
  372. }
  373. #endif
  374.  
  375.  
  376. /*
  377.  * Add a real number to a complex number.
  378.  */
  379. COMPLEX *
  380. caddq(c, q)
  381.     COMPLEX *c;
  382.     NUMBER *q;
  383. {
  384.     COMPLEX *r;
  385.  
  386.     if (qiszero(q))
  387.         return clink(c);
  388.     r = comalloc();
  389.     r->real = qadd(c->real, q);
  390.     r->imag = qlink(c->imag);
  391.     return r;
  392. }
  393.  
  394.  
  395. /*
  396.  * Subtract a real number from a complex number.
  397.  */
  398. COMPLEX *
  399. csubq(c, q)
  400.     COMPLEX *c;
  401.     NUMBER *q;
  402. {
  403.     COMPLEX *r;
  404.  
  405.     if (qiszero(q))
  406.         return clink(c);
  407.     r = comalloc();
  408.     r->real = qsub(c->real, q);
  409.     r->imag = qlink(c->imag);
  410.     return r;
  411. }
  412.  
  413.  
  414. /*
  415.  * Shift the components of a complex number left by the specified
  416.  * number of bits.  Negative values shift to the right.
  417.  */
  418. COMPLEX *
  419. cshift(c, n)
  420.     COMPLEX *c;
  421.     long n;
  422. {
  423.     COMPLEX *r;
  424.  
  425.     if (ciszero(c) || (n == 0))
  426.         return clink(c);
  427.     r = comalloc();
  428.     r->real = qshift(c->real, n);
  429.     r->imag = qshift(c->imag, n);
  430.     return r;
  431. }
  432.  
  433.  
  434. /*
  435.  * Scale a complex number by a power of two.
  436.  */
  437. COMPLEX *
  438. cscale(c, n)
  439.     COMPLEX *c;
  440.     long n;
  441. {
  442.     COMPLEX *r;
  443.  
  444.     if (ciszero(c) || (n == 0))
  445.         return clink(c);
  446.     r = comalloc();
  447.     r->real = qscale(c->real, n);
  448.     r->imag = qscale(c->imag, n);
  449.     return r;
  450. }
  451.  
  452.  
  453. /*
  454.  * Multiply a complex number by a real number.
  455.  */
  456. COMPLEX *
  457. cmulq(c, q)
  458.     COMPLEX *c;
  459.     NUMBER *q;
  460. {
  461.     COMPLEX *r;
  462.  
  463.     if (qiszero(q))
  464.         return clink(&_czero_);
  465.     if (qisone(q))
  466.         return clink(c);
  467.     if (qisnegone(q))
  468.         return cneg(c);
  469.     r = comalloc();
  470.     r->real = qmul(c->real, q);
  471.     r->imag = qmul(c->imag, q);
  472.     return r;
  473. }
  474.  
  475.  
  476. /*
  477.  * Divide a complex number by a real number.
  478.  */
  479. COMPLEX *
  480. cdivq(c, q)
  481.     COMPLEX *c;
  482.     NUMBER *q;
  483. {
  484.     COMPLEX *r;
  485.  
  486.     if (qiszero(q))
  487.         error("Division by zero");
  488.     if (qisone(q))
  489.         return clink(c);
  490.     if (qisnegone(q))
  491.         return cneg(c);
  492.     r = comalloc();
  493.     r->real = qdiv(c->real, q);
  494.     r->imag = qdiv(c->imag, q);
  495.     return r;
  496. }
  497.  
  498.  
  499. /*
  500.  * Take the integer quotient of a complex number by a real number.
  501.  * This is defined to be the result of doing the quotient for each component.
  502.  */
  503. COMPLEX *
  504. cquoq(c, q)
  505.     COMPLEX *c;
  506.     NUMBER *q;
  507. {
  508.     COMPLEX *r;
  509.  
  510.     if (qiszero(q))
  511.         error("Division by zero");
  512.     r = comalloc();
  513.     r->real = qquo(c->real, q);
  514.     r->imag = qquo(c->imag, q);
  515.     return r;
  516. }
  517.  
  518.  
  519. /*
  520.  * Take the modulus of a complex number by a real number.
  521.  * This is defined to be the result of doing the modulo for each component.
  522.  */
  523. COMPLEX *
  524. cmodq(c, q)
  525.     COMPLEX *c;
  526.     NUMBER *q;
  527. {
  528.     COMPLEX *r;
  529.  
  530.     if (qiszero(q))
  531.         error("Division by zero");
  532.     r = comalloc();
  533.     r->real = qmod(c->real, q);
  534.     r->imag = qmod(c->imag, q);
  535.     return r;
  536. }
  537.  
  538.  
  539. #if 0
  540. /*
  541.  * Construct a complex number given the real and imaginary components.
  542.  */
  543. COMPLEX *
  544. qqtoc(q1, q2)
  545.     NUMBER *q1, *q2;
  546. {
  547.     COMPLEX *r;
  548.  
  549.     if (qiszero(q1) && qiszero(q2))
  550.         return clink(&_czero_);
  551.     r = comalloc();
  552.     if (!qiszero(q1))
  553.         r->real = qlink(q1);
  554.     if (!qiszero(q2))
  555.         r->imag = qlink(q2);
  556.     return r;
  557. }
  558. #endif
  559.  
  560.  
  561. /*
  562.  * Compare two complex numbers for equality, returning FALSE if they are equal,
  563.  * and TRUE if they differ.
  564.  */
  565. BOOL
  566. ccmp(c1, c2)
  567.     COMPLEX *c1, *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.  * Allocate a new complex number.
  580.  */
  581. COMPLEX *
  582. comalloc()
  583. {
  584.     COMPLEX *r;
  585.  
  586.     r = (COMPLEX *) allocitem(&freelist);
  587.     if (r == NULL)
  588.         error("Cannot allocate complex number");
  589.     r->links = 1;
  590.     r->real = qlink(&_qzero_);
  591.     r->imag = qlink(&_qzero_);
  592.     return r;
  593. }
  594.  
  595.  
  596. /*
  597.  * Free a complex number.
  598.  */
  599. void
  600. comfree(c)
  601.     COMPLEX *c;
  602. {
  603.     if (--(c->links) > 0)
  604.         return;
  605.     qfree(c->real);
  606.     qfree(c->imag);
  607.     freeitem(&freelist, (FREEITEM *) c);
  608. }
  609.  
  610. /* END CODE */
  611.