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