C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
cinterval.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: cinterval.cpp,v 1.20 2014/01/30 17:23:44 cxsc Exp $ */
00025 
00026 #include "cinterval.hpp"
00027 #include "cidot.hpp"
00028 #include "dot.hpp"
00029 #include "rmath.hpp"
00030 #include "imath.hpp"
00031 
00032 namespace cxsc {
00033 
00034 #define CXSC_Zero 0.0
00035 
00036 cinterval::cinterval(const dotprecision &a) throw() : re(a),im(0) {}
00037 cinterval::cinterval(const idotprecision &a) throw() : re(a),im(0) {}
00038 cinterval::cinterval(const cdotprecision &a) throw() : re(Re(a)),im(Im(a)) {}
00039 cinterval::cinterval(const cidotprecision &a) throw() : 
00040    re(rnd(InfRe(a),RND_DOWN),rnd(SupRe(a),RND_UP)),
00041    im(rnd(InfIm(a),RND_DOWN),rnd(SupIm(a),RND_UP))  
00042 {
00043 }
00044 
00045 
00046 cinterval mult_operator(const cinterval & a,const cinterval & b) throw()
00047 {
00048    cidotprecision akku;
00049    akku=0.0;
00050    accumulate(akku,a,b);
00051    return rnd(akku);
00052 } 
00053  
00054 // In complex.cpp
00055 void product  (real, real, real, real, int&, real&, interval&);
00056 real quotient (real, interval, real, interval, int, int, int);
00057 
00058 // optimale komplexe Intervalldivision
00059 
00060 bool cxsc_complex_division_p[5];
00061 
00062 real cxsc_complex_division_f(real a, real b, real c, real d, int round)
00063 {
00064    int      zoverfl,noverfl;
00065    real     z1,n1;
00066    interval z2,n2;
00067 
00068    //  f:=(a*c+b*d)/(SQR(c)+SQR(d))
00069 
00070    product( a, c, b, d, zoverfl, z1,z2 );
00071    product( c, c, d, d, noverfl, n1,n2 );
00072 
00073    return quotient( z1,z2, n1,n2, round, zoverfl, noverfl );
00074 }
00075 
00076 static real minmax(int minimum, real a, real b, real y0,
00077                    interval x, int i, int j)
00078 // Calculates the inner minimum or maximum of f = (ac+bd)/(cc+dd)
00079 // on the interval x = [c.inf,c.sup] ( a,b,d=y0 fixated ).
00080 // If minimum = true the minimum will be calculated, otherwise the maximum.
00081 
00082 {
00083    real q1,q2,w1,w2,t1,t2,x1,x2,ay0, minmax;
00084    dotprecision akku;
00085    bool scaling = false;  // scaling = true <==> scaling is necessary
00086 
00087    if (minimum) 
00088       minmax = MaxReal;
00089    else         
00090       minmax = -MaxReal;
00091 
00092    if (Inf(x) == Sup(x)) 
00093    {
00094       if (cxsc_complex_division_p[i] && cxsc_complex_division_p[j])
00095          minmax = cxsc_complex_division_f( a, b, Inf(x), y0, 1-2*minimum );
00096       
00097       cxsc_complex_division_p[i] = false;
00098       cxsc_complex_division_p[j] = false;
00099   } else 
00100   if (a == 0.0) 
00101   {
00102       if ( b == CXSC_Zero  ||  y0 == CXSC_Zero ) 
00103       {
00104          minmax = 0.0;
00105          cxsc_complex_division_p[i]   = false;
00106          cxsc_complex_division_p[j]   = false;
00107       } else 
00108       {
00109          if (0.0 < x) {
00110             if (minimum  && sign(b) != sign(y0) ) 
00111             {
00112                minmax = divd(b, y0);
00113                cxsc_complex_division_p[i]   = false;
00114                cxsc_complex_division_p[j]   = false;
00115             } else 
00116             if (!minimum  &&  sign(b) == sign(y0) ) 
00117             {
00118                minmax =  divu(b, y0);
00119                cxsc_complex_division_p[i]   = false;
00120                cxsc_complex_division_p[j]   = false;
00121             }
00122          }   
00123       }
00124    } else 
00125    { 
00126       //  a != 0.0
00127       if (y0 == 0.0) 
00128       {
00129          if (minimum) 
00130          {
00131             if (a > 0.0) 
00132                minmax = divd(a, Sup(x));
00133             else         
00134                minmax = divd(a, Inf(x));
00135          } else 
00136          {
00137             if (a > 0.0) 
00138                minmax = divu(a, Inf(x));
00139             else         
00140                minmax = divu(a, Sup(x));
00141          }
00142          cxsc_complex_division_p[i] = false;
00143          cxsc_complex_division_p[j] = false;
00144       } else 
00145       { 
00146          // y0 != 0.0, Calculation of extrema points and minimum|maximum
00147          // values: 
00148          //  IF NOTBEKANNT THEN
00149          //  Calculation of: t = sign(a) * ( |b/a| + sqrt( 1+|b/a|^2 ) )
00150          //  in staggered presentation:  t ~ t1 + t2
00151 
00152          // Exponent over-/undeflow in |b/a| is now considered from 
00153          // Blomquist/Hofschuster 06.11.02.
00154 
00155          real invf2=1, a_skal;
00156          int exf1=0,  exf2=0,  exinf1=0,  exinf2=0;
00157 
00158          if (sign(b)==0) { t1 = 1;  t2 = 0; }
00159          else
00160          {  // b != 0, ---> expo(b) != -2147483647;
00161          if (a>0.0) 
00162             q2 =  abs(b);
00163          else       
00164             q2 = -abs(b); 
00165 
00166          // Skaling to avoid overflow by division b/a: 
00167          int expo_diff = expo(b)-expo(a), ex;
00168          if (expo_diff >= 512) 
00169          {
00170              int exdiff5 = expo_diff - 500;
00171              scaling = true;
00172              if (exdiff5 > 1024) // Two scaling factors necessary!
00173              {   
00174                  ex = exdiff5 / 2;
00175                  exf1 = ex-1;
00176                  exf2 = exdiff5 - ex;
00177                  exinf1 = -exf1;
00178                  invf2 = comp(0.5,1-exf2);
00179                  exinf2 = -exf2;
00180                  a_skal = a;
00181                  times2pown(a_skal,exf1);
00182                  times2pown(a_skal,exf2);
00183              }
00184              else // Scaling with only one factor!
00185              {   
00186                  exf2 = exdiff5-1;    // exf1 = 0;
00187                  invf2 = comp(0.5,2-exdiff5);  // invf1 = 1;
00188                  exinf2 = -exf2;  // exinf1 = 0;
00189                  a_skal = a;  
00190                  times2pown(a_skal,exf2);
00191              }
00192          }
00193          else // Scaling not necessary!
00194          { a_skal = a; }
00195 
00196          q1   = q2/a_skal;
00197          akku = q2;
00198          accumulate(akku, -q1, a_skal);
00199          q2   = rnd(akku) / a_skal;
00200 
00201          akku = 0.0;
00202          if (exinf1 == 0) accumulate(akku, invf2, invf2);        
00203          accumulate(akku, q1, q1);
00204          accumulate(akku, q1, q2);
00205          accumulate(akku, q1, q2);
00206          accumulate(akku, q2, q2);
00207 
00208          w1   = sqrt(rnd(akku));
00209 
00210          accumulate(akku, -w1, w1);
00211          w2 = rnd(akku) / (2.0*w1);
00212 
00213          akku  = q1;
00214          akku += q2;
00215          akku += w1;
00216          akku += w2;
00217 
00218          t1 = rnd(akku);
00219 
00220          akku -= t1;
00221          t2 = rnd(akku);
00222          }
00223          if (a<0.0)  // if (a_skale<0.0)
00224          {
00225             t1 = -t1;
00226             t2 = -t2;
00227          }
00228 
00229          // Fall differentiation for min-,max- calculation:
00230          ay0 = abs(y0);
00231          if (( sign(b) == sign(y0) ) == minimum) 
00232          {
00233             //   Calculation of:  x1 + x2  =  |y0| * ( t1 + t2 )
00234             akku  = 0.0;
00235             accumulate(akku,ay0,t1);
00236             accumulate(akku,ay0,t2);
00237             x1    = rnd(akku);
00238             if (expo(x1) == 2147483647) goto Ende;
00239             akku -= x1;
00240             x2    = rnd(akku);
00241             if (scaling) 
00242             { 
00243                 if (expo(x1)+exf1 > 1023) goto Ende;
00244                 times2pown(x1,exf1);
00245                 if (expo(x1)+exf2 > 1023) goto Ende;
00246                 times2pown(x1,exf2); 
00247                 times2pown(x2,exf1);
00248                 times2pown(x2,exf2);
00249             }
00250          } else 
00251          {
00252             //  Calculation of:  x1 + x2  =  |y0| / ( t1 + t2 )
00253             x1 = ay0 / t1;
00254             akku = ay0;
00255             accumulate(akku, -t1, x1);
00256             accumulate(akku, -t2, x1);
00257             x2 = rnd(akku) / t1;
00258             if (scaling) 
00259             {  
00260                 if (expo(x1)+exinf1 > 1023) goto Ende;
00261                 times2pown(x1,exinf1);
00262                 if (expo(x1)+exinf2 > 1023) goto Ende;
00263                 times2pown(x1,exinf2);
00264                 times2pown(x2,exinf1);
00265                 times2pown(x2,exinf2);
00266             }
00267          }
00268          if (minimum) 
00269          {
00270             x1 = -x1;
00271             x2 = -x2;
00272          }
00273 
00274          if (x1 < x) 
00275          {
00276             //  Calculation of:  a / ( 2*(x1+x2) )
00277             q1   = a/(2*x1);
00278             akku = 0.0;
00279             accumulate(akku, -x1, q1);  // vorher: accumulate(akku, x1, q1);
00280             accumulate(akku, -x2, q1);
00281 
00282             // exact calculation of (a + akku + akku) in new variable akku:
00283             akku += akku;  
00284             akku += a; 
00285             q2 = rnd(akku) / (2.0*x1);
00286             if (minimum) 
00287             {
00288                 if (sign(q2)==0 && sign(akku)!=0) minmax = pred(q1);
00289                 else minmax = addd(q1, q2);
00290             }
00291             else 
00292             {   
00293                 if (sign(q2)==0 && sign(akku)!=0) minmax = succ(q1);
00294                 else minmax = addu(q1, q2);
00295             }
00296         
00297             cxsc_complex_division_p[i] = false;
00298             cxsc_complex_division_p[j] = false;
00299          }
00300       Ende:;
00301       }  // y0 != 0.0
00302    }
00303    return minmax;
00304 } // *** minmax ***
00305 
00306 cinterval cidiv(const cinterval& A, const cinterval& B)
00307 {
00308    real     realteilINF, realteilSUP,
00309             imagteilINF, imagteilSUP;
00310    // da sonst eventuell zwischendurch illegale Intervalle entstehen
00311    real     a0,b0;
00312    bool     a_repeat,b_repeat; 
00313    int      i, rep, j;
00314    real     AREINF, ARESUP, AIMINF, AIMSUP,
00315             BREINF, BRESUP, BIMINF, BIMSUP;
00316    interval ARE, AIM, BRE, BIM;
00317 
00318    // keine Fehlerabfrage -> ist schon in CINTVAL.CPP
00319    //  IF ( 0.0 IN B.RE ) AND ( 0.0 IN B.IM ) THEN
00320    //   CIDIVISION:= COMPL( 1.0 / INTVAL(-1.0,1.0), INTVAL(0.0) );
00321    //   Fehlerabbruch erzwingen: Intervall enthaelt 0
00322 
00323    //  *** Berechnung des Realteils ***
00324 
00325    AREINF = Inf(Re(A));
00326    ARESUP = Sup(Re(A));
00327    AIMINF = Inf(Im(A));
00328    AIMSUP = Sup(Im(A));
00329    BREINF = Inf(Re(B));
00330    BRESUP = Sup(Re(B));
00331    BIMINF = Inf(Im(B));
00332    BIMSUP = Sup(Im(B));
00333    ARE    = Re(A);
00334    AIM    = Im(A);
00335    BRE    = Re(B);
00336    BIM    = Im(B);
00337 
00338    a_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
00339    b_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
00340 
00341    if (a_repeat || b_repeat) 
00342       rep = 2;
00343    else                      
00344       rep = 1;
00345 
00346    if (BREINF >= 0.0) 
00347       a0 = ARESUP;
00348    else               
00349       a0 = AREINF;
00350   
00351    if (BIMINF >= 0.0) 
00352       b0 = AIMSUP;
00353    else               
00354       b0 = AIMINF;
00355 
00356    realteilSUP = -MaxReal;
00357 
00358    for (j=1; j<=rep; j++) 
00359    {
00360       for (i=1; i<=4; cxsc_complex_division_p[i++] = true);
00361     
00362       realteilSUP =
00363              max( realteilSUP,
00364                    max( max( minmax( false, a0, b0, BIMINF, BRE, 1,2 ),
00365                              minmax( false, a0, b0, BIMSUP, BRE, 3,4 ) ),
00366                         max( minmax( false, b0, a0, BREINF, BIM, 1,3 ),
00367                              minmax( false, b0, a0, BRESUP, BIM, 2,4 ) ) )
00368 
00369                 );
00370                 
00371       if (cxsc_complex_division_p[1])
00372          realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BREINF, BIMINF, +1 ) );
00373       if (cxsc_complex_division_p[2])
00374          realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BRESUP, BIMINF, +1 ) );
00375       if (cxsc_complex_division_p[3])
00376          realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BREINF, BIMSUP, +1 ) );
00377       if (cxsc_complex_division_p[4])
00378          realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BRESUP, BIMSUP, +1 ) );
00379 
00380       if (a_repeat) 
00381          a0 = ARESUP; 
00382       else if (b_repeat) 
00383          b0 = AIMSUP;
00384    }
00385 
00386    if (BREINF >= 0.0) 
00387       a0 = AREINF;
00388    else               
00389       a0 = ARESUP;
00390    if (BIMINF >= 0.0) 
00391       b0 = AIMINF;
00392    else               
00393       b0 = AIMSUP;
00394 
00395    realteilINF = MaxReal;
00396 
00397    for (j=1; j<=rep; j++) 
00398    {
00399       for (i=1; i<=4; cxsc_complex_division_p[i++] = true);
00400       
00401       realteilINF =
00402               min( realteilINF,
00403                    min( min( minmax( true, a0, b0, BIMINF, BRE, 1,2 ),
00404                              minmax( true, a0, b0, BIMSUP, BRE, 3,4 ) ),
00405                         min( minmax( true, b0, a0, BREINF, BIM, 1,3 ),
00406                              minmax( true, b0, a0, BRESUP, BIM, 2,4 ) ) )
00407                  );
00408       if (cxsc_complex_division_p[1])
00409          realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BREINF, BIMINF, -1 ) );
00410       if (cxsc_complex_division_p[2])
00411          realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BRESUP, BIMINF, -1 ) );
00412       if (cxsc_complex_division_p[3])
00413          realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BREINF, BIMSUP, -1 ) );
00414       if (cxsc_complex_division_p[4])
00415          realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BRESUP, BIMSUP, -1 ) );
00416 
00417       if (a_repeat) 
00418          a0 = AREINF; 
00419       else if (b_repeat) 
00420          b0 = AIMINF;
00421    }
00422 
00423 
00424    //  Calculation of the img. part: 
00425    //  g(a, b, c, d) = cxsc_complex_division_f(b, -a, c, d) 
00426 
00427    a_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
00428    b_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
00429 
00430    //  IF a_repeat OR b_repeat THEN rep:= 2 ELSE rep:= 1;  
00431 
00432    if (BREINF >= 0.0) 
00433       b0 = AIMSUP;
00434    else 
00435       b0 = AIMINF;
00436    
00437    if (BIMINF >= 0.0) 
00438       a0 = AREINF;
00439    else 
00440       a0 = ARESUP;
00441 
00442    imagteilSUP = -MaxReal;
00443 
00444    for (j=1; j<=rep; j++) 
00445    {
00446       for (i=1; i<=4; cxsc_complex_division_p[i++] = true) ;
00447     
00448       imagteilSUP =
00449               max( imagteilSUP,
00450                    max( max( minmax( false,  b0, -a0, BIMINF, BRE, 1,2 ),
00451                              minmax( false,  b0, -a0, BIMSUP, BRE, 3,4 ) ),
00452                         max( minmax( false, -a0,  b0, BREINF, BIM, 1,3 ),
00453                              minmax( false, -a0,  b0, BRESUP, BIM, 2,4 ) ) )
00454                  );
00455     
00456       if (cxsc_complex_division_p[1])
00457          imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BREINF, BIMINF, +1 ) );
00458       if (cxsc_complex_division_p[2])
00459          imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BRESUP, BIMINF, +1 ) );
00460       if (cxsc_complex_division_p[3])
00461          imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BREINF, BIMSUP, +1 ) );
00462       if (cxsc_complex_division_p[4])
00463          imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BRESUP, BIMSUP, +1 ) );
00464 
00465       if (b_repeat) 
00466          b0 = AIMSUP;  
00467       else if (a_repeat) 
00468          a0 = AREINF  ;
00469    }
00470 
00471    if (BREINF >= 0.0) 
00472       b0 = AIMINF;
00473    else 
00474       b0 = AIMSUP;
00475   
00476    if (BIMINF >= 0.0) 
00477       a0 = ARESUP;
00478    else 
00479       a0 = AREINF;
00480 
00481    imagteilINF = MaxReal;
00482 
00483    for (j=1; j<=rep; j++) 
00484    {
00485       for (i=1; i<=4; cxsc_complex_division_p[i++] = true) ;
00486     
00487       imagteilINF =
00488               min( imagteilINF,
00489                    min( min( minmax( true,  b0, -a0, BIMINF, BRE, 1,2 ),
00490                              minmax( true,  b0, -a0, BIMSUP, BRE, 3,4 ) ),
00491                         min( minmax( true, -a0,  b0, BREINF, BIM, 1,3 ),
00492                              minmax( true, -a0,  b0, BRESUP, BIM, 2,4 ) ) )
00493                  );
00494     
00495       if (cxsc_complex_division_p[1])
00496          imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BREINF, BIMINF, -1 ) );
00497       if (cxsc_complex_division_p[2])
00498          imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BRESUP, BIMINF, -1 ) );
00499       if (cxsc_complex_division_p[3])
00500          imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BREINF, BIMSUP, -1 ) );
00501       if (cxsc_complex_division_p[4])
00502          imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BRESUP, BIMSUP, -1 ) );
00503 
00504       if (b_repeat) 
00505          b0 = AIMINF;  
00506       else if (a_repeat) 
00507          a0 = ARESUP;
00508    }
00509 
00510    return cinterval(interval(realteilINF, realteilSUP),
00511                      interval(imagteilINF, imagteilSUP));
00512 }  //    CIDIVISION
00513 
00514 cinterval C_point_div(const cinterval& z, const cinterval& n)
00515 // Division of complex point intervals; 
00516 // z,n must be point intervals!!  Blomquist, 07,11.02
00517 // This function only for internal use!
00518 {
00519     complex a,b,q1,q2;
00520     a = complex(InfRe(z),InfIm(z));
00521     b = complex(InfRe(n),InfIm(n));
00522     q1 = divd(a,b);
00523     q2 = divu(a,b);
00524 
00525     interval re, im;
00526     re = interval( Re(q1),Re(q2) );
00527     im = interval( Im(q1),Im(q2) );
00528 
00529     return cinterval(re,im);
00530 }  // C_point_div
00531 
00532 cinterval div_operator (const cinterval & a, const cinterval & b) throw(DIV_BY_ZERO)
00533 {
00534     bool a_point, b_point;
00535     a_point = InfRe(a)==SupRe(a) && InfIm(a)==SupIm(a);
00536     b_point = InfRe(b)==SupRe(b) && InfIm(b)==SupIm(b);
00537     if(a_point && b_point) return C_point_div(a,b); // a,b are point intervals
00538     else return cidiv(a,b);
00539 }
00540 
00541 interval abs(const cinterval &a) throw()
00542 {
00543 //    idotakku[2]=0;
00544 //    accumulate(idotakku[2],a.re,a.re);
00545 //    accumulate(idotakku[2],a.im,a.im);
00546 //    return sqrt(rnd(idotakku[2]));
00547     return sqrtx2y2(a.re,a.im);
00548 }
00549 
00550 
00551 // ---- Ausgabefunkt. ---------------------------------------
00552 
00553 std::ostream & operator << (std::ostream &s, const cinterval& a) throw()
00554 {
00555    s << '('          
00556      << a.re << ','  
00557      << a.im       
00558      << ')';
00559    return s;
00560 }
00561 std::string & operator << (std::string &s, const cinterval &a) throw()
00562 {
00563    s+='(';
00564    s << a.re;
00565    s+=',';
00566    s << a.im; 
00567    s+=')';
00568    return s;
00569 }
00570 
00571 std::istream & operator >> (std::istream &s, cinterval &a) throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
00572 { // New version for cinterval input; Blomquist, 27.10.02;
00573     char c;
00574     skipeolnflag = inpdotflag = true;
00575     c = skipwhitespacessinglechar (s, '(');
00576     if (inpdotflag) s.putback(c);
00577     c = skipwhitespacessinglechar (s, '[');
00578     if (inpdotflag) s.putback(c);
00579     s >> SaveOpt >> RndDown >> Inf(a.re);
00580     skipeolnflag = inpdotflag = true; 
00581     c = skipwhitespacessinglechar (s, ','); 
00582     if (inpdotflag) s.putback(c);
00583     s  >> RndUp >> Sup(a.re);
00584     c = skipwhitespacessinglechar (s, ']');
00585     if (inpdotflag) s.putback(c);
00586     c = skipwhitespacessinglechar (s, ',');
00587     if (inpdotflag) s.putback(c);
00588     
00589     c = skipwhitespacessinglechar (s, '[');
00590     if (inpdotflag) s.putback(c);
00591     s >> RndDown >> Inf(a.im);
00592     skipeolnflag = inpdotflag = true; 
00593     c = skipwhitespacessinglechar (s, ','); 
00594     if (inpdotflag) s.putback(c);
00595     s  >> RndUp >> Sup(a.im) >> RestoreOpt;
00596 
00597    if (!waseolnflag) 
00598    {
00599       skipeolnflag = false, inpdotflag = true;
00600       c = skipwhitespaces (s);
00601       if (inpdotflag && c != ']') 
00602          s.putback(c);
00603    }
00604    if (!waseolnflag) 
00605    {
00606       skipeolnflag = false, inpdotflag = true;
00607       c = skipwhitespaces (s);
00608       if (inpdotflag && c != ')') 
00609          s.putback(c);
00610          }
00611    if (Inf(a.re) > Sup(a.re) || Inf(a.im) > Sup(a.im)) 
00612    cxscthrow(ERROR_CINTERVAL_EMPTY_INTERVAL("std::istream & operator >> (std::istream &s, cinterval &a)"));
00613       
00614    return s;
00615 }
00616 
00617 std::string & operator >> (std::string &s, cinterval &a) throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
00618 {
00619    s = skipwhitespacessinglechar (s, '(');
00620    s = skipwhitespacessinglechar (s, '[');
00621    s = s >> SaveOpt >> RndDown >> Inf(a.re);
00622    s = skipwhitespacessinglechar (s, ',');
00623    s = s >> RndUp >> Sup(a.re);
00624    s = skipwhitespacessinglechar (s, ']');
00625    s = skipwhitespacessinglechar (s, ',');
00626    s = skipwhitespacessinglechar (s, '[');
00627    s = s >> RndDown >> Inf(a.im);
00628    s = skipwhitespacessinglechar (s, ',');
00629    s = s >> RndUp >> Sup(a.im) >> RestoreOpt;
00630    s = skipwhitespaces (s);
00631    if (s[0] == ']') 
00632       s.erase(0,1);
00633    s = skipwhitespaces (s);
00634    if (s[0] == ')') 
00635       s.erase(0,1);
00636 
00637    if (Inf(a.re) > Sup(a.re) || Inf(a.im) > Sup(a.im)) 
00638       cxscthrow(ERROR_CINTERVAL_EMPTY_INTERVAL("std::string & operator >> (std::string &s, cinterval &a)"));
00639 
00640    return s;
00641 }
00642 
00643 void operator >>(const std::string &s,cinterval &a) throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
00644 {
00645    std::string r(s);
00646    r>>a;
00647 }
00648 void operator >>(const char *s,cinterval &a) throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
00649 {
00650    std::string r(s);
00651    r>>a;
00652 }
00653 
00654 int in ( const cinterval& x, const cinterval& y )    // Contained-in-the-interior relation
00655 {                                        //-----------------------------------
00656   return ( in(Re(x),Re(y)) && in(Im(x),Im(y)) );
00657 }
00665 cinterval Blow ( cinterval x, const real& eps )           // Epsilon inflation
00666 {                                                         //------------------
00667   return cinterval(Blow(Re(x),eps),Blow(Im(x),eps));
00668 }
00669 } // namespace cxsc
00670 
00671 
00672 
00673 
00674 
00675 
00676 
00677 
00678 
00679 
00680 
00681 
00682 
00683 
00684