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))
60 bool cxsc_complex_division_p[5];
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 );
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;
308 real realteilINF, realteilSUP,
309 imagteilINF, imagteilSUP;
312 bool a_repeat,b_repeat;
314 real AREINF, ARESUP, AIMINF, AIMSUP,
315 BREINF, BRESUP, BIMINF, BIMSUP;
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 ) );
511 interval(imagteilINF, imagteilSUP));
520 a =
complex(InfRe(z),InfIm(z));
521 b =
complex(InfRe(n),InfIm(n));
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);
553 std::ostream & operator << (std::ostream &s,
const cinterval& a)
throw()
561 std::string & operator << (std::string &s,
const cinterval &a)
throw()
571 std::istream & operator >> (std::istream &s,
cinterval &a)
throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
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)"));
617 std::string & operator >> (std::string &s,
cinterval &a)
throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
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)"));
643 void operator >>(
const std::string &s,
cinterval &a)
throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
648 void operator >>(
const char *s,
cinterval &a)
throw(ERROR_CINTERVAL_EMPTY_INTERVAL)
656 return (
in(Re(x),Re(y)) &&
in(Im(x),Im(y)) );
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.
cinterval(void)
Constructor of class cinterval.
The Data Type idotprecision.
The Data Type dotprecision.
interval sqrtx2y2(const interval &x, const interval &y)
Calculates .
The namespace cxsc, providing all functionality of the class library C-XSC.
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
The Scalar Type interval.
The Data Type cidotprecision.
const real MaxReal
Greatest representable floating-point number.
The Data Type cdotprecision.
The Scalar Type cinterval.
cinterval sqrt(const cinterval &z)
Calculates .
void times2pown(cinterval &x, int n)
Fast multiplication of reference parameter [z] with .
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.