C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
l_cinterval.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_cinterval.cpp,v 1.18 2014/01/30 17:23:46 cxsc Exp $ */
25
26#include "l_cinterval.hpp"
27#include "cidot.hpp"
28#include "dot.hpp"
29#include "l_rmath.hpp"
30#include "l_imath.hpp"
31
32namespace cxsc {
33
34
35#define CXSC_Zero 0.0
36
38{
39 interval u,v;
40 u = a.re;
41 v = a.im;
42 *this = cinterval(u,v);
43}
44
46{
47 interval u,v;
48 u = a.re;
49 v = a.im;
50return *this = cinterval(u,v);
51}
52
53l_cinterval::l_cinterval(const dotprecision &a) noexcept : re(a),im(0) {}
54l_cinterval::l_cinterval(const idotprecision &a) noexcept : re(a),im(0) {}
56 noexcept : re(Re(a)),im(Im(a)) {}
58 re( l_interval(Re(a))),im(l_interval(Im(a)) ) {}
59
60l_cinterval operator * (const l_cinterval & a, const l_cinterval & b) noexcept
61{
62 idotprecision akku;
63 l_cinterval res;
64 l_interval u,v;
65 akku = 0.0;
66 accumulate(akku,a.re,b.re);
67 accumulate(akku,-a.im,b.im);
68 u = akku;
69 if (Inf(u)>Sup(u))
70 { std::cerr << "Error in l_cinterval * l_cinterval" << std::endl;
71 exit(1);
72 }
73 akku = 0.0;
74 accumulate(akku,a.im,b.re);
75 accumulate(akku,a.re,b.im);
76 v = akku; // v: Imaginaerteil
77 if (Inf(v)>Sup(v))
78 { std::cerr << "Error in l_cinterval * l_cinterval" << std::endl;
79 exit(1);
80 }
81 res = l_cinterval(u,v);
82 return res;
83}
84
85
86// *********************************************************************
87// In l_complex.cpp implemented: (In l_complex.hpp not declared!)
88void product(const l_real& a, const l_real& b, const l_real& c,
89 const l_real& d, int& ex, l_interval& res);
90void product(const l_real& c, const l_real& d, int& ex, l_interval& res);
91l_real quotient(const l_interval& z, const l_interval& n, int round,
92 int ex_z, int ex_n);
93//void Times2pown(l_interval& a, int n) noexcept;
94
95// *********************************************************************
96
97static const int max_expo = 1020, max_expo1 = 1023;
98
99// optimale komplexe Intervalldivision
100
101bool cxsc_l_complex_division_p[5];
102
103l_real cxsc_complex_division_f(l_real a, l_real b, l_real c, l_real d,
104 int round)
105{
106 int ex1, ex2;
107 l_interval z,n;
108
109 // f:=(a*c+b*d)/(SQR(c)+SQR(d))
110
111 product(a, c, b, d, ex1, z);
112 product(c, d, ex2, n);
113 return quotient(z, n, round, ex1, ex2);
114}
115
116// *************************************************************************
117// *************************************************************************
118
119static l_real minmax(int minimum, l_real a, l_real b, l_real y0,
120 l_interval x, int i, int j)
121// Calculates the inner minimum or maximum of f = (ac+bd)/(cc+dd)
122// on the interval x = [c.inf,c.sup] ( a,b,d=y0 fixated ).
123// If minimum = true the minimum will be calculated, otherwise the maximum.
124
125{
126 l_real ay0, minmax;
127 l_interval t,q,x0,two_Da;
128
129 int Da(0);
130
131 a += 0.0; b += 0.0; y0 += 0.0;
132
133 if (minimum)
134 minmax = MaxReal;
135 else
136 minmax = -MaxReal;
137
138 if (Inf(x) == Sup(x))
139 {
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 );
142
143 cxsc_l_complex_division_p[i] = false;
144 cxsc_l_complex_division_p[j] = false;
145 } else
146 if (a == 0.0)
147 {
148 if ( b == CXSC_Zero || y0 == CXSC_Zero )
149 {
150 minmax = 0.0;
151 cxsc_l_complex_division_p[i] = false;
152 cxsc_l_complex_division_p[j] = false;
153 } else
154 { // b*y0 <> 0:
155 if (0.0 < x) {
156 if (minimum && sign(b) != sign(y0) )
157 {
158 // minmax = divd(b, y0);
159 minmax = Inf(l_interval(b)/y0);
160 cxsc_l_complex_division_p[i] = false;
161 cxsc_l_complex_division_p[j] = false;
162 } else
163 if (!minimum && sign(b) == sign(y0) )
164 {
165 // minmax = divu(b, y0);
166 minmax = Sup(l_interval(b)/y0);
167 cxsc_l_complex_division_p[i] = false;
168 cxsc_l_complex_division_p[j] = false;
169 }
170 }
171 }
172 } else
173 {
174 // a != 0.0
175 if (y0 == 0.0)
176 {
177 if (minimum)
178 {
179 if (a > 0.0)
180 // minmax = divd(a, Sup(x));
181 minmax = Inf(l_interval(a)/l_interval(Sup(x)));
182 else
183 // minmax = divd(a, Inf(x));
184 minmax = Inf(l_interval(a)/l_interval(Inf(x)));
185 } else
186 {
187 if (a > 0.0)
188 // minmax = divu(a, Inf(x));
189 minmax = Sup(l_interval(a)/l_interval(Inf(x)));
190 else
191 // minmax = divu(a, Sup(x));
192 minmax = Sup(l_interval(a)/l_interval(Sup(x)));
193 }
194 cxsc_l_complex_division_p[i] = false;
195 cxsc_l_complex_division_p[j] = false;
196 } else
197 { // a != 0.0 and
198 // y0 != 0.0, Calculation of extrema points and minimum|maximum
199 // values:
200 // Calculation of: t = sign(a) * ( |b/a| + sqrt( 1+|b/a|^2 ) )
201
202 l_real invf2(1.0), a_skal;
203 // int exf1=0, exf2=0, exinf1=0, exinf2=0; // unused variable
204
205 // We first calculate: t = |b/a| + sqrt( 1+|b/a|^2 );
206
207 if (sign(b)==0) t = 1.0;
208 else
209 { // a != 0.0 and b != 0;
210 // To avoid overflow by calculating |b/a| + sqrt( 1+|b/a|^2 )
211 // we must multiply a with 2^Da:
212 int expo_diff = expo(b[1]) - expo(a[1]), ex;
213 a_skal = a;
214 if (expo_diff > max_expo)
215 {
216 Da = expo_diff-max_expo; // Da > 0;
217 // a must be multiplied with 2^Da to avoid overflow
218 // by calculating |b/a| + sqrt( 1+|b/a|^2 ) :
219 if (Da>max_expo1)
220 {
221 times2pown(a_skal,max_expo1);
222 ex = Da - max_expo1;
223 times2pown(a_skal,ex);
224 }
225 else times2pown(a_skal,Da);
226
227 // Now calculating an inclusion t of 2^(-Da):
228 if (Da>1022)
229 {
230 two_Da = l_interval( comp(0.5,-1021) );
231 times2pown(two_Da,1022-Da);
232 }
233 else two_Da = l_interval( comp(0.5,1-Da) );
234 // Now two_Da is am inclusion of 2^(-Da);
235 }
236 q = l_interval(b)/a_skal;
237 if (sign(q[1])<0) q = -q;
238 // q: Inclusion of |b/(a*2^Da)|;
239
240 t = (Da > 0)? q + sqrtx2y2(two_Da,q) : q + sqrt1px2(q);
241 }
242
243 if (a<0.0) t = -t;
244
245// if (Da > 0) the value t from the last line must additionally be
246// multiplied with 2^Da, to get an inclusion of the expression:
247// sign(a) * ( |b/a| + sqrt( 1+|b/a|^2 ) );
248
249// Now to a fall differentiation for min-,max- calculation:
250// First we will calculate an inclusion x0 of the point of the
251// relative minimum or maximum:
252
253 ay0 = abs(y0);
254
255 if ( (sign(b) == sign(y0)) == minimum )
256 { // Calculation of x0 = |y0| * t :
257 if (expo(ay0[1]) + expo(t[1]) + Da > max_expo1) goto Ende;
258 else x0 = ay0 * t;
259 if (Da>0) Times2pown(x0,Da);
260 }
261 else // Calculation of x0 = |y0| / t :
262 {
263 if (expo(ay0[1]) - expo(t[1]) - Da > max_expo1) goto Ende;
264 else x0 = ay0 / t;
265 if (Da>0) Times2pown(x0,-Da);
266 }
267
268 if (minimum) x0 = -x0;
269
270 if (x0 < x) // The minimum or maximum point lies in
271 { // the interior of x.
272 // Calculation of: a / ( 2*x0 )
273 q = a/x0;
274 times2pown(q,-1); // q: inclusion of a / ( 2*x0 );
275
276 if (minimum) minmax = Inf(q);
277 else minmax = Sup(q);
278
279 cxsc_l_complex_division_p[i] = false;
280 cxsc_l_complex_division_p[j] = false;
281 }
282 Ende:;
283 } // y0 != 0.0
284 }
285 return minmax;
286} // *** minmax ***
287
288l_real max(const l_real& u, const l_real& v)
289{
290 l_real res(u);
291 if (v>u) res = v;
292 return res;
293}
294
295l_real min(const l_real& u, const l_real& v)
296{
297 l_real res(u);
298 if (v<u) res = v;
299 return res;
300}
301
302l_cinterval cidiv(const l_cinterval& A, const l_cinterval& B)
303{
304 l_real realteilINF, realteilSUP,
305 imagteilINF, imagteilSUP;
306 // da sonst eventuell zwischendurch illegale Intervalle entstehen
307 l_real a0,b0;
308 bool a_repeat,b_repeat;
309 int i, rep, j;
310 l_real AREINF, ARESUP, AIMINF, AIMSUP,
311 BREINF, BRESUP, BIMINF, BIMSUP;
312 l_interval ARE, AIM, BRE, BIM;
313
314 // keine Fehlerabfrage -> ist schon in CINTVAL.CPP
315 // IF ( 0.0 IN B.RE ) AND ( 0.0 IN B.IM ) THEN
316 // CIDIVISION:= COMPL( 1.0 / INTVAL(-1.0,1.0), INTVAL(0.0) );
317 // Fehlerabbruch erzwingen: Intervall enthaelt 0
318
319 // *** Berechnung des Realteils ***
320
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;
329 ARE = Re(A);
330 AIM = Im(A);
331 BRE = Re(B);
332 BIM = Im(B);
333
334 a_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
335 b_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
336
337 if (a_repeat || b_repeat)
338 rep = 2;
339 else
340 rep = 1;
341
342 if (BREINF >= 0.0)
343 a0 = ARESUP;
344 else
345 a0 = AREINF;
346
347 if (BIMINF >= 0.0)
348 b0 = AIMSUP;
349 else
350 b0 = AIMINF;
351
352 realteilSUP = -MaxReal;
353
354 for (j=1; j<=rep; j++)
355 {
356 for (i=1; i<=4; cxsc_l_complex_division_p[i++] = true);
357
358 realteilSUP =
359 max( realteilSUP,
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 ) ) )
364
365 );
366
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 ) );
375
376 if (a_repeat)
377 a0 = ARESUP;
378 else if (b_repeat)
379 b0 = AIMSUP;
380 }
381
382 if (BREINF >= 0.0)
383 a0 = AREINF;
384 else
385 a0 = ARESUP;
386 if (BIMINF >= 0.0)
387 b0 = AIMINF;
388 else
389 b0 = AIMSUP;
390
391 realteilINF = MaxReal;
392
393 for (j=1; j<=rep; j++)
394 {
395 for (i=1; i<=4; cxsc_l_complex_division_p[i++] = true);
396
397 realteilINF =
398 min( realteilINF,
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 ) ) )
403 );
404
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 ) );
413
414 if (a_repeat)
415 a0 = AREINF;
416 else if (b_repeat)
417 b0 = AIMINF;
418 }
419
420
421 // Calculation of the img. part:
422 // g(a, b, c, d) = cxsc_complex_division_f(b, -a, c, d)
423
424 a_repeat = ( BIMINF < CXSC_Zero ) && ( CXSC_Zero < BIMSUP );
425 b_repeat = ( BREINF < CXSC_Zero ) && ( CXSC_Zero < BRESUP );
426
427 // IF a_repeat OR b_repeat THEN rep:= 2 ELSE rep:= 1;
428
429 if (BREINF >= 0.0)
430 b0 = AIMSUP;
431 else
432 b0 = AIMINF;
433
434 if (BIMINF >= 0.0)
435 a0 = AREINF;
436 else
437 a0 = ARESUP;
438
439 imagteilSUP = -MaxReal;
440
441 for (j=1; j<=rep; j++)
442 {
443 for (i=1; i<=4; cxsc_l_complex_division_p[i++] = true) ;
444
445 imagteilSUP =
446 max( imagteilSUP,
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 ) ) )
451 );
452
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 ) );
461
462 if (b_repeat)
463 b0 = AIMSUP;
464 else if (a_repeat)
465 a0 = AREINF ;
466 }
467
468 if (BREINF >= 0.0)
469 b0 = AIMINF;
470 else
471 b0 = AIMSUP;
472
473 if (BIMINF >= 0.0)
474 a0 = ARESUP;
475 else
476 a0 = AREINF;
477
478 imagteilINF = MaxReal;
479
480 for (j=1; j<=rep; j++)
481 {
482 for (i=1; i<=4; cxsc_l_complex_division_p[i++] = true) ;
483
484 imagteilINF =
485 min( imagteilINF,
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 ) ) )
490 );
491
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 ) );
500
501 if (b_repeat)
502 b0 = AIMINF;
503 else if (a_repeat)
504 a0 = ARESUP;
505 }
506
507 return l_cinterval(l_interval(realteilINF, realteilSUP),
508 l_interval(imagteilINF, imagteilSUP));
509} // cidiv
510
511l_cinterval C_point_div(const l_cinterval& z, const l_cinterval& n)
512// Division of complex point intervals;
513// z,n must be point intervals!! Blomquist, 07,11.02
514// This function only for internal use!
515{
516 l_complex a,b,q1,q2;
517 a = l_complex(InfRe(z),InfIm(z));
518 b = l_complex(InfRe(n),InfIm(n));
519 q1 = divd(a,b);
520 q2 = divu(a,b);
521
522 l_interval re, im;
523 re = l_interval( Re(q1),Re(q2) );
524 im = l_interval( Im(q1),Im(q2) );
525
526 return l_cinterval(re,im);
527} // C_point_div
528
529
531
532{
533 if (0.0 <= b.re && 0.0 <= b.im ) {
534// if (0.0 <= (sqr(b.re) + sqr(b.im))) {
535 cxscthrow(DIV_BY_ZERO("l_cinterval operator / (const l_cinterval&, const l_cinterval&)"));
536 return a; // dummy result
537 }
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); // a,b are point intervals
542 else return cidiv(a,b);
543}
544
545l_interval abs(const l_cinterval &a) noexcept
546{
547 return sqrtx2y2(a.re,a.im);
548}
549
550
551// ---- Ausgabefunkt. ---------------------------------------
552
553std::ostream & operator << (std::ostream &s, const l_cinterval& a) noexcept
554{
555 s << '('
556 << a.re << ','
557 << a.im
558 << ')';
559 return s;
560}
561
562std::string & operator << (std::string &s, const l_cinterval& a) noexcept
563{
564// string s; l_cinterval a;
565// s << a; s delivers the string of the value a in the form:
566// ([Inf(real-part(a)),Sup(real-part(a))],[Inf(img-part(a)),Sup(img-part(a))])
567 s+='(';
568 s << a.re;
569 s+=',';
570 s << a.im;
571 s+=')';
572 return s;
573}
574
575std::string & operator >> (std::string &s, l_cinterval &a)
576
577// With:
578// l_cinterval a;
579// string("([1.234,1.234],[2.567,2.567])") >> a;
580// the value a will be an inclusion of the above string.
581// The actual precisions of the staggered intervals a.re and a.im
582// will not be affected by the operator >> !
583// The above braces, brackets and commas must not be used in the
584// string, however the four numbers must then be seperated by spaces!
585// Thus, the following string will produce the same inclusion a:
586// string("1.234 1.234 2.567 2.567 ") >> a;
587// Blomquist, 15.11.2006;
588{
589 l_real Iar,Sar,Iai,Sai;
590 l_interval lr,li;
591 int stagprec_old(stagprec);
592 dotprecision dot;
593
594 s = skipwhitespacessinglechar (s, '(');
595 s = skipwhitespacessinglechar (s, '[');
596 s = s >> dot;
597 stagprec = StagPrec(a.re);
598 lr = l_interval(dot);
599 Iar = Inf(lr);
600 s = skipwhitespacessinglechar (s, ',');
601 s = s >> dot;
602 lr = l_interval(dot);
603 Sar = Sup(lr);
604 lr = l_interval(Iar,Sar);
605
606 stagprec = StagPrec(a.im);
607 s = skipwhitespacessinglechar (s, ']');
608 s = skipwhitespacessinglechar (s, ',');
609 s = skipwhitespacessinglechar (s, '[');
610 s = s >> dot;
611
612 li = l_interval(dot);
613 Iai = Inf(li);
614 s = skipwhitespacessinglechar (s, ',');
615 s = s >> dot;
616 li = l_interval(dot);
617 Sai = Sup(li);
618 li = l_interval(Iai,Sai);
619
620 a = l_cinterval(lr,li );
621 s = skipwhitespaces (s);
622 if (s[0] == ']')
623 s.erase(0,1);
624 s = skipwhitespaces (s);
625 if (s[0] == ')')
626 s.erase(0,1);
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)"));
631
632 return s;
633}
634
635std::istream & operator >> (std::istream & s, l_cinterval& a)
636
637// With:
638// l_cinterval lc;
639// cout << "([a,b],[c,d]) = ?" << endl;
640// cin >> lc;
641// the input string ([1.23,1.23],[3.45,3.45]) will be included by lc.
642// The actual precisions of the staggered intervals lc.re and lc.im
643// will not be affected by the operator >> !
644// The above braces, brackets and commas must not be used in the
645// string, however the four numbers a,b,c,d must then be seperated by
646// spaces! Thus, the following input string 1.23 1.23 3.45 3.45
647// will produce the same inclusion lc:
648// Blomquist, 15.11.2006;
649{
650 l_real Iar,Sar,Iai,Sai;
651 l_interval lr,li;
652 dotprecision dot;
653 // int stagprec_old(stagprec); // unused variable
654
655 char c;
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);
662 s >> dot;
663 lr = l_interval(dot);
664 Iar = Inf(lr);
665 skipeolnflag = inpdotflag = true;
666 c = skipwhitespacessinglechar (s, ',');
667 if (inpdotflag) s.putback(c);
668 s >> dot;
669 lr = l_interval(dot);
670 Sar = Sup(lr);
671 lr = l_interval(Iar,Sar);
672 c = skipwhitespacessinglechar (s, ']');
673 if (inpdotflag) s.putback(c);
674 c = skipwhitespacessinglechar (s, ',');
675 if (inpdotflag) s.putback(c);
676
677 c = skipwhitespacessinglechar (s, '[');
678 if (inpdotflag) s.putback(c);
679// s >> RndDown >> Inf(a.im);
680 stagprec = StagPrec(a.im);
681 s >> dot;
682 li = l_interval(dot);
683 Iai = Inf(li);
684 skipeolnflag = inpdotflag = true;
685 c = skipwhitespacessinglechar (s, ',');
686 if (inpdotflag) s.putback(c);
687 s >> dot;
688 li = l_interval(dot);
689 Sai = Sup(li);
690 li = l_interval(Iai,Sai);
691
692 a = l_cinterval(lr,li);
693
694 if (!waseolnflag)
695 {
696 skipeolnflag = false, inpdotflag = true;
697 c = skipwhitespaces (s);
698 if (inpdotflag && c != ']')
699 s.putback(c);
700 }
701 if (!waseolnflag)
702 {
703 skipeolnflag = false, inpdotflag = true;
704 c = skipwhitespaces (s);
705 if (inpdotflag && c != ')')
706 s.putback(c);
707 }
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)"));
711
712 return s;
713}
714
715void operator >> (const std::string &s, l_cinterval &a)
716{
717 std::string r(s);
718 r >> a;
719}
720
721void operator >> (const char *s, l_cinterval &a)
722{
723 std::string r(s);
724 r >> a;
725}
726
727} // namespace cxsc
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
The Data Type cdotprecision.
Definition cdot.hpp:61
The Data Type cidotprecision.
Definition cidot.hpp:58
The Scalar Type cinterval.
Definition cinterval.hpp:55
cinterval(void) noexcept
Constructor of class cinterval.
Definition cinterval.hpp:64
cinterval & operator=(const real &) noexcept
Implementation of standard assigning operator.
Definition cinterval.inl:53
The Data Type dotprecision.
Definition dot.hpp:112
The Data Type idotprecision.
Definition idot.hpp:48
The Scalar Type interval.
Definition interval.hpp:55
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.
Definition l_real.hpp:78
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
const real MaxReal
Greatest representable floating-point number.
Definition real.cpp:65
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
Definition cimatrix.inl:730
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1071
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
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
Definition cimatrix.inl:731
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition imath.cpp:80