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: interval.inl,v 1.31 2014/01/30 17:23:45 cxsc Exp $ */ 00025 00026 namespace cxsc { 00027 00028 // ---- Konstruktoren ---- 00029 00030 inline interval::interval(const real &a,const real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00031 : inf(a), sup(b) 00032 { 00033 if (a > b) 00034 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval::interval(const real &a,const real &b)")); 00035 } 00036 00037 /*inline interval::interval(int a,int b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00038 : inf(a), sup(b) 00039 { 00040 if (a > b) 00041 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval::interval(int a,int b)")); 00042 } 00043 00044 inline interval::interval(const double & a,const double & b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00045 : inf(a), sup(b) 00046 { 00047 if (a > b) 00048 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval::interval(const double & a,const double & b)")); 00049 } 00050 */ 00051 00052 inline interval& interval::operator= (const real& a) 00053 { 00054 inf = sup = a; 00055 return *this; 00056 } 00057 00058 00059 // ---- Typwandlungen ---- 00060 00066 inline interval _unchecked_interval(const real& a, const real& b) 00067 { 00068 interval tmp; 00069 tmp.inf = a; 00070 tmp.sup = b; 00071 return tmp; 00072 } 00073 00074 00075 // ---- Standardfunkt ---- (arithmetische Operatoren) 00076 00077 inline interval operator-(const interval &a) throw() { return interval (-a.sup, -a.inf); } 00078 inline interval operator+(const interval &a) throw() { return a; } 00079 00080 inline interval & operator +=(interval &a,const interval &b) throw() { return a=a+b; } 00081 inline interval & operator +=(interval &a,const real &b) throw() { return a=a+b; } 00082 inline interval & operator -=(interval &a,const interval &b) throw() { return a=a-b; } 00083 inline interval & operator -=(interval &a,const real &b) throw() { return a=a-b; } 00084 inline interval & operator *=(interval &a,const interval &b) throw() { return a=a*b; } 00085 inline interval & operator *=(interval &a,const real &b) throw() { return a=a*b; } 00086 inline interval & operator /=(interval &a,const interval &b) throw() { return a=a/b; } 00087 inline interval & operator /=(interval &a,const real &b) throw() { return a=a/b; } 00088 00089 00090 inline interval operator |(const interval &a,const interval &b) throw() 00091 { 00092 return interval((a.inf<b.inf)?a.inf:b.inf,(a.sup>b.sup)?a.sup:b.sup); 00093 } 00094 inline interval operator &(const interval &a,const interval &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00095 { 00096 return interval((a.inf>b.inf)?a.inf:b.inf,(a.sup<b.sup)?a.sup:b.sup); 00097 } 00098 inline interval operator |(const real &a,const interval &b) throw() 00099 { 00100 return interval((a<b.inf)?a:b.inf,(a>b.sup)?a:b.sup); 00101 } 00102 inline interval operator &(const real &a,const interval &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00103 { 00104 return interval((a>b.inf)?a:b.inf,(a<b.sup)?a:b.sup); 00105 } 00106 inline interval operator |(const interval &a,const real &b) throw() 00107 { 00108 return interval((a.inf<b)?a.inf:b,(a.sup>b)?a.sup:b); 00109 } 00110 inline interval operator |(const real &a,const real &b) throw() 00111 { 00112 if(a>b) return interval(b,a); 00113 else return interval(a,b); 00114 } 00115 inline interval operator &(const interval &a,const real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00116 { 00117 return interval((a.inf>b)?a.inf:b,(a.sup<b)?a.sup:b); 00118 } 00119 inline interval & operator |=(interval &a,const interval &b) throw() 00120 { 00121 a.inf=(a.inf<b.inf)?a.inf:b.inf,a.sup=(a.sup>b.sup)?a.sup:b.sup; 00122 return a; 00123 } 00124 inline interval & operator &=(interval &a,const interval &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00125 { 00126 a.inf=(a.inf>b.inf)?a.inf:b.inf,a.sup=(a.sup<b.sup)?a.sup:b.sup; 00127 if(a.inf>a.sup) 00128 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval & operator &=(interval &a,const interval &b)")); 00129 return a; 00130 } 00131 inline interval & operator |=(interval &a,const real &b) throw() 00132 { 00133 a.inf=(a.inf<b)?a.inf:b,a.sup=(a.sup>b)?a.sup:b; 00134 return a; 00135 } 00136 inline interval & operator &=(interval &a,const real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00137 { 00138 a.inf=(a.inf>b)?a.inf:b,a.sup=(a.sup<b)?a.sup:b; 00139 if(a.inf>a.sup) 00140 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval & operator &=(interval &a,const real &b)")); 00141 return a; 00142 } 00143 00144 // --- Vergleichsoperationen ---- 00145 inline bool operator ==(const interval &a,const interval &b) throw() { return(a.inf==b.inf && a.sup==b.sup); } 00146 inline bool operator !=(const interval &a,const interval &b) throw() { return(a.inf!=b.inf || a.sup!=b.sup); } 00147 inline bool operator ==(const real &r,const interval &a) throw() { return(r==a.inf && r==a.sup); } 00148 inline bool operator !=(const real &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); } 00149 inline bool operator ==(const interval &a,const real &r) throw() { return(r==a.inf && r==a.sup); } 00150 inline bool operator !=(const interval &a,const real &r) throw() { return(r!=a.inf || r!=a.sup); } 00151 00152 inline bool operator ==(const int &r,const interval &a) throw() { return(r==a.inf && r==a.sup); } 00153 inline bool operator !=(const int &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); } 00154 inline bool operator ==(const interval &a,const int &r) throw() { return(r==a.inf && r==a.sup); } 00155 inline bool operator !=(const interval &a,const int &r) throw() { return(r!=a.inf || r!=a.sup); } 00156 00157 inline bool operator ==(const long &r,const interval &a) throw() { return(r==a.inf && r==a.sup); } 00158 inline bool operator !=(const long &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); } 00159 inline bool operator ==(const interval &a,const long &r) throw() { return(r==a.inf && r==a.sup); } 00160 inline bool operator !=(const interval &a,const long &r) throw() { return(r!=a.inf || r!=a.sup); } 00161 00162 inline bool operator ==(const double &r,const interval &a) throw() { return(r==a.inf && r==a.sup); } 00163 inline bool operator !=(const double &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); } 00164 inline bool operator ==(const interval &a,const double &r) throw() { return(r==a.inf && r==a.sup); } 00165 inline bool operator !=(const interval &a,const double &r) throw() { return(r!=a.inf || r!=a.sup); } 00166 00167 // --- Mengenvergleiche --- 00168 // <,>,... 00169 inline bool operator <=(const interval &a,const interval &b) throw() 00170 { 00171 return(a.inf>=b.inf && a.sup<=b.sup); 00172 } 00173 inline bool operator >=(const interval &a,const interval &b) throw() 00174 { 00175 return(a.inf<=b.inf && a.sup>=b.sup); 00176 } 00177 inline bool operator <(const interval &a,const interval &b) throw() 00178 { 00179 return(a.inf>b.inf && a.sup<b.sup); 00180 } 00181 inline bool operator >(const interval &a,const interval &b) throw() 00182 { 00183 return(a.inf<b.inf && a.sup>b.sup); 00184 } 00185 00186 inline bool operator <=(const real &a,const interval &b) throw() 00187 { 00188 return(a>=b.inf && a<=b.sup); 00189 } 00190 inline bool operator >=(const real &a,const interval &b) throw() 00191 { 00192 return(a<=b.inf && a>=b.sup); 00193 } 00194 inline bool operator <(const real &a,const interval &b) throw() 00195 { 00196 return(a>b.inf && a<b.sup); 00197 } 00198 00199 inline bool operator <=(const interval &a,const real &b) throw() 00200 { 00201 return(a.inf>=b && a.sup<=b); 00202 } 00203 inline bool operator >=(const interval &a,const real &b) throw() 00204 { 00205 return(a.inf<=b && a.sup>=b); 00206 } 00207 inline bool operator >(const interval &a,const real &b) throw() 00208 { 00209 return(a.inf<b && a.sup>b); 00210 } 00211 00212 inline bool operator !(const interval &a) throw() { return (a.inf <= 0.0 && a.sup >= 0.0); } 00213 00214 inline real & Inf (interval& a) throw() { return a.inf; } 00215 inline const real & Inf (const interval &a) throw() { return a.inf; } 00216 inline real & Sup (interval& a) throw() { return a.sup; } 00217 inline const real & Sup (const interval &a) throw() { return a.sup; } 00218 00219 inline interval& SetInf (interval& a, const real& b) throw() {a.inf=b; return a;} 00220 inline interval& SetSup (interval& a, const real& b) throw() {a.sup=b; return a;} 00221 inline interval& UncheckedSetInf (interval& a, const real& b) throw() { a.inf=b; return a;} 00222 inline interval& UncheckedSetSup (interval& a, const real& b) throw() { a.sup=b; return a;} 00223 00224 inline bool IsEmpty(const interval &a) throw() { return (a.inf>a.sup); } 00225 00226 inline interval abs(const interval &a) throw() 00227 { 00228 real h1 = abs(a.inf); 00229 real h2 = abs(a.sup); 00230 00231 if (IsEmpty(a)) return a; 00232 if (!a) 00233 return interval(0.0, (h1 > h2) ? h1 : h2); 00234 if (h1 > h2) 00235 return interval(h2, h1); 00236 00237 return interval(h1, h2); 00238 } 00239 00240 inline real diam(const interval & a) throw() 00241 { 00242 return subup(a.sup,a.inf); 00243 } 00244 00245 inline real Mid(const interval & a) throw() 00246 { 00247 return addd( a.inf, subd(0.5*a.sup,0.5*a.inf) ); 00248 } 00249 00250 inline void times2pown(interval& x, const int& n) throw() 00251 { 00252 real r1,r2; 00253 int j; 00254 r1 = x.inf; r2 = x.sup; 00255 j = expo(r1) + n; 00256 if (j >= -1021) r1 = comp(mant(r1),j); 00257 else 00258 { 00259 j += 1021; 00260 r1 = comp(mant(r1), -1021); 00261 if (j<-53) 00262 { 00263 if (sign(r1)>=0) r1 = 0; 00264 else r1 = -minreal; 00265 } else 00266 r1 = muld(r1,comp(0.5,j+1)); 00267 } 00268 j = expo(r2) + n; 00269 if (j >= -1021) r2 = comp(mant(r2),j); 00270 else 00271 { 00272 j += 1021; 00273 r2 = comp(mant(r2), -1021); 00274 if (j<-53) 00275 { 00276 if (sign(r2)>0) r2 = minreal; 00277 else r2 = 0; 00278 } else r2 = mulu(r2,comp(0.5,j+1)); 00279 } 00280 x = _interval(r1,r2); 00281 } // times2pown(...) 00282 00283 interval operator+ (const interval& a, const interval& b) throw() 00284 { return interval (adddown(a.inf,b.inf),addup(a.sup,b.sup)); } 00285 interval operator+ (const interval& a, const real& b) throw() 00286 { return interval (adddown(a.inf,b),addup(a.sup,b)); } 00287 interval operator+ (const real& a, const interval& b) throw() 00288 { return interval (adddown(a,b.inf),addup(a,b.sup)); } 00289 00290 interval operator- (const interval& a, const interval& b) throw() 00291 { return interval ( subdown(a.inf,b.sup), subup(a.sup,b.inf)); } 00292 interval operator- (const interval& a, const real& b) throw() 00293 { return interval ( subdown(a.inf,b), subup(a.sup,b)); } 00294 interval operator- (const real& a, const interval& b) throw() 00295 { return interval ( subdown(a,b.sup), subup(a,b.inf)); } 00296 00297 //------------------------------------------------------------------------ 00298 // a * b = ueber Entscheidungstabelle : 00299 // 00300 // bi,bs >= 0 bi < 0, bs >=0 bi,bs < 0 00301 // ----------------+------------------+------------------+---------------- 00302 // ai,as >= 0 I ai*bi, as*bs I as*bi, as*bs I as*bi, ai*bs 00303 // ----------------+------------------+------------------+---------------- 00304 // ai < 0, as >= 0 I ai*bs, as*bs I .... I as*bi, ai*bi 00305 // ----------------+------------------+------------------+---------------- 00306 // ai,as < 0 I ai*bs, as*bi I ai*bs, ai*bi I as*bs, ai*bi 00307 // ----------------+------------------+------------------+---------------- 00308 // 00309 // .... : min(ai*bs, as*bi), max(ai*ai, as*as) 00310 // 00311 interval operator *(const interval& a, const interval& b) throw() 00312 { 00313 interval tmp; 00314 00315 if (sign(a.inf) >= 0) 00316 { // 1. Zeile der Entscheidungstabelle 00317 if (sign(b.inf) >= 0) 00318 { // 1. Spalte: [ai*bi, as*bs] 00319 tmp.inf = multdown(a.inf,b.inf); 00320 tmp.sup = multup(a.sup,b.sup); 00321 } else if (sign(b.sup) >= 0) 00322 { // 2. Spalte: [as*bi, as*bs] 00323 tmp.inf = multdown(a.sup,b.inf); 00324 tmp.sup = multup(a.sup,b.sup); 00325 } else 00326 { // 3. Spalte: [as*bi, ai*bs] 00327 tmp.inf = multdown(a.sup,b.inf); 00328 tmp.sup = multup(a.inf,b.sup); 00329 } 00330 } else if (sign(a.sup) >= 0) 00331 { // 2. Zeile der Entscheidungstabelle 00332 if (sign(b.inf) >= 0) 00333 { // 1. Spalte: [ai*bs, as*bs] 00334 tmp.inf = multdown(a.inf,b.sup); 00335 tmp.sup = multup(a.sup,b.sup); 00336 } else if (sign(b.sup) >= 0) 00337 { // 2. Spalte: [min(ai*bs, as*bi), 00338 real hlp; // max(ai*ai, as*as)] 00339 00340 tmp.inf = multdown(a.inf,b.sup); 00341 hlp = multdown(a.sup,b.inf); 00342 00343 if (hlp < tmp.inf) 00344 tmp.inf = hlp; 00345 00346 tmp.sup = multup(a.inf,b.inf); 00347 hlp = multup(a.sup,b.sup); 00348 00349 if (hlp > tmp.sup) 00350 tmp.sup = hlp; 00351 } else 00352 { // 3. Spalte: [as*bi, ai*bi] 00353 tmp.inf = multdown(a.sup,b.inf); 00354 tmp.sup = multup(a.inf,b.inf); 00355 } 00356 } else 00357 { // 3. Zeile der Entscheidungstabelle 00358 if (sign(b.inf) >= 0) 00359 { // 1. Spalte: [ai*bs, as*bi] 00360 tmp.inf = multdown(a.inf,b.sup); 00361 tmp.sup = multup(a.sup,b.inf); 00362 } else if (sign(b.sup) >= 0) 00363 { // 2. Spalte: [ai*bs, ai*bi] 00364 tmp.inf = multdown(a.inf,b.sup); 00365 tmp.sup = multup(a.inf,b.inf); 00366 } else 00367 { // 3. Spalte: [as*bs, ai*bi] 00368 tmp.inf = multdown(a.sup,b.sup); 00369 tmp.sup = multup(a.inf,b.inf); 00370 } 00371 } 00372 return tmp; 00373 } 00374 00375 //------------------------------------------------------------------------ 00376 // a / b = ueber Entscheidungstabelle : 00377 // 00378 // bi,bs > 0 bi,bs < 0 00379 // ----------------+------------------+------------------ 00380 // ai,as >= 0 I ai/bs, as/bi I as/bs, ai/bi 00381 // ----------------+------------------+------------------ 00382 // ai < 0, as >= 0 I ai/bi, as/bi I as/bs, ai/bs 00383 // ----------------+------------------+------------------ 00384 // ai,as < 0 I ai/bi, as/bs I as/bi, ai/bs 00385 // ----------------+------------------+------------------ 00386 // 00387 interval operator/ (const interval& a, const interval& b) throw(DIV_BY_ZERO) 00388 { 00389 interval tmp; 00390 00391 if ((sign(b.inf) <= 0) && (sign(b.sup) >= 0)) 00392 cxscthrow(DIV_BY_ZERO("interval::interval operator/(const interval&,const interval&)")); 00393 00394 if (sign(a.inf) >= 0) 00395 { // 1. Zeile der Entscheidungstabelle 00396 if (sign(b.inf) > 0) 00397 { // 1. Spalte: [ai/bs, as/bi] 00398 tmp.inf = divdown(a.inf,b.sup); 00399 tmp.sup = divup(a.sup,b.inf); 00400 } else 00401 { // 2. Spalte: [as/bs, ai/bi] 00402 tmp.inf = divdown(a.sup,b.sup); 00403 tmp.sup = divup(a.inf,b.inf); 00404 } 00405 } else if (sign(a.sup) >= 0) 00406 { // 2. Zeile der Entscheidungstabelle 00407 if (sign(b.inf) > 0) 00408 { // 1. Spalte: [ai/bi, as/bi] 00409 tmp.inf = divdown(a.inf,b.inf); 00410 tmp.sup = divup(a.sup,b.inf); 00411 } else 00412 { // 2. Spalte: [as/bs, ai/bs] 00413 tmp.inf = divdown(a.sup,b.sup); 00414 tmp.sup = divup(a.inf,b.sup); 00415 } 00416 } else 00417 { // 3. Zeile der Entscheidungstabelle 00418 if (sign(b.inf) > 0) 00419 { // 1. Spalte: [ai/bi, as/bs] 00420 tmp.inf = divdown(a.inf,b.inf); 00421 tmp.sup = divup(a.sup,b.sup); 00422 } else 00423 { // 2. Spalte: [as/bi, ai/bs] 00424 tmp.inf = divdown(a.sup,b.inf); 00425 tmp.sup = divup(a.inf,b.sup); 00426 } 00427 } 00428 00429 return tmp; 00430 } 00431 00432 interval operator* (const real& a, const interval& b) throw() 00433 { 00434 interval tmp; 00435 if (sign(a) == 0) 00436 { 00437 tmp.inf = 0.0; 00438 tmp.sup = 0.0; 00439 } else if (sign(a) > 0) 00440 { 00441 tmp.inf = multdown(a,b.inf); 00442 tmp.sup = multup(a,b.sup); 00443 } else // if (sign(a) < 0) 00444 { 00445 tmp.inf = multdown(a,b.sup); 00446 tmp.sup = multup(a,b.inf); 00447 } 00448 return tmp; 00449 } 00450 00451 interval operator* (const interval& a, const real& b) throw() 00452 { 00453 interval tmp; 00454 if (sign(b) == 0) 00455 { 00456 tmp.inf = 0.0; 00457 tmp.sup = 0.0; 00458 } else if (sign(b) > 0) 00459 { 00460 tmp.inf = multdown(a.inf,b); 00461 tmp.sup = multup(a.sup,b); 00462 } else // if (sign(b) < 0) 00463 { 00464 tmp.inf = multdown(a.sup,b); 00465 tmp.sup = multup(a.inf,b); 00466 } 00467 return tmp; 00468 } 00469 00470 interval operator/ (const real& a, const interval& b) throw() { return (interval(a) / b); } 00471 interval operator/ (const interval& a, const real& b) throw() { return (a / interval(b)); } 00472 00473 00474 } // namespace cxsc 00475