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
26namespace cxsc {
27
28// ---- Konstruktoren ----
29
30inline interval::interval(const real &a,const real &b)
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)
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
44inline interval::interval(const double & a,const double & b)
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
66inline 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
77inline interval operator-(const interval &a) noexcept { return interval (-a.sup, -a.inf); }
78inline interval operator+(const interval &a) noexcept { return a; }
79
80inline interval & operator +=(interval &a,const interval &b) noexcept { return a=a+b; }
81inline interval & operator +=(interval &a,const real &b) noexcept { return a=a+b; }
82inline interval & operator -=(interval &a,const interval &b) noexcept { return a=a-b; }
83inline interval & operator -=(interval &a,const real &b) noexcept { return a=a-b; }
84inline interval & operator *=(interval &a,const interval &b) noexcept { return a=a*b; }
85inline interval & operator *=(interval &a,const real &b) noexcept { return a=a*b; }
86inline interval & operator /=(interval &a,const interval &b) noexcept { return a=a/b; }
87inline interval & operator /=(interval &a,const real &b) noexcept { return a=a/b; }
88
89
90inline interval operator |(const interval &a,const interval &b) noexcept
91{
92 return interval((a.inf<b.inf)?a.inf:b.inf,(a.sup>b.sup)?a.sup:b.sup);
93}
94inline interval operator &(const interval &a,const interval &b)
95{
96 return interval((a.inf>b.inf)?a.inf:b.inf,(a.sup<b.sup)?a.sup:b.sup);
97}
98inline interval operator |(const real &a,const interval &b) noexcept
99{
100 return interval((a<b.inf)?a:b.inf,(a>b.sup)?a:b.sup);
101}
102inline interval operator &(const real &a,const interval &b)
103{
104 return interval((a>b.inf)?a:b.inf,(a<b.sup)?a:b.sup);
105}
106inline interval operator |(const interval &a,const real &b) noexcept
107{
108 return interval((a.inf<b)?a.inf:b,(a.sup>b)?a.sup:b);
109}
110inline interval operator |(const real &a,const real &b) noexcept
111{
112 if(a>b) return interval(b,a);
113 else return interval(a,b);
114}
115inline interval operator &(const interval &a,const real &b)
116{
117 return interval((a.inf>b)?a.inf:b,(a.sup<b)?a.sup:b);
118}
119inline interval & operator |=(interval &a,const interval &b) noexcept
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}
124inline interval & operator &=(interval &a,const interval &b)
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}
131inline interval & operator |=(interval &a,const real &b) noexcept
132{
133 a.inf=(a.inf<b)?a.inf:b,a.sup=(a.sup>b)?a.sup:b;
134 return a;
135}
136inline interval & operator &=(interval &a,const real &b)
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 ----
145inline bool operator ==(const interval &a,const interval &b) noexcept { return(a.inf==b.inf && a.sup==b.sup); }
146inline bool operator !=(const interval &a,const interval &b) noexcept { return(a.inf!=b.inf || a.sup!=b.sup); }
147inline bool operator ==(const real &r,const interval &a) noexcept { return(r==a.inf && r==a.sup); }
148inline bool operator !=(const real &r,const interval &a) noexcept { return(r!=a.inf || r!=a.sup); }
149inline bool operator ==(const interval &a,const real &r) noexcept { return(r==a.inf && r==a.sup); }
150inline bool operator !=(const interval &a,const real &r) noexcept { return(r!=a.inf || r!=a.sup); }
151
152inline bool operator ==(const int &r,const interval &a) noexcept { return(r==a.inf && r==a.sup); }
153inline bool operator !=(const int &r,const interval &a) noexcept { return(r!=a.inf || r!=a.sup); }
154inline bool operator ==(const interval &a,const int &r) noexcept { return(r==a.inf && r==a.sup); }
155inline bool operator !=(const interval &a,const int &r) noexcept { return(r!=a.inf || r!=a.sup); }
156
157inline bool operator ==(const long &r,const interval &a) noexcept { return(r==a.inf && r==a.sup); }
158inline bool operator !=(const long &r,const interval &a) noexcept { return(r!=a.inf || r!=a.sup); }
159inline bool operator ==(const interval &a,const long &r) noexcept { return(r==a.inf && r==a.sup); }
160inline bool operator !=(const interval &a,const long &r) noexcept { return(r!=a.inf || r!=a.sup); }
161
162inline bool operator ==(const double &r,const interval &a) noexcept { return(r==a.inf && r==a.sup); }
163inline bool operator !=(const double &r,const interval &a) noexcept { return(r!=a.inf || r!=a.sup); }
164inline bool operator ==(const interval &a,const double &r) noexcept { return(r==a.inf && r==a.sup); }
165inline bool operator !=(const interval &a,const double &r) noexcept { return(r!=a.inf || r!=a.sup); }
166
167// --- Mengenvergleiche ---
168// <,>,...
169inline bool operator <=(const interval &a,const interval &b) noexcept
170{
171 return(a.inf>=b.inf && a.sup<=b.sup);
172}
173inline bool operator >=(const interval &a,const interval &b) noexcept
174{
175 return(a.inf<=b.inf && a.sup>=b.sup);
176}
177inline bool operator <(const interval &a,const interval &b) noexcept
178{
179 return(a.inf>b.inf && a.sup<b.sup);
180}
181inline bool operator >(const interval &a,const interval &b) noexcept
182{
183 return(a.inf<b.inf && a.sup>b.sup);
184}
185
186inline bool operator <=(const real &a,const interval &b) noexcept
187{
188 return(a>=b.inf && a<=b.sup);
189}
190inline bool operator >=(const real &a,const interval &b) noexcept
191{
192 return(a<=b.inf && a>=b.sup);
193}
194inline bool operator <(const real &a,const interval &b) noexcept
195{
196 return(a>b.inf && a<b.sup);
197}
198
199inline bool operator <=(const interval &a,const real &b) noexcept
200{
201 return(a.inf>=b && a.sup<=b);
202}
203inline bool operator >=(const interval &a,const real &b) noexcept
204{
205 return(a.inf<=b && a.sup>=b);
206}
207inline bool operator >(const interval &a,const real &b) noexcept
208{
209 return(a.inf<b && a.sup>b);
210}
211
212inline bool operator !(const interval &a) noexcept { return (a.inf <= 0.0 && a.sup >= 0.0); }
213
214inline real & Inf (interval& a) noexcept { return a.inf; }
215inline const real & Inf (const interval &a) noexcept { return a.inf; }
216inline real & Sup (interval& a) noexcept { return a.sup; }
217inline const real & Sup (const interval &a) noexcept { return a.sup; }
218
219inline interval& SetInf (interval& a, const real& b) noexcept {a.inf=b; return a;}
220inline interval& SetSup (interval& a, const real& b) noexcept {a.sup=b; return a;}
221inline interval& UncheckedSetInf (interval& a, const real& b) noexcept { a.inf=b; return a;}
222inline interval& UncheckedSetSup (interval& a, const real& b) noexcept { a.sup=b; return a;}
223
224inline bool IsEmpty(const interval &a) noexcept { return (a.inf>a.sup); }
225
226inline interval abs(const interval &a) noexcept
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
240inline real diam(const interval & a) noexcept
241{
242 return subup(a.sup,a.inf);
243}
244
245inline real Mid(const interval & a) noexcept
246{
247 return addd( a.inf, subd(0.5*a.sup,0.5*a.inf) );
248}
249
250inline void times2pown(interval& x, const int& n) noexcept
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
283interval operator+ (const interval& a, const interval& b) noexcept
284 { return interval (adddown(a.inf,b.inf),addup(a.sup,b.sup)); }
285interval operator+ (const interval& a, const real& b) noexcept
286 { return interval (adddown(a.inf,b),addup(a.sup,b)); }
287interval operator+ (const real& a, const interval& b) noexcept
288 { return interval (adddown(a,b.inf),addup(a,b.sup)); }
289
290interval operator- (const interval& a, const interval& b) noexcept
291 { return interval ( subdown(a.inf,b.sup), subup(a.sup,b.inf)); }
292interval operator- (const interval& a, const real& b) noexcept
293 { return interval ( subdown(a.inf,b), subup(a.sup,b)); }
294interval operator- (const real& a, const interval& b) noexcept
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 //
311interval operator *(const interval& a, const interval& b) noexcept
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 //
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
432interval operator* (const real& a, const interval& b) noexcept
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
451interval operator* (const interval& a, const real& b) noexcept
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
470interval operator/ (const real& a, const interval& b) noexcept { return (interval(a) / b); }
471interval operator/ (const interval& a, const real& b) noexcept { return (a / interval(b)); }
472
473
474} // namespace cxsc
475
The Scalar Type interval.
Definition interval.hpp:55
interval & operator=(const real &a)
Implementation of standard assigning operator.
Definition interval.inl:52
interval()
Constructor of class interval.
Definition interval.hpp:64
The Scalar Type real.
Definition real.hpp:114
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc) noexcept
Implementation of standard algebraic addition and allocation operation.
Definition cdot.inl:251
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
Definition cimatrix.inl:730
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
Definition cimatrix.inl:738
const real minreal
Smallest positive denormalized representable floating-point number.
Definition real.cpp:63
interval _unchecked_interval(const real &a, const real &b)
Definition interval.inl:66
cimatrix & operator*=(cimatrix &m, const cinterval &c) noexcept
Implementation of multiplication and allocation operation.
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
cimatrix & operator/=(cimatrix &m, const cinterval &c) noexcept
Implementation of division and allocation operation.