C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
l_interval.cpp
1/*
2** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3**
4** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5** Universitaet Karlsruhe, Germany
6** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7** Universitaet Wuppertal, Germany
8**
9** This library is free software; you can redistribute it and/or
10** modify it under the terms of the GNU Library General Public
11** License as published by the Free Software Foundation; either
12** version 2 of the License, or (at your option) any later version.
13**
14** This library is distributed in the hope that it will be useful,
15** but WITHOUT ANY WARRANTY; without even the implied warranty of
16** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17** Library General Public License for more details.
18**
19** You should have received a copy of the GNU Library General Public
20** License along with this library; if not, write to the Free
21** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22*/
23
24/* CVS $Id: l_interval.cpp,v 1.21 2014/01/30 17:23:46 cxsc Exp $ */
25
26#include <math.h>
27#include "l_interval.hpp"
28#include "idot.hpp"
29#include "RtsFunc.h"
30
31namespace cxsc {
32
33#define CXSC_Zero 0.
34
35l_interval & l_interval::operator= (const l_interval & a) noexcept
36{
37 real *newdata=new real[a.prec+1];
38 int i;
39 for(i=0;i<=a.prec;i++)
40 newdata[i]=a.data[i];
41 delete [] data;
42 data=newdata;
43 prec=a.prec;
44 return *this;
45}
46
47// ---- Typwandlungen ----
48
49interval::interval(const l_interval & a) noexcept
50{
51 idotprecision idot(0.0);
52 a._akku_add(idot);
53 *this=rnd(idot);
54}
55
57{
58 idotprecision idot(0.0);
59 a._akku_add(idot);
60 return *this=rnd(idot);
61}
62
63
64interval _interval(const l_interval & a) noexcept
65{
66 idotprecision idot(0.0);
67 a._akku_add(idot);
68 return rnd(idot);
69}
70
72#if (CXSC_INDEX_CHECK)
73
74#else
75 noexcept
76#endif
77{
78 _allo(stagprec);
79 idotprecision idot(a);
80 _akku_out(idot);
81}
82
84#if (CXSC_INDEX_CHECK)
85
86#else
87
88#endif
89{
90 if(a>b)
91 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const dotprecision&,const dotprecision&)"));
92 _allo(stagprec);
93 idotprecision idot;
94 UncheckedSetInf(idot,a);
95 UncheckedSetSup(idot,b);
96 _akku_out(idot);
97}
98
99l_interval & l_interval::operator =(const dotprecision & a)
100#if (CXSC_INDEX_CHECK)
101
102#else
103 noexcept
104#endif
105{
106 if(stagprec!=prec)
107 {
108 delete [] data;
109 _allo(stagprec);
110 }
111 idotprecision idot(a);
112 _akku_out(idot);
113 return *this;
114}
115
117#if (CXSC_INDEX_CHECK)
118
119#else
120 noexcept
121#endif
122{
123 _allo(stagprec);
124 idotprecision idot(a);
125 _akku_out(idot);
126}
127
128l_interval & l_interval::operator =(const idotprecision & a)
129#if (CXSC_INDEX_CHECK)
130
131#else
132 noexcept
133#endif
134{
135 if(stagprec!=prec)
136 {
137 delete [] data;
138 _allo(stagprec);
139 }
140 idotprecision idot(a);
141 _akku_out(idot);
142 return *this;
143}
144
146#if (CXSC_INDEX_CHECK)
147
148#else
149
150#endif
151{
152 _allo(stagprec);
153 if(a>b)
154 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const l_real &a, const l_real &b)"));
155 dotprecision dot1, dot2;
156 dot1 = a;
157 dot2 = b;
158 idotprecision idot(dot1,dot2);
159 _akku_out(idot);
160}
161
163#if (CXSC_INDEX_CHECK)
164
165#else
166
167#endif
168{
169 _allo(stagprec);
170 if(a>b)
171 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const l_real &a, const l_real &b)"));
172 dotprecision dot1, dot2;
173 dot1 = a;
174 dot2 = b;
175 idotprecision idot(dot1,dot2);
176 _akku_out(idot);
177}
178
179void l_realz2l_interval(const l_real& lr, const interval& z,
180 l_interval& li) noexcept
181{ // converts lr+z to li of type l_interval; Blomquist, 12.10.06;
182 // lr+z is included by li.
183 // prec(lr) <= prec(li) must be realized!
184 int p = StagPrec(lr); int q = StagPrec(li);
185 if(p>q)
186 {
187 std::cerr << "l_realz2l_interval(const l_real& lr, const interval& z, l_interval& li): incorrect precisions of lr,li !"
188 << std:: endl;
189 exit (1);
190 }
191
192 for (int i=1; i<=q-1; i++) li[i] = 0;
193 li[q] = Inf(z);
194 li[q+1] = Sup(z);
195 if(p<q) for (int i=1; i<=p; i++) li[i] = lr[i];
196 else // p=q
197 {
198 p--;
199 for (int i=1; i<=p; i++) li[i] = lr[i];
200 li[q] = addd(lr[p+1],Inf(z));
201 li[q+1] = addu(lr[p+1],Sup(z));
202 }
203} // l_realz2l_interval
204
206#if (CXSC_INDEX_CHECK)
207
208#else
209
210#endif
211{
212 _allo(stagprec);
213 if(a>b)
214 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval::l_interval(const l_real &a, const l_real &b)"));
215 dotprecision dot1, dot2;
216 dot1 = a;
217 dot2 = b;
218 idotprecision idot(dot1,dot2);
219 _akku_out(idot);
220}
221
222
223/*
224l_interval _l_interval(const l_real & a, const l_real & b)
225{
226 if(a>b)
227 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("l_interval _l_interval(const l_real & a, const l_real & b)"));
228
229 dotakku[0]=0.0;
230 dotakku[1]=0.0;
231 a._akku_add(dotakku[0]);
232 b._akku_add(dotakku[1]);
233 idotakku[0]=_idotprecision(dotakku[0],dotakku[1]);
234 l_interval tmp;
235 tmp._akku_out();
236 return tmp;
237}
238*/
239
240/*
241l_interval _unchecked_l_interval(const l_real & lr1, const l_real & lr2) noexcept
242{
243 real inf, sup, tmp; // fuer Naeherungen, entspricht Interval z
244 int i=1;
245 l_interval li;
246 dotprecision dot1(0.0);
247 dotprecision dot2(0.0);
248 lr1._akku_add(dot1);
249 lr2._akku_add(dot2);
250 inf = rnd(dot1,RND_DOWN);
251 sup = rnd(dot2,RND_UP);
252 while (i<li.prec && !(inf<=0. && sup>=0.))
253 {
254 tmp = inf + (sup-inf)/2.0; // middle(interval)
255 li.elem(i) = tmp;
256 dot1 -= tmp;
257 dot2 -= tmp;
258 inf = rnd(dot1,RND_DOWN);
259 sup = rnd(dot2,RND_UP);
260 i++;
261 }
262 li.elem(li.prec)=inf; // Intervall in die letzten Stellen
263 li.elem(li.prec+1)=sup; // schreiben
264
265 return li;
266}
267*/
268
269l_interval _unchecked_l_interval(const l_real & lr1, const l_real & lr2) noexcept
270{
271 dotprecision dot1, dot2;
272 dot1 = lr1;
273 dot2 = lr2;
274 idotprecision idot(dot1,dot2);
275 l_interval li;
276 li._akku_out(idot);
277 return li;
278}
279
280
281idotprecision _idotprecision(const l_interval & a) noexcept
282{
283 return idotprecision(a);
284}
285
286idotprecision::idotprecision(const l_interval & a) noexcept : inf(0),
287 sup(0)
288{
289 a._akku_add(*this);
290}
291
293{
294 *this=0;
295 a._akku_add(*this);
296 return *this;
297}
298
299// ---- Standardfunkt ---- (arithmetische Operatoren)
300
301l_interval operator-(const l_interval & a) noexcept
302{
303 int precsave=stagprec;
304 stagprec=a.prec;
305
306 l_interval tmp;
307
308 stagprec=precsave;
309
310 int i;
311 for(i=0;i<a.prec-1;i++)
312 tmp.data[i]=-a.data[i];
313
314 tmp.data[a.prec-1]=-a.data[a.prec];
315 tmp.data[a.prec]=-a.data[a.prec-1];
316
317 return tmp;
318}
319
320// LI-LI
321
322l_interval operator+(const l_interval & li1, const l_interval & li2) noexcept
323{
324 l_interval li3;
325 idotprecision idot(0.0);
326 li1._akku_add(idot);
327 li2._akku_add(idot);
328 li3._akku_out(idot);
329 return li3;
330}
331
332l_interval operator-(const l_interval & li1, const l_interval & li2) noexcept
333{
334 l_interval li3;
335 idotprecision idot(0.0);
336 li1._akku_add(idot);
337 li2._akku_sub(idot);
338 li3._akku_out(idot);
339 return li3;
340}
341
342/* l_interval operator*(const l_interval & li1, const l_interval & li2) noexcept
343{
344 l_interval li3;
345 interval stdmul;
346
347 stdmul = _interval(li1)*_interval(li2);
348
349 if (abs(Inf(stdmul)) < MinReal)
350 li3 = _l_interval(stdmul) + 0.0; // Was das +0.0 soll ist mir ein Raetsel SR
351 else
352 { // if eingef�gt am 30.07.92 von Frederic Toussaint
353 idotakku[0]=0.0;
354 accumulate(idotakku[0], li1, li2);
355 li3._akku_out();
356 li3 = li3 & _l_interval(stdmul); // Nachkorrektur
357 // Wie bitte? (SR)
358 }
359 return li3;
360} */
361
362l_interval operator * (const l_interval& li1, const l_interval& li2) noexcept
363{ // Blomquist, Neue Version, 21.11.02;
364 l_interval li3;
365 interval stdmul;
366 stdmul = _interval(li1) * _interval(li2); // grobe Einschlie�ung
367 idotprecision idot(0.0);
368 accumulate(idot,li1,li2);
369 li3._akku_out(idot); // li3: meist feinere Einschlie�ung! Aber bei zu
370 // breiten li1,li2 ist dieses Ergebnis li3 breiter als stdmul!
371 li3 = li3 & _l_interval(stdmul); // Das optimale Ergebnis liegt daher
372 // im Durchschnitt der groben und der feineren Einschlie�ung!
373 return li3;
374}
375
376l_interval operator/(const l_interval & li1, const l_interval & li2)
377{
378// ge�ndert am 29.07.92 von Frederic Toussaint
379// 13.5.93 AW: neuer Algorithmus nach der Beschreibung von W. Kraemer
380
381 int i, j, k, Z_sign;
382 real dzmitte;
383 interval dn = _interval(li2),
384 dz = _interval(li1),
385 stddiv = dz/dn;
386 idotprecision idot;
387 dotprecision dot1,dot2;
388 l_interval li3;
389
390 li3._clear(1);
391
392 if(!li2)
393 {
394 cxscthrow(ERROR_LINTERVAL_DIV_BY_ZERO("l_interval operator/(const l_interval & li1, const l_interval & li2)"));
395 } else
396 {
397 dz = dz/dn;
398 if (stagprec > 1)
399 { // bei == 1 passiert nix!
400 idot = 0.0; // X, also Zaehler
401 li1._akku_add(idot);
402 k = 1;
403 Z_sign = (Inf(dz) > 0.0);
404 do {
405 dzmitte = (Inf(dz) + Sup(dz)) / 2.0; // mid(interval); grob!
406 if (!!dz)
407 if (Sup(abs(dz)) > MinReal)
408 if (diam(dz) < 1e-14 * abs(dzmitte) )
409 {
410 li3.elem(k) = dzmitte;
411 // Bestimme Rest in dotakku[2], [3] (Inf, Sup):
412 //
413 dot1 = 0.0;
414 for (i=1; i<=k; i++) // Rumpf
415 for (j=1; j<=li2.prec-1; j++)
416 accumulate(dot1, -li2.elem(j), li3.elem(i));
417
418 dot2 = dot1; // Rumpf ist fuer Inf, Sup gleich
419 dot1 += Inf(idot); // Inf += Inf(X)
420 dot2 += Sup(idot); // Sup += Sup(X)
421 if (Z_sign) // INCL(Z) > 0.0
422 for (i=1; i<=k; i++)
423 {
424 accumulate(dot1, -li2.elem(li2.prec+1), li3.elem(i));
425 accumulate(dot2, -li2.elem(li2.prec), li3.elem(i));
426 } else // INCL(Z) < 0.0
427 for (i=1; i<=k; i++)
428 {
429 accumulate(dot1, -li2.elem(li2.prec), li3.elem(i));
430 accumulate(dot2, -li2.elem(li2.prec+1), li3.elem(i));
431 }
432 //
433 // Rausrunden in dz
434 dz = _interval(rnd(dot1, RND_DOWN),
435 rnd(dot2, RND_UP));
436 dz = dz / dn;
437 }
438 k++;
439 } while (k < stagprec);
440 } // if (stagprec > 1)
441 li3.elem(stagprec) = Inf(dz); // Intervall in die letzten Stelle
442 li3.elem(stagprec+1) = Sup(dz);
443 li3 = li3 & _l_interval(stddiv); // Nachkorrektur
444 } // if keine Null im Nenner
445 return li3;
446}
447
448// ---- Vergleichsop. ----
449bool operator!(const l_interval & li) noexcept
450{
451 idotprecision idot(0.0);
452 li._akku_add(idot);
453 return (!idot);
454}
455
456/*l_interval::operator void *(void) noexcept
457{
458 idotakku[0]=0.0;
459 _akku_add(idotakku[0]);
460 return (idotakku[0]);
461}*/
462
463bool operator==(const l_interval & li1, const l_interval & li2) noexcept
464{
465 idotprecision idot1(0.0);
466 idotprecision idot2(0.0);
467 li1._akku_add(idot1);
468 li2._akku_add(idot2);
469 return (idot1==idot2);
470}
471
472// ---- Mengenvergleiche ----
473
474bool operator<(const l_interval & li1, const l_interval & li2) noexcept
475{
476 idotprecision idot1(0.0);
477 idotprecision idot2(0.0);
478 li1._akku_add(idot1);
479 li2._akku_add(idot2);
480 return (idot1<idot2);
481}
482
483bool operator>(const l_interval & li1, const l_interval & li2) noexcept
484{
485 idotprecision idot1(0.0);
486 idotprecision idot2(0.0);
487 li1._akku_add(idot1);
488 li2._akku_add(idot2);
489 return (idot1>idot2);
490}
491
492bool operator<=(const l_interval & li1, const l_interval & li2) noexcept
493{
494 idotprecision idot1(0.0);
495 idotprecision idot2(0.0);
496 li1._akku_add(idot1);
497 li2._akku_add(idot2);
498 return (idot1<=idot2);
499}
500
501bool operator>=(const l_interval & li1, const l_interval & li2) noexcept
502{
503 idotprecision idot1(0.0);
504 idotprecision idot2(0.0);
505 li1._akku_add(idot1);
506 li2._akku_add(idot2);
507 return (idot1>=idot2);
508}
509
510void ConvexHull(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4) noexcept
511{
512 if(li1<=li2)
513 { // Trivialfall 1
514 li3=li2;
515 li4=li2;
516 } else if(li2<=li1)
517 { // Trivialfall 2
518 li3=li1;
519 li4=li1;
520 } else
521 { // rechne
522 idotprecision idot1(0.0);
523 idotprecision idot2(0.0);
524 li1._akku_add(idot1);
525 li2._akku_add(idot2);
526
527 idot1 |= idot2;
528 // nun steht das richtige Ergebnis in idotakku[0]
529 idot2 = idot1; // sichern
530 li3._akku_out_inn(idot1); // Rundung nach innen
531 idot1 = idot2; // und wiederherstellen
532 li4._akku_out(idot1); // ... nach aussen
533 }
534}
535
536void Intersection(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4)
537{
538 if(li1<=li2)
539 { // Trivialfall 1
540 li3=li1;
541 li4=li1;
542 } else if(li2<=li1)
543 { // Trivialfall 2
544 li3=li2;
545 li4=li2;
546 } else
547 { // rechne
548 idotprecision idot1(0.0);
549 idotprecision idot2(0.0);
550 idotprecision idot;
551 li1._akku_add(idot1);
552 li2._akku_add(idot2);
553 try
554 {
555 idot=(idot1 &= idot2);
556 }
557 catch(const EMPTY_INTERVAL &)
558 {
559 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL("void Intersection(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4)"));
560 }
561 // nun steht das richtige Ergebnis in idotakku[0]
562 li3._akku_out_inn(idot); // Rundung nach innen
563 idot = idot1;
564 li4._akku_out(idot); // ... nach aussen
565 }
566}
567
568l_real mid (const l_interval & li) noexcept
569{
570 l_real lr;
571
572 // --------------------------------------------------------
573 // dotakku[0] = Inf(li) + Sup(li)
574
575 dotprecision dot(0.0);
576 for (int j=1; j<=li.prec-1; j++)
577 dot += li.elem(j);
578
579 dot += dot; // mal 2
580 dot += li.elem(li.prec);
581 dot += li.elem(li.prec+1);
582
583 // ------------------------------------------------------
584 // Division nur bei ungleich 0
585
586 if (dot != 0.0)
587 {
588 Dotprecision ptr = *dot.ptr();
589
590 // --------------------------------------------------------
591 // Dividiere dotakku[0] durch 2, mittels 1 Rechtsshift
592
593 ptr[(a_intg)++ptr[A_END]] = ZERO;
594
595 b_shr1 (&ptr[(a_intg)ptr[A_BEGIN]], (a_intg)(ptr[A_END]-ptr[A_BEGIN]+1));
596
597 if (ptr[(a_intg)ptr[A_END]] == ZERO)
598 --ptr[A_END];
599 if (ptr[(a_intg)ptr[A_BEGIN]] == ZERO)
600 ++ptr[A_BEGIN];
601
602 // --------------------------------------------------------
603 }
604
605 lr._akku_out(dot);
606
607 return lr;
608}
609
610/* void accumulate(idotprecision & d, const l_interval & li1, const l_interval & li2) noexcept
611{
612 // �nderungen am 24.9.92 von F. Toussaint wegen Underflow-Fehlern
613
614 int i, j, n = 0;
615 real tmp;
616
617 for (i=1; i<=li1.prec-1; i++)
618 for (j=1; j<=li2.prec-1; j++)
619 {
620 tmp = abs(li1.elem(i)*li2.elem(j));
621 if (tmp < MinReal && tmp != CXSC_Zero)
622 n++;
623 else
624 accumulate(d, _interval(li1.elem(i)), _interval(li2.elem(j)));
625 }
626 for (i=1; i<=li2.prec-1; i++)
627 {
628 tmp = abs(Inf(_interval(li1.elem(li1.prec), li1.elem(li1.prec+1)) * _interval(li2.elem(i))));
629 if (tmp < MinReal && tmp != CXSC_Zero)
630 n++;
631 else
632 accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
633 _interval(li2.elem(i)));
634 }
635 for (i=1; i<=li1.prec-1; i++)
636 {
637 tmp = abs(Inf(_interval(li2.elem(li2.prec), li2.elem(li2.prec+1)) * _interval(li1.elem(i))));
638 if (tmp < MinReal && tmp != CXSC_Zero)
639 n++;
640 else
641 accumulate(d, _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)),
642 _interval(li1.elem(i)));
643 }
644 tmp = abs(Inf(_interval(li1.elem(li1.prec), li1.elem(li1.prec+1)) *
645 _interval(li2.elem(li2.prec), li2.elem(li2.prec+1))));
646 if (tmp < MinReal && tmp != CXSC_Zero)
647 n++;
648 else
649 accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
650 _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)));
651 if (n > 0)
652 accumulate(d, _interval(_real(n)),
653 _interval(-MinReal, MinReal));
654} */
655
656void accumulate(idotprecision & d, const l_interval & li1,
657 const l_interval & li2) noexcept
658{ // Blomquist, Neue Version vom 21.11.02;
659 // Die �nderungen von Toussaint wurden r�ckg�ngig gemacht.
660 int i,j;
661
662 for (i=1; i<=li1.prec-1; i++)
663 for (j=1; j<=li2.prec-1; j++)
664 {
665 accumulate(d, _interval(li1.elem(i)), _interval(li2.elem(j)));
666 }
667 for (i=1; i<=li2.prec-1; i++)
668 {
669 accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
670 _interval(li2.elem(i)));
671 }
672 for (i=1; i<=li1.prec-1; i++)
673 {
674 accumulate(d, _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)),
675 _interval(li1.elem(i)));
676 }
677 accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
678 _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)));
679}
680
681
682/* void l_interval::_akku_out() noexcept
683{
684 // ein im idotakku[0] liegendes Zwischenergebnis
685 // wird in der entsprechenden Praezision in das aufrufende l_interval
686 // gerundet. Rundung, also Einschluss von aussen!!
687 // Alle MinReal wurden ge�ndert in minreal, Blomquist 19.11.02;
688 _clear(1);
689 interval z;
690 real tmp, r, s;
691 int i = 1,
692
693 n = 0;
694 z = rnd(idotakku[0]);
695
696 while (i<prec && !(CXSC_Zero <= z))
697 {
698 // z ist Intervall, deshalb <=
699 if (succ(succ(Inf(z))) < Sup(z))
700 break; // bei zu grobem Einschluss:
701 // Intervall direkt in
702 // Einschlusskomponente
703 // schreiben
704 tmp = Inf(z) + (Sup(z)-Inf(z))/2.0; // middle(interval)
705 if (abs(Inf(z)) < minreal)
706 {
707 if (abs(Sup(z)) < minreal)
708 {
709 if (tmp != 0.0)
710 {
711 tmp = 0.0;
712 n++;
713 }
714 i = prec;
715 }
716 }
717
718 this->elem(i) = tmp;
719 idotakku[0] -= tmp;
720 z = rnd(idotakku[0]);
721 i++;
722 } // while
723 r = Inf(z);
724 if (abs(r) < minreal)
725 {
726 if (r < 0.0)
727 r = -minreal;
728 else
729 r = 0.0;
730 }
731 s = Sup(z);
732 if (abs(s) < minreal)
733 {
734 if (s <= 0.0)
735 s = 0.0;
736 else
737 s = minreal;
738 }
739 if (n > 0)
740 {
741 this->elem(prec) = r-_real(n+1)*minreal; // Intervall in die letzten
742 this->elem(prec+1) = s+_real(n+1)*minreal; // Stellen schreiben
743 // (n+1), da Rundungsrichtung nicht bestimmt werden kann!
744 } else
745 {
746 this->elem(prec) = r; // Intervall in die letzten Stellen
747 this->elem(prec+1) = s; // schreiben
748 }
749} */
750
751void l_interval::_akku_out(idotprecision& idot) noexcept
752{
753 // ein im idotakku[0] liegendes Zwischenergebnis
754 // wird in der entsprechenden Praezision in das aufrufende l_interval
755 // so gerundet, dass idotakku[0] eingeschlossen wird.
756 // Neueste Version: Blomquist 21.07.03;
757 _clear(1);
758 interval z;
759 real tmp;
760 int i = 1;
761
762 z = rnd(idot);
763 while (i<prec && !(CXSC_Zero <= z))
764 {
765 // z ist Intervall, deshalb <=
766 if (succ(succ(Inf(z))) < Sup(z))
767 break; // bei zu grobem Einschluss: Intervall direkt in
768 // Einschlusskomponente schreiben.
769 tmp = mid(z); // middle(interval)
770 this->elem(i) = tmp;
771 idot -= tmp;
772 z = rnd(idot);
773 i++;
774 } // while
775 this->elem(prec) = Inf(z); // Intervall in die letzten Stellen
776 this->elem(prec+1) = Sup(z); // schreiben
777} // _akku_out()
778
779void l_interval::_akku_out_inn(idotprecision& idot) noexcept
780{
781 // ein im idotakku[0] liegendes Zwischenergebnis
782 // wird in der entsprechenden Praezision in die aufrufende l_interval Zahl
783 // li so gerundet, dass gilt: li enthalten in idotakku[0].
784 _clear(1);
785 real inf, sup, tmp; // fuer Naeherungen, entspricht Interval z
786 int i=1;
787 inf = rnd(Inf(idot),RND_UP);
788 sup = rnd(Sup(idot),RND_DOWN);
789
790 if (inf>sup)
791 inf = sup; // Vorsichtsmassnahme
792
793 while (i<prec && !(inf<=CXSC_Zero&&sup>=CXSC_Zero))
794 {
795 tmp = inf + (sup-inf)/2.0; // middle(interval)
796 this->elem(i) = tmp;
797 idot -= tmp;
798 inf = rnd(Inf(idot),RND_UP);
799 sup = rnd(Sup(idot),RND_DOWN);
800 if (inf>sup)
801 inf = sup; // Vorsichtsmassnahme
802 i++;
803 }
804 this->elem(prec)=inf; // Intervall in die letzten Stellen
805 this->elem(prec+1)=sup; // schreiben
806}
807
808/* void l_interval::_akku_add(idotprecision& d) const noexcept
809{
810 // addiert aufrufenden l_interval auf iakku d.
811 // �nderung am 23.9.92 von F. Toussaint, da Fehler im Unterlaufbereich
812 // Ausgeblendet von Blomquist, 20.11.02, da Fehler im Unterlaufbereich
813 // nicht erkennbar. Meine neue Version direkt anschliessend!
814 int n = 0;
815 real r, s;
816
817 for (int i=1; i<=prec-1; i++)
818 {
819 r = this->elem(i);
820 if (abs(r) < MinReal)
821 {
822 if (r != 0.0)
823 n++;
824 } else
825 d += _interval(r);
826 }
827 r = this->elem(prec);
828 if (abs(r) < MinReal)
829 {
830 if (r < 0.0)
831 r = -MinReal;
832 else
833 r = 0.0;
834 }
835 s = this->elem(prec+1);
836 if (abs(s) < MinReal)
837 {
838 if (s <= 0.0)
839 s = 0.0;
840 else
841 s = MinReal;
842 }
843 if (r != CXSC_Zero || s != CXSC_Zero)
844 d += _interval(r, s);
845 if (n > 0)
846 {
847 d += _interval(-_real(n+1)*MinReal, _real(n+1)*MinReal);
848 // (n+1), da Rundungsrichtung nicht bestimmt werden kann!
849 }
850} */
851
852void l_interval::_akku_add(idotprecision& d) const noexcept
853{
854 // Addiert aufrufendes Intervall vom Typ l_interval auf d.
855 // Meine neue Version; Blomquist, 20.11.02;
856 real r,s;
857 for (int i=1; i<=prec-1; i++)
858 {
859 r = this->elem(i);
860 if (sign(r) != CXSC_Zero) // Addition nur, falls n�tig!
861 d += _interval(r);
862 }
863 r = this->elem(prec);
864 s = this->elem(prec+1);
865 if (r != CXSC_Zero || s != CXSC_Zero) // Addition nur, falls n�tig!
866 d += _interval(r, s);
867}
868
869/* void l_interval::_akku_sub(idotprecision& d) const noexcept
870{
871 // Subtrahiert aufrufendes Intervall vom Typ l_interval von d.
872 // Intervallsubtraktion!!
873 // �nderung am 23.9.92 von F. Toussaint, da Fehler im Unterlaufbereich
874 // Blomquist: Mir sind keine solchen Fehler bekannt, daher wurde diese
875 // Version von mir ausgeklammert, 20.11.02; Neue Version anschliessend!
876
877 int n = 0;
878 real r, s;
879
880 for (int i=1; i<=prec-1; i++)
881 {
882 r = this->elem(i);
883 if (abs(r) < MinReal)
884 {
885 if (r != 0.0)
886 n++;
887 } else
888 d -= _interval(r);
889 }
890 r = this->elem(prec);
891 if (abs(r) < MinReal)
892 {
893 if (r < 0.0)
894 r = -MinReal;
895 else
896 r = 0.0;
897 }
898 s = this->elem(prec+1);
899 if (abs(s) < MinReal)
900 {
901 if (s <= 0.0)
902 s = 0.0;
903 else
904 s = MinReal;
905 }
906 d -= _interval(r, s);
907 if (n > 0)
908 {
909 d -= _interval(-_real(n+1)*MinReal, _real(n+1)*MinReal);
910 // (n+1), da Rundungsrichtung nicht bestimmt werden kann!
911 }
912} */
913
914void l_interval::_akku_sub(idotprecision& d) const noexcept
915{
916 // Subtrahiert aufrufendes Intervall vom Typ l_interval von d.
917 // Meine neue Version; Blomquist, 20.11.02;
918 real r,s;
919
920 for (int i=1; i<=prec-1; i++)
921 {
922 r = this->elem(i);
923 if (sign(r) != CXSC_Zero) // Subtraktion nur, falls n�tig!
924 d -= _interval(r);
925 }
926 r = this->elem(prec);
927 s = this->elem(prec+1);
928 if (r != CXSC_Zero || s != CXSC_Zero) // Subtraktion nur, falls n�tig!
929 d -= _interval(r, s);
930}
931
932// ---- Ausgabefunkt. ---------------------------------------
933
934std::ostream & operator << (std::ostream &s, const l_interval & a) noexcept
935{
936 idotprecision idot(0.0);
937 a._akku_add(idot);
938 s << idot;
939 return s;
940}
941std::string & operator << (std::string &s, const l_interval & a) noexcept
942{
943 idotprecision idot(0.0);
944 a._akku_add(idot);
945 s << idot;
946 return s;
947}
948
949std::istream & operator >> (std::istream &s, l_interval & a) noexcept
950{
951 idotprecision idot;
952 s >> idot;
953 a._akku_out(idot);
954 return s;
955}
956
957std::string & operator >> (std::string &s, l_interval & a) noexcept
958{
959 idotprecision idot;
960 s >> idot;
961 a._akku_out(idot);
962 return s;
963
964}
965
966void operator >>(const std::string &s,l_interval & a) noexcept
967{
968 std::string r(s);
969 idotprecision idot;
970 r >> idot;
971 a._akku_out(idot);
972}
973void operator >>(const char *s,l_interval & a) noexcept
974{
975 std::string r(s);
976 idotprecision idot;
977 r>>idot;
978 a._akku_out(idot);
979}
980
981int in ( const real& x, const l_interval& y ) // Contained-in relation
982 { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //----------------------
983
984int in ( const l_real& x, const l_interval& y ) // Contained-in relation
985 { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //----------------------
986
987int in ( const interval& x, const l_interval& y ) // Contained-in relation
988{ //----------------------
989 return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
990}
991
992int in ( const l_interval& x, const l_interval& y ) // Contained-in relation
993{ //----------------------
994 return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
995}
996
997l_interval Blow (const l_interval& x, const real& eps )
998{
999 int n;
1000 l_interval y;
1001 l_real lr1,lr2;
1002 y = x + interval(-eps,eps) * diam(x);
1003 lr1 = Inf(y);
1004 n = StagPrec(lr1);
1005 lr1[n] = pred(lr1[n]);
1006 lr2 = Sup(y);
1007 lr2[n] = succ(lr2[n]);
1008 return l_interval(lr1,lr2);
1009}
1010
1011int Disjoint (const l_interval& a, const l_interval& b )
1012 // Test for disjointedness
1013{ //------------------------
1014 return (Inf(a) > Sup(b)) || (Inf(b) > Sup(a));
1015}
1016
1017l_real AbsMin ( const l_interval& x ) // Absolute minimum of
1018{ // an interval
1019 if ( in(0.0,x) ) return l_real(0.0);
1020 else
1021 {
1022 l_real y(Inf(x));
1023 if (y > 0.0) return y;
1024 else return -Sup(x);
1025 }
1026
1027}
1028
1029l_real AbsMax (const l_interval& x ) // Absolute maximum of
1030{ // an interval
1031 l_real a, b; //--------------------
1032
1033 a = abs(Inf(x));
1034 b = abs(Sup(x));
1035
1036 if (a > b)
1037 return a;
1038 else
1039 return b;
1040}
1041
1042l_real RelDiam ( const l_interval x ) // Relative diameter
1043{ // of an interval
1044 if ( in(0.0,x) ) //------------------
1045 return diam(x);
1046 else
1047 return( Sup( l_interval(diam(x)) / l_interval(AbsMin(x)) ) );
1048}
1049
1050inline void times2pown(l_interval& x, int n) noexcept
1051{ // x = x * 2^n; Blomquist, 28.11.02;
1052 real mt,t;
1053 interval z;
1054 int p = StagPrec(x);
1055 if ( n<LI_Min_Exp_ || n>LI_maxexpo1 )
1056 { std::cerr << "Error in: "
1057 << "times2pown(l_interval& x, const int& n): " << std::endl
1058 << " -1074 <= n <= +1023 not fulfilled" << std::endl; exit(0);
1059 }
1060 real F = comp(0.5,n+1);
1061 z = _interval(x[p],x[p+1]);
1062 times2pown(z,n); // Scaling the original interval;
1063
1064 for (int i=1; i<=p-1; i++)
1065 {
1066 mt = mant(x[i]);
1067 t = x[i];
1068 times2pown(x[i],n);
1069 if ( mt != mant(x[i]) )
1070 {
1071 x[i] = 0;
1072 z += _interval(t) * F;
1073 }
1074 }
1075 x[p] = Inf(z);
1076 x[p+1] = Sup(z);
1077} // times2pown(...)
1078
1079
1080void Times2pown(l_interval& a, const real& p) noexcept
1081// The first parameter delivers an inclusion of a * 2^p;
1082// For p in [-2100,+2100] p must be an integer value.
1083// This condition is NOT tested in this function!
1084// For p outside [-2100,+2100] an inclusion of a * 2^p is
1085// calculated for any p of type real, unless an overflow occurs.
1086// If the function is activated with the second parameter of type int,
1087// then the first parameter delivers correct inclusions of a * 2^p,
1088// unless an overflow occurs.
1089// Blomquist, 28.01.2008;
1090{
1091 const int c2 = 2100;
1092 int ex(expo_gr(a)),fac,rest,n;
1093 double dbl;
1094
1095 if (ex > -1000000)
1096 {
1097 if (p>=0)
1098 if (p>c2)
1099 times2pown(a,c2); // Produces an error
1100 else // 0 <= p <= 2100
1101 {
1102 dbl = _double(p);
1103 n = (int) dbl;
1104 fac = n/LI_maxexpo1;
1105 rest = n%LI_maxexpo1;
1106 for (int k=1; k<=fac; k++)
1107 times2pown(a,LI_maxexpo1);
1108 times2pown(a,rest);
1109 }
1110 else // p<0; No overflow or underflow!
1111 if (p<-c2)
1112 {
1113 if (0<a)
1115 else
1116 if (Inf(a)>=0)
1117 a = l_interval(0,minreal);
1118 else a = l_interval(-minreal,0);
1119 }
1120 else // -2100 <= p < 0
1121 {
1122 dbl = _double(p);
1123 n = (int) dbl;
1124 fac = n/LI_Min_Exp_;
1125 rest = n%LI_Min_Exp_;
1126 for (int k=1; k<=fac; k++)
1127 times2pown(a,LI_Min_Exp_);
1128 times2pown(a,rest);
1129 }
1130 }
1131} // Times2pown(...)
1132
1133} // namespace cxsc
1134
The Data Type dotprecision.
Definition dot.hpp:112
The Data Type idotprecision.
Definition idot.hpp:48
idotprecision & operator=(const real &a)
Implementation of standard assigning operator.
Definition idot.hpp:96
idotprecision()
Constructor of class idotprecision.
Definition idot.hpp:57
The Scalar Type interval.
Definition interval.hpp:55
interval & operator=(const real &a)
Implementation of standard assigning operator.
Definition interval.inl:52
interval()
Constructor of class interval.
Definition interval.hpp:64
The Multiple-Precision Data Type l_interval.
l_interval() noexcept
Constructor of class l_interval.
The Multiple-Precision Data Type l_real.
Definition l_real.hpp:78
The Scalar Type real.
Definition real.hpp:114
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
const real MinReal
Smallest normalized representable floating-point number.
Definition real.cpp:62
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
Definition interval.cpp:288
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
Definition cimatrix.inl:730
int expo_gr(const l_interval &x)
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
Definition cimatrix.inl:738
const real minreal
Smallest positive denormalized representable floating-point number.
Definition real.cpp:63
l_interval _l_interval(const real &a) noexcept
real RelDiam(const interval &x)
Computes the relative diameter .
Definition interval.cpp:316
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
idotprecision _idotprecision(const real &a)
Definition idot.inl:51
real AbsMax(const interval &x)
Computes the greatest absolute value .
Definition interval.cpp:303
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition cimath.cpp:2059
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
Definition cimatrix.inl:731
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.
real AbsMin(const interval &x)
Computes the smallest absolute value .
Definition interval.cpp:293