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
39namespace cxsc {
40
42{
43 *this=rnd(a);
44}
45
47{
48 return *this=rnd(a);
49}
50
52{
53 rnd(a,inf,sup);
54}
55
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
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/*
84interval operator+(const interval &a, const interval &b) noexcept { return interval(adddown(a.inf,b.inf),addup(a.sup,b.sup)); }
85interval operator+(const interval &a, const real &b) noexcept { return interval(adddown(a.inf,b),addup(a.sup,b)); }
86interval operator+(const real &a, const interval &b) noexcept { return interval(adddown(a,b.inf),addup(a,b.sup)); }
87
88interval operator-(const interval &a, const interval &b) noexcept { return interval(subdown(a.inf,b.inf),subup(a.sup,b.sup)); }
89interval operator-(const interval &a, const real &b) noexcept { return interval(subdown(a.inf,b),subup(a.sup,b)); }
90interval operator-(const real &a, const interval &b) noexcept { return interval(subdown(a,b.sup),subup(a,b.inf)); }
91
92interval operator*(const interval &a, const interval &b) noexcept { return mul_ii (a, b); }
93interval operator*(const interval &a, const real &b) noexcept { return mul_id (a, *(double*)&b); }
94interval operator*(const real &a, const interval &b) noexcept { return mul_di (*(double*)&a, b); }
95
96interval operator/(const interval &a, const interval &b) noexcept { return div_ii (a, b); }
97interval operator/(const interval &a, const real &b) noexcept { return div_id (a, *(double*)&b); }
98interval operator/(const real &a, const interval &b) noexcept { return div_di (*(double*)&a, b); }
99*/
100
101
102bool operator ==(const interval &a,const dotprecision &r) noexcept { return(r==a.inf && r==a.sup); }
103bool operator !=(const interval &a,const dotprecision &r) noexcept { return(r!=a.inf || r!=a.sup); }
104bool operator ==(const dotprecision &r,const interval &a) noexcept { return(r==a.inf && r==a.sup); }
105bool operator !=(const dotprecision &r,const interval &a) noexcept { return(r!=a.inf || r!=a.sup); }
106
107bool operator <=(const dotprecision &a,const interval &b) noexcept
108{
109 return(a>=b.inf && a<=b.sup);
110}
111bool operator >=(const dotprecision &a,const interval &b) noexcept
112{
113 return(a<=b.inf && a>=b.sup);
114}
115bool operator <(const dotprecision &a,const interval &b) noexcept
116{
117 return(a>b.inf && a<b.sup);
118}
119
120bool operator <=(const interval &a,const dotprecision &b) noexcept
121{
122 return(a.inf>=b && a.sup<=b);
123}
124bool operator >=(const interval &a,const dotprecision &b) noexcept
125{
126 return(a.inf<=b && a.sup>=b);
127}
128bool operator >(const interval &a,const dotprecision &b) noexcept
129{
130 return(a.inf<b && a.sup>b);
131}
132
133
134
135// ---- Ausgabefunkt. ---------------------------------------
136
137std::ostream & operator << (std::ostream &s, const interval& a) noexcept
138{
139 s << '[' << SaveOpt << RndDown
140 << a.inf << ',' << RndUp
141 << a.sup << RestoreOpt
142 << ']';
143 return s;
144}
145std::string & operator << (std::string &s, const interval &a) noexcept
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
157std::istream & operator >> (std::istream &s, interval &a) noexcept
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
188std::string & operator >> (std::string &s, interval &a) noexcept
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
205void operator >>(const std::string &s,interval &a) noexcept
206{
207 std::string r(s);
208 r>>a;
209}
210void operator >>(const char *s,interval &a) noexcept
211{
212 std::string r(s);
213 r>>a;
214}
215
216real mid(const interval & a) noexcept
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
246int in (const real& x, const interval& y ) // Contained-in relation
247 { return ( (Inf(y) <= x) && (x <= Sup(y)) ); } //----------------------
248
249int 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
281interval 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
288int Disjoint ( const interval& a, const interval& b ) // Test for disjointedness
289{ //------------------------
290 return (Inf(a) > Sup(b)) || (Inf(b) > Sup(a));
291}
292
293real 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
303real 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
316real 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//----------------------------------------------------------------------------
335int 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
The Data Type dotprecision.
Definition dot.hpp:112
The Data Type idotprecision.
Definition idot.hpp:48
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 Multiple-Precision Data Type l_real.
Definition l_real.hpp:78
The Scalar Type real.
Definition real.hpp:114
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
const interval Sqrt2Pi_interval
Enclosure-Interval for .
Definition interval.cpp:380
const interval EpPi_interval
Enclosure-Interval for .
Definition interval.cpp:452
const interval Sqrt5_interval
Enclosure-Interval for .
Definition interval.cpp:396
const interval Er_interval
Enclosure-Interval for .
Definition interval.cpp:443
const interval Sqrt3d2_interval
Enclosure-Interval for .
Definition interval.cpp:412
const interval Pid4_interval
Enclosure-Interval for .
Definition interval.cpp:363
const interval Pi2r_interval
Enclosure-Interval for .
Definition interval.cpp:369
const interval Sqrt2_interval
Enclosure-Interval for .
Definition interval.cpp:392
const interval Ln2Pi_interval
Enclosure-Interval for .
Definition interval.cpp:436
const interval Sqrt3r_interval
Enclosure-Interval for .
Definition interval.cpp:416
const interval Ln2r_interval
Enclosure-Interval for .
Definition interval.cpp:423
const interval EpPid4_interval
Enclosure-Interval for .
Definition interval.cpp:463
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
Definition interval.cpp:288
const interval E_interval
Enclosure-Interval for .
Definition interval.cpp:440
const interval LnPi_interval
Enclosure-Interval for .
Definition interval.cpp:433
const interval Pid2_interval
Enclosure-Interval for .
Definition interval.cpp:357
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
Definition cimatrix.inl:738
int UlpAcc(const interval &x, int n)
Checks if the diameter of the interval is ulps.
Definition interval.cpp:335
const interval Ep2_interval
Enclosure-Interval for .
Definition interval.cpp:446
real RelDiam(const interval &x)
Computes the relative diameter .
Definition interval.cpp:316
const interval Sqrt7_interval
Enclosure-Interval for .
Definition interval.cpp:400
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
const interval Pir_interval
Enclosure-Interval for .
Definition interval.cpp:366
const interval Pip2_interval
Enclosure-Interval for .
Definition interval.cpp:373
const interval Ln2_interval
Enclosure-Interval for .
Definition interval.cpp:420
const interval Pid3_interval
Enclosure-Interval for .
Definition interval.cpp:360
const interval SqrtPir_interval
Enclosure-Interval for .
Definition interval.cpp:384
const interval Sqrt2r_interval
Enclosure-Interval for .
Definition interval.cpp:404
const interval Pi2_interval
Enclosure-Interval for .
Definition interval.cpp:351
const interval Ep2Pi_interval
Enclosure-Interval for .
Definition interval.cpp:455
const interval Ln10_interval
Enclosure-Interval for .
Definition interval.cpp:426
const interval Ep2r_interval
Enclosure-Interval for .
Definition interval.cpp:449
real AbsMax(const interval &x)
Computes the greatest absolute value .
Definition interval.cpp:303
const interval Ln10r_interval
Enclosure-Interval for .
Definition interval.cpp:429
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739
const interval EpPid2_interval
Enclosure-Interval for .
Definition interval.cpp:459
const interval Sqrt3_interval
Enclosure-Interval for .
Definition interval.cpp:408
const interval Sqrt2Pir_interval
Enclosure-Interval for .
Definition interval.cpp:388
const interval Pi3_interval
Enclosure-Interval for .
Definition interval.cpp:354
const interval SqrtPi_interval
Enclosure-Interval for .
Definition interval.cpp:376
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.
const interval Pi_interval
Enclosure-Interval for .
Definition interval.cpp:348
real AbsMin(const interval &x)
Computes the smallest absolute value .
Definition interval.cpp:293