C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_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: l_interval.cpp,v 1.21 2014/01/30 17:23:46 cxsc Exp $ */
00025 
00026 #include <math.h>
00027 #include "l_interval.hpp"
00028 #include "idot.hpp"
00029 #include "RtsFunc.h"
00030 
00031 namespace cxsc {
00032 
00033 #define CXSC_Zero 0.
00034 
00035 l_interval & l_interval::operator= (const l_interval & a)    throw()
00036 {
00037    real *newdata=new real[a.prec+1];
00038    int i;
00039    for(i=0;i<=a.prec;i++)
00040       newdata[i]=a.data[i];
00041    delete [] data;
00042    data=newdata;
00043    prec=a.prec;
00044    return *this;
00045 }   
00046 
00047 // ---- Typwandlungen ----
00048 
00049 interval::interval(const l_interval & a) throw()
00050 {
00051    idotprecision idot(0.0);
00052    a._akku_add(idot);
00053    *this=rnd(idot);
00054 }
00055 
00056 interval & interval::operator =(const l_interval & a) throw()
00057 {
00058    idotprecision idot(0.0);
00059    a._akku_add(idot);
00060    return *this=rnd(idot);
00061 }
00062 
00063 
00064 interval _interval(const l_interval & a) throw()
00065 {
00066    idotprecision idot(0.0);
00067    a._akku_add(idot);
00068    return rnd(idot);
00069 }
00070 
00071 l_interval::l_interval(const dotprecision & a)
00072 #if (CXSC_INDEX_CHECK)
00073    throw(ERROR_LINTERVAL_WRONG_STAGPREC)
00074 #else
00075    throw()
00076 #endif
00077 {
00078    _allo(stagprec);
00079    idotprecision idot(a);
00080    _akku_out(idot);   
00081 }
00082 
00083 l_interval::l_interval(const dotprecision & a,const dotprecision & b)
00084 #if (CXSC_INDEX_CHECK)
00085    throw(ERROR_LINTERVAL_WRONG_STAGPREC,ERROR_LINTERVAL_EMPTY_INTERVAL)
00086 #else
00087    throw(ERROR_LINTERVAL_EMPTY_INTERVAL)
00088 #endif
00089 {
00090    if(a>b)
00091       cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const dotprecision&,const dotprecision&)"));
00092    _allo(stagprec);
00093    idotprecision idot;
00094    UncheckedSetInf(idot,a);
00095    UncheckedSetSup(idot,b);
00096    _akku_out(idot);   
00097 }
00098 
00099 l_interval & l_interval::operator =(const dotprecision & a)
00100 #if (CXSC_INDEX_CHECK)
00101    throw(ERROR_LINTERVAL_WRONG_STAGPREC)
00102 #else
00103    throw()
00104 #endif
00105 {
00106    if(stagprec!=prec)
00107    {
00108       delete [] data;
00109       _allo(stagprec);
00110    }
00111    idotprecision idot(a);
00112    _akku_out(idot);   
00113    return *this;
00114 }
00115 
00116 l_interval::l_interval(const idotprecision & a)
00117 #if (CXSC_INDEX_CHECK)
00118    throw(ERROR_LINTERVAL_WRONG_STAGPREC)
00119 #else
00120    throw()
00121 #endif
00122 {
00123    _allo(stagprec);
00124    idotprecision idot(a);
00125    _akku_out(idot);   
00126 }
00127 
00128 l_interval & l_interval::operator =(const idotprecision & a)
00129 #if (CXSC_INDEX_CHECK)
00130    throw(ERROR_LINTERVAL_WRONG_STAGPREC)
00131 #else
00132    throw()
00133 #endif
00134 {
00135    if(stagprec!=prec)
00136    {
00137       delete [] data;
00138       _allo(stagprec);
00139    }
00140    idotprecision idot(a);
00141    _akku_out(idot);   
00142    return *this;
00143 }
00144 
00145 l_interval::l_interval(const l_real &a, const l_real &b) 
00146 #if (CXSC_INDEX_CHECK)
00147    throw(ERROR_LINTERVAL_WRONG_STAGPREC,ERROR_LINTERVAL_EMPTY_INTERVAL)
00148 #else
00149    throw(ERROR_LINTERVAL_EMPTY_INTERVAL)
00150 #endif
00151 {
00152    _allo(stagprec);
00153    if(a>b)
00154       cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const l_real &a, const l_real &b)"));
00155    dotprecision dot1, dot2;
00156    dot1 = a;
00157    dot2 = b;
00158    idotprecision idot(dot1,dot2);
00159    _akku_out(idot);
00160 }
00161 
00162 l_interval::l_interval(const real &a, const l_real &b) 
00163 #if (CXSC_INDEX_CHECK)
00164    throw(ERROR_LINTERVAL_WRONG_STAGPREC,ERROR_LINTERVAL_EMPTY_INTERVAL)
00165 #else
00166    throw(ERROR_LINTERVAL_EMPTY_INTERVAL)
00167 #endif
00168 {
00169    _allo(stagprec);
00170    if(a>b)
00171       cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const l_real &a, const l_real &b)"));
00172    dotprecision dot1, dot2;
00173    dot1 = a;
00174    dot2 = b;
00175    idotprecision idot(dot1,dot2);
00176    _akku_out(idot);
00177 }
00178 
00179 void l_realz2l_interval(const l_real& lr, const interval& z,
00180                                                 l_interval& li) throw()
00181 {  // converts lr+z to li of type l_interval; Blomquist, 12.10.06;
00182    // lr+z is included by li.
00183    // prec(lr) <= prec(li)  must be realized!
00184     int p = StagPrec(lr);  int q = StagPrec(li);
00185     if(p>q) 
00186     { 
00187         std::cerr << "l_realz2l_interval(const l_real& lr, const interval& z, l_interval& li): incorrect precisions of lr,li !" 
00188                   << std:: endl; 
00189         exit (1); 
00190     }
00191 
00192     for (int i=1; i<=q-1; i++) li[i] = 0; 
00193     li[q]   = Inf(z);
00194     li[q+1] = Sup(z);
00195     if(p<q) for (int i=1; i<=p; i++) li[i] = lr[i];
00196     else // p=q
00197     {
00198         p--;
00199         for (int i=1; i<=p; i++) li[i] = lr[i];
00200         li[q]   = addd(lr[p+1],Inf(z)); 
00201         li[q+1] = addu(lr[p+1],Sup(z));
00202     }
00203 } // l_realz2l_interval
00204 
00205 l_interval::l_interval(const l_real &a, const real &b) 
00206 #if (CXSC_INDEX_CHECK)
00207    throw(ERROR_LINTERVAL_WRONG_STAGPREC,ERROR_LINTERVAL_EMPTY_INTERVAL)
00208 #else
00209    throw(ERROR_LINTERVAL_EMPTY_INTERVAL)
00210 #endif
00211 {
00212    _allo(stagprec);
00213    if(a>b)
00214       cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const l_real &a, const l_real &b)"));
00215    dotprecision dot1, dot2;
00216    dot1 = a;
00217    dot2 = b;
00218    idotprecision idot(dot1,dot2);
00219    _akku_out(idot);
00220 }
00221 
00222 
00223 /*
00224 l_interval _l_interval(const l_real & a, const l_real & b) throw(ERROR_LINTERVAL_EMPTY_INTERVAL)
00225 {
00226    if(a>b)
00227       cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval _l_interval(const l_real & a, const l_real & b)"));
00228 
00229    dotakku[0]=0.0;
00230    dotakku[1]=0.0;
00231    a._akku_add(dotakku[0]);
00232    b._akku_add(dotakku[1]);
00233    idotakku[0]=_idotprecision(dotakku[0],dotakku[1]);
00234    l_interval tmp;
00235    tmp._akku_out();
00236    return tmp;
00237 }   
00238 */
00239 
00240 /*
00241 l_interval _unchecked_l_interval(const l_real & lr1, const l_real & lr2) throw()
00242 {
00243    real inf, sup, tmp;          // fuer Naeherungen, entspricht Interval z
00244    int i=1;
00245    l_interval li;
00246    dotprecision dot1(0.0);
00247    dotprecision dot2(0.0);
00248    lr1._akku_add(dot1);
00249    lr2._akku_add(dot2);
00250    inf = rnd(dot1,RND_DOWN);
00251    sup = rnd(dot2,RND_UP);
00252    while (i<li.prec && !(inf<=0. && sup>=0.)) 
00253    {
00254       tmp = inf + (sup-inf)/2.0; // middle(interval)
00255       li.elem(i) = tmp;
00256       dot1 -= tmp;
00257       dot2 -= tmp;
00258       inf = rnd(dot1,RND_DOWN);
00259       sup = rnd(dot2,RND_UP);
00260       i++;
00261    }
00262    li.elem(li.prec)=inf;           // Intervall in die letzten Stellen
00263    li.elem(li.prec+1)=sup;         // schreiben
00264 
00265    return li;
00266 }
00267 */
00268 
00269 l_interval _unchecked_l_interval(const l_real & lr1, const l_real & lr2) throw()
00270 {
00271    dotprecision dot1, dot2;
00272    dot1 = lr1;
00273    dot2 = lr2;
00274    idotprecision idot(dot1,dot2);
00275    l_interval li;
00276    li._akku_out(idot);  
00277    return li;
00278 }
00279 
00280 
00281 idotprecision _idotprecision(const l_interval & a) throw()
00282 {
00283    return idotprecision(a);
00284 }
00285 
00286 idotprecision::idotprecision(const l_interval & a) throw() : inf(0),
00287                                                              sup(0)
00288 { 
00289    a._akku_add(*this);
00290 }
00291 
00292 idotprecision & idotprecision::operator =(const l_interval & a) throw()
00293 {
00294    *this=0;
00295    a._akku_add(*this);
00296    return *this;
00297 }
00298 
00299 // ---- Standardfunkt ---- (arithmetische Operatoren)
00300 
00301 l_interval operator-(const l_interval & a) throw()
00302 {
00303    int precsave=stagprec;
00304    stagprec=a.prec;
00305 
00306    l_interval tmp;
00307 
00308    stagprec=precsave;
00309 
00310    int i;
00311    for(i=0;i<a.prec-1;i++)
00312       tmp.data[i]=-a.data[i];
00313 
00314    tmp.data[a.prec-1]=-a.data[a.prec];
00315    tmp.data[a.prec]=-a.data[a.prec-1];
00316 
00317    return tmp;
00318 } 
00319 
00320 // LI-LI
00321 
00322 l_interval operator+(const l_interval & li1, const l_interval & li2) throw()
00323 {
00324    l_interval li3;
00325    idotprecision idot(0.0);
00326    li1._akku_add(idot);
00327    li2._akku_add(idot);
00328    li3._akku_out(idot);
00329    return li3;
00330 }
00331 
00332 l_interval operator-(const l_interval & li1, const l_interval & li2) throw()
00333 {
00334    l_interval li3;
00335    idotprecision idot(0.0);
00336    li1._akku_add(idot);
00337    li2._akku_sub(idot);
00338    li3._akku_out(idot);
00339    return li3;
00340 } 
00341  
00342 /* l_interval operator*(const l_interval & li1, const l_interval & li2) throw()
00343 {
00344    l_interval li3;
00345    interval stdmul;
00346    
00347    stdmul = _interval(li1)*_interval(li2);
00348 
00349    if (abs(Inf(stdmul)) < MinReal) 
00350       li3 = _l_interval(stdmul) + 0.0; // Was das +0.0 soll ist mir ein Raetsel SR
00351    else 
00352    {               // if eingef�gt am 30.07.92 von Frederic Toussaint
00353       idotakku[0]=0.0;
00354       accumulate(idotakku[0], li1, li2);
00355       li3._akku_out();
00356       li3 = li3 & _l_interval(stdmul);  // Nachkorrektur
00357       // Wie bitte? (SR)
00358    }
00359    return li3;
00360 }  */
00361 
00362 l_interval operator * (const l_interval& li1, const l_interval& li2) throw()
00363 {  // Blomquist, Neue Version, 21.11.02;
00364     l_interval li3;
00365     interval stdmul;
00366     stdmul = _interval(li1) * _interval(li2); // grobe Einschlie�ung
00367     idotprecision idot(0.0);
00368     accumulate(idot,li1,li2);
00369     li3._akku_out(idot);       // li3: meist feinere Einschlie�ung! Aber bei zu
00370              // breiten li1,li2 ist dieses Ergebnis li3 breiter als stdmul!
00371     li3 = li3 & _l_interval(stdmul);   // Das optimale Ergebnis liegt daher 
00372               // im Durchschnitt der groben und der feineren Einschlie�ung!
00373     return li3;
00374 }    
00375 
00376 l_interval operator/(const l_interval & li1, const l_interval & li2) throw(ERROR_LINTERVAL_DIV_BY_ZERO)
00377 { 
00378 // ge�ndert am 29.07.92 von Frederic Toussaint
00379 // 13.5.93 AW: neuer Algorithmus nach der Beschreibung von W. Kraemer
00380 
00381    int         i, j, k, Z_sign;
00382    real        dzmitte;
00383    interval    dn = _interval(li2),
00384                dz = _interval(li1),
00385                stddiv = dz/dn;
00386    idotprecision idot;
00387    dotprecision dot1,dot2;
00388    l_interval  li3;
00389 
00390    li3._clear(1);
00391 
00392    if(!li2) 
00393    {
00394       cxscthrow(ERROR_LINTERVAL_DIV_BY_ZERO("l_interval operator/(const l_interval & li1, const l_interval & li2)"));
00395    } else 
00396    {
00397       dz = dz/dn;
00398       if (stagprec > 1) 
00399       {         // bei == 1 passiert nix!
00400          idot = 0.0;          // X, also Zaehler
00401          li1._akku_add(idot);
00402          k = 1;
00403          Z_sign = (Inf(dz) > 0.0);
00404          do {
00405             dzmitte = (Inf(dz) + Sup(dz)) / 2.0;  // mid(interval); grob!
00406             if (!!dz)                           
00407                if (Sup(abs(dz)) > MinReal)
00408                   if (diam(dz) < 1e-14 * abs(dzmitte) ) 
00409                   {
00410                      li3.elem(k) = dzmitte;
00411                      // Bestimme Rest in dotakku[2], [3] (Inf, Sup):
00412                      //
00413                      dot1 = 0.0;
00414                      for (i=1; i<=k; i++)                   // Rumpf
00415                         for (j=1; j<=li2.prec-1; j++)
00416                            accumulate(dot1, -li2.elem(j), li3.elem(i));
00417 
00418                      dot2 = dot1;        // Rumpf ist fuer Inf, Sup gleich
00419                      dot1 += Inf(idot);  // Inf += Inf(X)
00420                      dot2 += Sup(idot);  // Sup += Sup(X)
00421                      if (Z_sign)                      // INCL(Z) > 0.0
00422                      for (i=1; i<=k; i++) 
00423                      {
00424                         accumulate(dot1, -li2.elem(li2.prec+1), li3.elem(i));
00425                         accumulate(dot2, -li2.elem(li2.prec),   li3.elem(i));
00426                      } else                             // INCL(Z) < 0.0
00427                         for (i=1; i<=k; i++) 
00428                         {
00429                            accumulate(dot1, -li2.elem(li2.prec),   li3.elem(i));
00430                            accumulate(dot2, -li2.elem(li2.prec+1), li3.elem(i));
00431                         }
00432                      //
00433                      // Rausrunden in dz
00434                      dz = _interval(rnd(dot1, RND_DOWN),
00435                                     rnd(dot2, RND_UP));
00436                      dz = dz / dn;
00437                   }
00438             k++;
00439          } while (k < stagprec);
00440       } // if (stagprec > 1)
00441       li3.elem(stagprec) = Inf(dz);     // Intervall in die letzten Stelle
00442       li3.elem(stagprec+1) = Sup(dz);
00443       li3 = li3 & _l_interval(stddiv);  // Nachkorrektur
00444    } // if keine Null im Nenner
00445    return li3;
00446 }
00447 
00448 // ---- Vergleichsop. ----
00449 bool operator!(const l_interval & li) throw()
00450 {
00451    idotprecision idot(0.0);
00452    li._akku_add(idot);
00453    return (!idot);
00454 }
00455 
00456 /*l_interval::operator void *(void) throw()
00457 {
00458    idotakku[0]=0.0;
00459    _akku_add(idotakku[0]);
00460    return (idotakku[0]);
00461 }*/
00462 
00463 bool operator==(const l_interval & li1, const l_interval & li2) throw()
00464 {
00465    idotprecision idot1(0.0);
00466    idotprecision idot2(0.0);
00467    li1._akku_add(idot1);
00468    li2._akku_add(idot2);
00469    return (idot1==idot2);
00470 }   
00471 
00472 // ---- Mengenvergleiche ----
00473 
00474 bool operator<(const l_interval & li1, const l_interval & li2) throw()
00475 {
00476    idotprecision idot1(0.0);
00477    idotprecision idot2(0.0);
00478    li1._akku_add(idot1);
00479    li2._akku_add(idot2);
00480    return (idot1<idot2);
00481 } 
00482 
00483 bool operator>(const l_interval & li1, const l_interval & li2) throw()
00484 {
00485    idotprecision idot1(0.0);
00486    idotprecision idot2(0.0);
00487    li1._akku_add(idot1);
00488    li2._akku_add(idot2);
00489    return (idot1>idot2);
00490 } 
00491 
00492 bool operator<=(const l_interval & li1, const l_interval & li2) throw()
00493 {
00494    idotprecision idot1(0.0);
00495    idotprecision idot2(0.0);
00496    li1._akku_add(idot1);
00497    li2._akku_add(idot2);
00498    return (idot1<=idot2);
00499 } 
00500 
00501 bool operator>=(const l_interval & li1, const l_interval & li2) throw()
00502 { 
00503    idotprecision idot1(0.0);
00504    idotprecision idot2(0.0);
00505    li1._akku_add(idot1);
00506    li2._akku_add(idot2);
00507    return (idot1>=idot2);
00508 }
00509 
00510 void ConvexHull(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4) throw()
00511 {
00512    if(li1<=li2) 
00513    {                                      // Trivialfall 1
00514       li3=li2;
00515       li4=li2;
00516    } else if(li2<=li1) 
00517    {                                 // Trivialfall 2
00518       li3=li1;
00519       li4=li1;
00520    } else 
00521    {                                              // rechne
00522       idotprecision idot1(0.0);
00523       idotprecision idot2(0.0);
00524       li1._akku_add(idot1);
00525       li2._akku_add(idot2);
00526 
00527       idot1 |= idot2;
00528       // nun steht das richtige Ergebnis in idotakku[0]
00529       idot2 = idot1;                              // sichern
00530       li3._akku_out_inn(idot1);                              // Rundung nach innen
00531       idot1 = idot2;                              // und wiederherstellen
00532       li4._akku_out(idot1);                                  // ... nach aussen
00533    }
00534 }
00535 
00536 void Intersection(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4) throw(ERROR_LINTERVAL_EMPTY_INTERVAL)
00537 {
00538    if(li1<=li2) 
00539    {                                      // Trivialfall 1
00540       li3=li1;
00541       li4=li1;
00542    } else if(li2<=li1) 
00543    {                                 // Trivialfall 2
00544       li3=li2;
00545       li4=li2;
00546    } else 
00547    {                                              // rechne
00548       idotprecision idot1(0.0);
00549       idotprecision idot2(0.0);
00550       idotprecision idot;
00551       li1._akku_add(idot1);
00552       li2._akku_add(idot2);
00553       try
00554       {
00555          idot=(idot1 &= idot2);
00556       }
00557       catch(const EMPTY_INTERVAL &)
00558       {
00559          cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("void Intersection(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4)"));
00560       }
00561       // nun steht das richtige Ergebnis in idotakku[0]
00562       li3._akku_out_inn(idot);                              // Rundung nach innen
00563       idot = idot1;
00564       li4._akku_out(idot);                                  // ... nach aussen
00565    }
00566 }
00567 
00568 l_real     mid  (const l_interval & li) throw()
00569 {
00570    l_real lr;
00571 
00572    // --------------------------------------------------------
00573    //   dotakku[0] = Inf(li) + Sup(li)
00574 
00575    dotprecision dot(0.0);
00576    for (int j=1; j<=li.prec-1; j++)
00577       dot += li.elem(j);
00578 
00579    dot += dot; // mal 2
00580    dot += li.elem(li.prec);
00581    dot += li.elem(li.prec+1);
00582 
00583    //  ------------------------------------------------------
00584    //  Division nur bei ungleich 0
00585 
00586    if (dot != 0.0) 
00587    {
00588       Dotprecision ptr = *dot.ptr();
00589 
00590       // --------------------------------------------------------
00591       //  Dividiere dotakku[0] durch 2, mittels 1 Rechtsshift
00592 
00593       ptr[(a_intg)++ptr[A_END]] = ZERO;
00594           
00595       b_shr1 (&ptr[(a_intg)ptr[A_BEGIN]], (a_intg)(ptr[A_END]-ptr[A_BEGIN]+1));
00596           
00597       if (ptr[(a_intg)ptr[A_END]]   == ZERO) 
00598          --ptr[A_END];
00599       if (ptr[(a_intg)ptr[A_BEGIN]] == ZERO) 
00600          ++ptr[A_BEGIN];
00601 
00602       // --------------------------------------------------------
00603    }
00604 
00605    lr._akku_out(dot);
00606 
00607    return lr;
00608 }
00609  
00610 /* void accumulate(idotprecision & d, const l_interval & li1, const l_interval & li2) throw()
00611 {
00612    // �nderungen am 24.9.92 von F. Toussaint wegen Underflow-Fehlern
00613 
00614    int  i, j, n = 0;
00615    real tmp;
00616 
00617    for (i=1; i<=li1.prec-1; i++)
00618       for (j=1; j<=li2.prec-1; j++) 
00619       {
00620          tmp = abs(li1.elem(i)*li2.elem(j));
00621          if (tmp < MinReal && tmp != CXSC_Zero) 
00622             n++;
00623          else 
00624             accumulate(d, _interval(li1.elem(i)), _interval(li2.elem(j)));
00625       }
00626       for (i=1; i<=li2.prec-1; i++) 
00627       {
00628          tmp = abs(Inf(_interval(li1.elem(li1.prec), li1.elem(li1.prec+1)) * _interval(li2.elem(i))));
00629          if (tmp < MinReal && tmp != CXSC_Zero) 
00630             n++;
00631          else 
00632             accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
00633                           _interval(li2.elem(i)));
00634       }
00635       for (i=1; i<=li1.prec-1; i++) 
00636       {
00637          tmp = abs(Inf(_interval(li2.elem(li2.prec), li2.elem(li2.prec+1)) * _interval(li1.elem(i))));
00638          if (tmp < MinReal && tmp != CXSC_Zero) 
00639             n++;
00640          else 
00641             accumulate(d, _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)),
00642                           _interval(li1.elem(i)));
00643       }
00644       tmp = abs(Inf(_interval(li1.elem(li1.prec), li1.elem(li1.prec+1)) *
00645                     _interval(li2.elem(li2.prec), li2.elem(li2.prec+1))));
00646       if (tmp < MinReal && tmp != CXSC_Zero) 
00647          n++;
00648       else 
00649          accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
00650                        _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)));
00651       if (n > 0) 
00652          accumulate(d, _interval(_real(n)), 
00653                        _interval(-MinReal, MinReal));
00654 } */
00655 
00656 void accumulate(idotprecision & d, const l_interval & li1, 
00657                                    const l_interval & li2) throw()
00658 {   // Blomquist, Neue Version vom 21.11.02;
00659     // Die �nderungen von Toussaint wurden r�ckg�ngig gemacht.
00660    int  i,j;
00661 
00662    for (i=1; i<=li1.prec-1; i++)
00663       for (j=1; j<=li2.prec-1; j++) 
00664       {
00665             accumulate(d, _interval(li1.elem(i)), _interval(li2.elem(j)));
00666       }
00667       for (i=1; i<=li2.prec-1; i++) 
00668       {
00669             accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
00670                           _interval(li2.elem(i)));
00671       }
00672       for (i=1; i<=li1.prec-1; i++) 
00673       {
00674             accumulate(d, _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)),
00675                           _interval(li1.elem(i)));
00676       }
00677       accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
00678                     _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)));
00679 }
00680 
00681 
00682 /* void l_interval::_akku_out() throw()
00683 {  
00684    // ein im idotakku[0] liegendes Zwischenergebnis
00685    // wird in der entsprechenden Praezision in das aufrufende l_interval
00686    // gerundet. Rundung, also Einschluss von aussen!!
00687    // Alle MinReal wurden ge�ndert in minreal, Blomquist 19.11.02;
00688    _clear(1);
00689    interval z;
00690    real tmp, r, s;
00691    int i = 1,
00692    
00693    n = 0;
00694    z = rnd(idotakku[0]);
00695 
00696    while (i<prec && !(CXSC_Zero <= z)) 
00697    {       
00698       // z ist Intervall, deshalb <=
00699       if (succ(succ(Inf(z))) < Sup(z))
00700          break;                            // bei zu grobem Einschluss:
00701                                            // Intervall direkt in
00702                                            // Einschlusskomponente
00703                                            // schreiben
00704       tmp = Inf(z) + (Sup(z)-Inf(z))/2.0; // middle(interval)
00705       if (abs(Inf(z)) < minreal) 
00706       {
00707          if (abs(Sup(z)) < minreal) 
00708          {
00709             if (tmp != 0.0) 
00710             {
00711                tmp = 0.0;
00712                n++;
00713             }
00714             i = prec;
00715          }
00716       }
00717           
00718       this->elem(i) = tmp;
00719       idotakku[0] -= tmp;
00720       z = rnd(idotakku[0]);
00721       i++;
00722    } // while
00723    r = Inf(z);
00724    if (abs(r) < minreal) 
00725    {
00726       if (r < 0.0) 
00727          r = -minreal;
00728       else         
00729          r = 0.0;
00730    }
00731    s = Sup(z);
00732    if (abs(s) < minreal) 
00733    {
00734       if (s <= 0.0) 
00735          s = 0.0;
00736       else          
00737          s = minreal;
00738    }
00739    if (n > 0) 
00740    {  
00741       this->elem(prec) = r-_real(n+1)*minreal;   // Intervall in die letzten
00742       this->elem(prec+1) = s+_real(n+1)*minreal; // Stellen schreiben
00743       // (n+1), da Rundungsrichtung nicht bestimmt werden kann!
00744    } else 
00745    {
00746       this->elem(prec) = r;           // Intervall in die letzten Stellen
00747       this->elem(prec+1) = s;         // schreiben
00748    }
00749 }  */
00750 
00751 void l_interval::_akku_out(idotprecision& idot) throw()
00752 {  
00753    // ein im idotakku[0] liegendes Zwischenergebnis
00754    // wird in der entsprechenden Praezision in das aufrufende l_interval
00755    // so gerundet, dass idotakku[0] eingeschlossen wird.
00756    // Neueste Version: Blomquist 21.07.03;
00757    _clear(1);
00758    interval z;
00759    real tmp;
00760    int i = 1;
00761    
00762    z = rnd(idot);
00763    while (i<prec && !(CXSC_Zero <= z)) 
00764    {       
00765       // z ist Intervall, deshalb <=
00766       if (succ(succ(Inf(z))) < Sup(z))
00767          break;  // bei zu grobem Einschluss: Intervall direkt in
00768                  // Einschlusskomponente schreiben.
00769       tmp = mid(z);  // middle(interval)
00770       this->elem(i) = tmp;
00771       idot -= tmp;
00772       z = rnd(idot);
00773       i++;
00774    } // while
00775    this->elem(prec) = Inf(z);      // Intervall in die letzten Stellen
00776    this->elem(prec+1) = Sup(z);    // schreiben
00777 } // _akku_out()
00778 
00779 void l_interval::_akku_out_inn(idotprecision& idot) throw()
00780 { 
00781    // ein im idotakku[0] liegendes Zwischenergebnis
00782    // wird in der entsprechenden Praezision in die aufrufende l_interval Zahl
00783    // li so gerundet, dass gilt:  li enthalten in idotakku[0].
00784    _clear(1);
00785    real inf, sup, tmp;          // fuer Naeherungen, entspricht Interval z
00786    int i=1;
00787    inf = rnd(Inf(idot),RND_UP);
00788    sup = rnd(Sup(idot),RND_DOWN);
00789    
00790    if (inf>sup) 
00791       inf = sup;      // Vorsichtsmassnahme
00792         
00793    while (i<prec && !(inf<=CXSC_Zero&&sup>=CXSC_Zero)) 
00794    {
00795       tmp = inf + (sup-inf)/2.0; // middle(interval)
00796       this->elem(i) = tmp;
00797       idot -= tmp;
00798       inf = rnd(Inf(idot),RND_UP);
00799       sup = rnd(Sup(idot),RND_DOWN);
00800       if (inf>sup) 
00801          inf = sup;    // Vorsichtsmassnahme
00802       i++;
00803    }
00804    this->elem(prec)=inf;           // Intervall in die letzten Stellen
00805    this->elem(prec+1)=sup;         // schreiben
00806 }
00807 
00808 /* void l_interval::_akku_add(idotprecision& d) const throw()
00809 { 
00810    // addiert aufrufenden l_interval auf iakku d.
00811    // �nderung am 23.9.92 von F. Toussaint, da Fehler im Unterlaufbereich
00812    // Ausgeblendet von Blomquist, 20.11.02, da Fehler im Unterlaufbereich
00813    // nicht erkennbar. Meine neue Version direkt anschliessend!
00814    int        n = 0;
00815    real       r, s;
00816 
00817    for (int i=1; i<=prec-1; i++) 
00818    {
00819       r = this->elem(i);
00820       if (abs(r) < MinReal) 
00821       {
00822          if (r != 0.0) 
00823             n++;
00824       } else 
00825          d += _interval(r);
00826    }
00827    r = this->elem(prec);
00828    if (abs(r) < MinReal) 
00829    {
00830       if (r < 0.0) 
00831          r = -MinReal;
00832       else         
00833          r = 0.0;
00834    }
00835    s = this->elem(prec+1);
00836    if (abs(s) < MinReal) 
00837    {
00838       if (s <= 0.0) 
00839          s = 0.0;
00840       else          
00841          s = MinReal;
00842    }
00843    if (r != CXSC_Zero || s != CXSC_Zero) 
00844       d += _interval(r, s);
00845    if (n > 0) 
00846    {
00847       d += _interval(-_real(n+1)*MinReal, _real(n+1)*MinReal);
00848       // (n+1), da Rundungsrichtung nicht bestimmt werden kann!
00849    }
00850 } */
00851 
00852 void l_interval::_akku_add(idotprecision& d) const throw()
00853 { 
00854     // Addiert aufrufendes Intervall vom Typ l_interval auf d.
00855     // Meine neue Version; Blomquist, 20.11.02;
00856    real r,s;
00857    for (int i=1; i<=prec-1; i++)
00858    {   
00859        r = this->elem(i);
00860        if (sign(r) != CXSC_Zero) // Addition nur, falls n�tig!
00861            d += _interval(r); 
00862    }
00863    r = this->elem(prec);
00864    s = this->elem(prec+1);
00865    if (r != CXSC_Zero || s != CXSC_Zero) // Addition nur, falls n�tig!
00866       d += _interval(r, s);
00867 }
00868 
00869 /* void l_interval::_akku_sub(idotprecision& d) const throw()
00870 { 
00871    // Subtrahiert aufrufendes Intervall vom Typ l_interval von d.
00872    // Intervallsubtraktion!!
00873    // �nderung am 23.9.92 von F. Toussaint, da Fehler im Unterlaufbereich
00874    // Blomquist: Mir sind keine solchen Fehler bekannt, daher wurde diese
00875    // Version von mir ausgeklammert, 20.11.02; Neue Version anschliessend!
00876 
00877    int        n = 0;
00878    real       r, s;
00879 
00880    for (int i=1; i<=prec-1; i++) 
00881    {
00882       r = this->elem(i);
00883       if (abs(r) < MinReal) 
00884       {
00885          if (r != 0.0) 
00886             n++;
00887       } else 
00888          d -= _interval(r);
00889    }
00890    r = this->elem(prec);
00891    if (abs(r) < MinReal) 
00892    {
00893       if (r < 0.0) 
00894          r = -MinReal;
00895       else         
00896          r = 0.0;
00897    }
00898    s = this->elem(prec+1);
00899    if (abs(s) < MinReal) 
00900    {
00901       if (s <= 0.0) 
00902          s = 0.0;
00903       else          
00904          s = MinReal;
00905    }
00906    d -= _interval(r, s);
00907    if (n > 0) 
00908    {
00909       d -= _interval(-_real(n+1)*MinReal, _real(n+1)*MinReal);
00910       // (n+1), da Rundungsrichtung nicht bestimmt werden kann!
00911    }
00912 } */
00913 
00914 void l_interval::_akku_sub(idotprecision& d) const throw()
00915 { 
00916     // Subtrahiert aufrufendes Intervall vom Typ l_interval von d.
00917     // Meine neue Version; Blomquist, 20.11.02;
00918    real r,s;
00919 
00920    for (int i=1; i<=prec-1; i++)
00921    { 
00922        r = this->elem(i);
00923        if (sign(r) != CXSC_Zero) // Subtraktion nur, falls n�tig!
00924            d -= _interval(r); 
00925    }
00926    r = this->elem(prec);
00927    s = this->elem(prec+1);
00928    if (r != CXSC_Zero || s != CXSC_Zero) // Subtraktion nur, falls n�tig!
00929       d -= _interval(r, s);
00930 }
00931 
00932 // ---- Ausgabefunkt. ---------------------------------------
00933 
00934 std::ostream & operator << (std::ostream &s, const l_interval & a) throw()
00935 {
00936    idotprecision idot(0.0);
00937    a._akku_add(idot);
00938    s << idot;
00939    return s;
00940 }
00941 std::string & operator << (std::string &s, const l_interval & a) throw()
00942 {
00943    idotprecision idot(0.0);
00944    a._akku_add(idot);
00945    s << idot;
00946    return s;
00947 }
00948 
00949 std::istream & operator >> (std::istream &s, l_interval & a) throw()
00950 {
00951    idotprecision idot;
00952    s >> idot;
00953    a._akku_out(idot);
00954    return s;
00955 }
00956 
00957 std::string & operator >> (std::string &s, l_interval & a) throw()
00958 {
00959    idotprecision idot;
00960    s >> idot;
00961    a._akku_out(idot);
00962    return s;
00963 
00964 }
00965 
00966 void operator >>(const std::string &s,l_interval & a) throw()
00967 {
00968    std::string r(s);
00969    idotprecision idot;
00970    r >> idot;
00971    a._akku_out(idot);
00972 }
00973 void operator >>(const char *s,l_interval & a) throw()
00974 {
00975    std::string r(s);
00976    idotprecision idot;
00977    r>>idot;
00978    a._akku_out(idot);
00979 }
00980 
00981 int in ( const real& x, const l_interval& y )         // Contained-in relation
00982   { return ( (Inf(y) <= x) && (x <= Sup(y)) ); }      //----------------------
00983 
00984 int in ( const l_real& x, const l_interval& y )       // Contained-in relation
00985   { return ( (Inf(y) <= x) && (x <= Sup(y)) ); }      //----------------------
00986 
00987 int in ( const interval& x, const l_interval& y )     // Contained-in relation
00988 {                                                     //----------------------
00989   return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
00990 }
00991 
00992 int in ( const l_interval& x, const l_interval& y )   // Contained-in relation
00993 {                                                     //----------------------
00994   return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
00995 }
00996 
00997 l_interval Blow (const l_interval& x, const real& eps )
00998 {  
00999     int n;
01000     l_interval y;
01001     l_real lr1,lr2;
01002     y = x + interval(-eps,eps) * diam(x);
01003     lr1 = Inf(y);
01004     n = StagPrec(lr1);
01005     lr1[n] = pred(lr1[n]);
01006     lr2 = Sup(y);
01007     lr2[n] = succ(lr2[n]);
01008     return l_interval(lr1,lr2);
01009 }
01010 
01011 int Disjoint (const l_interval& a, const l_interval& b ) 
01012                                                  // Test for disjointedness
01013 {                                                //------------------------
01014     return (Inf(a) > Sup(b)) || (Inf(b) > Sup(a));
01015 }
01016 
01017 l_real AbsMin ( const l_interval& x )              // Absolute minimum of
01018 {                                                  // an interval
01019    if ( in(0.0,x) ) return l_real(0.0);
01020    else 
01021    { 
01022        l_real y(Inf(x));
01023        if (y > 0.0) return y;
01024        else return -Sup(x);
01025    }
01026 
01027 }
01028 
01029 l_real AbsMax (const l_interval& x )            // Absolute maximum of
01030 {                                               // an interval
01031     l_real a, b;                                //--------------------
01032 
01033     a = abs(Inf(x));
01034     b = abs(Sup(x));
01035 
01036     if (a > b)
01037         return a;
01038     else
01039         return b;
01040 }
01041 
01042 l_real RelDiam ( const l_interval x )                     // Relative diameter
01043 {                                                         // of an interval
01044   if ( in(0.0,x) )                                        //------------------
01045     return diam(x);
01046   else
01047     return( Sup( l_interval(diam(x)) / l_interval(AbsMin(x)) ) );
01048 }
01049 
01050 inline void times2pown(l_interval& x, int n) throw()
01051 { // x = x * 2^n;  Blomquist, 28.11.02;
01052     real mt,t;
01053     interval z;
01054     int p = StagPrec(x);
01055     if ( n<LI_Min_Exp_ || n>LI_maxexpo1 ) 
01056     { std::cerr << "Error in:  " 
01057            << "times2pown(l_interval& x, const int& n): " << std::endl
01058            << " -1074 <= n <= +1023 not fulfilled" << std::endl; exit(0); 
01059     }
01060     real F = comp(0.5,n+1);
01061     z = _interval(x[p],x[p+1]);
01062     times2pown(z,n);  // Scaling the original interval;
01063 
01064     for (int i=1; i<=p-1; i++)
01065     {
01066         mt = mant(x[i]);
01067         t = x[i];
01068         times2pown(x[i],n);
01069         if ( mt != mant(x[i]) ) 
01070         {
01071             x[i] = 0;
01072             z += _interval(t) * F;
01073         }
01074     }
01075     x[p]   = Inf(z);
01076     x[p+1] = Sup(z);
01077 } // times2pown(...)
01078 
01079 
01080 void Times2pown(l_interval& a, const real& p) throw()
01081 // The first parameter delivers an inclusion of a * 2^p;
01082 // For p in [-2100,+2100] p must be an integer value.
01083 // This condition is NOT tested in this function!
01084 // For p outside [-2100,+2100] an inclusion of a * 2^p is
01085 // calculated for any p of type real, unless an overflow occurs.
01086 // If the function is activated with the second parameter of type int, 
01087 // then the first parameter delivers correct inclusions of a * 2^p, 
01088 // unless an overflow occurs.
01089 // Blomquist, 28.01.2008;
01090 {
01091     const int c2 = 2100;
01092     int ex(expo_gr(a)),fac,rest,n;
01093     double dbl;
01094 
01095     if (ex > -1000000)
01096     {
01097         if (p>=0)
01098             if (p>c2)
01099                 times2pown(a,c2); // Produces an error
01100             else // 0 <= p <= 2100
01101             {
01102                 dbl = _double(p);
01103                 n = (int) dbl;
01104                 fac = n/LI_maxexpo1;
01105                 rest = n%LI_maxexpo1;
01106                 for (int k=1; k<=fac; k++)
01107                     times2pown(a,LI_maxexpo1);
01108                 times2pown(a,rest);
01109             }
01110         else // p<0; No overflow or underflow!
01111             if (p<-c2)
01112             {
01113                 if (0<a)
01114                     a = l_interval(-minreal,minreal);
01115                 else
01116                     if (Inf(a)>=0)
01117                         a = l_interval(0,minreal);
01118                     else a = l_interval(-minreal,0);
01119             }
01120             else // -2100 <= p < 0
01121             {
01122                 dbl = _double(p);
01123                 n = (int) dbl;
01124                 fac = n/LI_Min_Exp_;
01125                 rest = n%LI_Min_Exp_;
01126                 for (int k=1; k<=fac; k++)
01127                     times2pown(a,LI_Min_Exp_);
01128                 times2pown(a,rest);
01129             }
01130     }
01131 } // Times2pown(...)
01132 
01133 } // namespace cxsc
01134