C-XSC - A C++ Class Library for Extended Scientific Computing: interval.inl Source File

00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2011 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.27 2011/06/07 15:17:39 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