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 
32 namespace cxsc {
33 
34 
35 #define CXSC_Zero 0.0
36 
37 cinterval::cinterval(const l_cinterval & a) noexcept
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;
50 return *this = cinterval(u,v);
51 }
52 
53 l_cinterval::l_cinterval(const dotprecision &a) noexcept : re(a),im(0) {}
54 l_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 
60 l_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!)
88 void product(const l_real& a, const l_real& b, const l_real& c,
89  const l_real& d, int& ex, l_interval& res);
90 void product(const l_real& c, const l_real& d, int& ex, l_interval& res);
91 l_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 
97 static const int max_expo = 1020, max_expo1 = 1023;
98 
99 // optimale komplexe Intervalldivision
100 
101 bool cxsc_l_complex_division_p[5];
102 
103 l_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 
119 static 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 
288 l_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 
295 l_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 
302 l_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 
511 l_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 
545 l_interval abs(const l_cinterval &a) noexcept
546 {
547  return sqrtx2y2(a.re,a.im);
548 }
549 
550 
551 // ---- Ausgabefunkt. ---------------------------------------
552 
553 std::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 
562 std::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 
575 std::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 
635 std::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 
715 void operator >> (const std::string &s, l_cinterval &a)
716 {
717  std::string r(s);
718  r >> a;
719 }
720 
721 void 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.
Definition: l_cinterval.hpp:54
l_cinterval(void) noexcept
Constructor of class l_cinterval.
Definition: l_cinterval.hpp:63
The Multiple-Precision Data Type l_interval.
Definition: l_interval.hpp:72
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