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 
31 namespace cxsc {
32 
33 #define CXSC_Zero 0.
34 
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 
49 interval::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 
64 interval _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 
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 
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 
162 l_interval::l_interval(const real &a, const l_real &b)
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 
179 void 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 
205 l_interval::l_interval(const l_real &a, const real &b)
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 /*
224 l_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 /*
241 l_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 
269 l_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 
281 idotprecision _idotprecision(const l_interval & a) noexcept
282 {
283  return idotprecision(a);
284 }
285 
286 idotprecision::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 
301 l_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 
322 l_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 
332 l_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 
362 l_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 
376 l_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. ----
449 bool 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 
463 bool 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 
474 bool 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 
483 bool 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 
492 bool 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 
501 bool 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 
510 void 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 
536 void 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 
568 l_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 
656 void 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 
751 void 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 
779 void 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 
852 void 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 
914 void 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 
934 std::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 }
941 std::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 
949 std::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 
957 std::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 
966 void 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 }
973 void 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 
981 int in ( const real& x, const l_interval& y ) // Contained-in relation
982  { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //----------------------
983 
984 int in ( const l_real& x, const l_interval& y ) // Contained-in relation
985  { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //----------------------
986 
987 int in ( const interval& x, const l_interval& y ) // Contained-in relation
988 { //----------------------
989  return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
990 }
991 
992 int 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 
997 l_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 
1011 int 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 
1017 l_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 
1029 l_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 
1042 l_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 
1050 inline 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 
1080 void 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)
1114  a = l_interval(-minreal,minreal);
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.
Definition: l_interval.hpp:72
l_interval() noexcept
Constructor of class l_interval.
Definition: l_interval.inl:45
l_interval & operator=(const real &a) noexcept
Implementation of standard assigning operator.
Definition: l_interval.inl:89
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)
Definition: l_interval.inl:522
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
Definition: l_interval.hpp:879
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.
Definition: cinterval.cpp:654
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.
Definition: cinterval.cpp:665
real AbsMin(const interval &x)
Computes the smallest absolute value .
Definition: interval.cpp:293