C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
complex.inl
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: complex.inl,v 1.26 2014/01/30 18:13:52 cxsc Exp $ */
00025 #include "cinterval.hpp"
00026 #include "idot.hpp"
00027 
00028 namespace cxsc {
00029 // ---- Constructors ----------------------------------------------
00030 
00031 inline complex & complex::operator= (const real & r) throw()
00032 {
00033    re=r;im=0;
00034    return *this;
00035 }
00036 
00037       // ---- Std.Operators ---------------------------------------
00038 inline complex operator -(const complex &a) throw () 
00039 {
00040    return complex(-a.re,-a.im);
00041 }
00042 
00043 inline complex operator +(const complex &a) throw ()
00044 {
00045    return a;
00046 }
00047 
00048 inline complex operator +(const complex &a,const complex &b) throw()
00049 {
00050    return complex(a.re+b.re,a.im+b.im);
00051 }
00052 
00053 inline complex operator -(const complex &a,const complex &b) throw()
00054 {
00055    return complex(a.re-b.re,a.im-b.im);
00056 }
00057 
00058 inline complex & operator +=(complex &a, const complex &b) throw() { return a=a+b; }
00059 inline complex & operator -=(complex &a, const complex &b) throw() { return a=a-b; }
00060 inline complex & operator *=(complex &a, const complex &b) throw() { return a=a*b; }
00061 inline complex & operator /=(complex &a, const complex &b) throw() { return a=a/b; }
00062 
00063 inline complex operator +(const complex & a,const real & b) throw() 
00064 { 
00065    return complex(a.re+b,a.im);
00066 }
00067 
00068 inline complex operator +(const real & a,const complex & b) throw()
00069 {
00070    return complex(a+b.re,b.im);
00071 }
00072 
00073 inline complex operator -(const complex & a,const real & b) throw()
00074 {
00075    return complex(a.re-b,a.im);
00076 }
00077 
00078 inline complex operator -(const real & a,const complex & b) throw()
00079 {
00080    return complex(a-b.re,-b.im);
00081 }
00082 
00083 inline complex operator *(const complex & a,const real & b) throw()
00084 {
00085 //   return a*_complex(b);
00086      return complex(a.re*b,a.im*b);  // Blomquist, 07.11.02;
00087 }
00088 
00089 inline complex operator *(const real & a,const complex & b) throw()
00090 {
00091 //   return _complex(a)*b;
00092      return complex(a*b.re, a*b.im);  // Blomquist, 07.11.02;
00093 }
00094 
00095 inline complex operator /(const complex & a,const real & b) throw()
00096 {
00097 //   return a/_complex(b);
00098      return complex(a.re/b, a.im/b);  // Blomquist, 07.11.02;
00099 }
00100 
00101 inline complex operator /(const real & a,const complex & b) throw()
00102 {
00103    return _complex(a)/b;
00104 }
00105 
00106 inline complex & operator +=(complex & a, const real & b) throw() { return a=a+b; }
00107 inline complex & operator -=(complex & a, const real & b) throw() { return a=a-b; }
00108 inline complex & operator *=(complex & a, const real & b) throw() { return a=a*b; }
00109 inline complex & operator /=(complex & a, const real & b) throw() { return a=a/b; }
00110 
00111       // ---- Comp.Operat.  ---------------------------------------
00112 inline bool operator!  (const complex & a)                    throw() { return !a.re && !a.im; }
00113 inline bool operator== (const complex & a, const complex & b) throw() { return a.re==b.re && a.im==b.im; }
00114 inline bool operator!= (const complex & a, const complex & b) throw() { return a.re!=b.re || a.im!=b.im; }
00115 inline bool operator== (const complex & a, const real & b)    throw() { return !a.im && a.re==b; }
00116 inline bool operator== (const real & a, const complex & b)    throw() { return !b.im && a==b.re; }
00117 inline bool operator!= (const complex & a, const real & b)    throw() { return !!a.im || a.re!=b; }
00118 inline bool operator!= (const real & a, const complex & b)    throw() { return !!b.im || a!=b.re; }
00119 
00120       // ---- Others   -------------------------------------------
00121 
00122 inline complex conj(const complex & a) throw() { return complex(a.re,-a.im); }
00123 
00124 
00125 // ----------- Directed Rounding, Blomquist -------------------------------
00126 // ------------------------------------------------------------------------
00127 
00128    // -------------------- addition --------------------------------
00129 
00130 inline complex addd(const complex& a, const complex& b) throw()
00131 { return complex(addd(a.re,b.re), addd(a.im,b.im)); }
00132 
00133 inline complex addu(const complex& a, const complex& b) throw()
00134 { return complex(addu(a.re,b.re), addu(a.im,b.im)); }
00135 
00136 inline complex addd(const complex& a, const real& b) throw()
00137 { return complex(addd(a.re,b), a.im); }
00138 
00139 inline complex addu(const complex& a, const real& b) throw()
00140 { return complex(addu(a.re,b), a.im); }
00141 
00142 inline complex addd(const real& a, const complex& b) throw()
00143 { return complex(addd(a,b.re), b.im); }
00144 
00145 inline complex addu(const real& a, const complex& b) throw()
00146 { return complex(addu(a,b.re), b.im); }
00147    // ----------------- subtraction: ----------------------------
00148 
00149 inline complex subd(const complex& a, const complex& b) throw()
00150 { return complex(subd(a.re,b.re), subd(a.im,b.im)); }
00151 
00152 inline complex subu(const complex& a, const complex& b) throw()
00153 { return complex(subu(a.re,b.re), subu(a.im,b.im)); }
00154 
00155 inline complex subd(const complex& a, const real& b) throw()
00156 { return complex(subd(a.re,b), a.im); }
00157 
00158 inline complex subu(const complex& a, const real& b) throw()
00159 { return complex(subu(a.re,b), a.im); }
00160 
00161 inline complex subd(const real& a, const complex& b) throw()
00162 { return complex(subd(a,b.re), -b.im); }
00163 
00164 inline complex subu(const real& a, const complex& b) throw()
00165 { return complex(subu(a,b.re), -b.im); }
00166 
00167    // --------------- multiplikation ------------------------
00168 
00169 inline complex muld(const complex &a, const real &b) throw()
00170 { return complex( muld(a.re,b), muld(a.im,b) ); }
00171 
00172 inline complex mulu(const complex &a, const real &b) throw()
00173 { return complex( mulu(a.re,b), mulu(a.im,b) ); }
00174 
00175 inline complex muld(const real &a, const complex &b) throw()
00176 { return complex( muld(a,b.re), muld(a,b.im) ); }
00177 
00178 inline complex mulu(const real &a, const complex &b) throw()
00179 { return complex( mulu(a,b.re), mulu(a,b.im) ); }
00180 
00181    // -------------- division ---------------------------------
00182 
00183 inline complex divd(const complex &a, const real &b) throw()
00184 { return complex( divd(a.re,b), divd(a.im,b) ); }
00185 
00186 inline complex divu(const complex &a, const real &b) throw()
00187 { return complex( divu(a.re,b), divu(a.im,b) ); }
00188 
00189 inline complex divd(const real &a, const complex &b) throw()
00190 { return divd(_complex(a),b); }
00191 
00192 inline complex divu(const real &a, const complex &b) throw()
00193 { return divu(_complex(a),b); }
00194 
00195 inline complex operator *(const complex &a,const complex &b) throw()
00196 {
00197 #ifdef CXSC_FAST_COMPLEX_OPERATIONS
00198    return complex(Re(a)*Re(b)-Im(a)*Im(b), Re(a)*Im(b)+Im(a)*Re(b));
00199 #else
00200    complex tmp;
00201    dotprecision dot(0.0);
00202     
00203    accumulate (dot,  a.re, b.re);
00204    accumulate (dot, -a.im, b.im);
00205    rnd (dot, tmp.re, RND_NEXT);
00206 
00207    dot = 0.0;
00208    accumulate (dot, a.re, b.im);
00209    accumulate (dot, a.im, b.re);
00210    rnd (dot, tmp.im, RND_NEXT);
00211 
00212    return tmp;
00213 #endif
00214 }
00215 
00216 inline complex muld(const complex &a, const complex &b) throw()
00217 {  // Blomquist 07.11.02;
00218    complex tmp;
00219    dotprecision dot(0.0);
00220     
00221    accumulate (dot,  a.re, b.re);
00222    accumulate (dot, -a.im, b.im);
00223    rnd (dot, tmp.re, RND_DOWN);
00224 
00225    dot = 0.0;
00226    accumulate (dot, a.re, b.im);
00227    accumulate (dot, a.im, b.re);
00228    rnd (dot, tmp.im, RND_DOWN);
00229 
00230    return tmp;
00231 }
00232 
00233 inline complex mulu(const complex &a, const complex &b) throw()
00234 {  // Blomquist 07.11.02;
00235    complex tmp;
00236    dotprecision dot(0.0);
00237     
00238    accumulate (dot,  a.re, b.re);
00239    accumulate (dot, -a.im, b.im);
00240    rnd (dot, tmp.re, RND_UP);
00241 
00242    dot = 0.0;
00243    accumulate (dot, a.re, b.im);
00244    accumulate (dot, a.im, b.re);
00245    rnd (dot, tmp.im, RND_UP);
00246 
00247    return tmp;
00248 }
00249 
00250 
00251 static const int Min_Exp_ = 1074, minexpom = -914, 
00252                  maxexpo1 = 1022, MANT_W   = 52;
00253 
00254 inline void product(real a, real b, real c, real d,
00255              int& overfl, real& p1, interval& p2)
00256 // New version of function product(...) from Blomquist, 26.10.02;
00257 // Input data: a,b,c,d;  Output data: overfl, p1, p2;
00258 // In case of overfl=0 the interval p1+p2 is an inclusion of a*b + c*d;
00259 // overfl=1 (overflow) signalizes that p1+p2 must be multiplied with 2^1074
00260 // to be an inclusion of a*b + c*d;
00261 // overfl=-1 (underflow) signalizes that p1+p2 must be multiplied with 
00262 // 2^-1074 to be an inclusion of a*b + c*d; 
00263 
00264 {  
00265    int exa, exb, exc, exd;  // Exp. von a-d
00266    dotprecision dot;
00267    int inexact;
00268 
00269    overfl  = 0;
00270    inexact = 0; // false
00271 
00272    dot = 0.0;
00273    exa = expo(a);
00274    exb = expo(b);
00275    exc = expo(c);
00276    exd = expo(d);
00277 
00278    if ( sign(a) == 0  ||   sign(b) == 0 )    //  a * b == 0.0
00279       if ( sign(c) == 0  ||  sign(d) == 0 )  //  a * b == c * d == 0
00280                                              //  dot := #(0);
00281          ;    // No Operation necessary
00282       else 
00283       {
00284          //  a * b == 0;  c * d != 0;
00285          if (exc+exd > maxexpo1) 
00286          {
00287             //  overflow !
00288              if ( exc > exd ) c = comp( mant(c), exc-Min_Exp_ );
00289              else d = comp( mant(d), exd-Min_Exp_ );
00290              overfl = 1;
00291          } else 
00292          if  ( exc+exd < minexpom ) 
00293          { // undeflow; 
00294            // c = comp( mant(c), exc+Min_Exp_ ); kann Overflow erzeugen!
00295              if (exc < exd) c = comp( mant(c), exc+Min_Exp_ );
00296              else d = comp( mant(d), exd+Min_Exp_ );
00297              overfl = -1;
00298          }
00299          accumulate(dot,c,d);
00300       } else //  a,b != 0
00301       if ( sign(c) == 0  ||  sign(d) == 0 ) 
00302       {
00303          //  a*b != 0, c * d == 0
00304          if (exa+exb > maxexpo1) 
00305          {
00306             //  overflow !
00307             if ( exa > exb ) a = comp( mant(a), exa-Min_Exp_ );
00308             else b = comp( mant(b), exb-Min_Exp_ );
00309             overfl = 1;
00310          } else 
00311          if (exa+exb < minexpom) 
00312          { // undeflow; 
00313            // a = comp( mant(a), exa+Min_Exp_ ); kann Overflow erzeugen!
00314              if (exa < exb) a = comp( mant(a), exa+Min_Exp_ );
00315              else b = comp( mant(b), exb+Min_Exp_ );
00316             overfl = -1;
00317          }
00318          accumulate(dot,a,b);
00319       } else 
00320       {
00321          // a,b,c,d != 0
00322          if (exa+exb > maxexpo1) 
00323          {  //  overflow bei a*b
00324             if ( exa > exb ) a = comp( mant(a), exa-Min_Exp_ );
00325             else b = comp( mant(b), exb-Min_Exp_ );
00326             if (exc > MANT_W) c = comp( mant(c), exc-Min_Exp_ );
00327             else if (exd > MANT_W)
00328                d = comp( mant(d), exd-Min_Exp_ );
00329             else 
00330             {
00331                // underflow wegen Skalierung bei c*d
00332                 c = 0.0;
00333                 inexact = 1; // true
00334             }
00335             overfl = 1;  // Hat vorher gefehlt!! Blomquist, 24.10.02;
00336          } else if (exc+exd > maxexpo1) 
00337          {
00338             // overflow bei c*d
00339             if ( exc > exd ) c = comp( mant(c), exc-Min_Exp_ );
00340             else d = comp( mant(d), exd-Min_Exp_ );
00341             if (exa > MANT_W) a = comp( mant(a), exa-Min_Exp_ );
00342             else if (exb > MANT_W)
00343                b = comp( mant(b), exb-Min_Exp_ );
00344             else 
00345             {
00346                // underflow wegen Skalierung bei a*b
00347                a = 0.0;
00348                inexact = 1; // true
00349             }
00350             overfl = 1;
00351          } else 
00352          if ( exa+exb < minexpom  &&  exc+exd < minexpom ) 
00353          {
00354             //  underflow bei a*b und bei c*d
00355              if (exa < exb) a = comp( mant(a), exa+Min_Exp_ );
00356              else b = comp( mant(b), exb+Min_Exp_ );
00357              if (exc < exd) c = comp( mant(c), exc+Min_Exp_ );
00358              else d = comp( mant(d), exd+Min_Exp_ );
00359              overfl = -1;
00360          }
00361          accumulate(dot, a, b);
00362          accumulate(dot, c, d);
00363       }
00364 
00365    p1 = rnd(dot);
00366    dot -= p1;
00367    rnd(dot,p2);  // Blomquist, 07.11.02;
00368 
00369    if (inexact)
00370        p2 = interval( pred(Inf(p2)), succ(Sup(p2)) );
00371 
00372 } // end product
00373 
00374 inline real quotient (real z1, interval z2, real n1, 
00375                interval n2, int round, int zoverfl, int noverfl)
00376 // z1+z2 is an inclusion of a numerator.
00377 // n1+n2 is an inclusion of a denominator.
00378 // quotient(...) calculates with q1 an approximation of (z1+z2)/(n1+n2)
00379 // using staggered arithmetic.
00380 // Rounding with round (-1,0,+1) is considered.
00381 // zoverfl and noverfl are considered by scaling back with 2^1074 or 2^-1074
00382 // n1+n2 > 0 is assumed.
00383 {
00384    real q1=0, scale;
00385    interval q2, nh;
00386    idotprecision id;
00387    int vorz, anz_scale, ex = 0;
00388 
00389    vorz = sign(z1); // n1,n2 > 0 is assumed!!
00390                     // so the sign of q1 ist sign of z1.
00391    if ( zoverfl == -1  &&  noverfl == 1 ) 
00392    {
00393       //  result in the undeflow range:
00394       switch (round) 
00395       {
00396          case RND_DOWN:
00397             if (vorz >= 0) 
00398                 q1 = 0.0;
00399             else         
00400                 q1 = -minreal; // Blomquist: MinReal --> minreal;
00401             break;
00402          case RND_NEXT:
00403             q1 = 0.0;
00404             break;
00405          case RND_UP:
00406             if (vorz <= 0) 
00407                 q1 = 0.0;
00408             else         
00409                 q1 = minreal;  // Blomquist: MinReal --> minreal;
00410             break;
00411       } // switch
00412    } else if ( zoverfl==1  &&  noverfl==-1 ) 
00413    {  
00414       //  result in the overflow range:
00415       if (vorz >= 0) q1 = MaxReal+MaxReal;  // Overflow
00416       else q1 = -MaxReal-MaxReal;           // Overflow
00417    } else 
00418    {  
00419       q1 = divd(z1, n1);  // down, to get q2 >= 0
00420       nh = interval(addd(n1, Inf(n2)),
00421                    addu(n1, Sup(n2)));
00422 
00423       // q2:= ##( z1 - q1*n1 + z2 - q1*n2 );
00424       id = z1;
00425       accumulate(id, -q1, n1);
00426       id += z2;
00427       accumulate(id, -q1, n2);
00428       q2 = rnd(id);
00429 
00430       switch (round) // Considering the rounding before scaling
00431       {
00432          case RND_DOWN:
00433             q1 = adddown(q1, divd(Inf(q2), Sup(nh)));
00434             break;
00435          case RND_NEXT:
00436             q1 = q1 + (Inf(q2)+Sup(q2))*0.5/n1;
00437             break;
00438          case RND_UP:
00439             q1 = addup(q1, divu(Sup(q2), Inf(nh)));
00440             break;
00441       } // switch
00442 
00443 
00444       //  scaling back, if numerator|denominator - over-|underflow :
00445 
00446       //  actuell as follows:
00447       //  q1:= comp( mant(q1), expo(q1) + (zoverfl-noverfl)*1074 );
00448       //  The scaling with 2^1074 must be done in two steps with 2^ex
00449       //  and with the factor scale:
00450 
00451       anz_scale = zoverfl - noverfl; // |anz_scale| <= 1; 
00452       if (anz_scale > 0)
00453       {
00454           scale = comp(0.5, +1024);
00455           ex = MANT_W-1;
00456       }
00457       else if (anz_scale < 0)
00458       {
00459           scale = comp(0.5, -1022);
00460           ex = -MANT_W+1;
00461       }
00462       else scale = 1.0;  // ex = 0 is already initialized for this case 
00463 
00464 
00465       if (ex) times2pown(q1,ex); // EXACT part scaling, if ex != 0.
00466       switch (round) // correct rounding with the second factor scale:
00467       {
00468           case RND_DOWN:
00469               q1 = multdown(q1, scale);
00470               break;
00471           case RND_NEXT:
00472               q1 = q1 * scale;
00473               break;
00474           case RND_UP:
00475               q1 = multup(q1, scale);
00476               break;
00477       }  // switch
00478    }
00479    return q1;
00480 } // end of quotient
00481 
00482 inline complex _c_division(complex a, complex b, int round) throw(DIV_BY_ZERO)
00483 {
00484     if (0.0 == (sqr(Re(b))+sqr(Im(b)))) {
00485       cxscthrow(DIV_BY_ZERO("complex operator / (const complex&,const complex&)"));
00486     }      
00487     
00488    int zoverflow, noverflow;
00489    real z1, n1;
00490    interval z2, n2;
00491    complex tmp;
00492 
00493    product (Re(b), Re(b), Im(b), Im(b), noverflow, n1, n2);
00494    product (Re(a), Re(b), Im(a), Im(b), zoverflow, z1, z2);
00495    SetRe (tmp, quotient (z1, z2, n1, n2, round, zoverflow, noverflow));
00496    product (Im(a), Re(b), -Re(a), Im(b), zoverflow, z1, z2);
00497    SetIm (tmp, quotient (z1, z2, n1, n2, round, zoverflow, noverflow));
00498    return tmp;
00499 }
00500 
00501 inline complex divn (const complex & a, const complex & b)  
00502 {  // Blomquist: vorher c_divd(...), 07.11.02;
00503    return _c_division(a, b, RND_NEXT);
00504 }
00505 
00506 inline complex divd (const complex & a, const complex & b) 
00507 {  // Blomquist: vorher c_divd(...), 07.11.02;
00508    return _c_division(a, b, RND_DOWN);
00509 }
00510 
00511 inline complex divu (const complex & a, const complex & b)  
00512 {  // Blomquist: vorher c_divu(...), 07.11.02;
00513    return _c_division(a, b, RND_UP);
00514 }
00515 
00516 inline complex operator / (const complex &a,const complex &b) throw()
00517 {
00518 #ifdef CXSC_FAST_COMPLEX_OPERATIONS
00519    real q = Re(b)*Re(b) + Im(b)*Im(b);
00520    return complex((Re(a)*Re(b)+Im(a)*Im(b))/q, (Im(a)*Re(b)-Re(a)*Im(b))/q);
00521 #else
00522    return divn(a,b);
00523 #endif
00524 }
00525 
00526 inline real abs2(const complex &a) throw()
00527 {
00528    dotprecision dot(0.0);
00529    accumulate(dot,a.re,a.re);
00530    accumulate(dot,a.im,a.im);
00531    return rnd(dot);
00532 }
00533 
00534 inline real abs (complex z) throw()
00535 {   //  calculation of |z|; Blomquist 06.12.02;
00536 #ifdef CXSC_FAST_COMPLEX_OPERATIONS
00537     return sqrt(Re(z)*Re(z)+Im(z)*Im(z));
00538 #else
00539     return sqrtx2y2(Re(z),Im(z));
00540 #endif
00541 }
00542 
00543 
00544 } // namespace cxsc