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 
31 namespace cxsc {
32 extern "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 
39 inline double _double(const real &a) throw()
40 {
41  return a.w;
42 }
43 
49 inline real _real(const double &a) throw()
50 {
51  return real(a);
52 }
53 
54 inline real pred(const real& r) throw()
55 {
56  real ret;
57  ret = fi_lib::q_pred(r.w);
58  return ret;
59 }
60 
61 inline real succ(const real& r) throw()
62 {
63  real ret;
64  ret = fi_lib::q_succ(r.w);
65  return ret;
66 }
67 
68 inline a_intg expo(const real &r) throw()
69 {
70  return r_expo(*(a_real*)&r.w);
71 }
72 
73 inline real comp(const real &r, a_intg i) throw()
74 {
75  return real(r_comp(*(a_real*)&r.w,i));
76 }
77 
78 inline real mant(const real &r) throw()
79 {
80  return real(r_mant(*(a_real*)&r.w));
81 }
82 
83 inline real operator -(const real &a) throw() { return -a.w; }
84 inline real operator +(const real &a) throw() { return a; }
85 
86 inline real operator +(const real &a,const real &b) throw () { return a.w+b.w; }
87 inline real operator -(const real &a,const real &b) throw () { return a.w-b.w; }
88 inline real operator *(const real &a,const real &b) throw () { return a.w*b.w; }
89 inline real operator /(const real &a,const real &b) throw () { return a.w/b.w; }
90 
91 inline real& operator +=(real &a, const real &b) throw () { a.w+=b.w; return a; }
92 inline real& operator -=(real &a, const real &b) throw () { a.w-=b.w; return a; }
93 inline real& operator *=(real &a, const real &b) throw () { a.w*=b.w; return a; }
94 inline real& operator /=(real &a, const real &b) throw () { a.w/=b.w; return a; }
95 
96 inline bool operator! (const real& a) throw () { return !a.w; }
97 inline bool operator== (const real& a, const real& b) throw () { return a.w==b.w; }
98 inline bool operator!= (const real& a, const real& b) throw () { return a.w!=b.w; }
99 inline bool operator< (const real& a, const real& b) throw () { return a.w<b.w; }
100 inline bool operator<= (const real& a, const real& b) throw () { return a.w<=b.w; }
101 inline bool operator>= (const real& a, const real& b) throw () { return a.w>=b.w; }
102 inline bool operator> (const real& a, const real& b) throw () { return a.w>b.w; }
103 
104 inline bool operator== (const real& a, const int & b) throw () { return a.w==b; }
105 inline bool operator!= (const real& a, const int & b) throw () { return a.w!=b; }
106 inline bool operator== (const int & a, const real& b) throw () { return a==b.w; }
107 inline bool operator!= (const int & a, const real& b) throw () { return a!=b.w; }
108 inline bool operator== (const real& a, const long & b) throw () { return a.w==b; }
109 inline bool operator!= (const real& a, const long & b) throw () { return a.w!=b; }
110 inline bool operator== (const long & a, const real& b) throw () { return a==b.w; }
111 inline bool operator!= (const long & a, const real& b) throw () { return a!=b.w; }
112 inline bool operator== (const real& a, const float & b) throw () { return a.w==b; }
113 inline bool operator!= (const real& a, const float & b) throw () { return a.w!=b; }
114 inline bool operator== (const float & a, const real& b) throw () { return a==b.w; }
115 inline bool operator!= (const float & a, const real& b) throw () { return a!=b.w; }
116 inline bool operator== (const real& a, const double & b) throw () { return a.w==b; }
117 inline bool operator!= (const real& a, const double & b) throw () { return a.w!=b; }
118 inline bool operator== (const double & a, const real& b) throw () { return a==b.w; }
119 inline bool operator!= (const double & a, const real& b) throw () { return a!=b.w; }
120 
121 
122 inline real abs(const real &a) throw()
123 {
124  /* if(a.w>=0)
125  return a.w;
126  return -a.w;*/
127  return fabs(a.w);
128 }
129 
130 inline int sign(const real &a) throw()
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 
141 inline 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 
156 inline real pow2n(const int n) throw() // Blomquist 03.10.02.
157 { // Fast and exact calculation of 2^n; -1074 <= n <= +1023;
158  return comp(0.5,n+1);
159 }
160 
161 inline real max(const real & a, const real & b) { return a>b?a:b; }
162 inline real min(const real & a, const real & b) { return a<b?a:b; }
163 inline 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 
214 namespace cxsc {
215 
216 inline 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 
249 inline 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 
282 inline 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 
315 inline 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 
348 inline 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 
381 inline 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 
414 inline 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 
447 inline 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 //
486 inline 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 //
504 inline 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 //
526 inline 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
cimatrix & operator/=(cimatrix &m, const cinterval &c)
Implementation of division and allocation operation.
Definition: cimatrix.inl:1623
bool IsInfinity(const real &a)
Returns if the given real value represents the value infinity.
bool IsQuietNaN(const real &a)
Returns if the given real value represents the value of a quiet NaN.
bool IsSignalingNaN(const real &a)
Returns if the given real value represents the value of a signaling NaN.
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
real Max(const real &a, const real &b)
Returns the greater value of two real values (for Compatibility with former r_util....
real pow2n(const int n)
Returns the value of .
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
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
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737