"Fossies" - the Fresh Open Source Software Archive

Member "grpn-1.5.2/src/complex.c" (24 Nov 2018, 14615 Bytes) of package /linux/misc/grpn-1.5.2.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "complex.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.5.1_vs_1.5.2.

    1 /*
    2 
    3 Copyright (C) 2002  Paul Wilkins
    4 
    5 This program is free software; you can redistribute it and/or
    6 modify it under the terms of the GNU General Public License
    7 as published by the Free Software Foundation; either version 2
    8 of the License, or (at your option) any later version.
    9 
   10 This program is distributed in the hope that it will be useful,
   11 but WITHOUT ANY WARRANTY; without even the implied warranty of
   12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13 GNU General Public License for more details.
   14 
   15 You should have received a copy of the GNU General Public License
   16 along with this program; if not, write to the Free Software
   17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
   18 
   19 */
   20 /* complex.c  by Paul Wilkins */
   21 
   22 #include <stdio.h>
   23 #include <stdlib.h>
   24 #include <math.h>
   25 #include <string.h>
   26 
   27 #include "complex.h"
   28 #include "real.h"
   29 #include "number.h"
   30 #include "constant.h"
   31 #include "mode.h"
   32 
   33 #include "lcd.h"  /* we need to know the display size */
   34 
   35 
   36 
   37 /* create a new complex number */
   38 Cmplx * newCmplx(){
   39    Cmplx *p;
   40    if(NULL == (p = (Cmplx *)malloc(sizeof(Cmplx)))){
   41       perror("Malloc");
   42       exit(0);
   43    }
   44    p->re = NULL;
   45    p->im = NULL;
   46 
   47    return p;
   48 }
   49 
   50 void freeCmplx(Cmplx *a){
   51    if(a){
   52       if(a->re) freeReal(a->re);
   53       if(a->im) freeReal(a->im);
   54       free((char *)a);
   55       a = NULL; /* make sure we die... */
   56    }
   57 }
   58 
   59 /* assumes data in in rectangular form */
   60 Cmplx * setCmplxReal(Cmplx *a, Real *rp, Real *ip){
   61    if(a == NULL || rp == NULL || ip == NULL)
   62       { fprintf(stderr, "setCmplxReal(NULL)\n"); exit(0); }
   63 
   64    if(a->re) free((char *)a->re);
   65    if(a->im) free((char *)a->im);
   66 
   67    a->re = setRealReal(newReal(), rp);
   68    a->im = setRealReal(newReal(), ip);
   69 
   70    return a;
   71 }
   72 
   73 /* assumes data in in form defined by getRadixMode() and getPolarMode() */
   74 Cmplx * inputCmplxReal(Cmplx *a, Real *rp, Real *ip){
   75    Cmplx *c1, *c2;
   76    if(a == NULL || rp == NULL || ip == NULL)
   77       { fprintf(stderr, "setCmplxReal(NULL)\n"); exit(0); }
   78 
   79    if(a->re) free((char *)a->re);
   80    if(a->im) free((char *)a->im);
   81 
   82    if(getPolarMode() == POLAR){
   83       c1 = newCmplx();
   84       c1->re = setRealReal(newReal(), rp);
   85 
   86       if(getRadixMode() == DEGREES) c1->im = divReal(ip, real180Pi);
   87       else c1->im = setRealReal(newReal(), ip);
   88 
   89       c2 = rectCmplx(c1);
   90       setCmplxCmplx(a, c2);
   91       freeCmplx(c1);
   92       freeCmplx(c2);
   93    } else {
   94       a->re = setRealReal(newReal(), rp);
   95       a->im = setRealReal(newReal(), ip);
   96    }
   97 
   98    return a;
   99 }
  100 
  101 Cmplx * setCmplxCmplx(Cmplx *a, Cmplx *b){
  102    if(a == NULL || b == NULL)
  103       { fprintf(stderr, "setCmplxCmplx(NULL)\n"); exit(0); }
  104 
  105    if(a->re) free((char *)a->re);
  106    if(a->im) free((char *)a->im);
  107 
  108    a->re = setRealReal(newReal(), b->re);
  109    a->im = setRealReal(newReal(), b->im);
  110 
  111    return a;
  112 }
  113 
  114 #define COMPLEX_PRINT_SIZE 90
  115 
  116 char * printCmplxShort(Cmplx *a){
  117    char *c;
  118    char *p1;
  119    char *p2;
  120    Cmplx *c1;
  121 
  122    if(NULL == (c=(char *)malloc(COMPLEX_PRINT_SIZE)))
  123       { perror("Malloc"); exit(0); }
  124 
  125    if(getPolarMode() == POLAR){
  126       c1 = polarCmplx(a);
  127       p1 = printReal(c1->re);
  128       if(getRadixMode() == DEGREES){
  129          mulEqReal(c1->im, real180Pi);
  130          p2 = printReal(c1->im);
  131       } else {
  132          p2 = printReal(c1->im);
  133       }
  134       freeCmplx(c1);
  135       sprintf(c, "(%s< %s)", p1, p2);
  136    } else {
  137       p1 = printReal(a->re);
  138       p2 = printReal(a->im);
  139       sprintf(c, "(%s; %s)", p1, p2);
  140    }
  141 
  142    free(p1);
  143    free(p2);
  144 
  145    return c;
  146 }
  147 
  148 char * printCmplx(Cmplx *a){
  149    char *c;
  150    char *p1;
  151    char *p2;
  152    Cmplx *c1;
  153 
  154    if(NULL == (c=(char *)malloc(COMPLEX_PRINT_SIZE)))
  155       { perror("Malloc"); exit(0); }
  156    *c = '\0';
  157 
  158    if(getPolarMode() == POLAR){
  159       c1 = polarCmplx(a);
  160       p1 = printReal(c1->re);
  161       if(getRadixMode() == DEGREES){
  162          mulEqReal(c1->im, real180Pi);
  163          p2 = printReal(c1->im);
  164       } else {
  165          p2 = printReal(c1->im);
  166       }
  167       freeCmplx(c1);
  168       sprintf(c, "(%s< %s)", p1, p2);
  169    } else {
  170       p1 = printReal(a->re);
  171       p2 = printReal(a->im);
  172       sprintf(c, "(%s; %s)", p1, p2);
  173    }
  174 
  175 
  176    if(strlen(p1)+strlen(p2)+4 > lcdWidth-4) *(c+strlen(p1)+2) = '\n';
  177    free(p1);
  178    free(p2);
  179 
  180    return c;
  181 }
  182 
  183 Cmplx * negCmplx(Cmplx *a){
  184    Cmplx *p = newCmplx();
  185 
  186    p->re = negReal(a->re);
  187    p->im = negReal(a->im);
  188 
  189    return p;
  190 }
  191 
  192 Cmplx * negEqCmplx(Cmplx *a){
  193    negEqReal(a->re);
  194    negEqReal(a->im);
  195 
  196    return a;
  197 }
  198 
  199 /* calculate 1/(re+im) */
  200 Cmplx * invCmplx(Cmplx *a){
  201    Real *r1, *r2, *r3;
  202    Cmplx *p = newCmplx();
  203 
  204    /* r1 = re^2 + im^2 */
  205    r1 = mulReal(a->re, a->re);
  206    r2 = mulReal(a->im, a->im);
  207    addEqReal(r1, r2);
  208 
  209    p->re = divReal(a->re, r1);
  210 
  211    r3 = divReal(a->im, r1);
  212    p->im = negEqReal(r3);
  213 
  214    freeReal(r1);
  215    freeReal(r2);
  216 
  217    return p;
  218 }
  219 
  220 /* calculate 1/(re+im) */
  221 Cmplx * invEqCmplx(Cmplx *a){
  222    Real *r1, *r2;
  223 
  224    /* r1 = re^2 + im^2 */
  225    r1 = mulReal(a->re, a->re);
  226    r2 = mulReal(a->im, a->im);
  227    addEqReal(r1, r2);
  228 
  229    divEqReal(a->re, r1);
  230 
  231    divEqReal(a->im, r1);
  232    negEqReal(a->im);
  233 
  234    freeReal(r1);
  235    freeReal(r2);
  236 
  237    return a;
  238 }
  239 
  240 Cmplx * powCmplxInt(Cmplx *a, int b){
  241    int i;
  242    Cmplx *c1, *c2;
  243 
  244    c1 = setCmplxCmplx(newCmplx(), a);;
  245 
  246    for(i=1; i<b; i++){
  247       c2 = mulCmplx(c1, a);
  248       freeCmplx(c1);
  249       c1 = c2;
  250    }
  251 
  252    return c1;
  253 }
  254 
  255 Cmplx * powCmplx(Cmplx *a, Cmplx *b){
  256    Cmplx *c1;
  257    Cmplx *p;
  258 
  259    c1 = lnCmplx(a);
  260 
  261    p = expEqCmplx(mulCmplx(c1, b));
  262 
  263    freeCmplx(c1);
  264 
  265    return p;
  266 }
  267 
  268 Cmplx * powCmplxReal(Cmplx *a, Real *b){
  269    Cmplx *c1, *c2;
  270    Cmplx *p;
  271 
  272    c1 = lnCmplx(a);
  273    c2 = mulCmplxReal(c1, b);
  274    freeCmplx(c1);
  275 
  276    p = expCmplx(c2);
  277    freeCmplx(c2);
  278 
  279    return p;
  280 }
  281 
  282 Cmplx * powRealCmplx(Real *a, Cmplx *b){
  283    Real *r1;
  284    Cmplx *c1, *c2, *c3;
  285    Cmplx *p;
  286 
  287    if(-1 == cmpReal(a, realZero)){
  288       /* we can't take lnReal(<0) */
  289       c2 = setCmplxReal(newCmplx(), a, realZero);
  290       c3 = lnCmplx(c2);
  291       c1 = mulCmplx(b, c3);
  292       freeCmplx(c2);
  293       freeCmplx(c3);
  294    } else {
  295       r1 = lnReal(a);
  296       c1 = mulCmplxReal(b, r1);
  297       freeReal(r1);
  298    }
  299 
  300    p = expCmplx(c1);
  301    freeCmplx(c1);
  302 
  303    return p;
  304 }
  305 
  306 Cmplx * polarCmplx(Cmplx *a){
  307    Cmplx *p = newCmplx();
  308 
  309    p->re = absCmplx(a);
  310    p->im = thetaCmplx(a);
  311 
  312    return p;
  313 }
  314 
  315 Real * absCmplx(Cmplx *a){
  316    Real *re, *ri;
  317    Real *r1, *r2;
  318 
  319    /* Implements re = re^2 + im^2 but without the range problem */
  320 
  321    re = absReal(a->re);
  322    ri = absReal(a->im);
  323 
  324    if(1 == cmpReal(re, ri)){
  325       r2 = re;  /* bigger */
  326       r1 = ri;  /* smaller */
  327    } else {
  328       r2 = ri;
  329       r1 = re;
  330    }
  331 
  332    divEqReal(r1, r2);
  333    mulEqReal(r1, r1);
  334    addEqReal(r1, realOne);
  335    powEqReal(r1, realHalf);
  336    mulEqReal(r1, r2);
  337 
  338    freeReal(r2);
  339 
  340    return r1;
  341 
  342 }
  343 
  344 Real * thetaCmplx(Cmplx *a){
  345    int sign_re;
  346    Real *r1;
  347    Real *theta;
  348 
  349    r1 = atanEqReal(divReal(a->im, a->re));
  350 
  351    sign_re = cmpReal(a->re, realZero);
  352 
  353    if(0 == sign_re){
  354 
  355       /* -90 deg */
  356       if(-1 == cmpReal(a->im, realZero))
  357          theta = setRealDouble(newReal(), -M_PI/2.0);
  358       else /* 90 deg */
  359          theta = setRealDouble(newReal(), M_PI/2.0);
  360 
  361       freeReal(r1);
  362    }
  363    /* quadrant 2 and 3 */
  364    else if(-1 == sign_re){
  365 
  366       /* quad 3 */
  367       if(-1 == cmpReal(a->im, realZero))
  368          theta = subEqReal(r1, realPi);
  369       else /* quad 2 */
  370          theta = addEqReal(r1, realPi);
  371 
  372    }
  373    /* quadrant 1 and 4 */
  374    else{
  375       theta = r1;
  376    }
  377 
  378    return theta;
  379 }
  380 
  381 Cmplx * rectCmplx(Cmplx *a){
  382    Real *r1, *r2;
  383    Cmplx *p = newCmplx();
  384 
  385    r1 = cosReal(a->im);
  386    r2 = sinReal(a->im);
  387    p->re = mulEqReal(r1, a->re);
  388    p->im = mulEqReal(r2, a->re);
  389 
  390    return p;
  391 }
  392 
  393 Cmplx * lnCmplx(Cmplx *a){
  394    Cmplx *p;
  395 
  396    p = polarCmplx(a);
  397 
  398    lnEqReal(p->re);
  399 
  400    return p;
  401 
  402 }
  403 
  404 /* Note: instead of doing lnCmplx(a)/lnReal(10.0) I did it 
  405  *   like this.  Hopefully this way is a little faster???
  406  */
  407 Cmplx * logCmplx(Cmplx *a){
  408    Real *r1;
  409    Cmplx *p;
  410 
  411    p = polarCmplx(a);
  412 
  413    r1 = lnReal(realTen);
  414 
  415    logEqReal(p->re);
  416    divEqReal(p->im, r1);
  417 
  418    freeReal(r1);
  419 
  420    return p;
  421 
  422 }
  423 
  424 Cmplx * expCmplx(Cmplx *a){
  425    Real *rr, *ri, *re;
  426    Cmplx *p = newCmplx();
  427 
  428    rr = cosReal(a->im);
  429    ri = sinReal(a->im);
  430 
  431    re = expReal(a->re);
  432 
  433    p->re = mulEqReal(rr, re);
  434    p->im = mulEqReal(ri, re);
  435 
  436    freeReal(re);
  437 
  438    return p;
  439 }
  440 
  441 Cmplx * expEqCmplx(Cmplx *a){
  442    Real *r1, *r2;
  443 
  444    r1 = expReal(a->re);
  445 
  446    /* res.re = e^re * cos(im) */
  447    mulEqReal(expEqReal(a->re), (r2=cosReal(a->im)));
  448 
  449    /* res.im = e^re * sin(im) */
  450    mulEqReal(sinEqReal(a->im), r1);
  451 
  452    freeReal(r1);
  453    freeReal(r2);
  454 
  455    return a;
  456 }
  457 
  458 
  459 /***************** TRIG *************************/
  460 
  461 Cmplx * sinCmplx(Cmplx *a){
  462    Real *r1;
  463    Cmplx *p;
  464    Cmplx *c1, *c2, *c4, *c5;
  465 
  466    c1 = mulCmplx(cmplxI, a);
  467    c2 = expCmplx(c1);
  468    negEqCmplx(c1);
  469    c4 = expCmplx(c1);
  470 
  471    c5 = subCmplx(c2, c4);
  472 
  473    freeCmplx(c1);
  474    freeCmplx(c2);
  475    freeCmplx(c4);
  476 
  477    r1 = setRealDouble(newReal(), 2.0);
  478    c1 = setCmplxReal(newCmplx(), realZero, r1);
  479    freeReal(r1);
  480 
  481    p = divCmplx(c5, c1);
  482 
  483    freeCmplx(c1);
  484    freeCmplx(c5);
  485    
  486    return p;
  487 }
  488 
  489 Cmplx * cosCmplx(Cmplx *a){
  490    Cmplx *p;
  491    Cmplx *c1, *c2, *c4;
  492 
  493    c1 = mulCmplx(cmplxI, a);
  494    c2 = expCmplx(c1);
  495    negEqCmplx(c1);
  496    c4 = expCmplx(c1);
  497 
  498    p = addCmplx(c2, c4);
  499 
  500    freeCmplx(c1);
  501    freeCmplx(c2);
  502    freeCmplx(c4);
  503 
  504    mulEqReal(p->re, realHalf);
  505    mulEqReal(p->im, realHalf);
  506 
  507    return p;
  508 }
  509 
  510 Cmplx * tanCmplx(Cmplx *a){
  511    Cmplx *p;
  512    Cmplx *c1, *c2;
  513 
  514    c1 = sinCmplx(a);
  515    c2 = cosCmplx(a);
  516 
  517    p = divCmplx(c1, c2);
  518 
  519    freeCmplx(c1);
  520    freeCmplx(c2);
  521 
  522    return p;
  523 }
  524 
  525 Cmplx * asinCmplx(Cmplx *a){
  526    Real *r1;
  527    Cmplx *z, *sqzp1, *sqzm1;
  528    Cmplx *c1, *c2, *c3, *c4, *c5, *c6, *c7;
  529 
  530    z = newCmplx();
  531    setCmplxCmplx(z, a);
  532    if(1 == cmpReal(realZero, a->re)) negEqCmplx(z);
  533 
  534    /* sqrt(z + 1) */
  535    c1 = addCmplxReal(z, realOne);
  536    sqzp1 = powCmplxReal(c1, realHalf);
  537    freeCmplx(c1);
  538 
  539    /* sqrt(z - 1) */
  540    c1 = subCmplxReal(z, realOne);
  541    sqzm1 = powCmplxReal(c1, realHalf);
  542    freeCmplx(c1);
  543 
  544    /* if imag_part(sqzp1) < 0.0 then sqzp1 = -sqzp1 */
  545    if(1 == cmpReal(realZero, sqzp1->im)) negEqCmplx(sqzp1);
  546 
  547    c2 = mulCmplx(sqzp1, sqzm1);
  548    c3 = addCmplx(c2, z);
  549    c4 = lnCmplx(c3);
  550    c5 = mulCmplx(c4, cmplxI);
  551 
  552    c6 = subRealCmplx(realPi2, c5);
  553 
  554    if(1 == cmpReal(c6->re, realPi2)){
  555       c7 = subRealCmplx(realPi, c6);
  556       freeCmplx(c6);
  557       c6 = c7;
  558    }
  559    r1 = negReal(realPi);
  560    if(-1 == cmpReal(c6->re, r1)){
  561       c7 = subRealCmplx(r1, c6);
  562       freeCmplx(c6);
  563       c6 = c7;
  564    }
  565    freeReal(r1);
  566    if(1 == cmpReal(realZero, a->re)) negEqCmplx(c6);
  567       
  568 
  569    freeCmplx(sqzp1);
  570    freeCmplx(sqzm1);
  571    freeCmplx(z);
  572    freeCmplx(c2);
  573    freeCmplx(c3);
  574    freeCmplx(c4);
  575    freeCmplx(c5);
  576    
  577    return c6;
  578 }
  579 
  580 Cmplx * acosCmplx(Cmplx *a){
  581    Cmplx *p, *c1;
  582 
  583    c1 = asinCmplx(a);
  584    p = subRealCmplx(realPi2, c1);
  585    freeCmplx(c1);
  586 
  587    return p;
  588 }
  589 
  590 Cmplx * atanCmplx(Cmplx *a){
  591    Real *r, *x, *y, *rsq;
  592    Real *r1, *r2, *r3, *r4, *r5, *r6, *r7;
  593    Cmplx *p = newCmplx();
  594 
  595    r = absCmplx(a);
  596    x = a->re;
  597    y = a->im;
  598    rsq = mulReal(r, r);
  599    freeReal(r);
  600 
  601    r1 = mulReal(x, realTwo);
  602    r2 = subReal(realOne, rsq);
  603    r3 = atan2Real(r1, r2);
  604    p->re = mulReal(r3, realHalf);
  605    freeReal(r1);
  606    freeReal(r2);
  607    freeReal(r3);
  608 
  609    r1 = mulReal(y, realTwo);
  610    r2 = addReal(rsq, realOne);
  611    r3 = addReal(r2, r1);
  612    r4 = subReal(r2, r1);
  613    r5 = divReal(r3, r4);
  614    r6 = lnReal(r5);
  615    r7 = mulReal(realHalf, realHalf);
  616    p->im = mulReal(r6, r7);
  617    freeReal(r1);
  618    freeReal(r2);
  619    freeReal(r3);
  620    freeReal(r4);
  621    freeReal(r5);
  622    freeReal(r6);
  623    freeReal(r7);
  624    freeReal(rsq);
  625 
  626 
  627    return p; 
  628 }
  629 
  630 
  631 /*************** MULTIPLY **********************/
  632 
  633 /* multiply 2 Cmplx numbers */
  634 Cmplx * mulCmplx(Cmplx *a, Cmplx *b){
  635    Real *r1, *r2, *r3, *r4;
  636    Cmplx *p = newCmplx();
  637 
  638    r1 = mulReal(a->re, b->re);
  639    r2 = mulReal(a->im, b->im);
  640 
  641    r3 = mulReal(a->re, b->im);
  642    r4 = mulReal(a->im, b->re);
  643 
  644    p->re = subEqReal(r1, r2);
  645    p->im = addEqReal(r3, r4);
  646 
  647    freeReal(r2);
  648    freeReal(r4);
  649 
  650    return p;
  651 }
  652 
  653 /* multiply a Cmplx by a Real number */
  654 Cmplx * mulCmplxReal(Cmplx *a, Real *b){
  655    Cmplx *p = newCmplx();
  656 
  657    p->re = mulReal(a->re, b);
  658    p->im = mulReal(a->im, b);
  659 
  660    return p;
  661 }
  662 
  663 /***************** DIVIDE ***********************/
  664    
  665 
  666 /* divide 2 Cmplx numbers */
  667 Cmplx * divCmplx(Cmplx *a, Cmplx *b){
  668    Real *r2, *r3;
  669    Real *r5, *r6;
  670    Real *r8, *r9;
  671    Cmplx *p = newCmplx();
  672 
  673    /* r3 = bre^2 + bim^2 */
  674    r3 = mulReal(b->re, b->re);
  675    r2 = mulReal(b->im, b->im);
  676    addEqReal(r3, r2);
  677 
  678    /* r6 = bre*are + bim*aim */
  679    r6 = mulReal(b->re, a->re);
  680    r5 = mulReal(b->im, a->im);
  681    addEqReal(r6, r5);
  682 
  683    /* r9 = bre*aim - bim*are */
  684    r9 = mulReal(b->re, a->im);
  685    r8 = mulReal(b->im, a->re);
  686    subEqReal(r9, r8);
  687 
  688    p->re = divEqReal(r6, r3);
  689    p->im = divEqReal(r9, r3);
  690 
  691    freeReal(r2); freeReal(r3);
  692    freeReal(r5); 
  693    freeReal(r8); 
  694 
  695    return p;
  696 }
  697 
  698 /* divide a Cmplx by a Real number */
  699 Cmplx * divCmplxReal(Cmplx *a, Real *b){
  700    Cmplx *p = newCmplx();
  701 
  702    p->re = divReal(a->re, b);
  703    p->im = divReal(a->im, b);
  704 
  705    return p;
  706 }
  707 
  708 /* divide a Real by a Cmplx number */
  709 Cmplx * divRealCmplx(Real *a, Cmplx *b){
  710    Real *r2, *r3, *r4, *r5, *r6;
  711    Cmplx *p = newCmplx();
  712 
  713    /* r3 = bre^2 + bim^2 */
  714    r3 = mulReal(b->re, b->re);
  715    r2 = mulReal(b->im, b->im);
  716    addEqReal(r3, r2);
  717 
  718    r4 = divReal(b->re, r3);
  719    r5 = divReal(b->im, r3);
  720 
  721    r6 = negReal(a);
  722 
  723    p->re = mulEqReal(r4, a);
  724    p->im = mulEqReal(r5, r6);
  725 
  726    freeReal(r2);
  727    freeReal(r3);
  728    freeReal(r6);
  729 
  730    return p;
  731 }
  732 
  733 
  734 
  735 
  736 
  737 
  738 /***************** ADD ***************************/
  739 
  740 
  741 /* add 2 Cmplx numbers */
  742 Cmplx * addCmplx(Cmplx *a, Cmplx *b){
  743    Cmplx *p = newCmplx();
  744 
  745    p->re = addReal(a->re, b->re);
  746    p->im = addReal(a->im, b->im);
  747 
  748    return p;
  749 }
  750 
  751 /* add a Cmplx and a Real number */
  752 Cmplx * addCmplxReal(Cmplx *a, Real *b){
  753    Cmplx *p = newCmplx();
  754 
  755    p->re = addReal(a->re, b);
  756    p->im = setRealReal(newReal(), a->im);
  757 
  758    return p;
  759 }
  760 
  761 
  762 /******************* SUBTRACT *******************/
  763 
  764 
  765 /* subtract 2 Cmplx numbers */
  766 Cmplx * subCmplx(Cmplx *a, Cmplx *b){
  767    Cmplx *p = newCmplx();
  768 
  769    p->re = subReal(a->re, b->re);
  770    p->im = subReal(a->im, b->im);
  771 
  772    return p;
  773 }
  774 
  775 /* subtract a Real from a Cmplx */
  776 Cmplx * subCmplxReal(Cmplx *a, Real *b){
  777    Cmplx *p = newCmplx();
  778 
  779    p->re = subReal(a->re, b);
  780    p->im = setRealReal(newReal(), a->im);
  781 
  782    return p;
  783 }
  784 
  785 /* subtract a Cmplx from a Real */
  786 Cmplx * subRealCmplx(Real *a, Cmplx *b){
  787    Cmplx *p = newCmplx();
  788 
  789    p->re = subReal(a, b->re);
  790    p->im = negReal(b->im);
  791 
  792    return p;
  793 }
  794