C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
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