26#include "cinterval.hpp"
40 re(rnd(InfRe(a),RND_DOWN),rnd(SupRe(a),RND_UP)),
41 im(rnd(InfIm(a),RND_DOWN),rnd(SupIm(a),RND_UP))
55void product (real, real, real, real,
int&, real&, interval&);
56real quotient (real, interval, real, interval,
int,
int,
int);
60bool cxsc_complex_division_p[5];
62real cxsc_complex_division_f(real a, real b, real c, real d,
int round)
70 product( a, c, b, d, zoverfl, z1,z2 );
71 product( c, c, d, d, noverfl, n1,n2 );
73 return quotient( z1,z2, n1,n2, round, zoverfl, noverfl );
76static real minmax(
int minimum, real a, real b, real y0,
77 interval x,
int i,
int j)
83 real q1,q2,w1,w2,t1,t2,x1,x2,ay0, minmax;
94 if (cxsc_complex_division_p[i] && cxsc_complex_division_p[j])
95 minmax = cxsc_complex_division_f( a, b, Inf(x), y0, 1-2*minimum );
97 cxsc_complex_division_p[i] =
false;
98 cxsc_complex_division_p[j] =
false;
102 if ( b == CXSC_Zero || y0 == CXSC_Zero )
105 cxsc_complex_division_p[i] =
false;
106 cxsc_complex_division_p[j] =
false;
110 if (minimum && sign(b) != sign(y0) )
112 minmax = divd(b, y0);
113 cxsc_complex_division_p[i] =
false;
114 cxsc_complex_division_p[j] =
false;
116 if (!minimum && sign(b) == sign(y0) )
118 minmax = divu(b, y0);
119 cxsc_complex_division_p[i] =
false;
120 cxsc_complex_division_p[j] =
false;
132 minmax = divd(a, Sup(x));
134 minmax = divd(a, Inf(x));
138 minmax = divu(a, Inf(x));
140 minmax = divu(a, Sup(x));
142 cxsc_complex_division_p[i] =
false;
143 cxsc_complex_division_p[j] =
false;
155 real invf2=1, a_skal;
156 int exf1=0, exf2=0, exinf1=0, exinf2=0;
158 if (sign(b)==0) { t1 = 1; t2 = 0; }
167 int expo_diff = expo(b)-expo(a), ex;
168 if (expo_diff >= 512)
170 int exdiff5 = expo_diff - 500;
178 invf2 = comp(0.5,1-exf2);
187 invf2 = comp(0.5,2-exdiff5);
198 accumulate(akku, -q1, a_skal);
199 q2 = rnd(akku) / a_skal;
202 if (exinf1 == 0) accumulate(akku, invf2, invf2);
203 accumulate(akku, q1, q1);
204 accumulate(akku, q1, q2);
205 accumulate(akku, q1, q2);
206 accumulate(akku, q2, q2);
208 w1 =
sqrt(rnd(akku));
210 accumulate(akku, -w1, w1);
211 w2 = rnd(akku) / (2.0*w1);
231 if (( sign(b) == sign(y0) ) == minimum)
235 accumulate(akku,ay0,t1);
236 accumulate(akku,ay0,t2);
238 if (expo(x1) == 2147483647)
goto Ende;
243 if (expo(x1)+exf1 > 1023)
goto Ende;
245 if (expo(x1)+exf2 > 1023)
goto Ende;
255 accumulate(akku, -t1, x1);
256 accumulate(akku, -t2, x1);
260 if (expo(x1)+exinf1 > 1023)
goto Ende;
262 if (expo(x1)+exinf2 > 1023)
goto Ende;
279 accumulate(akku, -x1, q1);
280 accumulate(akku, -x2, q1);
285 q2 = rnd(akku) / (2.0*x1);
288 if (sign(q2)==0 && sign(akku)!=0) minmax = pred(q1);
289 else minmax = addd(q1, q2);
293 if (sign(q2)==0 && sign(akku)!=0) minmax = succ(q1);
294 else minmax = addu(q1, q2);
297 cxsc_complex_division_p[i] =
false;
298 cxsc_complex_division_p[j] =
false;
306cinterval cidiv(
const cinterval& A,
const cinterval& B)
308 real realteilINF, realteilSUP,
309 imagteilINF, imagteilSUP;
312 bool a_repeat,b_repeat;
314 real AREINF, ARESUP, AIMINF, AIMSUP,
315 BREINF, BRESUP, BIMINF, BIMSUP;
316 interval ARE, AIM, BRE, BIM;
338 a_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
339 b_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
341 if (a_repeat || b_repeat)
358 for (j=1; j<=rep; j++)
360 for (i=1; i<=4; cxsc_complex_division_p[i++] =
true);
364 max( max( minmax(
false, a0, b0, BIMINF, BRE, 1,2 ),
365 minmax(
false, a0, b0, BIMSUP, BRE, 3,4 ) ),
366 max( minmax(
false, b0, a0, BREINF, BIM, 1,3 ),
367 minmax(
false, b0, a0, BRESUP, BIM, 2,4 ) ) )
371 if (cxsc_complex_division_p[1])
372 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BREINF, BIMINF, +1 ) );
373 if (cxsc_complex_division_p[2])
374 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BRESUP, BIMINF, +1 ) );
375 if (cxsc_complex_division_p[3])
376 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BREINF, BIMSUP, +1 ) );
377 if (cxsc_complex_division_p[4])
378 realteilSUP = max( realteilSUP, cxsc_complex_division_f( a0, b0, BRESUP, BIMSUP, +1 ) );
397 for (j=1; j<=rep; j++)
399 for (i=1; i<=4; cxsc_complex_division_p[i++] =
true);
403 min( min( minmax(
true, a0, b0, BIMINF, BRE, 1,2 ),
404 minmax(
true, a0, b0, BIMSUP, BRE, 3,4 ) ),
405 min( minmax(
true, b0, a0, BREINF, BIM, 1,3 ),
406 minmax(
true, b0, a0, BRESUP, BIM, 2,4 ) ) )
408 if (cxsc_complex_division_p[1])
409 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BREINF, BIMINF, -1 ) );
410 if (cxsc_complex_division_p[2])
411 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BRESUP, BIMINF, -1 ) );
412 if (cxsc_complex_division_p[3])
413 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BREINF, BIMSUP, -1 ) );
414 if (cxsc_complex_division_p[4])
415 realteilINF = min( realteilINF, cxsc_complex_division_f( a0, b0, BRESUP, BIMSUP, -1 ) );
427 a_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
428 b_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
444 for (j=1; j<=rep; j++)
446 for (i=1; i<=4; cxsc_complex_division_p[i++] =
true) ;
450 max( max( minmax(
false, b0, -a0, BIMINF, BRE, 1,2 ),
451 minmax(
false, b0, -a0, BIMSUP, BRE, 3,4 ) ),
452 max( minmax(
false, -a0, b0, BREINF, BIM, 1,3 ),
453 minmax(
false, -a0, b0, BRESUP, BIM, 2,4 ) ) )
456 if (cxsc_complex_division_p[1])
457 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BREINF, BIMINF, +1 ) );
458 if (cxsc_complex_division_p[2])
459 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BRESUP, BIMINF, +1 ) );
460 if (cxsc_complex_division_p[3])
461 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BREINF, BIMSUP, +1 ) );
462 if (cxsc_complex_division_p[4])
463 imagteilSUP = max( imagteilSUP, cxsc_complex_division_f( b0, -a0, BRESUP, BIMSUP, +1 ) );
483 for (j=1; j<=rep; j++)
485 for (i=1; i<=4; cxsc_complex_division_p[i++] =
true) ;
489 min( min( minmax(
true, b0, -a0, BIMINF, BRE, 1,2 ),
490 minmax(
true, b0, -a0, BIMSUP, BRE, 3,4 ) ),
491 min( minmax(
true, -a0, b0, BREINF, BIM, 1,3 ),
492 minmax(
true, -a0, b0, BRESUP, BIM, 2,4 ) ) )
495 if (cxsc_complex_division_p[1])
496 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BREINF, BIMINF, -1 ) );
497 if (cxsc_complex_division_p[2])
498 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BRESUP, BIMINF, -1 ) );
499 if (cxsc_complex_division_p[3])
500 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BREINF, BIMSUP, -1 ) );
501 if (cxsc_complex_division_p[4])
502 imagteilINF = min( imagteilINF, cxsc_complex_division_f( b0, -a0, BRESUP, BIMSUP, -1 ) );
510 return cinterval(interval(realteilINF, realteilSUP),
511 interval(imagteilINF, imagteilSUP));
514cinterval C_point_div(
const cinterval& z,
const cinterval& n)
520 a = complex(InfRe(z),InfIm(z));
521 b = complex(InfRe(n),InfIm(n));
526 re = interval( Re(q1),Re(q2) );
527 im = interval( Im(q1),Im(q2) );
529 return cinterval(re,im);
532cinterval div_operator (
const cinterval & a,
const cinterval & b)
534 bool a_point, b_point;
535 a_point = InfRe(a)==SupRe(a) && InfIm(a)==SupIm(a);
536 b_point = InfRe(b)==SupRe(b) && InfIm(b)==SupIm(b);
537 if(a_point && b_point)
return C_point_div(a,b);
538 else return cidiv(a,b);
553std::ostream & operator << (std::ostream &s,
const cinterval& a)
noexcept
561std::string & operator << (std::string &s,
const cinterval &a)
noexcept
571std::istream & operator >> (std::istream &s,
cinterval &a)
574 skipeolnflag = inpdotflag =
true;
575 c = skipwhitespacessinglechar (s,
'(');
576 if (inpdotflag) s.putback(c);
577 c = skipwhitespacessinglechar (s,
'[');
578 if (inpdotflag) s.putback(c);
579 s >> SaveOpt >> RndDown >> Inf(a.re);
580 skipeolnflag = inpdotflag =
true;
581 c = skipwhitespacessinglechar (s,
',');
582 if (inpdotflag) s.putback(c);
583 s >> RndUp >> Sup(a.re);
584 c = skipwhitespacessinglechar (s,
']');
585 if (inpdotflag) s.putback(c);
586 c = skipwhitespacessinglechar (s,
',');
587 if (inpdotflag) s.putback(c);
589 c = skipwhitespacessinglechar (s,
'[');
590 if (inpdotflag) s.putback(c);
591 s >> RndDown >> Inf(a.im);
592 skipeolnflag = inpdotflag =
true;
593 c = skipwhitespacessinglechar (s,
',');
594 if (inpdotflag) s.putback(c);
595 s >> RndUp >> Sup(a.im) >> RestoreOpt;
599 skipeolnflag =
false, inpdotflag =
true;
600 c = skipwhitespaces (s);
601 if (inpdotflag && c !=
']')
606 skipeolnflag =
false, inpdotflag =
true;
607 c = skipwhitespaces (s);
608 if (inpdotflag && c !=
')')
611 if (Inf(a.re) > Sup(a.re) || Inf(a.im) > Sup(a.im))
612 cxscthrow(ERROR_CINTERVAL_EMPTY_INTERVAL(
"std::istream & operator >> (std::istream &s, cinterval &a)"));
619 s = skipwhitespacessinglechar (s,
'(');
620 s = skipwhitespacessinglechar (s,
'[');
621 s = s >> SaveOpt >> RndDown >> Inf(a.re);
622 s = skipwhitespacessinglechar (s,
',');
623 s = s >> RndUp >> Sup(a.re);
624 s = skipwhitespacessinglechar (s,
']');
625 s = skipwhitespacessinglechar (s,
',');
626 s = skipwhitespacessinglechar (s,
'[');
627 s = s >> RndDown >> Inf(a.im);
628 s = skipwhitespacessinglechar (s,
',');
629 s = s >> RndUp >> Sup(a.im) >> RestoreOpt;
630 s = skipwhitespaces (s);
633 s = skipwhitespaces (s);
637 if (Inf(a.re) > Sup(a.re) || Inf(a.im) > Sup(a.im))
638 cxscthrow(ERROR_CINTERVAL_EMPTY_INTERVAL(
"std::string & operator >> (std::string &s, cinterval &a)"));
656 return (
in(Re(x),Re(y)) &&
in(Im(x),Im(y)) );
The Data Type cdotprecision.
The Data Type cidotprecision.
The Scalar Type cinterval.
cinterval(void) noexcept
Constructor of class cinterval.
The Data Type dotprecision.
The Data Type idotprecision.
The Scalar Type interval.
The namespace cxsc, providing all functionality of the class library C-XSC.
const real MaxReal
Greatest representable floating-point number.
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
cinterval sqrt(const cinterval &z) noexcept
Calculates .
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.