27#include "l_interval.hpp"
39 for(i=0;i<=a.prec;i++)
60 return *
this=rnd(idot);
91 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL(
"l_interval::l_interval(const dotprecision&,const dotprecision&)"));
94 UncheckedSetInf(idot,a);
95 UncheckedSetSup(idot,b);
100#if (CXSC_INDEX_CHECK)
117#if (CXSC_INDEX_CHECK)
129#if (CXSC_INDEX_CHECK)
146#if (CXSC_INDEX_CHECK)
154 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL(
"l_interval::l_interval(const l_real &a, const l_real &b)"));
163#if (CXSC_INDEX_CHECK)
171 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL(
"l_interval::l_interval(const l_real &a, const l_real &b)"));
184 int p = StagPrec(lr);
int q = StagPrec(li);
187 std::cerr <<
"l_realz2l_interval(const l_real& lr, const interval& z, l_interval& li): incorrect precisions of lr,li !"
192 for (
int i=1; i<=q-1; i++) li[i] = 0;
195 if(p<q)
for (
int i=1; i<=p; i++) li[i] = lr[i];
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));
206#if (CXSC_INDEX_CHECK)
214 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL(
"l_interval::l_interval(const l_real &a, const l_real &b)"));
283 return idotprecision(a);
303 int precsave=stagprec;
311 for(i=0;i<a.prec-1;i++)
312 tmp.data[i]=-a.data[i];
314 tmp.data[a.prec-1]=-a.data[a.prec];
315 tmp.data[a.prec]=-a.data[a.prec-1];
366 stdmul = _interval(li1) * _interval(li2);
368 accumulate(idot,li1,li2);
394 cxscthrow(ERROR_LINTERVAL_DIV_BY_ZERO(
"l_interval operator/(const l_interval & li1, const l_interval & li2)"));
403 Z_sign = (Inf(dz) > 0.0);
405 dzmitte = (Inf(dz) + Sup(dz)) / 2.0;
408 if (
diam(dz) < 1e-14 *
abs(dzmitte) )
410 li3.elem(k) = dzmitte;
415 for (j=1; j<=li2.prec-1; j++)
416 accumulate(dot1, -li2.elem(j), li3.elem(i));
424 accumulate(dot1, -li2.elem(li2.prec+1), li3.elem(i));
425 accumulate(dot2, -li2.elem(li2.prec), li3.elem(i));
429 accumulate(dot1, -li2.elem(li2.prec), li3.elem(i));
430 accumulate(dot2, -li2.elem(li2.prec+1), li3.elem(i));
434 dz = _interval(rnd(dot1, RND_DOWN),
439 }
while (k < stagprec);
441 li3.elem(stagprec) = Inf(dz);
442 li3.elem(stagprec+1) = Sup(dz);
467 li1._akku_add(idot1);
468 li2._akku_add(idot2);
469 return (idot1==idot2);
478 li1._akku_add(idot1);
479 li2._akku_add(idot2);
480 return (idot1<idot2);
487 li1._akku_add(idot1);
488 li2._akku_add(idot2);
489 return (idot1>idot2);
496 li1._akku_add(idot1);
497 li2._akku_add(idot2);
498 return (idot1<=idot2);
505 li1._akku_add(idot1);
506 li2._akku_add(idot2);
507 return (idot1>=idot2);
524 li1._akku_add(idot1);
525 li2._akku_add(idot2);
530 li3._akku_out_inn(idot1);
532 li4._akku_out(idot1);
551 li1._akku_add(idot1);
552 li2._akku_add(idot2);
555 idot=(idot1 &= idot2);
557 catch(
const EMPTY_INTERVAL &)
559 cxscthrow(ERROR_LINTERVAL_EMPTY_INTERVAL(
"void Intersection(const l_interval & li1, const l_interval & li2, l_interval & li3, l_interval & li4)"));
562 li3._akku_out_inn(idot);
576 for (
int j=1; j<=li.prec-1; j++)
580 dot += li.elem(li.prec);
581 dot += li.elem(li.prec+1);
588 Dotprecision ptr = *dot.ptr();
593 ptr[(a_intg)++ptr[A_END]] = ZERO;
595 b_shr1 (&ptr[(a_intg)ptr[A_BEGIN]], (a_intg)(ptr[A_END]-ptr[A_BEGIN]+1));
597 if (ptr[(a_intg)ptr[A_END]] == ZERO)
599 if (ptr[(a_intg)ptr[A_BEGIN]] == ZERO)
662 for (i=1; i<=li1.prec-1; i++)
663 for (j=1; j<=li2.prec-1; j++)
665 accumulate(d, _interval(li1.elem(i)), _interval(li2.elem(j)));
667 for (i=1; i<=li2.prec-1; i++)
669 accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
670 _interval(li2.elem(i)));
672 for (i=1; i<=li1.prec-1; i++)
674 accumulate(d, _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)),
675 _interval(li1.elem(i)));
677 accumulate(d, _interval(li1.elem(li1.prec), li1.elem(li1.prec+1)),
678 _interval(li2.elem(li2.prec), li2.elem(li2.prec+1)));
763 while (i<prec && !(CXSC_Zero <= z))
766 if (succ(succ(Inf(z))) < Sup(z))
775 this->elem(prec) = Inf(z);
776 this->elem(prec+1) = Sup(z);
779void l_interval::_akku_out_inn(idotprecision& idot)
noexcept
787 inf = rnd(Inf(idot),RND_UP);
788 sup = rnd(Sup(idot),RND_DOWN);
793 while (i<prec && !(inf<=CXSC_Zero&&sup>=CXSC_Zero))
795 tmp = inf + (sup-inf)/2.0;
798 inf = rnd(Inf(idot),RND_UP);
799 sup = rnd(Sup(idot),RND_DOWN);
804 this->elem(prec)=inf;
805 this->elem(prec+1)=sup;
852void l_interval::_akku_add(idotprecision& d)
const noexcept
857 for (
int i=1; i<=prec-1; i++)
860 if (sign(r) != CXSC_Zero)
863 r = this->elem(prec);
864 s = this->elem(prec+1);
865 if (r != CXSC_Zero || s != CXSC_Zero)
866 d += _interval(r, s);
914void l_interval::_akku_sub(idotprecision& d)
const noexcept
920 for (
int i=1; i<=prec-1; i++)
923 if (sign(r) != CXSC_Zero)
926 r = this->elem(prec);
927 s = this->elem(prec+1);
928 if (r != CXSC_Zero || s != CXSC_Zero)
929 d -= _interval(r, s);
934std::ostream & operator << (std::ostream &s,
const l_interval & a)
noexcept
941std::string & operator << (std::string &s,
const l_interval & a)
noexcept
949std::istream & operator >> (std::istream &s,
l_interval & a)
noexcept
957std::string & operator >> (std::string &s,
l_interval & a)
noexcept
966void operator >>(
const std::string &s,
l_interval & a)
noexcept
982 {
return ( (Inf(y) <= x) && (x <= Sup(y)) ); }
984int in (
const l_real& x,
const l_interval& y )
985 {
return ( (Inf(y) <= x) && (x <= Sup(y)) ); }
987int in (
const interval& x,
const l_interval& y )
989 return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
992int in (
const l_interval& x,
const l_interval& y )
994 return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
997l_interval
Blow (
const l_interval& x,
const real& eps )
1002 y = x + interval(-eps,eps) *
diam(x);
1005 lr1[n] = pred(lr1[n]);
1007 lr2[n] = succ(lr2[n]);
1008 return l_interval(lr1,lr2);
1011int Disjoint (
const l_interval& a,
const l_interval& b )
1014 return (Inf(a) > Sup(b)) || (Inf(b) > Sup(a));
1017l_real
AbsMin (
const l_interval& x )
1019 if (
in(0.0,x) )
return l_real(0.0);
1023 if (y > 0.0)
return y;
1024 else return -Sup(x);
1029l_real
AbsMax (
const l_interval& x )
1042l_real
RelDiam (
const l_interval x )
1047 return( Sup( l_interval(
diam(x)) / l_interval(
AbsMin(x)) ) );
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);
1060 real F = comp(0.5,n+1);
1061 z = _interval(x[p],x[p+1]);
1064 for (
int i=1; i<=p-1; i++)
1069 if ( mt != mant(x[i]) )
1072 z += _interval(t) * F;
1091 const int c2 = 2100;
1092 int ex(
expo_gr(a)),fac,rest,n;
1104 fac = n/LI_maxexpo1;
1105 rest = n%LI_maxexpo1;
1106 for (
int k=1; k<=fac; k++)
1124 fac = n/LI_Min_Exp_;
1125 rest = n%LI_Min_Exp_;
1126 for (
int k=1; k<=fac; k++)
The Data Type dotprecision.
The Data Type idotprecision.
idotprecision & operator=(const real &a)
Implementation of standard assigning operator.
idotprecision()
Constructor of class idotprecision.
The Scalar Type interval.
interval & operator=(const real &a)
Implementation of standard assigning operator.
interval()
Constructor of class interval.
The Multiple-Precision Data Type l_interval.
l_interval() noexcept
Constructor of class l_interval.
The Multiple-Precision Data Type l_real.
The namespace cxsc, providing all functionality of the class library C-XSC.
const real MinReal
Smallest normalized representable floating-point number.
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
int expo_gr(const l_interval &x)
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
const real minreal
Smallest positive denormalized representable floating-point number.
l_interval _l_interval(const real &a) noexcept
real RelDiam(const interval &x)
Computes the relative diameter .
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
idotprecision _idotprecision(const real &a)
real AbsMax(const interval &x)
Computes the greatest absolute value .
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.
real AbsMin(const interval &x)
Computes the smallest absolute value .