26#include "l_cinterval.hpp"
56 noexcept : re(Re(a)),im(Im(a)) {}
66 accumulate(akku,a.re,b.re);
67 accumulate(akku,-a.im,b.im);
70 { std::cerr <<
"Error in l_cinterval * l_cinterval" << std::endl;
74 accumulate(akku,a.im,b.re);
75 accumulate(akku,a.re,b.im);
78 { std::cerr <<
"Error in l_cinterval * l_cinterval" << std::endl;
97static const int max_expo = 1020, max_expo1 = 1023;
101bool cxsc_l_complex_division_p[5];
111 product(a, c, b, d, ex1, z);
112 product(c, d, ex2, n);
113 return quotient(z, n, round, ex1, ex2);
119static l_real minmax(
int minimum, l_real a, l_real b, l_real y0,
120 l_interval x,
int i,
int j)
127 l_interval t,q,x0,two_Da;
131 a += 0.0; b += 0.0; y0 += 0.0;
138 if (Inf(x) == Sup(x))
140 if (cxsc_l_complex_division_p[i] && cxsc_l_complex_division_p[j])
141 minmax = cxsc_complex_division_f( a, b, Inf(x), y0, 1-2*minimum );
143 cxsc_l_complex_division_p[i] =
false;
144 cxsc_l_complex_division_p[j] =
false;
148 if ( b == CXSC_Zero || y0 == CXSC_Zero )
151 cxsc_l_complex_division_p[i] =
false;
152 cxsc_l_complex_division_p[j] =
false;
156 if (minimum && sign(b) != sign(y0) )
159 minmax = Inf(l_interval(b)/y0);
160 cxsc_l_complex_division_p[i] =
false;
161 cxsc_l_complex_division_p[j] =
false;
163 if (!minimum && sign(b) == sign(y0) )
166 minmax = Sup(l_interval(b)/y0);
167 cxsc_l_complex_division_p[i] =
false;
168 cxsc_l_complex_division_p[j] =
false;
181 minmax = Inf(l_interval(a)/l_interval(Sup(x)));
184 minmax = Inf(l_interval(a)/l_interval(Inf(x)));
189 minmax = Sup(l_interval(a)/l_interval(Inf(x)));
192 minmax = Sup(l_interval(a)/l_interval(Sup(x)));
194 cxsc_l_complex_division_p[i] =
false;
195 cxsc_l_complex_division_p[j] =
false;
202 l_real invf2(1.0), a_skal;
207 if (sign(b)==0) t = 1.0;
212 int expo_diff = expo(b[1]) - expo(a[1]), ex;
214 if (expo_diff > max_expo)
216 Da = expo_diff-max_expo;
230 two_Da = l_interval( comp(0.5,-1021) );
233 else two_Da = l_interval( comp(0.5,1-Da) );
236 q = l_interval(b)/a_skal;
237 if (sign(q[1])<0) q = -q;
255 if ( (sign(b) == sign(y0)) == minimum )
257 if (expo(ay0[1]) + expo(t[1]) + Da > max_expo1)
goto Ende;
259 if (Da>0) Times2pown(x0,Da);
263 if (expo(ay0[1]) - expo(t[1]) - Da > max_expo1)
goto Ende;
265 if (Da>0) Times2pown(x0,-Da);
268 if (minimum) x0 = -x0;
276 if (minimum) minmax = Inf(q);
277 else minmax = Sup(q);
279 cxsc_l_complex_division_p[i] =
false;
280 cxsc_l_complex_division_p[j] =
false;
288l_real max(
const l_real& u,
const l_real& v)
295l_real min(
const l_real& u,
const l_real& v)
302l_cinterval cidiv(
const l_cinterval& A,
const l_cinterval& B)
304 l_real realteilINF, realteilSUP,
305 imagteilINF, imagteilSUP;
308 bool a_repeat,b_repeat;
310 l_real AREINF, ARESUP, AIMINF, AIMSUP,
311 BREINF, BRESUP, BIMINF, BIMSUP;
312 l_interval ARE, AIM, BRE, BIM;
321 AREINF = Inf(Re(A)); AREINF += 0.0;
322 ARESUP = Sup(Re(A)); ARESUP += 0.0;
323 AIMINF = Inf(Im(A)); AIMINF += 0.0;
324 AIMSUP = Sup(Im(A)); AIMSUP += 0.0;
325 BREINF = Inf(Re(B)); BREINF += 0.0;
326 BRESUP = Sup(Re(B)); BRESUP += 0.0;
327 BIMINF = Inf(Im(B)); BIMINF += 0.0;
328 BIMSUP = Sup(Im(B)); BIMSUP += 0.0;
334 a_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
335 b_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
337 if (a_repeat || b_repeat)
354 for (j=1; j<=rep; j++)
356 for (i=1; i<=4; cxsc_l_complex_division_p[i++] =
true);
360 max( max( minmax(
false, a0, b0, BIMINF, BRE, 1,2 ),
361 minmax(
false, a0, b0, BIMSUP, BRE, 3,4 ) ),
362 max( minmax(
false, b0, a0, BREINF, BIM, 1,3 ),
363 minmax(
false, b0, a0, BRESUP, BIM, 2,4 ) ) )
367 if (cxsc_l_complex_division_p[1])
368 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BREINF, BIMINF, +1 ) );
369 if (cxsc_l_complex_division_p[2])
370 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BRESUP, BIMINF, +1 ) );
371 if (cxsc_l_complex_division_p[3])
372 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BREINF, BIMSUP, +1 ) );
373 if (cxsc_l_complex_division_p[4])
374 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BRESUP, BIMSUP, +1 ) );
393 for (j=1; j<=rep; j++)
395 for (i=1; i<=4; cxsc_l_complex_division_p[i++] =
true);
399 min( min( minmax(
true, a0, b0, BIMINF, BRE, 1,2 ),
400 minmax(
true, a0, b0, BIMSUP, BRE, 3,4 ) ),
401 min( minmax(
true, b0, a0, BREINF, BIM, 1,3 ),
402 minmax(
true, b0, a0, BRESUP, BIM, 2,4 ) ) )
405 if (cxsc_l_complex_division_p[1])
406 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BREINF, BIMINF, -1 ) );
407 if (cxsc_l_complex_division_p[2])
408 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BRESUP, BIMINF, -1 ) );
409 if (cxsc_l_complex_division_p[3])
410 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BREINF, BIMSUP, -1 ) );
411 if (cxsc_l_complex_division_p[4])
412 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BRESUP, BIMSUP, -1 ) );
424 a_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
425 b_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
441 for (j=1; j<=rep; j++)
443 for (i=1; i<=4; cxsc_l_complex_division_p[i++] =
true) ;
447 max( max( minmax(
false, b0, -a0, BIMINF, BRE, 1,2 ),
448 minmax(
false, b0, -a0, BIMSUP, BRE, 3,4 ) ),
449 max( minmax(
false, -a0, b0, BREINF, BIM, 1,3 ),
450 minmax(
false, -a0, b0, BRESUP, BIM, 2,4 ) ) )
453 if (cxsc_l_complex_division_p[1])
454 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BREINF, BIMINF, +1 ) );
455 if (cxsc_l_complex_division_p[2])
456 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BRESUP, BIMINF, +1 ) );
457 if (cxsc_l_complex_division_p[3])
458 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BREINF, BIMSUP, +1 ) );
459 if (cxsc_l_complex_division_p[4])
460 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BRESUP, BIMSUP, +1 ) );
480 for (j=1; j<=rep; j++)
482 for (i=1; i<=4; cxsc_l_complex_division_p[i++] =
true) ;
486 min( min( minmax(
true, b0, -a0, BIMINF, BRE, 1,2 ),
487 minmax(
true, b0, -a0, BIMSUP, BRE, 3,4 ) ),
488 min( minmax(
true, -a0, b0, BREINF, BIM, 1,3 ),
489 minmax(
true, -a0, b0, BRESUP, BIM, 2,4 ) ) )
492 if (cxsc_l_complex_division_p[1])
493 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BREINF, BIMINF, -1 ) );
494 if (cxsc_l_complex_division_p[2])
495 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BRESUP, BIMINF, -1 ) );
496 if (cxsc_l_complex_division_p[3])
497 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BREINF, BIMSUP, -1 ) );
498 if (cxsc_l_complex_division_p[4])
499 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BRESUP, BIMSUP, -1 ) );
507 return l_cinterval(l_interval(realteilINF, realteilSUP),
508 l_interval(imagteilINF, imagteilSUP));
511l_cinterval C_point_div(
const l_cinterval& z,
const l_cinterval& n)
517 a = l_complex(InfRe(z),InfIm(z));
518 b = l_complex(InfRe(n),InfIm(n));
523 re = l_interval( Re(q1),Re(q2) );
524 im = l_interval( Im(q1),Im(q2) );
526 return l_cinterval(re,im);
533 if (0.0 <= b.re && 0.0 <= b.im ) {
535 cxscthrow(DIV_BY_ZERO(
"l_cinterval operator / (const l_cinterval&, const l_cinterval&)"));
538 bool a_point, b_point;
539 a_point = InfRe(a)==SupRe(a) && InfIm(a)==SupIm(a);
540 b_point = InfRe(b)==SupRe(b) && InfIm(b)==SupIm(b);
541 if(a_point && b_point)
return C_point_div(a,b);
542 else return cidiv(a,b);
553std::ostream & operator << (std::ostream &s,
const l_cinterval& a)
noexcept
562std::string & operator << (std::string &s,
const l_cinterval& a)
noexcept
591 int stagprec_old(stagprec);
594 s = skipwhitespacessinglechar (s,
'(');
595 s = skipwhitespacessinglechar (s,
'[');
597 stagprec = StagPrec(a.re);
600 s = skipwhitespacessinglechar (s,
',');
606 stagprec = StagPrec(a.im);
607 s = skipwhitespacessinglechar (s,
']');
608 s = skipwhitespacessinglechar (s,
',');
609 s = skipwhitespacessinglechar (s,
'[');
614 s = skipwhitespacessinglechar (s,
',');
621 s = skipwhitespaces (s);
624 s = skipwhitespaces (s);
627 stagprec = stagprec_old;
628 if (Inf(a.re) > Sup(a.re) || Inf(a.im) > Sup(a.im))
629 cxscthrow(EMPTY_INTERVAL
630 (
"std::string & operator >> (std::string &s, cinterval &a)"));
656 skipeolnflag = inpdotflag =
true;
657 stagprec = StagPrec(a.re);
658 c = skipwhitespacessinglechar (s,
'(');
659 if (inpdotflag) s.putback(c);
660 c = skipwhitespacessinglechar (s,
'[');
661 if (inpdotflag) s.putback(c);
665 skipeolnflag = inpdotflag =
true;
666 c = skipwhitespacessinglechar (s,
',');
667 if (inpdotflag) s.putback(c);
672 c = skipwhitespacessinglechar (s,
']');
673 if (inpdotflag) s.putback(c);
674 c = skipwhitespacessinglechar (s,
',');
675 if (inpdotflag) s.putback(c);
677 c = skipwhitespacessinglechar (s,
'[');
678 if (inpdotflag) s.putback(c);
680 stagprec = StagPrec(a.im);
684 skipeolnflag = inpdotflag =
true;
685 c = skipwhitespacessinglechar (s,
',');
686 if (inpdotflag) s.putback(c);
696 skipeolnflag =
false, inpdotflag =
true;
697 c = skipwhitespaces (s);
698 if (inpdotflag && c !=
']')
703 skipeolnflag =
false, inpdotflag =
true;
704 c = skipwhitespaces (s);
705 if (inpdotflag && c !=
')')
708 if (Inf(a.re) > Sup(a.re) || Inf(a.im) > Sup(a.im))
709 cxscthrow(EMPTY_INTERVAL
710 (
"std::istream & operator >> (std::istream &s, cinterval &a)"));
The Data Type cdotprecision.
The Data Type cidotprecision.
The Scalar Type cinterval.
cinterval(void) noexcept
Constructor of class cinterval.
cinterval & operator=(const real &) noexcept
Implementation of standard assigning operator.
The Data Type dotprecision.
The Data Type idotprecision.
The Scalar Type interval.
The Multiple-Precision Data Type l_cinterval.
l_cinterval(void) noexcept
Constructor of class l_cinterval.
The Multiple-Precision Data Type l_interval.
The Multiple-Precision Data Type l_real.
The namespace cxsc, providing all functionality of the class library C-XSC.
const real MaxReal
Greatest representable floating-point number.
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
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 .
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .