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.cpp,v 1.38 2014/01/30 17:23:45 cxsc Exp $ */ 00025 00026 #include "interval.hpp" 00027 #include "ioflags.hpp" 00028 #include "idot.hpp" 00029 #include "dot.hpp" 00030 #include "RtsFunc.h" // b_shr1 in mid() 00031 00032 // Auch fi_lib hat ein eigenes interval (wie RTS sein dotprecision) 00033 //typedef struct fi_interval { double INF, SUP;} fi_interval; 00034 // class interval; 00035 #undef LINT_ARGS 00036 #define CXSC_INCLUDE 00037 #include <fi_lib.hpp> 00038 00039 namespace cxsc { 00040 00041 interval::interval(const idotprecision & a) throw() 00042 { 00043 *this=rnd(a); 00044 } 00045 00046 interval & interval::operator =(const idotprecision & a) throw() 00047 { 00048 return *this=rnd(a); 00049 } 00050 00051 interval::interval(const dotprecision & a) throw() 00052 { 00053 rnd(a,inf,sup); 00054 } 00055 00056 interval::interval(const dotprecision &a,const dotprecision &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00057 { 00058 inf=rnd(a,RND_DOWN); 00059 sup=rnd(b,RND_UP); 00060 if(inf>sup) 00061 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("interval::interval(const dotprecision &,const dotprecision &)")); 00062 } 00063 00064 interval::interval(const l_real &a,const l_real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL) 00065 { 00066 dotprecision dot(a); 00067 inf=rnd(dot,RND_DOWN); 00068 dot = b; 00069 sup = rnd(dot,RND_UP); 00070 if(inf>sup) 00071 cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("interval::interval(const l_real &,const l_real &)")); 00072 } 00073 00074 00075 interval & interval::operator =(const dotprecision & a) throw() 00076 { 00077 rnd(a,inf,sup); 00078 return *this; 00079 } 00080 00081 00082 // ---- Standardfunkt ---- (arithmetische Operatoren) 00083 /* 00084 interval operator+(const interval &a, const interval &b) throw() { return interval(adddown(a.inf,b.inf),addup(a.sup,b.sup)); } 00085 interval operator+(const interval &a, const real &b) throw() { return interval(adddown(a.inf,b),addup(a.sup,b)); } 00086 interval operator+(const real &a, const interval &b) throw() { return interval(adddown(a,b.inf),addup(a,b.sup)); } 00087 00088 interval operator-(const interval &a, const interval &b) throw() { return interval(subdown(a.inf,b.inf),subup(a.sup,b.sup)); } 00089 interval operator-(const interval &a, const real &b) throw() { return interval(subdown(a.inf,b),subup(a.sup,b)); } 00090 interval operator-(const real &a, const interval &b) throw() { return interval(subdown(a,b.sup),subup(a,b.inf)); } 00091 00092 interval operator*(const interval &a, const interval &b) throw() { return mul_ii (a, b); } 00093 interval operator*(const interval &a, const real &b) throw() { return mul_id (a, *(double*)&b); } 00094 interval operator*(const real &a, const interval &b) throw() { return mul_di (*(double*)&a, b); } 00095 00096 interval operator/(const interval &a, const interval &b) throw() { return div_ii (a, b); } 00097 interval operator/(const interval &a, const real &b) throw() { return div_id (a, *(double*)&b); } 00098 interval operator/(const real &a, const interval &b) throw() { return div_di (*(double*)&a, b); } 00099 */ 00100 00101 00102 bool operator ==(const interval &a,const dotprecision &r) throw() { return(r==a.inf && r==a.sup); } 00103 bool operator !=(const interval &a,const dotprecision &r) throw() { return(r!=a.inf || r!=a.sup); } 00104 bool operator ==(const dotprecision &r,const interval &a) throw() { return(r==a.inf && r==a.sup); } 00105 bool operator !=(const dotprecision &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); } 00106 00107 bool operator <=(const dotprecision &a,const interval &b) throw() 00108 { 00109 return(a>=b.inf && a<=b.sup); 00110 } 00111 bool operator >=(const dotprecision &a,const interval &b) throw() 00112 { 00113 return(a<=b.inf && a>=b.sup); 00114 } 00115 bool operator <(const dotprecision &a,const interval &b) throw() 00116 { 00117 return(a>b.inf && a<b.sup); 00118 } 00119 00120 bool operator <=(const interval &a,const dotprecision &b) throw() 00121 { 00122 return(a.inf>=b && a.sup<=b); 00123 } 00124 bool operator >=(const interval &a,const dotprecision &b) throw() 00125 { 00126 return(a.inf<=b && a.sup>=b); 00127 } 00128 bool operator >(const interval &a,const dotprecision &b) throw() 00129 { 00130 return(a.inf<b && a.sup>b); 00131 } 00132 00133 00134 00135 // ---- Ausgabefunkt. --------------------------------------- 00136 00137 std::ostream & operator << (std::ostream &s, const interval& a) throw() 00138 { 00139 s << '[' << SaveOpt << RndDown 00140 << a.inf << ',' << RndUp 00141 << a.sup << RestoreOpt 00142 << ']'; 00143 return s; 00144 } 00145 std::string & operator << (std::string &s, const interval &a) throw() 00146 { 00147 s+='['; 00148 s << SaveOpt << RndDown 00149 << a.inf; 00150 s+=','; 00151 s << RndUp 00152 << a.sup << RestoreOpt; 00153 s+=']'; 00154 return s; 00155 } 00156 00157 std::istream & operator >> (std::istream &s, interval &a) throw() 00158 { 00159 char c; 00160 00161 skipeolnflag = inpdotflag = true; 00162 c = skipwhitespacessinglechar (s, '['); 00163 if (inpdotflag) 00164 s.putback(c); 00165 00166 s >> SaveOpt >> RndDown >> a.inf; 00167 00168 skipeolnflag = inpdotflag = true; 00169 c = skipwhitespacessinglechar (s, ','); 00170 if (inpdotflag) s.putback(c); 00171 00172 s >> RndUp >> a.sup >> RestoreOpt; 00173 00174 if (!waseolnflag) 00175 { 00176 skipeolnflag = false, inpdotflag = true; 00177 c = skipwhitespaces (s); 00178 if (inpdotflag && c != ']') 00179 s.putback(c); 00180 } 00181 00182 /*if (a.inf > a.sup) { 00183 errmon (ERR_INTERVAL(EMPTY)); 00184 } */ 00185 return s; 00186 } 00187 00188 std::string & operator >> (std::string &s, interval &a) throw() 00189 { 00190 s = skipwhitespacessinglechar (s, '['); 00191 s >> SaveOpt >> RndDown >> a.inf; 00192 s = skipwhitespacessinglechar (s, ','); 00193 s >> RndUp >> a.sup >> RestoreOpt; 00194 s = skipwhitespaces (s); 00195 00196 if (s[0] == ']') 00197 s.erase(0,1); 00198 00199 /*if (a.inf > a.sup) { 00200 errmon (ERR_INTERVAL(EMPTY)); 00201 } */ 00202 return s; 00203 } 00204 00205 void operator >>(const std::string &s,interval &a) throw() 00206 { 00207 std::string r(s); 00208 r>>a; 00209 } 00210 void operator >>(const char *s,interval &a) throw() 00211 { 00212 std::string r(s); 00213 r>>a; 00214 } 00215 00216 real mid(const interval & a) throw() 00217 { 00218 dotprecision dot(a.inf); 00219 dot += a.sup; 00220 00221 if (dot != 0.0) 00222 { 00223 // Division nur bei ungleich 0.0 00224 Dotprecision ptr = *dot.ptr(); 00225 00226 // -------------------------------------------------------- 00227 // Dividiere dot durch 2, mittels 1 Rechtsshift 00228 00229 ptr[(a_intg)++ptr[A_END]] = ZERO; 00230 00231 b_shr1 (&ptr[(a_intg)ptr[A_BEGIN]], (a_intg)(ptr[A_END]-ptr[A_BEGIN]+1)); 00232 00233 if (ptr[(a_intg)ptr[A_END]] == ZERO) 00234 --ptr[A_END]; 00235 00236 if (ptr[(a_intg)ptr[A_BEGIN]] == ZERO) 00237 ++ptr[A_BEGIN]; 00238 00239 // -------------------------------------------------------- 00240 } 00241 00242 return rnd(dot); 00243 } 00244 00245 // for compatibility with CToolbox library - from former i_util.hpp 00246 int in (const real& x, const interval& y ) // Contained-in relation 00247 { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //---------------------- 00248 00249 int in ( const interval& x, const interval& y ) // Contained-in relation 00250 { //---------------------- 00251 return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) ); 00252 } 00253 00254 // function in dot.hpp and dot.cpp 00255 //void rnd ( const dotprecision& d, interval& x ) // Round dotprecision to interval 00256 //{ //------------------------------- 00257 // real Lower, Upper; 00258 // 00259 // rnd(d,Lower,Upper); // Rounding to interval bounds 00260 // x = interval(Lower,Upper); 00261 //} 00262 00281 interval Blow ( const interval& x, const real& eps ) // Epsilon inflation 00282 { //------------------ 00283 interval tmp; 00284 tmp = (1.0+eps)*x - eps*x; 00285 return interval(pred(Inf(tmp)),succ(Sup(tmp))); 00286 } 00287 00288 int Disjoint ( const interval& a, const interval& b ) // Test for disjointedness 00289 { //------------------------ 00290 return (Inf(a) > Sup(b)) || (Inf(b) > Sup(a)); 00291 } 00292 00293 real AbsMin ( const interval& x ) // Absolute minimum of 00294 { // an interval 00295 if ( in(0.0,x) ) //-------------------- 00296 return 0.0; 00297 else if (Inf(x) > 0.0) 00298 return Inf(x); 00299 else 00300 return -Sup(x); 00301 } 00302 00303 real AbsMax ( const interval& x ) // Absolute maximum of 00304 { // an interval 00305 real a, b; //-------------------- 00306 00307 a = abs(Inf(x)); 00308 b = abs(Sup(x)); 00309 00310 if (a > b) 00311 return a; 00312 else 00313 return b; 00314 } 00315 00316 real RelDiam ( const interval& x ) // Relative diameter 00317 { // of an interval 00318 if ( in(0.0,x) ) //------------------ 00319 return diam(x); 00320 else 00321 return( divu(diam(x),AbsMin(x)) ); 00322 } 00323 00324 //---------------------------------------------------------------------------- 00325 // Checks if the diameter of the interval 'x' is less or equal to 'n' ulps. 00326 // Ulp is an abbreviation for units in the last place. 00327 //---------------------------------------------------------------------------- 00335 int UlpAcc ( const interval& x, int n ) 00336 { 00337 int i; 00338 real Infimum; 00339 00340 Infimum = Inf(x); 00341 for (i = 1; i <= n; i++) Infimum = succ(Infimum); 00342 return (Infimum >= Sup(x)); 00343 } 00344 00345 // The following interval constants are optimal inclusions: 00346 00347 const real Pi_Inf = 7074237752028440.0 / 2251799813685248.0; 00348 const interval Pi_interval = interval(Pi_Inf,succ(Pi_Inf)); // Pi 00349 00350 const real Pi2_Inf = 7074237752028440.0 / 1125899906842624.0; 00351 const interval Pi2_interval = interval(Pi2_Inf,succ(Pi2_Inf)); // 2*Pi 00352 00353 const real Pi3_Inf = 5305678314021330.0 / 562949953421312.0; 00354 const interval Pi3_interval = interval(Pi3_Inf,succ(Pi3_Inf)); // 3*Pi 00355 00356 const real Pid2_Inf = 7074237752028440.0 / 4503599627370496.0; 00357 const interval Pid2_interval = interval(Pid2_Inf,succ(Pid2_Inf)); // Pi/2 00358 00359 const real Pid3_Inf = 4716158501352293.0 / 4503599627370496.0; 00360 const interval Pid3_interval = interval(Pid3_Inf,succ(Pid3_Inf)); // Pi/3 00361 00362 const real Pid4_Inf = 7074237752028440.0 / 9007199254740992.0; 00363 const interval Pid4_interval = interval(Pid4_Inf,succ(Pid4_Inf)); // Pi/4 00364 00365 const real Pir_Inf = 5734161139222658.0 / 18014398509481984.0; 00366 const interval Pir_interval = interval(Pir_Inf,succ(Pir_Inf)); // 1/Pi 00367 00368 const real Pi2r_Inf = 5734161139222658.0 / 36028797018963968.0; 00369 const interval Pi2r_interval = interval(Pi2r_Inf,succ(Pi2r_Inf)); 00370 // 1/(2*Pi) 00371 00372 const real Pip2_Inf = 5556093337880030.0 / 562949953421312.0; 00373 const interval Pip2_interval = interval(Pip2_Inf,succ(Pip2_Inf)); // Pi^2 00374 00375 const real SqrtPi_Inf = 7982422502469482.0 / 4503599627370496.0; 00376 const interval SqrtPi_interval = interval(SqrtPi_Inf,succ(SqrtPi_Inf)); 00377 // sqrt(Pi) 00378 00379 const real Sqrt2Pi_Inf = 5644425081792261.0 / 2251799813685248.0; 00380 const interval Sqrt2Pi_interval = interval(Sqrt2Pi_Inf,succ(Sqrt2Pi_Inf)); 00381 // sqrt(2Pi) 00382 00383 const real SqrtPir_Inf = 5081767996463981.0 / 9007199254740992.0; 00384 const interval SqrtPir_interval = interval(SqrtPir_Inf,succ(SqrtPir_Inf)); 00385 // 1/sqrt(Pi) 00386 00387 const real Sqrt2Pir_Inf = 7186705221432912.0 / 18014398509481984.0; 00388 const interval Sqrt2Pir_interval=interval(Sqrt2Pir_Inf,succ(Sqrt2Pir_Inf)); 00389 // 1/sqrt(2Pi) 00390 00391 const real Sqrt2_Inf = 6369051672525772.0 / 4503599627370496.0; 00392 const interval Sqrt2_interval=interval(Sqrt2_Inf,succ(Sqrt2_Inf)); 00393 // sqrt(2) 00394 00395 const real Sqrt5_Inf = 5035177455121575.0 / 2251799813685248.0; 00396 const interval Sqrt5_interval=interval(Sqrt5_Inf,succ(Sqrt5_Inf)); 00397 // sqrt(5) 00398 00399 const real Sqrt7_Inf = 5957702309312745.0 / 2251799813685248.0; 00400 const interval Sqrt7_interval=interval(Sqrt7_Inf,succ(Sqrt7_Inf)); 00401 // sqrt(7) 00402 00403 const real Sqrt2r_Inf = 6369051672525772.0 / 9007199254740992.0; 00404 const interval Sqrt2r_interval=interval(Sqrt2r_Inf,succ(Sqrt2r_Inf)); 00405 // 1/sqrt(2) 00406 00407 const real Sqrt3_Inf = 7800463371553962.0 / 4503599627370496.0; 00408 const interval Sqrt3_interval=interval(Sqrt3_Inf,succ(Sqrt3_Inf)); 00409 // sqrt(3) 00410 00411 const real Sqrt3d2_Inf = 7800463371553962.0 / 9007199254740992.0; 00412 const interval Sqrt3d2_interval=interval(Sqrt3d2_Inf,succ(Sqrt3d2_Inf)); 00413 // sqrt(3)/2 00414 00415 const real Sqrt3r_Inf = 5200308914369308.0 / 9007199254740992.0; 00416 const interval Sqrt3r_interval=interval(Sqrt3r_Inf,succ(Sqrt3r_Inf)); 00417 // 1/sqrt(3) 00418 00419 const real Ln2_Inf = 6243314768165359.0 / 9007199254740992.0; 00420 const interval Ln2_interval=interval(Ln2_Inf,succ(Ln2_Inf)); // ln(2) 00421 00422 const real Ln2r_Inf = 6497320848556798.0 / 4503599627370496.0; 00423 const interval Ln2r_interval=interval(Ln2r_Inf,succ(Ln2r_Inf)); // 1/ln(2) 00424 00425 const real Ln10_Inf = 5184960683398421.0 / 2251799813685248.0; 00426 const interval Ln10_interval=interval(Ln10_Inf,succ(Ln10_Inf)); // ln(10) 00427 00428 const real Ln10r_Inf = 7823553867474190.0 / 18014398509481984.0; 00429 const interval Ln10r_interval=interval(Ln10r_Inf,succ(Ln10r_Inf)); 00430 // 1/ln(10) 00431 00432 const real LnPi_Inf = 5155405087351229.0 / 4503599627370496.0; 00433 const interval LnPi_interval=interval(LnPi_Inf,succ(LnPi_Inf)); // ln(Pi) 00434 00435 const real Ln2Pi_Inf = 8277062471433908.0 / 4503599627370496.0; 00436 const interval Ln2Pi_interval=interval(Ln2Pi_Inf,succ(Ln2Pi_Inf)); 00437 // ln(2Pi) 00438 00439 const real E_Inf = 6121026514868073.0 / 2251799813685248.0; 00440 const interval E_interval=interval(E_Inf,succ(E_Inf)); // e 00441 00442 const real Er_Inf = 6627126856707895.0 / 18014398509481984.0; 00443 const interval Er_interval=interval(Er_Inf,succ(Er_Inf)); // 1/e 00444 00445 const real Ep2_Inf = 8319337573440941.0 / 1125899906842624.0; 00446 const interval Ep2_interval=interval(Ep2_Inf,succ(Ep2_Inf)); // e^2 00447 00448 const real Ep2r_Inf = 4875967449235915.0 / 36028797018963968.0; 00449 const interval Ep2r_interval=interval(Ep2r_Inf,succ(Ep2r_Inf)); // 1/e^2 00450 00451 const real EpPi_Inf = 6513525919879993.0 / 281474976710656.0; 00452 const interval EpPi_interval=interval(EpPi_Inf,succ(EpPi_Inf)); // e^(Pi) 00453 00454 const real Ep2Pi_Inf = 4710234414611993.0 / 8796093022208.0; 00455 const interval Ep2Pi_interval=interval(Ep2Pi_Inf,succ(Ep2Pi_Inf)); 00456 // e^(2Pi) 00457 00458 const real EpPid2_Inf = 5416116035097439.0 / 1125899906842624.0; 00459 const interval EpPid2_interval=interval(EpPid2_Inf,succ(EpPid2_Inf)); 00460 // e^(Pi/2) 00461 00462 const real EpPid4_Inf = 4938827609611434.0 / 2251799813685248.0; 00463 const interval EpPid4_interval=interval(EpPid4_Inf,succ(EpPid4_Inf)); 00464 // e^(Pi/4) 00465 00466 00467 } // namespace cxsc 00468