C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
interval.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: interval.cpp,v 1.38 2014/01/30 17:23:45 cxsc Exp $ */
25 
26 #include "interval.hpp"
27 #include "ioflags.hpp"
28 #include "idot.hpp"
29 #include "dot.hpp"
30 #include "RtsFunc.h" // b_shr1 in mid()
31 
32 // Auch fi_lib hat ein eigenes interval (wie RTS sein dotprecision)
33 //typedef struct fi_interval { double INF, SUP;} fi_interval;
34 // class interval;
35 #undef LINT_ARGS
36 #define CXSC_INCLUDE
37 #include <fi_lib.hpp>
38 
39 namespace cxsc {
40 
42 {
43  *this=rnd(a);
44 }
45 
47 {
48  return *this=rnd(a);
49 }
50 
51 interval::interval(const dotprecision & a) throw()
52 {
53  rnd(a,inf,sup);
54 }
55 
56 interval::interval(const dotprecision &a,const dotprecision &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
57 {
58  inf=rnd(a,RND_DOWN);
59  sup=rnd(b,RND_UP);
60  if(inf>sup)
61  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("interval::interval(const dotprecision &,const dotprecision &)"));
62 }
63 
64 interval::interval(const l_real &a,const l_real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
65 {
66  dotprecision dot(a);
67  inf=rnd(dot,RND_DOWN);
68  dot = b;
69  sup = rnd(dot,RND_UP);
70  if(inf>sup)
71  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("interval::interval(const l_real &,const l_real &)"));
72 }
73 
74 
76 {
77  rnd(a,inf,sup);
78  return *this;
79 }
80 
81 
82 // ---- Standardfunkt ---- (arithmetische Operatoren)
83 /*
84 interval operator+(const interval &a, const interval &b) throw() { return interval(adddown(a.inf,b.inf),addup(a.sup,b.sup)); }
85 interval operator+(const interval &a, const real &b) throw() { return interval(adddown(a.inf,b),addup(a.sup,b)); }
86 interval operator+(const real &a, const interval &b) throw() { return interval(adddown(a,b.inf),addup(a,b.sup)); }
87 
88 interval operator-(const interval &a, const interval &b) throw() { return interval(subdown(a.inf,b.inf),subup(a.sup,b.sup)); }
89 interval operator-(const interval &a, const real &b) throw() { return interval(subdown(a.inf,b),subup(a.sup,b)); }
90 interval operator-(const real &a, const interval &b) throw() { return interval(subdown(a,b.sup),subup(a,b.inf)); }
91 
92 interval operator*(const interval &a, const interval &b) throw() { return mul_ii (a, b); }
93 interval operator*(const interval &a, const real &b) throw() { return mul_id (a, *(double*)&b); }
94 interval operator*(const real &a, const interval &b) throw() { return mul_di (*(double*)&a, b); }
95 
96 interval operator/(const interval &a, const interval &b) throw() { return div_ii (a, b); }
97 interval operator/(const interval &a, const real &b) throw() { return div_id (a, *(double*)&b); }
98 interval operator/(const real &a, const interval &b) throw() { return div_di (*(double*)&a, b); }
99 */
100 
101 
102 bool operator ==(const interval &a,const dotprecision &r) throw() { return(r==a.inf && r==a.sup); }
103 bool operator !=(const interval &a,const dotprecision &r) throw() { return(r!=a.inf || r!=a.sup); }
104 bool operator ==(const dotprecision &r,const interval &a) throw() { return(r==a.inf && r==a.sup); }
105 bool operator !=(const dotprecision &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); }
106 
107 bool operator <=(const dotprecision &a,const interval &b) throw()
108 {
109  return(a>=b.inf && a<=b.sup);
110 }
111 bool operator >=(const dotprecision &a,const interval &b) throw()
112 {
113  return(a<=b.inf && a>=b.sup);
114 }
115 bool operator <(const dotprecision &a,const interval &b) throw()
116 {
117  return(a>b.inf && a<b.sup);
118 }
119 
120 bool operator <=(const interval &a,const dotprecision &b) throw()
121 {
122  return(a.inf>=b && a.sup<=b);
123 }
124 bool operator >=(const interval &a,const dotprecision &b) throw()
125 {
126  return(a.inf<=b && a.sup>=b);
127 }
128 bool operator >(const interval &a,const dotprecision &b) throw()
129 {
130  return(a.inf<b && a.sup>b);
131 }
132 
133 
134 
135 // ---- Ausgabefunkt. ---------------------------------------
136 
137 std::ostream & operator << (std::ostream &s, const interval& a) throw()
138 {
139  s << '[' << SaveOpt << RndDown
140  << a.inf << ',' << RndUp
141  << a.sup << RestoreOpt
142  << ']';
143  return s;
144 }
145 std::string & operator << (std::string &s, const interval &a) throw()
146 {
147  s+='[';
148  s << SaveOpt << RndDown
149  << a.inf;
150  s+=',';
151  s << RndUp
152  << a.sup << RestoreOpt;
153  s+=']';
154  return s;
155 }
156 
157 std::istream & operator >> (std::istream &s, interval &a) throw()
158 {
159  char c;
160 
161  skipeolnflag = inpdotflag = true;
162  c = skipwhitespacessinglechar (s, '[');
163  if (inpdotflag)
164  s.putback(c);
165 
166  s >> SaveOpt >> RndDown >> a.inf;
167 
168  skipeolnflag = inpdotflag = true;
169  c = skipwhitespacessinglechar (s, ',');
170  if (inpdotflag) s.putback(c);
171 
172  s >> RndUp >> a.sup >> RestoreOpt;
173 
174  if (!waseolnflag)
175  {
176  skipeolnflag = false, inpdotflag = true;
177  c = skipwhitespaces (s);
178  if (inpdotflag && c != ']')
179  s.putback(c);
180  }
181 
182  /*if (a.inf > a.sup) {
183  errmon (ERR_INTERVAL(EMPTY));
184  } */
185  return s;
186 }
187 
188 std::string & operator >> (std::string &s, interval &a) throw()
189 {
190  s = skipwhitespacessinglechar (s, '[');
191  s >> SaveOpt >> RndDown >> a.inf;
192  s = skipwhitespacessinglechar (s, ',');
193  s >> RndUp >> a.sup >> RestoreOpt;
194  s = skipwhitespaces (s);
195 
196  if (s[0] == ']')
197  s.erase(0,1);
198 
199  /*if (a.inf > a.sup) {
200  errmon (ERR_INTERVAL(EMPTY));
201  } */
202  return s;
203 }
204 
205 void operator >>(const std::string &s,interval &a) throw()
206 {
207  std::string r(s);
208  r>>a;
209 }
210 void operator >>(const char *s,interval &a) throw()
211 {
212  std::string r(s);
213  r>>a;
214 }
215 
216 real mid(const interval & a) throw()
217 {
218  dotprecision dot(a.inf);
219  dot += a.sup;
220 
221  if (dot != 0.0)
222  {
223  // Division nur bei ungleich 0.0
224  Dotprecision ptr = *dot.ptr();
225 
226  // --------------------------------------------------------
227  // Dividiere dot durch 2, mittels 1 Rechtsshift
228 
229  ptr[(a_intg)++ptr[A_END]] = ZERO;
230 
231  b_shr1 (&ptr[(a_intg)ptr[A_BEGIN]], (a_intg)(ptr[A_END]-ptr[A_BEGIN]+1));
232 
233  if (ptr[(a_intg)ptr[A_END]] == ZERO)
234  --ptr[A_END];
235 
236  if (ptr[(a_intg)ptr[A_BEGIN]] == ZERO)
237  ++ptr[A_BEGIN];
238 
239  // --------------------------------------------------------
240  }
241 
242  return rnd(dot);
243 }
244 
245 // for compatibility with CToolbox library - from former i_util.hpp
246 int in (const real& x, const interval& y ) // Contained-in relation
247  { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //----------------------
248 
249 int in ( const interval& x, const interval& y ) // Contained-in relation
250 { //----------------------
251  return ( (Inf(y) < Inf(x)) && (Sup(x) < Sup(y)) );
252 }
253 
254 // function in dot.hpp and dot.cpp
255 //void rnd ( const dotprecision& d, interval& x ) // Round dotprecision to interval
256 //{ //-------------------------------
257 // real Lower, Upper;
258 //
259 // rnd(d,Lower,Upper); // Rounding to interval bounds
260 // x = interval(Lower,Upper);
261 //}
262 
281 interval Blow ( const interval& x, const real& eps ) // Epsilon inflation
282 { //------------------
283  interval tmp;
284  tmp = (1.0+eps)*x - eps*x;
285  return interval(pred(Inf(tmp)),succ(Sup(tmp)));
286 }
287 
288 int Disjoint ( const interval& a, const interval& b ) // Test for disjointedness
289 { //------------------------
290  return (Inf(a) > Sup(b)) || (Inf(b) > Sup(a));
291 }
292 
293 real AbsMin ( const interval& x ) // Absolute minimum of
294 { // an interval
295  if ( in(0.0,x) ) //--------------------
296  return 0.0;
297  else if (Inf(x) > 0.0)
298  return Inf(x);
299  else
300  return -Sup(x);
301 }
302 
303 real AbsMax ( const interval& x ) // Absolute maximum of
304 { // an interval
305  real a, b; //--------------------
306 
307  a = abs(Inf(x));
308  b = abs(Sup(x));
309 
310  if (a > b)
311  return a;
312  else
313  return b;
314 }
315 
316 real RelDiam ( const interval& x ) // Relative diameter
317 { // of an interval
318  if ( in(0.0,x) ) //------------------
319  return diam(x);
320  else
321  return( divu(diam(x),AbsMin(x)) );
322 }
323 
324 //----------------------------------------------------------------------------
325 // Checks if the diameter of the interval 'x' is less or equal to 'n' ulps.
326 // Ulp is an abbreviation for units in the last place.
327 //----------------------------------------------------------------------------
335 int UlpAcc ( const interval& x, int n )
336 {
337  int i;
338  real Infimum;
339 
340  Infimum = Inf(x);
341  for (i = 1; i <= n; i++) Infimum = succ(Infimum);
342  return (Infimum >= Sup(x));
343 }
344 
345 // The following interval constants are optimal inclusions:
346 
347  const real Pi_Inf = 7074237752028440.0 / 2251799813685248.0;
348  const interval Pi_interval = interval(Pi_Inf,succ(Pi_Inf)); // Pi
349 
350  const real Pi2_Inf = 7074237752028440.0 / 1125899906842624.0;
351  const interval Pi2_interval = interval(Pi2_Inf,succ(Pi2_Inf)); // 2*Pi
352 
353  const real Pi3_Inf = 5305678314021330.0 / 562949953421312.0;
354  const interval Pi3_interval = interval(Pi3_Inf,succ(Pi3_Inf)); // 3*Pi
355 
356  const real Pid2_Inf = 7074237752028440.0 / 4503599627370496.0;
357  const interval Pid2_interval = interval(Pid2_Inf,succ(Pid2_Inf)); // Pi/2
358 
359  const real Pid3_Inf = 4716158501352293.0 / 4503599627370496.0;
360  const interval Pid3_interval = interval(Pid3_Inf,succ(Pid3_Inf)); // Pi/3
361 
362  const real Pid4_Inf = 7074237752028440.0 / 9007199254740992.0;
363  const interval Pid4_interval = interval(Pid4_Inf,succ(Pid4_Inf)); // Pi/4
364 
365  const real Pir_Inf = 5734161139222658.0 / 18014398509481984.0;
366  const interval Pir_interval = interval(Pir_Inf,succ(Pir_Inf)); // 1/Pi
367 
368  const real Pi2r_Inf = 5734161139222658.0 / 36028797018963968.0;
369  const interval Pi2r_interval = interval(Pi2r_Inf,succ(Pi2r_Inf));
370  // 1/(2*Pi)
371 
372  const real Pip2_Inf = 5556093337880030.0 / 562949953421312.0;
373  const interval Pip2_interval = interval(Pip2_Inf,succ(Pip2_Inf)); // Pi^2
374 
375  const real SqrtPi_Inf = 7982422502469482.0 / 4503599627370496.0;
376  const interval SqrtPi_interval = interval(SqrtPi_Inf,succ(SqrtPi_Inf));
377  // sqrt(Pi)
378 
379  const real Sqrt2Pi_Inf = 5644425081792261.0 / 2251799813685248.0;
380  const interval Sqrt2Pi_interval = interval(Sqrt2Pi_Inf,succ(Sqrt2Pi_Inf));
381  // sqrt(2Pi)
382 
383  const real SqrtPir_Inf = 5081767996463981.0 / 9007199254740992.0;
384  const interval SqrtPir_interval = interval(SqrtPir_Inf,succ(SqrtPir_Inf));
385  // 1/sqrt(Pi)
386 
387  const real Sqrt2Pir_Inf = 7186705221432912.0 / 18014398509481984.0;
388  const interval Sqrt2Pir_interval=interval(Sqrt2Pir_Inf,succ(Sqrt2Pir_Inf));
389  // 1/sqrt(2Pi)
390 
391  const real Sqrt2_Inf = 6369051672525772.0 / 4503599627370496.0;
392  const interval Sqrt2_interval=interval(Sqrt2_Inf,succ(Sqrt2_Inf));
393  // sqrt(2)
394 
395  const real Sqrt5_Inf = 5035177455121575.0 / 2251799813685248.0;
396  const interval Sqrt5_interval=interval(Sqrt5_Inf,succ(Sqrt5_Inf));
397  // sqrt(5)
398 
399  const real Sqrt7_Inf = 5957702309312745.0 / 2251799813685248.0;
400  const interval Sqrt7_interval=interval(Sqrt7_Inf,succ(Sqrt7_Inf));
401  // sqrt(7)
402 
403  const real Sqrt2r_Inf = 6369051672525772.0 / 9007199254740992.0;
404  const interval Sqrt2r_interval=interval(Sqrt2r_Inf,succ(Sqrt2r_Inf));
405  // 1/sqrt(2)
406 
407  const real Sqrt3_Inf = 7800463371553962.0 / 4503599627370496.0;
408  const interval Sqrt3_interval=interval(Sqrt3_Inf,succ(Sqrt3_Inf));
409  // sqrt(3)
410 
411  const real Sqrt3d2_Inf = 7800463371553962.0 / 9007199254740992.0;
412  const interval Sqrt3d2_interval=interval(Sqrt3d2_Inf,succ(Sqrt3d2_Inf));
413  // sqrt(3)/2
414 
415  const real Sqrt3r_Inf = 5200308914369308.0 / 9007199254740992.0;
416  const interval Sqrt3r_interval=interval(Sqrt3r_Inf,succ(Sqrt3r_Inf));
417  // 1/sqrt(3)
418 
419  const real Ln2_Inf = 6243314768165359.0 / 9007199254740992.0;
420  const interval Ln2_interval=interval(Ln2_Inf,succ(Ln2_Inf)); // ln(2)
421 
422  const real Ln2r_Inf = 6497320848556798.0 / 4503599627370496.0;
423  const interval Ln2r_interval=interval(Ln2r_Inf,succ(Ln2r_Inf)); // 1/ln(2)
424 
425  const real Ln10_Inf = 5184960683398421.0 / 2251799813685248.0;
426  const interval Ln10_interval=interval(Ln10_Inf,succ(Ln10_Inf)); // ln(10)
427 
428  const real Ln10r_Inf = 7823553867474190.0 / 18014398509481984.0;
429  const interval Ln10r_interval=interval(Ln10r_Inf,succ(Ln10r_Inf));
430  // 1/ln(10)
431 
432  const real LnPi_Inf = 5155405087351229.0 / 4503599627370496.0;
433  const interval LnPi_interval=interval(LnPi_Inf,succ(LnPi_Inf)); // ln(Pi)
434 
435  const real Ln2Pi_Inf = 8277062471433908.0 / 4503599627370496.0;
436  const interval Ln2Pi_interval=interval(Ln2Pi_Inf,succ(Ln2Pi_Inf));
437  // ln(2Pi)
438 
439  const real E_Inf = 6121026514868073.0 / 2251799813685248.0;
440  const interval E_interval=interval(E_Inf,succ(E_Inf)); // e
441 
442  const real Er_Inf = 6627126856707895.0 / 18014398509481984.0;
443  const interval Er_interval=interval(Er_Inf,succ(Er_Inf)); // 1/e
444 
445  const real Ep2_Inf = 8319337573440941.0 / 1125899906842624.0;
446  const interval Ep2_interval=interval(Ep2_Inf,succ(Ep2_Inf)); // e^2
447 
448  const real Ep2r_Inf = 4875967449235915.0 / 36028797018963968.0;
449  const interval Ep2r_interval=interval(Ep2r_Inf,succ(Ep2r_Inf)); // 1/e^2
450 
451  const real EpPi_Inf = 6513525919879993.0 / 281474976710656.0;
452  const interval EpPi_interval=interval(EpPi_Inf,succ(EpPi_Inf)); // e^(Pi)
453 
454  const real Ep2Pi_Inf = 4710234414611993.0 / 8796093022208.0;
455  const interval Ep2Pi_interval=interval(Ep2Pi_Inf,succ(Ep2Pi_Inf));
456  // e^(2Pi)
457 
458  const real EpPid2_Inf = 5416116035097439.0 / 1125899906842624.0;
459  const interval EpPid2_interval=interval(EpPid2_Inf,succ(EpPid2_Inf));
460  // e^(Pi/2)
461 
462  const real EpPid4_Inf = 4938827609611434.0 / 2251799813685248.0;
463  const interval EpPid4_interval=interval(EpPid4_Inf,succ(EpPid4_Inf));
464  // e^(Pi/4)
465 
466 
467 } // namespace cxsc
468 
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.
Definition: cinterval.cpp:665
const interval Ln2_interval
Enclosure-Interval for .
Definition: interval.cpp:420
const interval Pi3_interval
Enclosure-Interval for .
Definition: interval.cpp:354
const interval Ln10_interval
Enclosure-Interval for .
Definition: interval.cpp:426
const interval Pi2_interval
Enclosure-Interval for .
Definition: interval.cpp:351
The Data Type idotprecision.
Definition: idot.hpp:47
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:77
The Data Type dotprecision.
Definition: dot.hpp:111
const interval Pip2_interval
Enclosure-Interval for .
Definition: interval.cpp:373
const interval Sqrt2_interval
Enclosure-Interval for .
Definition: interval.cpp:392
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
const interval Sqrt7_interval
Enclosure-Interval for .
Definition: interval.cpp:400
const interval EpPid4_interval
Enclosure-Interval for .
Definition: interval.cpp:463
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
Definition: cinterval.cpp:654
cvector mid(const cimatrix_subv &mv)
Returns the middle of the matrix.
Definition: cimatrix.inl:739
The Scalar Type interval.
Definition: interval.hpp:54
const interval EpPid2_interval
Enclosure-Interval for .
Definition: interval.cpp:459
real AbsMin(const interval &x)
Computes the smallest absolute value .
Definition: interval.cpp:293
const interval Sqrt5_interval
Enclosure-Interval for .
Definition: interval.cpp:396
int UlpAcc(const interval &x, int n)
Checks if the diameter of the interval is ulps.
Definition: interval.cpp:335
real RelDiam(const interval &x)
Computes the relative diameter .
Definition: interval.cpp:316
const interval Pid2_interval
Enclosure-Interval for .
Definition: interval.cpp:357
const interval Sqrt2Pir_interval
Enclosure-Interval for .
Definition: interval.cpp:388
const interval Er_interval
Enclosure-Interval for .
Definition: interval.cpp:443
const interval Pid3_interval
Enclosure-Interval for .
Definition: interval.cpp:360
const interval Sqrt3r_interval
Enclosure-Interval for .
Definition: interval.cpp:416
const interval Sqrt3_interval
Enclosure-Interval for .
Definition: interval.cpp:408
const interval LnPi_interval
Enclosure-Interval for .
Definition: interval.cpp:433
const interval Ln2r_interval
Enclosure-Interval for .
Definition: interval.cpp:423
const interval Pid4_interval
Enclosure-Interval for .
Definition: interval.cpp:363
const interval Pir_interval
Enclosure-Interval for .
Definition: interval.cpp:366
const interval Ep2Pi_interval
Enclosure-Interval for .
Definition: interval.cpp:455
const interval EpPi_interval
Enclosure-Interval for .
Definition: interval.cpp:452
const interval Pi2r_interval
Enclosure-Interval for .
Definition: interval.cpp:369
const interval SqrtPir_interval
Enclosure-Interval for .
Definition: interval.cpp:384
cvector diam(const cimatrix_subv &mv)
Returns the diameter of the matrix.
Definition: cimatrix.inl:738
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
Definition: interval.cpp:288
real AbsMax(const interval &x)
Computes the greatest absolute value .
Definition: interval.cpp:303
interval()
Constructor of class interval.
Definition: interval.hpp:64
interval & operator=(const real &a)
Implementation of standard assigning operator.
Definition: interval.inl:52
const interval SqrtPi_interval
Enclosure-Interval for .
Definition: interval.cpp:376
const interval Ep2_interval
Enclosure-Interval for .
Definition: interval.cpp:446
const interval Sqrt3d2_interval
Enclosure-Interval for .
Definition: interval.cpp:412
const interval Pi_interval
Enclosure-Interval for .
Definition: interval.cpp:348
const interval Sqrt2Pi_interval
Enclosure-Interval for .
Definition: interval.cpp:380
const interval Ln2Pi_interval
Enclosure-Interval for .
Definition: interval.cpp:436
const interval Sqrt2r_interval
Enclosure-Interval for .
Definition: interval.cpp:404
const interval E_interval
Enclosure-Interval for .
Definition: interval.cpp:440
const interval Ep2r_interval
Enclosure-Interval for .
Definition: interval.cpp:449
The Scalar Type real.
Definition: real.hpp:113
const interval Ln10r_interval
Enclosure-Interval for .
Definition: interval.cpp:429
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737