C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
real.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: real.inl,v 1.33 2014/01/30 17:23:48 cxsc Exp $ */
25
26#ifdef _CXSC_REAL_HPP_INCLUDED
27
28#include "fi_lib.hpp"
29#include <cmath>
30
31namespace cxsc {
32extern "C"
33{
34 a_real r_comp(a_real,a_intg);
35 a_real r_mant(a_real);
36 a_intg r_expo(a_real);
37}
38
39inline double _double(const real &a) noexcept
40{
41 return a.w;
42}
43
49inline real _real(const double &a) noexcept
50{
51 return real(a);
52}
53
54inline real pred(const real& r) noexcept
55{
56 real ret;
57 ret = fi_lib::q_pred(r.w);
58 return ret;
59}
60
61inline real succ(const real& r) noexcept
62{
63 real ret;
64 ret = fi_lib::q_succ(r.w);
65 return ret;
66}
67
68inline a_intg expo(const real &r) noexcept
69{
70 return r_expo(*(a_real*)&r.w);
71}
72
73inline real comp(const real &r, a_intg i) noexcept
74{
75 return real(r_comp(*(a_real*)&r.w,i));
76}
77
78inline real mant(const real &r) noexcept
79{
80 return real(r_mant(*(a_real*)&r.w));
81}
82
83inline real operator -(const real &a) noexcept { return -a.w; }
84inline real operator +(const real &a) noexcept { return a; }
85
86inline real operator +(const real &a,const real &b) noexcept { return a.w+b.w; }
87inline real operator -(const real &a,const real &b) noexcept { return a.w-b.w; }
88inline real operator *(const real &a,const real &b) noexcept { return a.w*b.w; }
89inline real operator /(const real &a,const real &b) noexcept { return a.w/b.w; }
90
91inline real& operator +=(real &a, const real &b) noexcept { a.w+=b.w; return a; }
92inline real& operator -=(real &a, const real &b) noexcept { a.w-=b.w; return a; }
93inline real& operator *=(real &a, const real &b) noexcept { a.w*=b.w; return a; }
94inline real& operator /=(real &a, const real &b) noexcept { a.w/=b.w; return a; }
95
96inline bool operator! (const real& a) noexcept { return !a.w; }
97inline bool operator== (const real& a, const real& b) noexcept { return a.w==b.w; }
98inline bool operator!= (const real& a, const real& b) noexcept { return a.w!=b.w; }
99inline bool operator< (const real& a, const real& b) noexcept { return a.w<b.w; }
100inline bool operator<= (const real& a, const real& b) noexcept { return a.w<=b.w; }
101inline bool operator>= (const real& a, const real& b) noexcept { return a.w>=b.w; }
102inline bool operator> (const real& a, const real& b) noexcept { return a.w>b.w; }
103
104inline bool operator== (const real& a, const int & b) noexcept { return a.w==b; }
105inline bool operator!= (const real& a, const int & b) noexcept { return a.w!=b; }
106inline bool operator== (const int & a, const real& b) noexcept { return a==b.w; }
107inline bool operator!= (const int & a, const real& b) noexcept { return a!=b.w; }
108inline bool operator== (const real& a, const long & b) noexcept { return a.w==b; }
109inline bool operator!= (const real& a, const long & b) noexcept { return a.w!=b; }
110inline bool operator== (const long & a, const real& b) noexcept { return a==b.w; }
111inline bool operator!= (const long & a, const real& b) noexcept { return a!=b.w; }
112inline bool operator== (const real& a, const float & b) noexcept { return a.w==b; }
113inline bool operator!= (const real& a, const float & b) noexcept { return a.w!=b; }
114inline bool operator== (const float & a, const real& b) noexcept { return a==b.w; }
115inline bool operator!= (const float & a, const real& b) noexcept { return a!=b.w; }
116inline bool operator== (const real& a, const double & b) noexcept { return a.w==b; }
117inline bool operator!= (const real& a, const double & b) noexcept { return a.w!=b; }
118inline bool operator== (const double & a, const real& b) noexcept { return a==b.w; }
119inline bool operator!= (const double & a, const real& b) noexcept { return a!=b.w; }
120
121
122inline real abs(const real &a) noexcept
123{
124 /* if(a.w>=0)
125 return a.w;
126 return -a.w;*/
127 return fabs(a.w);
128}
129
130inline int sign(const real &a) noexcept
131{
132 if(a.w>0)
133 return 1;
134 else
135 if(a.w)
136 return -1;
137 else
138 return 0;
139}
140
141inline void times2pown(real& r, const int n) // Blomquist 01.10.02.
142{ // times2pown(r,n): The reference parameter r gets the new value: r*2^n;
143 // r*2^n in the normalized range --> r*2^n is the exact value.
144 // r*2^n in the denormalized range --> r*2^n is in general not exact.
145 int j = expo(r) + n;
146 if (j >= -1021) { r = comp(mant(r),j); } // r is normalized up to now !
147 else
148 { // result now in the denormalized range:
149 j += 1021;
150 r = comp(mant(r), -1021); // r is still normalized
151 if (j<-53) r = 0;
152 else r = r * comp(0.5,j+1);
153 }
154} // end of times2pown(...)
155
156inline real pow2n(const int n) noexcept // Blomquist 03.10.02.
157{ // Fast and exact calculation of 2^n; -1074 <= n <= +1023;
158 return comp(0.5,n+1);
159}
160
161inline real max(const real & a, const real & b) { return a>b?a:b; }
162inline real min(const real & a, const real & b) { return a<b?a:b; }
163inline real Max(const real & a, const real & b) { return a>b?a:b; }
164
165
166// ------- Verknuepfungen mit nach oben bzw. unten gerichteter Rundung ------
167// ------- (Operators with rounding direction upwards or downwards) ---------
168
169} // namespace cxsc
170
171#include "rts_real.hpp"
172
173#if ROUND_ASM
174#include <r_ari.h>
175#endif
176
177#if ROUND_C99_SAVE+ROUND_C99_QUICK
178#include <fenv.h>
179#endif
180
181#ifdef _MSC_VER
182#include <float.h>
183#pragma fenv_access (on)
184 static inline void setround(int r) {
185 unsigned int control_word;
186 _controlfp_s(&control_word,0,0);
187 switch(r) {
188 case -1:
189 _controlfp_s(&control_word,_RC_DOWN,_MCW_RC);
190 break;
191 case 0:
192 _controlfp_s(&control_word,_RC_NEAR,_MCW_RC);
193 break;
194 case 1:
195 _controlfp_s(&control_word,_RC_UP,_MCW_RC);
196 break;
197 case 2:
198 _controlfp_s(&control_word,_RC_CHOP,_MCW_RC);
199 break;
200 default:
201 _controlfp_s(&control_word,_RC_NEAR,_MCW_RC);
202 }
203 }
204
205#endif
206
207#if ROUND_C96_SAVE+ROUND_C96_QUICK
208#include <fenv96.h>
209#define fegetround fegetround96
210#define fesetround fesetround96
211#endif
212
213
214namespace cxsc {
215
216inline real addup(const real &x, const real &y)
217{
218 #if ROUND_ASM
219 double ret;
220 ret = ra_addu(x.w,y.w);
221 return real(ret);
222 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
223 int mode;
224 mode=fegetround();
225 fesetround(FE_UPWARD);
226 double ret;
227 ret = x.w + y.w;
228 fesetround(mode);
229 return real(ret);
230 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
231 fesetround(FE_UPWARD);
232 return( x.w + y.w);
233 #elif _MSC_VER
234 unsigned int cw;
235 unsigned int mask = 0xFFFFF;
236 _controlfp_s(&cw,0,0);
237 setround(1);
238 double ret;
239 ret = x.w + y.w;
240 _controlfp_s(&cw,cw,mask);
241 return real(ret);
242 #else
243 a_real ret;
244 ret = r_addu(_a_real(x.w), _a_real(y.w));
245 return _real(ret);
246 #endif
247}
248
249inline real adddown(const real &x, const real &y)
250{
251 #if ROUND_ASM
252 double ret;
253 ret = ra_addd(x.w,y.w);
254 return real(ret);
255 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
256 int mode;
257 mode=fegetround();
258 fesetround(FE_DOWNWARD);
259 double ret;
260 ret = x.w + y.w;
261 fesetround(mode);
262 return real(ret);
263 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
264 fesetround(FE_DOWNWARD);
265 return( x.w + y.w);
266 #elif _MSC_VER
267 unsigned int cw;
268 unsigned int mask = 0xFFFFF;
269 _controlfp_s(&cw,0,0);
270 setround(-1);
271 double ret;
272 ret = x.w + y.w;
273 _controlfp_s(&cw,cw,mask);
274 return real(ret);
275 #else
276 a_real ret;
277 ret = r_addd(_a_real(x.w), _a_real(y.w));
278 return _real(ret);
279 #endif
280}
281
282inline real subup(const real &x, const real &y)
283{
284 #if ROUND_ASM
285 double ret;
286 ret = ra_subu(x.w,y.w);
287 return real(ret);
288 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
289 int mode;
290 mode=fegetround();
291 fesetround(FE_UPWARD);
292 double ret;
293 ret = x.w - y.w;
294 fesetround(mode);
295 return real(ret);
296 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
297 fesetround(FE_UPWARD);
298 return( x.w - y.w);
299 #elif _MSC_VER
300 unsigned int cw;
301 unsigned int mask = 0xFFFFF;
302 _controlfp_s(&cw,0,0);
303 setround(1);
304 double ret;
305 ret = x.w - y.w;
306 _controlfp_s(&cw,cw,mask);
307 return real(ret);
308 #else
309 a_real ret;
310 ret = r_subu(_a_real(x.w), _a_real(y.w));
311 return _real(ret);
312 #endif
313}
314
315inline real subdown(const real &x, const real &y)
316{
317 #if ROUND_ASM
318 double ret;
319 ret = ra_subd(x.w,y.w);
320 return real(ret);
321 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
322 int mode;
323 mode=fegetround();
324 fesetround(FE_DOWNWARD);
325 double ret;
326 ret = x.w - y.w;
327 fesetround(mode);
328 return real(ret);
329 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
330 fesetround(FE_DOWNWARD);
331 return( x.w - y.w);
332 #elif _MSC_VER
333 unsigned int cw;
334 unsigned int mask = 0xFFFFF;
335 _controlfp_s(&cw,0,0);
336 setround(-1);
337 double ret;
338 ret = x.w - y.w;
339 _controlfp_s(&cw,cw,mask);
340 return real(ret);
341 #else
342 a_real ret;
343 ret = r_subd(_a_real(x.w), _a_real(y.w));
344 return _real(ret);
345 #endif
346}
347
348inline real multup(const real &x, const real &y)
349{
350 #if ROUND_ASM
351 double ret;
352 ret = ra_mulu(x.w,y.w);
353 return real(ret);
354 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
355 int mode;
356 mode=fegetround();
357 fesetround(FE_UPWARD);
358 double ret;
359 ret = x.w * y.w;
360 fesetround(mode);
361 return real(ret);
362 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
363 fesetround(FE_UPWARD);
364 return( x.w * y.w);
365 #elif _MSC_VER
366 unsigned int cw;
367 unsigned int mask = 0xFFFFF;
368 _controlfp_s(&cw,0,0);
369 setround(1);
370 double ret;
371 ret = x.w * y.w;
372 _controlfp_s(&cw,cw,mask);
373 return real(ret);
374 #else
375 a_real ret;
376 ret = r_mulu(_a_real(x.w), _a_real(y.w));
377 return _real(ret);
378 #endif
379}
380
381inline real multdown(const real &x, const real &y)
382{
383 #if ROUND_ASM
384 double ret;
385 ret = ra_muld(x.w,y.w);
386 return real(ret);
387 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
388 int mode;
389 mode=fegetround();
390 fesetround(FE_DOWNWARD);
391 double ret;
392 ret = x.w * y.w;
393 fesetround(mode);
394 return real(ret);
395 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
396 fesetround(FE_DOWNWARD);
397 return( x.w * y.w);
398 #elif _MSC_VER
399 unsigned int cw;
400 unsigned int mask = 0xFFFFF;
401 _controlfp_s(&cw,0,0);
402 setround(-1);
403 double ret;
404 ret = x.w * y.w;
405 _controlfp_s(&cw,cw,mask);
406 return real(ret);
407 #else
408 a_real ret;
409 ret = r_muld(_a_real(x.w), _a_real(y.w));
410 return _real(ret);
411 #endif
412}
413
414inline real divup(const real &x, const real &y)
415{
416 #if ROUND_ASM
417 double ret;
418 ret = ra_divu(x.w,y.w);
419 return real(ret);
420 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
421 int mode;
422 mode=fegetround();
423 fesetround(FE_UPWARD);
424 double ret;
425 ret = x.w / y.w;
426 fesetround(mode);
427 return real(ret);
428 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
429 fesetround(FE_UPWARD);
430 return( x.w / y.w);
431 #elif _MSC_VER
432 unsigned int cw;
433 unsigned int mask = 0xFFFFF;
434 _controlfp_s(&cw,0,0);
435 setround(1);
436 double ret;
437 ret = x.w / y.w;
438 _controlfp_s(&cw,cw,mask);
439 return real(ret);
440 #else
441 a_real ret;
442 ret = r_divu(_a_real(x.w), _a_real(y.w));
443 return _real(ret);
444 #endif
445}
446
447inline real divdown(const real &x, const real &y)
448{
449 #if ROUND_ASM
450 double ret;
451 ret = ra_divd(x.w,y.w);
452 return real(ret);
453 #elif ROUND_C99_SAVE+ROUND_C96_SAVE
454 int mode;
455 mode=fegetround();
456 fesetround(FE_DOWNWARD);
457 double ret;
458 ret = x.w / y.w;
459 fesetround(mode);
460 return real(ret);
461 #elif ROUND_C99_QUICK+ROUND_C96_QUICK
462 fesetround(FE_DOWNWARD);
463 return( x.w / y.w);
464 #elif _MSC_VER
465 unsigned int cw;
466 unsigned int mask = 0xFFFFF;
467 _controlfp_s(&cw,0,0);
468 setround(-1);
469 double ret;
470 ret = x.w / y.w;
471 _controlfp_s(&cw,cw,mask);
472 return real(ret);
473 #else
474 a_real ret;
475 ret = r_divd(_a_real(x.w), _a_real(y.w));
476 return _real(ret);
477 #endif
478}
479
480//----------------------------------------------------------------------------
481// IsInfinity - prueft ob die uebergebene IEEE 64-Bit Gleitkommazahl
482// 'Unendlich' darstellt, dies ist der Fall, falls:
483// alle Bits des Exponenten 1 sind
484// alle Bits der Mantisse 0 sind
485//
486inline bool IsInfinity (const real &a)
487{
488 if ((((a_btyp*)&a)[HIGHREAL] & 0x7FF00000L) != 0x7FF00000L)
489 return false;
490 if ((((a_btyp*)&a)[HIGHREAL] & 0x000FFFFFL) != 0L)
491 return false;
492 if (((a_btyp*)&a)[LOWREAL] != 0L)
493 return false;
494 return true; // a ist 'unendlich'
495}
496
497//----------------------------------------------------------------------------
498// IsNan - prueft ob die uebergebene IEEE 64-Bit Gleitkommazahl
499// eine ungueltige Zahl (Not a number) darstellt,
500// dies ist der Fall, falls:
501// alle Bits des Exponenten 1 sind
502// und nicht alle Bits der Mantisse 0 sind
503//
504inline bool IsQuietNaN (const real &a)
505{
506 if ((((a_btyp*)&a)[HIGHREAL] & 0x7FF00000L) != 0x7FF00000L)
507 return false;
508
509 if ((((a_btyp*)&a)[HIGHREAL] & 0x000FFFFFL) == 0L)
510 {
511 if (((a_btyp*)&a)[LOWREAL] == 0L)
512 return false;
513 }
514 return true;
515}
516
517//----------------------------------------------------------------------------
518// IsSignalingNaN - prueft ob die uebergebene IEEE 64-Bit Gleitkommazahl
519// undefiniert ist (z.B. Berechnung der Wurzel einer negativen
520// Zahl), dies ist der Fall, falls
521// das Vorzeichenbit 1 ist
522// alle Bits des Exponenten 1 sind
523// das hoechste Bit der Mantisse 1 ist
524// und alle anderen Bits der Mantisse 0 sind
525//
526inline bool IsSignalingNaN (const real &a)
527{
528 if ((((a_btyp*)&a)[HIGHREAL] & 0xFFF00000L) != 0xFFF00000L)
529 return false;
530 if ((((a_btyp*)&a)[HIGHREAL] & 0x000FFFFFL) != 0x00080000L)
531 return false;
532 if (((a_btyp*)&a)[LOWREAL] != 0L)
533 return false;
534 return true;
535}
536
537} // namespace cxsc
538
539#endif // _CXSC_REAL_HPP_INCLUDED
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
bool IsQuietNaN(const real &a)
Returns if the given real value represents the value of a quiet NaN.
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
Definition cimatrix.inl:730
bool IsSignalingNaN(const real &a)
Returns if the given real value represents the value of a signaling NaN.
real pow2n(const int n) noexcept
Returns the value of .
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
real Max(const real &a, const real &b)
Returns the greater value of two real values (for Compatibility with former r_util....
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
Definition cimatrix.inl:731
bool IsInfinity(const real &a)
Returns if the given real value represents the value infinity.
cimatrix & operator/=(cimatrix &m, const cinterval &c) noexcept
Implementation of division and allocation operation.