C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
interval.inl
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.inl,v 1.31 2014/01/30 17:23:45 cxsc Exp $ */
25 
26 namespace cxsc {
27 
28 // ---- Konstruktoren ----
29 
30 inline interval::interval(const real &a,const real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
31  : inf(a), sup(b)
32 {
33  if (a > b)
34  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval::interval(const real &a,const real &b)"));
35 }
36 
37 /*inline interval::interval(int a,int b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
38  : inf(a), sup(b)
39 {
40  if (a > b)
41  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval::interval(int a,int b)"));
42 }
43 
44 inline interval::interval(const double & a,const double & b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
45  : inf(a), sup(b)
46 {
47  if (a > b)
48  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval::interval(const double & a,const double & b)"));
49 }
50 */
51 
53 {
54  inf = sup = a;
55  return *this;
56 }
57 
58 
59 // ---- Typwandlungen ----
60 
66 inline interval _unchecked_interval(const real& a, const real& b)
67 {
68  interval tmp;
69  tmp.inf = a;
70  tmp.sup = b;
71  return tmp;
72 }
73 
74 
75 // ---- Standardfunkt ---- (arithmetische Operatoren)
76 
77 inline interval operator-(const interval &a) throw() { return interval (-a.sup, -a.inf); }
78 inline interval operator+(const interval &a) throw() { return a; }
79 
80 inline interval & operator +=(interval &a,const interval &b) throw() { return a=a+b; }
81 inline interval & operator +=(interval &a,const real &b) throw() { return a=a+b; }
82 inline interval & operator -=(interval &a,const interval &b) throw() { return a=a-b; }
83 inline interval & operator -=(interval &a,const real &b) throw() { return a=a-b; }
84 inline interval & operator *=(interval &a,const interval &b) throw() { return a=a*b; }
85 inline interval & operator *=(interval &a,const real &b) throw() { return a=a*b; }
86 inline interval & operator /=(interval &a,const interval &b) throw() { return a=a/b; }
87 inline interval & operator /=(interval &a,const real &b) throw() { return a=a/b; }
88 
89 
90 inline interval operator |(const interval &a,const interval &b) throw()
91 {
92  return interval((a.inf<b.inf)?a.inf:b.inf,(a.sup>b.sup)?a.sup:b.sup);
93 }
94 inline interval operator &(const interval &a,const interval &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
95 {
96  return interval((a.inf>b.inf)?a.inf:b.inf,(a.sup<b.sup)?a.sup:b.sup);
97 }
98 inline interval operator |(const real &a,const interval &b) throw()
99 {
100  return interval((a<b.inf)?a:b.inf,(a>b.sup)?a:b.sup);
101 }
102 inline interval operator &(const real &a,const interval &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
103 {
104  return interval((a>b.inf)?a:b.inf,(a<b.sup)?a:b.sup);
105 }
106 inline interval operator |(const interval &a,const real &b) throw()
107 {
108  return interval((a.inf<b)?a.inf:b,(a.sup>b)?a.sup:b);
109 }
110 inline interval operator |(const real &a,const real &b) throw()
111 {
112  if(a>b) return interval(b,a);
113  else return interval(a,b);
114 }
115 inline interval operator &(const interval &a,const real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
116 {
117  return interval((a.inf>b)?a.inf:b,(a.sup<b)?a.sup:b);
118 }
119 inline interval & operator |=(interval &a,const interval &b) throw()
120 {
121  a.inf=(a.inf<b.inf)?a.inf:b.inf,a.sup=(a.sup>b.sup)?a.sup:b.sup;
122  return a;
123 }
124 inline interval & operator &=(interval &a,const interval &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
125 {
126  a.inf=(a.inf>b.inf)?a.inf:b.inf,a.sup=(a.sup<b.sup)?a.sup:b.sup;
127  if(a.inf>a.sup)
128  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval & operator &=(interval &a,const interval &b)"));
129  return a;
130 }
131 inline interval & operator |=(interval &a,const real &b) throw()
132 {
133  a.inf=(a.inf<b)?a.inf:b,a.sup=(a.sup>b)?a.sup:b;
134  return a;
135 }
136 inline interval & operator &=(interval &a,const real &b) throw(ERROR_INTERVAL_EMPTY_INTERVAL)
137 {
138  a.inf=(a.inf>b)?a.inf:b,a.sup=(a.sup<b)?a.sup:b;
139  if(a.inf>a.sup)
140  cxscthrow(ERROR_INTERVAL_EMPTY_INTERVAL("inline interval & operator &=(interval &a,const real &b)"));
141  return a;
142 }
143 
144 // --- Vergleichsoperationen ----
145 inline bool operator ==(const interval &a,const interval &b) throw() { return(a.inf==b.inf && a.sup==b.sup); }
146 inline bool operator !=(const interval &a,const interval &b) throw() { return(a.inf!=b.inf || a.sup!=b.sup); }
147 inline bool operator ==(const real &r,const interval &a) throw() { return(r==a.inf && r==a.sup); }
148 inline bool operator !=(const real &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); }
149 inline bool operator ==(const interval &a,const real &r) throw() { return(r==a.inf && r==a.sup); }
150 inline bool operator !=(const interval &a,const real &r) throw() { return(r!=a.inf || r!=a.sup); }
151 
152 inline bool operator ==(const int &r,const interval &a) throw() { return(r==a.inf && r==a.sup); }
153 inline bool operator !=(const int &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); }
154 inline bool operator ==(const interval &a,const int &r) throw() { return(r==a.inf && r==a.sup); }
155 inline bool operator !=(const interval &a,const int &r) throw() { return(r!=a.inf || r!=a.sup); }
156 
157 inline bool operator ==(const long &r,const interval &a) throw() { return(r==a.inf && r==a.sup); }
158 inline bool operator !=(const long &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); }
159 inline bool operator ==(const interval &a,const long &r) throw() { return(r==a.inf && r==a.sup); }
160 inline bool operator !=(const interval &a,const long &r) throw() { return(r!=a.inf || r!=a.sup); }
161 
162 inline bool operator ==(const double &r,const interval &a) throw() { return(r==a.inf && r==a.sup); }
163 inline bool operator !=(const double &r,const interval &a) throw() { return(r!=a.inf || r!=a.sup); }
164 inline bool operator ==(const interval &a,const double &r) throw() { return(r==a.inf && r==a.sup); }
165 inline bool operator !=(const interval &a,const double &r) throw() { return(r!=a.inf || r!=a.sup); }
166 
167 // --- Mengenvergleiche ---
168 // <,>,...
169 inline bool operator <=(const interval &a,const interval &b) throw()
170 {
171  return(a.inf>=b.inf && a.sup<=b.sup);
172 }
173 inline bool operator >=(const interval &a,const interval &b) throw()
174 {
175  return(a.inf<=b.inf && a.sup>=b.sup);
176 }
177 inline bool operator <(const interval &a,const interval &b) throw()
178 {
179  return(a.inf>b.inf && a.sup<b.sup);
180 }
181 inline bool operator >(const interval &a,const interval &b) throw()
182 {
183  return(a.inf<b.inf && a.sup>b.sup);
184 }
185 
186 inline bool operator <=(const real &a,const interval &b) throw()
187 {
188  return(a>=b.inf && a<=b.sup);
189 }
190 inline bool operator >=(const real &a,const interval &b) throw()
191 {
192  return(a<=b.inf && a>=b.sup);
193 }
194 inline bool operator <(const real &a,const interval &b) throw()
195 {
196  return(a>b.inf && a<b.sup);
197 }
198 
199 inline bool operator <=(const interval &a,const real &b) throw()
200 {
201  return(a.inf>=b && a.sup<=b);
202 }
203 inline bool operator >=(const interval &a,const real &b) throw()
204 {
205  return(a.inf<=b && a.sup>=b);
206 }
207 inline bool operator >(const interval &a,const real &b) throw()
208 {
209  return(a.inf<b && a.sup>b);
210 }
211 
212 inline bool operator !(const interval &a) throw() { return (a.inf <= 0.0 && a.sup >= 0.0); }
213 
214 inline real & Inf (interval& a) throw() { return a.inf; }
215 inline const real & Inf (const interval &a) throw() { return a.inf; }
216 inline real & Sup (interval& a) throw() { return a.sup; }
217 inline const real & Sup (const interval &a) throw() { return a.sup; }
218 
219 inline interval& SetInf (interval& a, const real& b) throw() {a.inf=b; return a;}
220 inline interval& SetSup (interval& a, const real& b) throw() {a.sup=b; return a;}
221 inline interval& UncheckedSetInf (interval& a, const real& b) throw() { a.inf=b; return a;}
222 inline interval& UncheckedSetSup (interval& a, const real& b) throw() { a.sup=b; return a;}
223 
224 inline bool IsEmpty(const interval &a) throw() { return (a.inf>a.sup); }
225 
226 inline interval abs(const interval &a) throw()
227 {
228  real h1 = abs(a.inf);
229  real h2 = abs(a.sup);
230 
231  if (IsEmpty(a)) return a;
232  if (!a)
233  return interval(0.0, (h1 > h2) ? h1 : h2);
234  if (h1 > h2)
235  return interval(h2, h1);
236 
237  return interval(h1, h2);
238 }
239 
240 inline real diam(const interval & a) throw()
241 {
242  return subup(a.sup,a.inf);
243 }
244 
245 inline real Mid(const interval & a) throw()
246 {
247  return addd( a.inf, subd(0.5*a.sup,0.5*a.inf) );
248 }
249 
250 inline void times2pown(interval& x, const int& n) throw()
251 {
252  real r1,r2;
253  int j;
254  r1 = x.inf; r2 = x.sup;
255  j = expo(r1) + n;
256  if (j >= -1021) r1 = comp(mant(r1),j);
257  else
258  {
259  j += 1021;
260  r1 = comp(mant(r1), -1021);
261  if (j<-53)
262  {
263  if (sign(r1)>=0) r1 = 0;
264  else r1 = -minreal;
265  } else
266  r1 = muld(r1,comp(0.5,j+1));
267  }
268  j = expo(r2) + n;
269  if (j >= -1021) r2 = comp(mant(r2),j);
270  else
271  {
272  j += 1021;
273  r2 = comp(mant(r2), -1021);
274  if (j<-53)
275  {
276  if (sign(r2)>0) r2 = minreal;
277  else r2 = 0;
278  } else r2 = mulu(r2,comp(0.5,j+1));
279  }
280  x = _interval(r1,r2);
281 } // times2pown(...)
282 
283 interval operator+ (const interval& a, const interval& b) throw()
284  { return interval (adddown(a.inf,b.inf),addup(a.sup,b.sup)); }
285 interval operator+ (const interval& a, const real& b) throw()
286  { return interval (adddown(a.inf,b),addup(a.sup,b)); }
287 interval operator+ (const real& a, const interval& b) throw()
288  { return interval (adddown(a,b.inf),addup(a,b.sup)); }
289 
290 interval operator- (const interval& a, const interval& b) throw()
291  { return interval ( subdown(a.inf,b.sup), subup(a.sup,b.inf)); }
292 interval operator- (const interval& a, const real& b) throw()
293  { return interval ( subdown(a.inf,b), subup(a.sup,b)); }
294 interval operator- (const real& a, const interval& b) throw()
295  { return interval ( subdown(a,b.sup), subup(a,b.inf)); }
296 
297  //------------------------------------------------------------------------
298  // a * b = ueber Entscheidungstabelle :
299  //
300  // bi,bs >= 0 bi < 0, bs >=0 bi,bs < 0
301  // ----------------+------------------+------------------+----------------
302  // ai,as >= 0 I ai*bi, as*bs I as*bi, as*bs I as*bi, ai*bs
303  // ----------------+------------------+------------------+----------------
304  // ai < 0, as >= 0 I ai*bs, as*bs I .... I as*bi, ai*bi
305  // ----------------+------------------+------------------+----------------
306  // ai,as < 0 I ai*bs, as*bi I ai*bs, ai*bi I as*bs, ai*bi
307  // ----------------+------------------+------------------+----------------
308  //
309  // .... : min(ai*bs, as*bi), max(ai*ai, as*as)
310  //
311 interval operator *(const interval& a, const interval& b) throw()
312 {
313  interval tmp;
314 
315  if (sign(a.inf) >= 0)
316  { // 1. Zeile der Entscheidungstabelle
317  if (sign(b.inf) >= 0)
318  { // 1. Spalte: [ai*bi, as*bs]
319  tmp.inf = multdown(a.inf,b.inf);
320  tmp.sup = multup(a.sup,b.sup);
321  } else if (sign(b.sup) >= 0)
322  { // 2. Spalte: [as*bi, as*bs]
323  tmp.inf = multdown(a.sup,b.inf);
324  tmp.sup = multup(a.sup,b.sup);
325  } else
326  { // 3. Spalte: [as*bi, ai*bs]
327  tmp.inf = multdown(a.sup,b.inf);
328  tmp.sup = multup(a.inf,b.sup);
329  }
330  } else if (sign(a.sup) >= 0)
331  { // 2. Zeile der Entscheidungstabelle
332  if (sign(b.inf) >= 0)
333  { // 1. Spalte: [ai*bs, as*bs]
334  tmp.inf = multdown(a.inf,b.sup);
335  tmp.sup = multup(a.sup,b.sup);
336  } else if (sign(b.sup) >= 0)
337  { // 2. Spalte: [min(ai*bs, as*bi),
338  real hlp; // max(ai*ai, as*as)]
339 
340  tmp.inf = multdown(a.inf,b.sup);
341  hlp = multdown(a.sup,b.inf);
342 
343  if (hlp < tmp.inf)
344  tmp.inf = hlp;
345 
346  tmp.sup = multup(a.inf,b.inf);
347  hlp = multup(a.sup,b.sup);
348 
349  if (hlp > tmp.sup)
350  tmp.sup = hlp;
351  } else
352  { // 3. Spalte: [as*bi, ai*bi]
353  tmp.inf = multdown(a.sup,b.inf);
354  tmp.sup = multup(a.inf,b.inf);
355  }
356  } else
357  { // 3. Zeile der Entscheidungstabelle
358  if (sign(b.inf) >= 0)
359  { // 1. Spalte: [ai*bs, as*bi]
360  tmp.inf = multdown(a.inf,b.sup);
361  tmp.sup = multup(a.sup,b.inf);
362  } else if (sign(b.sup) >= 0)
363  { // 2. Spalte: [ai*bs, ai*bi]
364  tmp.inf = multdown(a.inf,b.sup);
365  tmp.sup = multup(a.inf,b.inf);
366  } else
367  { // 3. Spalte: [as*bs, ai*bi]
368  tmp.inf = multdown(a.sup,b.sup);
369  tmp.sup = multup(a.inf,b.inf);
370  }
371  }
372  return tmp;
373 }
374 
375  //------------------------------------------------------------------------
376  // a / b = ueber Entscheidungstabelle :
377  //
378  // bi,bs > 0 bi,bs < 0
379  // ----------------+------------------+------------------
380  // ai,as >= 0 I ai/bs, as/bi I as/bs, ai/bi
381  // ----------------+------------------+------------------
382  // ai < 0, as >= 0 I ai/bi, as/bi I as/bs, ai/bs
383  // ----------------+------------------+------------------
384  // ai,as < 0 I ai/bi, as/bs I as/bi, ai/bs
385  // ----------------+------------------+------------------
386  //
387 interval operator/ (const interval& a, const interval& b) throw(DIV_BY_ZERO)
388 {
389  interval tmp;
390 
391  if ((sign(b.inf) <= 0) && (sign(b.sup) >= 0))
392  cxscthrow(DIV_BY_ZERO("interval::interval operator/(const interval&,const interval&)"));
393 
394  if (sign(a.inf) >= 0)
395  { // 1. Zeile der Entscheidungstabelle
396  if (sign(b.inf) > 0)
397  { // 1. Spalte: [ai/bs, as/bi]
398  tmp.inf = divdown(a.inf,b.sup);
399  tmp.sup = divup(a.sup,b.inf);
400  } else
401  { // 2. Spalte: [as/bs, ai/bi]
402  tmp.inf = divdown(a.sup,b.sup);
403  tmp.sup = divup(a.inf,b.inf);
404  }
405  } else if (sign(a.sup) >= 0)
406  { // 2. Zeile der Entscheidungstabelle
407  if (sign(b.inf) > 0)
408  { // 1. Spalte: [ai/bi, as/bi]
409  tmp.inf = divdown(a.inf,b.inf);
410  tmp.sup = divup(a.sup,b.inf);
411  } else
412  { // 2. Spalte: [as/bs, ai/bs]
413  tmp.inf = divdown(a.sup,b.sup);
414  tmp.sup = divup(a.inf,b.sup);
415  }
416  } else
417  { // 3. Zeile der Entscheidungstabelle
418  if (sign(b.inf) > 0)
419  { // 1. Spalte: [ai/bi, as/bs]
420  tmp.inf = divdown(a.inf,b.inf);
421  tmp.sup = divup(a.sup,b.sup);
422  } else
423  { // 2. Spalte: [as/bi, ai/bs]
424  tmp.inf = divdown(a.sup,b.inf);
425  tmp.sup = divup(a.inf,b.sup);
426  }
427  }
428 
429  return tmp;
430 }
431 
432 interval operator* (const real& a, const interval& b) throw()
433 {
434  interval tmp;
435  if (sign(a) == 0)
436  {
437  tmp.inf = 0.0;
438  tmp.sup = 0.0;
439  } else if (sign(a) > 0)
440  {
441  tmp.inf = multdown(a,b.inf);
442  tmp.sup = multup(a,b.sup);
443  } else // if (sign(a) < 0)
444  {
445  tmp.inf = multdown(a,b.sup);
446  tmp.sup = multup(a,b.inf);
447  }
448  return tmp;
449 }
450 
451 interval operator* (const interval& a, const real& b) throw()
452 {
453  interval tmp;
454  if (sign(b) == 0)
455  {
456  tmp.inf = 0.0;
457  tmp.sup = 0.0;
458  } else if (sign(b) > 0)
459  {
460  tmp.inf = multdown(a.inf,b);
461  tmp.sup = multup(a.sup,b);
462  } else // if (sign(b) < 0)
463  {
464  tmp.inf = multdown(a.sup,b);
465  tmp.sup = multup(a.inf,b);
466  }
467  return tmp;
468 }
469 
470 interval operator/ (const real& a, const interval& b) throw() { return (interval(a) / b); }
471 interval operator/ (const interval& a, const real& b) throw() { return (a / interval(b)); }
472 
473 
474 } // namespace cxsc
475 
cimatrix & operator/=(cimatrix &m, const cinterval &c)
Implementation of division and allocation operation.
Definition: cimatrix.inl:1623
interval _unchecked_interval(const real &a, const real &b)
Definition: interval.inl:66
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
civector operator*(const cimatrix_subv &rv, const cinterval &s)
Implementation of multiplication operation.
Definition: cimatrix.inl:731
The Scalar Type interval.
Definition: interval.hpp:54
cvector diam(const cimatrix_subv &mv)
Returns the diameter of the matrix.
Definition: cimatrix.inl:738
interval()
Constructor of class interval.
Definition: interval.hpp:64
interval & operator=(const real &a)
Implementation of standard assigning operator.
Definition: interval.inl:52
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc)
Implementation of standard algebraic addition and allocation operation.
Definition: cdot.inl:251
void times2pown(cinterval &x, int n)
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
const real minreal
Smallest positive denormalized representable floating-point number.
Definition: real.cpp:63
civector operator/(const cimatrix_subv &rv, const cinterval &s)
Implementation of division operation.
Definition: cimatrix.inl:730
cimatrix & operator*=(cimatrix &m, const cinterval &c)
Implementation of multiplication and allocation operation.
Definition: cimatrix.inl:1605
The Scalar Type real.
Definition: real.hpp:113
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737