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: 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