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 /* 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