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 
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) throw() : re(a),im(0) {}
54 l_cinterval::l_cinterval(const idotprecision &a) throw() : re(a),im(0) {}
56  throw() : 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) throw()
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) throw();
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  throw(DIV_BY_ZERO)
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) throw()
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) throw()
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) throw()
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  throw(EMPTY_INTERVAL)
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  throw(EMPTY_INTERVAL)
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) throw(EMPTY_INTERVAL)
716 {
717  std::string r(s);
718  r >> a;
719 }
720 
721 void operator >> (const char *s, l_cinterval &a) throw(EMPTY_INTERVAL)
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 
cxsc::cinterval::cinterval
cinterval(void)
Constructor of class cinterval.
Definition: cinterval.hpp:64
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:55
cxsc::sqrt1px2
cinterval sqrt1px2(const cinterval &z)
Calculates .
Definition: cimath.cpp:1071
cxsc::idotprecision
The Data Type idotprecision.
Definition: idot.hpp:48
cxsc::l_cinterval::l_cinterval
l_cinterval(void)
Constructor of class l_cinterval.
Definition: l_cinterval.hpp:63
cxsc::abs
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::operator*
civector operator*(const cimatrix_subv &rv, const cinterval &s)
Implementation of multiplication operation.
Definition: cimatrix.inl:731
cxsc::dotprecision
The Data Type dotprecision.
Definition: dot.hpp:112
cxsc::l_interval
The Multiple-Precision Data Type l_interval.
Definition: l_interval.hpp:72
cxsc::times2pown
void times2pown(cinterval &x, int n)
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cxsc::sqrtx2y2
interval sqrtx2y2(const interval &x, const interval &y)
Calculates .
Definition: imath.cpp:80
cxsc::cidotprecision
The Data Type cidotprecision.
Definition: cidot.hpp:58
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::MaxReal
const real MaxReal
Greatest representable floating-point number.
Definition: real.cpp:65
cxsc::l_real
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:78
cxsc::l_cinterval
The Multiple-Precision Data Type l_cinterval.
Definition: l_cinterval.hpp:54
cxsc::cdotprecision
The Data Type cdotprecision.
Definition: cdot.hpp:61
cxsc::cinterval::operator=
cinterval & operator=(const real &)
Implementation of standard assigning operator.
Definition: cinterval.inl:53
cxsc::operator/
civector operator/(const cimatrix_subv &rv, const cinterval &s)
Implementation of division operation.
Definition: cimatrix.inl:730
cxsc::cinterval
The Scalar Type cinterval.
Definition: cinterval.hpp:55