C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_complex.cpp
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 /* CVS $Id: l_complex.cpp,v 1.21 2014/01/30 17:23:46 cxsc Exp $ */
00024 
00025 #include "l_complex.hpp"  // corresponding header file
00026 #include "l_interval.hpp"
00027 
00028 namespace cxsc {
00029 
00030 // ---------------- Unary Operators -----------------------------------------
00031 l_complex operator-(const l_complex& x)
00032 {
00033     return l_complex(-x.re, -x.im);
00034 }
00035 
00036 l_complex operator+(const l_complex& x)
00037 {
00038     return x;
00039 }
00040 
00041 // ----------------- cdotprecision +/- l_complex ----------------------------
00042 
00043 cdotprecision operator+(const l_complex& lc, const cdotprecision& cd) throw() 
00044 { 
00045    return _cdotprecision(lc) + cd; 
00046 }
00047 
00048 cdotprecision operator+(const cdotprecision& cd, const l_complex& lc) throw() 
00049 { 
00050    return _cdotprecision(lc) + cd; 
00051 }   
00052 
00053 cdotprecision operator-(const l_complex& lc, const cdotprecision& cd) throw() 
00054 { 
00055    return _cdotprecision(lc) - cd; 
00056 }
00057 
00058 cdotprecision operator-(const cdotprecision& cd, const l_complex& lc) throw() 
00059 { 
00060    return cd - _cdotprecision(lc); 
00061 }
00062 
00063 // ------------------ dotprecision +/- l_complex ----------------------------
00064 
00065 cdotprecision operator+(const l_complex& lc, const dotprecision& cd) throw() 
00066 { 
00067    return _cdotprecision(lc) + cd; 
00068 }
00069 
00070 cdotprecision operator+(const dotprecision& cd, const l_complex& lc) throw() 
00071 { 
00072    return _cdotprecision(lc) + cd; 
00073 }   
00074 
00075 cdotprecision operator-(const l_complex& lc, const dotprecision& cd) throw() 
00076 { 
00077    return _cdotprecision(lc) - cd; 
00078 }
00079 
00080 cdotprecision operator-(const dotprecision& cd, const l_complex& lc) throw() 
00081 { 
00082    return cd - _cdotprecision(lc); 
00083 }
00084 
00085 // ---------------- l_complex * l_complex -----------------------------------
00086 l_complex operator * (const l_complex& a, const l_complex& b)
00087     throw()
00088 {
00089     l_real r, i; 
00090     dotprecision dot(0.0);
00091 
00092     accumulate(dot,a.re,b.re);
00093     accumulate(dot,-a.im,b.im);
00094     r = dot;
00095 
00096     dot = 0.0;
00097     accumulate(dot,a.im,b.re);
00098     accumulate(dot,a.re,b.im);
00099     i = dot;
00100 
00101     return( l_complex(r,i) );            
00102 }
00103 
00104 // ------------------ l_complex * complex -----------------------------------
00105 l_complex operator * (const l_complex& a, const complex& b)
00106     throw()
00107 {
00108     l_real r, i; 
00109     dotprecision dot(0.0);
00110     
00111     accumulate(dot,a.re,Re(b));
00112     accumulate(dot,-a.im,Im(b));
00113     r = dot;  // r is real-part
00114 
00115     dot = 0.0;
00116     accumulate(dot,a.im,Re(b));
00117     accumulate(dot,a.re,Im(b));
00118     i = dot;  // i is imaginary part
00119 
00120     return( l_complex(r,i) );
00121 }
00122 
00123 l_complex operator * ( const complex& b, const l_complex& a )
00124     throw()
00125 {
00126     l_real r, i; 
00127     dotprecision dot(0.0);
00128 
00129     accumulate(dot,a.re,Re(b));
00130     accumulate(dot,-a.im,Im(b));
00131     r = dot;  // r is real part
00132 
00133     dot = 0.0;
00134     accumulate(dot,a.im,Re(b));
00135     accumulate(dot,a.re,Im(b));
00136     i = dot;  // i is imaginary part
00137 
00138     return( l_complex(r,i) );
00139 }
00140 
00141 // ---------------- l_complex / l_complex -----------------------------------
00142 
00143 static const int maxexpo  = 1020;
00144 
00145 void skale_down_exp(int ex1, int ex2, int D, int& d1, int& d2)
00146 // l_interval x1,x2 sind Punktintervalle mit den Zweierexponenten
00147 // ex1,ex2 > -maxint, d.h. x1,x2 <> 0.0;
00148 // x1*x2 ist mit 2^D, D<0, so zu skalieren, dass fuer das Produkt
00149 // x1*x2 moeglichst wenig Stellen verloren gehen! x1 ist dazu mit
00150 // 2^d1 und x2 ist mit 2^d2 zu multiplizieren. Es muss also gelten:
00151 // d1+d2 = D<0;
00152 // Diese Funktion berechnet zu den gegebenen Zweierexponenten
00153 // ex1 und ex2 die Exponenten d1 und d2.
00154 // ex1 berechnet sich z.B. durch:  ex1 = expo(x1).
00155 // Blomquist, 25.10.2006;
00156   
00157 {
00158     bool change(false);
00159     int Diff, D1, c;
00160     d1 = 0;  d2 = 0;
00161 
00162     if (D<0)
00163     {
00164         if (ex2>ex1)
00165         { c = ex1;  ex1 = ex2;  ex2 = c;  change = true; }
00166         // ab jetzt gilt:  ex1 >= ex2;
00167         c = ex1 + D; 
00168         if (c<ex2)
00169         {
00170             Diff = ex2 - c;   D1 = Diff/2;
00171             d1 = D + D1;      d2 = D - d1; // d1+d2 == D<0;
00172         }
00173         else d1 = D;  // d2 = 0, s.o.
00174         if (change)
00175         { c = d1;  d1 = d2;  d2 = c; }
00176     }
00177 
00178 } // skale_down_exp(...)
00179 
00180 void skale_up_exp1(int ex1, int ex2, int& fillin, int& d1, int& d2)
00181 // l_interval x1,x2 sind Punktintervalle mit den Zweierexponenten
00182 // ex1,ex2 > -maxint, d.h. x1,x2 <> 0.0;
00183 // Das Maximum (ex1+ex2) der beiden gegebenen Exponentensummen
00184 // ist <= 1022; Die Punktintervalle x1,x2 sind so mit 2^d1 bzw. mit
00185 // 2^d2 zu multiplizieren, dass gilt.
00186 //     I.   2*|(x1*x2)| < MaxReal;
00187 //     II.  Bei dem Produkt x1*x2 duerfen moeglichst wenig Stellen 
00188 //          verloren gehen.
00189 // Blomquist, 25.10.2006;
00190 {
00191     bool change(false);
00192     int c, pot2;
00193     d1 = 0;  d2 = 0;
00194     fillin = maxexpo - (ex1+ex2);
00195     if (fillin>0)
00196     {
00197         if (ex2>ex1)
00198         { c = ex1;  ex1 = ex2;  ex2 = c;  change = true; }
00199         // ab jetzt gilt:  ex1 >= ex2;
00200         pot2 = maxexpo - ex2;  
00201         // Um pot2 kann x2 ohne Overflow hochskaliert werden.
00202         if (fillin <= pot2) d2 = fillin;
00203         else { d2 = pot2;  d1 = fillin - pot2; }
00204         if (change)
00205         { c = d1;  d1 = d2;  d2 = c; }
00206     }
00207 } // skale_up_exp1(...)
00208 
00209 void skale_up_exp2(int ex1, int ex2, int fillin, int& d1, int& d2)
00210 // l_interval x1,x2 sind Punktintervalle mit den Zweierexponenten
00211 // ex1,ex2 > -maxint, d.h. x1,x2 <> 0.0;
00212 // Das Minimum (ex1+ex2) der beiden gegebenen Exponentensummen
00213 // ist <= 1022; Die Punktintervalle x1,x2 sind so mit 2^d1 bzw. mit
00214 // 2^d2 zu multiplizieren, dass gilt.
00215 //     I.   2*|(x1*x2)| < MaxReal;
00216 //     II.  Bei dem Produkt x1*x2 duerfen moeglichst wenig Stellen 
00217 //          verloren gehen.
00218 // Blomquist, 25.10.2006;
00219 {
00220     bool change(false);
00221     int c, pot2;
00222     d1 = 0;  d2 = 0;
00223     if (fillin>0)
00224     {
00225         if (ex2>ex1)
00226         { c = ex1;  ex1 = ex2;  ex2 = c;  change = true; }
00227         // ab jetzt gilt:  ex1 >= ex2;
00228         pot2 = 1022 - ex2; // 1022 ??
00229         // Um pot2 kann x2 ohne Overflow hochskaliert werden.
00230         if (fillin <= pot2) d2 = fillin;
00231         else { d2 = pot2;  d1 = fillin - pot2; }
00232         if (change)
00233         { c = d1;  d1 = d2;  d2 = c; }
00234     }
00235 } // skale_up_exp2(...)
00236 
00237 void product(const l_real& a, const l_real& b, const l_real& c, 
00238               const l_real& d, int& ex, l_interval& res)
00239 // Calulation of an inclusion of: a*b + c*d
00240 {
00241     l_real a1(a), b1(b), c1(c), d1(d);
00242 
00243     a1 += 0.0;   b1 += 0.0;   c1 += 0.0;   d1 += 0.0;
00244     int ex_a1( expo(a1[1]) ), ex_b1( expo(b1[1]) ),
00245         ex_c1( expo(c1[1]) ), ex_d1( expo(d1[1]) ), m,p,D1,D2;
00246     l_interval a_i(a1),b_i(b1),c_i(c1),d_i(d1),li;  // point intervals!
00247     idotprecision Akku(0.0);
00248     ex = expo(0.0);  res = 0.0;  // Initialization for a*b + c*d == 0;
00249 
00250     if ( ex_a1 == ex || ex_b1 == ex ) // a*b == 0;
00251         if ( ex_c1 == ex || ex_d1 == ex ) ;  // a*b == c*d == 0.0;
00252         else 
00253         {// a*b == 0;  c*d != 0;
00254             m = ex_c1 + ex_d1;
00255             if (m > maxexpo) 
00256             {   
00257                 p = maxexpo - m;
00258                 skale_down_exp(ex_c1, ex_d1, p, D1, D2); 
00259                 Times2pown(c_i,D1);
00260                 Times2pown(d_i,D2);
00261                 Akku = 0.0;
00262                 accumulate(Akku,c_i,d_i);
00263                 res = Akku;
00264                 ex = -p;
00265             }
00266             else
00267             {
00268                 skale_up_exp1(ex_c1, ex_d1, p, D1, D2);
00269                 Times2pown(c_i,D1);
00270                 Times2pown(d_i,D2);
00271                 Akku = 0.0;
00272                 accumulate(Akku,c_i,d_i);
00273                 res = Akku;
00274                 ex = -p;
00275             }   
00276         }
00277     else // a*b != 0.0;
00278         if ( ex_c1 == ex || ex_d1 == ex )
00279         {// a*b != 0.0;  c*d == 0.0;
00280             m = ex_a1 + ex_b1;
00281             if (m > maxexpo)
00282             {
00283                 p = maxexpo - m;
00284                 skale_down_exp(ex_a1, ex_b1, p, D1, D2); 
00285                 Times2pown(a_i,D1);
00286                 Times2pown(b_i,D2);
00287                 Akku = 0.0;
00288                 accumulate(Akku,a_i,b_i);
00289                 res = Akku;
00290                 ex = -p;
00291             }
00292             else
00293             {
00294                 skale_up_exp1(ex_a1, ex_b1, p, D1, D2);
00295                 Times2pown(a_i,D1);
00296                 Times2pown(b_i,D2);
00297                 Akku = 0.0;
00298                 accumulate(Akku,a_i,b_i);
00299                 res = Akku;
00300                 ex = -p;
00301             }
00302         }
00303         else // a*b != 0.0  and  c*d != 0.0;  
00304         {
00305             if ( (ex_c1+ex_d1) > (ex_a1+ex_b1) )
00306             {
00307                 li = a_i;  a_i = c_i;  c_i = li;  // a_i <--> c_i
00308                 li = b_i;  b_i = d_i;  d_i = li;  // b_i <--> d_i
00309                 m = ex_a1;  ex_a1 = ex_c1;  ex_c1 = m;
00310                 m = ex_b1;  ex_b1 = ex_d1;  ex_d1 = m;
00311             }
00312             m = ex_a1 + ex_b1;
00313             // ab jetzt gilt:  m = (ex_a1+ex_b1) >= (ex_c1+ex_d1);
00314             if (m > maxexpo)
00315             {   
00316                 p = maxexpo - m;
00317                 skale_down_exp(ex_a1,ex_b1,p,D1,D2);
00318                 Times2pown(a_i,D1);
00319                 Times2pown(b_i,D2);
00320                 skale_down_exp(ex_c1,ex_d1,p,D1,D2);
00321                 Times2pown(c_i,D1);
00322                 Times2pown(d_i,D2);
00323                 Akku = 0.0;
00324                 accumulate(Akku,a_i,b_i);
00325                 accumulate(Akku,c_i,d_i);
00326                 res = Akku;
00327                 ex = -p; 
00328             }
00329             else
00330             {
00331                 skale_up_exp1(ex_a1, ex_b1, p, D1, D2);
00332                 Times2pown(a_i,D1);
00333                 Times2pown(b_i,D2);
00334                 skale_up_exp2(ex_c1, ex_d1, p, D1, D2);
00335                 Times2pown(c_i,D1);
00336                 Times2pown(d_i,D2);
00337                 Akku = 0.0;
00338                 accumulate(Akku,a_i,b_i);
00339                 accumulate(Akku,c_i,d_i);
00340                 res = Akku;
00341                 ex = -p; 
00342             }
00343         }
00344 } // product(...)
00345 
00346 void product(const l_real& c, const l_real& d, int& ex, l_interval& res)
00347 // Calulation of an inclusion of: c*c + d*d
00348 {
00349     l_real c1(c), d1(d);
00350 
00351     c1 += 0.0;   d1 += 0.0;
00352     int ex_c1( expo(c1[1]) ), ex_d1( expo(d1[1]) ), m;
00353 
00354     l_interval c_i(c1),d_i(d1),li;  // point intervals!
00355     idotprecision Akku(0.0);
00356     ex = expo(0.0);  res = 0.0;  // Initialization for c*c + d*d == 0;
00357 
00358     if (ex_c1 == ex) // c*c == 0;
00359         if (ex_d1 == ex) ; // c*c == d*d == 0; 
00360         else // c*c == 0;  d*d != 0;
00361         {
00362             times2pown(d_i,-ex_d1); // d_i about 1.0;
00363             Akku = 0.0;
00364             accumulate(Akku,d_i,d_i);
00365             res = Akku;
00366             ex = 2*ex_d1;
00367         }
00368     else // c*c != 0;
00369         if (ex_d1 == ex) 
00370         { // c*c != 0;   d*d == 0;
00371             times2pown(c_i,-ex_c1); // c_i about 1.0;
00372             Akku = 0.0;
00373             accumulate(Akku,c_i,c_i);
00374             res = Akku;
00375             ex = 2*ex_c1;
00376         }
00377         else // c*c != 0;   d*d != 0; 
00378         {
00379             if (ex_d1 > ex_c1)
00380             { 
00381                 li = c_i;  c_i = d_i;  d_i = li; // c_i <--> d_i
00382                 m = ex_c1; ex_c1 = ex_d1;  ex_d1 = m; 
00383             }  // c*c >= d*d:
00384             times2pown(c_i,-ex_c1); // c_i about 1.0;
00385             times2pown(d_i,-ex_c1);
00386             Akku = 0.0;
00387             accumulate(Akku,c_i,c_i);
00388             accumulate(Akku,d_i,d_i);
00389             res = Akku;
00390             ex = 2*ex_c1;
00391         }
00392 
00393 } // product(...)
00394 
00395 l_real quotient(const l_interval& z, const l_interval& n, int round, 
00396                 int ex_z, int ex_n)
00397 // z is an inclusion of a numerator.
00398 // n is an inclusion of a denominator.
00399 // quotient(...) calculates with q1 an approximation of z/n
00400 // using staggered arithmetic.
00401 // Rounding with round (-1,0,+1) is considered.
00402 
00403 { 
00404     l_real q1; // return value
00405     int ex_diff; 
00406     l_interval res;
00407 
00408     if (0.0<=n) // 0 in denominator n leads to error message:
00409     {
00410         std::cerr << "quotient1(const l_interval& z, const l_interval& n, int round, int ex_z, int ex_n):  Division by zero" << std::endl;
00411         exit(1);
00412     }
00413 
00414     if ( zero_(z) ) { q1=0;  return q1; };
00415 
00416     ex_diff = ex_z - ex_n;
00417     res = z/n;
00418     Times2pown(res,ex_diff);
00419     switch(round)
00420     {
00421         case RND_DOWN:
00422             q1 = Inf(res);
00423             break;
00424         case RND_NEXT:
00425             q1 = mid(res);
00426             break;
00427         case RND_UP:
00428             q1 = Sup(res);
00429             break;
00430     } // switch
00431     return q1;
00432 } // quotient
00433 
00434 l_complex _c_division(l_complex a, l_complex b, int round)
00435 {
00436     int ex1, ex2;
00437     l_interval z, n;
00438     l_complex tmp;
00439 
00440     product(Re(b),Im(b),ex2,n);
00441     product(Re(a),Re(b),Im(a),Im(b),ex1,z);
00442     SetRe(tmp, quotient(z,n,round,ex1,ex2));
00443     product(Im(a),Re(b),-Re(a),Im(b),ex1,z);
00444     SetIm(tmp, quotient(z,n,round,ex1,ex2));
00445     return tmp;
00446 } // _c_division
00447 
00448 l_complex divn (const l_complex & a, const l_complex & b)  
00449 {  
00450    return _c_division(a,b,RND_NEXT);
00451 }
00452 
00453 l_complex divd (const l_complex & a, const l_complex & b) 
00454 {  
00455    return _c_division(a,b,RND_DOWN);
00456 }
00457 
00458 l_complex divu (const l_complex & a, const l_complex & b)  
00459 {  
00460    return _c_division(a,b,RND_UP);
00461 }
00462 
00463 l_complex operator / (const l_complex &a, const l_complex &b) throw()
00464 {
00465    return divn(a,b);
00466 }
00467 
00468 int StagPrec(const l_complex& lc) throw()
00469 {
00470     return StagPrec(lc.re);
00471 }      
00472 
00473 void accumulate(cdotprecision& cd, const l_complex& lc1, 
00474                                    const l_complex& lc2) throw()
00475 {
00476     accumulate(Re(cd),lc1.re,lc2.re); 
00477     accumulate(Re(cd),-lc1.im,lc2.im);
00478     accumulate(Im(cd),lc1.im,lc2.re);
00479     accumulate(Im(cd),lc1.re,lc2.im);
00480 }
00481 
00482 void accumulate(cdotprecision& cd, const l_complex& lc, 
00483                                    const complex& c) throw()
00484 {
00485     accumulate(Re(cd),lc.re,Re(c)); 
00486     accumulate(Re(cd),-lc.im,Im(c));
00487     accumulate(Im(cd),lc.im,Re(c));
00488     accumulate(Im(cd),lc.re,Im(c));
00489 }
00490 
00491 void accumulate(cdotprecision& cd, const l_complex& lc, 
00492                                    const real& r) throw()
00493 {
00494     accumulate(Re(cd),lc.re,r); 
00495     accumulate(Im(cd),lc.im,r);
00496 }
00497 
00498 void accumulate(cdotprecision& cd, const l_complex& lc, 
00499                                    const l_real& lr) throw()
00500 {
00501     accumulate(Re(cd),lc.re,lr); 
00502     accumulate(Im(cd),lc.im,lr);
00503 }
00504 
00505 l_real abs2(const l_complex &a) throw()
00506 {
00507    dotprecision dot(0.0);
00508    accumulate(dot,a.re,a.re);
00509    accumulate(dot,a.im,a.im);
00510    return l_real(dot);
00511 }
00512 
00513 l_real abs(const l_complex& z) throw()
00514 // Calculation of an approximation of |z|.
00515 // In general the maximum precision is stagprec=19, predifined by the used
00516 // sqrt-function declared in l_rmath.hpp.
00517 // If the difference |exa-exb| of the exponents (base 2) is sufficient high,
00518 // precision and accuracy can be choosen greater 19.
00519 {
00520     return sqrtx2y2(Re(z),Im(z));
00521 } // abs(...)
00522 
00523 l_complex & SetIm(l_complex & a,const l_real & b) 
00524         { a.im=b; return a; } // See SetRe(...);
00525 
00526 l_complex & SetRe(l_complex & a,const l_real & b)
00527         { a.re=b; return a; } // The real part of a is substituted by b.
00528 
00529 l_real & Re(l_complex& a) { return a.re; }
00530 l_real   Re(const l_complex& a) { return a.re; }
00531 l_real & Im(l_complex& a) { return a.im; }
00532 l_real   Im(const l_complex& a) { return a.im; }
00533 
00534 complex & complex::operator = (const l_complex& a) throw()
00535 {
00536         real x(Re(a)), y(Im(a));
00537         return *this = complex(x,y);
00538 }
00539 
00540 } // end namespace cxsc
00541 
00542 
00543 
00544 
00545 
00546