- Generated by
1.7.5.1
|
C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: 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
1.7.5.1