- Generated by
1.7.5.1
|
C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: real.inl,v 1.33 2014/01/30 17:23:48 cxsc Exp $ */ 00025 00026 #ifdef _CXSC_REAL_HPP_INCLUDED 00027 00028 #include "fi_lib.hpp" 00029 #include <cmath> 00030 00031 namespace cxsc { 00032 extern "C" 00033 { 00034 a_real r_comp(a_real,a_intg); 00035 a_real r_mant(a_real); 00036 a_intg r_expo(a_real); 00037 } 00038 00039 inline double _double(const real &a) throw() 00040 { 00041 return a.w; 00042 } 00043 00049 inline real _real(const double &a) throw() 00050 { 00051 return real(a); 00052 } 00053 00054 inline real pred(const real& r) throw() 00055 { 00056 real ret; 00057 ret = fi_lib::q_pred(r.w); 00058 return ret; 00059 } 00060 00061 inline real succ(const real& r) throw() 00062 { 00063 real ret; 00064 ret = fi_lib::q_succ(r.w); 00065 return ret; 00066 } 00067 00068 inline a_intg expo(const real &r) throw() 00069 { 00070 return r_expo(*(a_real*)&r.w); 00071 } 00072 00073 inline real comp(const real &r, a_intg i) throw() 00074 { 00075 return real(r_comp(*(a_real*)&r.w,i)); 00076 } 00077 00078 inline real mant(const real &r) throw() 00079 { 00080 return real(r_mant(*(a_real*)&r.w)); 00081 } 00082 00083 inline real operator -(const real &a) throw() { return -a.w; } 00084 inline real operator +(const real &a) throw() { return a; } 00085 00086 inline real operator +(const real &a,const real &b) throw () { return a.w+b.w; } 00087 inline real operator -(const real &a,const real &b) throw () { return a.w-b.w; } 00088 inline real operator *(const real &a,const real &b) throw () { return a.w*b.w; } 00089 inline real operator /(const real &a,const real &b) throw () { return a.w/b.w; } 00090 00091 inline real& operator +=(real &a, const real &b) throw () { a.w+=b.w; return a; } 00092 inline real& operator -=(real &a, const real &b) throw () { a.w-=b.w; return a; } 00093 inline real& operator *=(real &a, const real &b) throw () { a.w*=b.w; return a; } 00094 inline real& operator /=(real &a, const real &b) throw () { a.w/=b.w; return a; } 00095 00096 inline bool operator! (const real& a) throw () { return !a.w; } 00097 inline bool operator== (const real& a, const real& b) throw () { return a.w==b.w; } 00098 inline bool operator!= (const real& a, const real& b) throw () { return a.w!=b.w; } 00099 inline bool operator< (const real& a, const real& b) throw () { return a.w<b.w; } 00100 inline bool operator<= (const real& a, const real& b) throw () { return a.w<=b.w; } 00101 inline bool operator>= (const real& a, const real& b) throw () { return a.w>=b.w; } 00102 inline bool operator> (const real& a, const real& b) throw () { return a.w>b.w; } 00103 00104 inline bool operator== (const real& a, const int & b) throw () { return a.w==b; } 00105 inline bool operator!= (const real& a, const int & b) throw () { return a.w!=b; } 00106 inline bool operator== (const int & a, const real& b) throw () { return a==b.w; } 00107 inline bool operator!= (const int & a, const real& b) throw () { return a!=b.w; } 00108 inline bool operator== (const real& a, const long & b) throw () { return a.w==b; } 00109 inline bool operator!= (const real& a, const long & b) throw () { return a.w!=b; } 00110 inline bool operator== (const long & a, const real& b) throw () { return a==b.w; } 00111 inline bool operator!= (const long & a, const real& b) throw () { return a!=b.w; } 00112 inline bool operator== (const real& a, const float & b) throw () { return a.w==b; } 00113 inline bool operator!= (const real& a, const float & b) throw () { return a.w!=b; } 00114 inline bool operator== (const float & a, const real& b) throw () { return a==b.w; } 00115 inline bool operator!= (const float & a, const real& b) throw () { return a!=b.w; } 00116 inline bool operator== (const real& a, const double & b) throw () { return a.w==b; } 00117 inline bool operator!= (const real& a, const double & b) throw () { return a.w!=b; } 00118 inline bool operator== (const double & a, const real& b) throw () { return a==b.w; } 00119 inline bool operator!= (const double & a, const real& b) throw () { return a!=b.w; } 00120 00121 00122 inline real abs(const real &a) throw() 00123 { 00124 /* if(a.w>=0) 00125 return a.w; 00126 return -a.w;*/ 00127 return fabs(a.w); 00128 } 00129 00130 inline int sign(const real &a) throw() 00131 { 00132 if(a.w>0) 00133 return 1; 00134 else 00135 if(a.w) 00136 return -1; 00137 else 00138 return 0; 00139 } 00140 00141 inline void times2pown(real& r, const int n) // Blomquist 01.10.02. 00142 { // times2pown(r,n): The reference parameter r gets the new value: r*2^n; 00143 // r*2^n in the normalized range --> r*2^n is the exact value. 00144 // r*2^n in the denormalized range --> r*2^n is in general not exact. 00145 int j = expo(r) + n; 00146 if (j >= -1021) { r = comp(mant(r),j); } // r is normalized up to now ! 00147 else 00148 { // result now in the denormalized range: 00149 j += 1021; 00150 r = comp(mant(r), -1021); // r is still normalized 00151 if (j<-53) r = 0; 00152 else r = r * comp(0.5,j+1); 00153 } 00154 } // end of times2pown(...) 00155 00156 inline real pow2n(const int n) throw() // Blomquist 03.10.02. 00157 { // Fast and exact calculation of 2^n; -1074 <= n <= +1023; 00158 return comp(0.5,n+1); 00159 } 00160 00161 inline real max(const real & a, const real & b) { return a>b?a:b; } 00162 inline real min(const real & a, const real & b) { return a<b?a:b; } 00163 inline real Max(const real & a, const real & b) { return a>b?a:b; } 00164 00165 00166 // ------- Verknuepfungen mit nach oben bzw. unten gerichteter Rundung ------ 00167 // ------- (Operators with rounding direction upwards or downwards) --------- 00168 00169 } // namespace cxsc 00170 00171 #include "rts_real.hpp" 00172 00173 #if ROUND_ASM 00174 #include <r_ari.h> 00175 #endif 00176 00177 #if ROUND_C99_SAVE+ROUND_C99_QUICK 00178 #include <fenv.h> 00179 #endif 00180 00181 #ifdef _MSC_VER 00182 #include <float.h> 00183 #pragma fenv_access (on) 00184 static inline void setround(int r) { 00185 unsigned int control_word; 00186 _controlfp_s(&control_word,0,0); 00187 switch(r) { 00188 case -1: 00189 _controlfp_s(&control_word,_RC_DOWN,_MCW_RC); 00190 break; 00191 case 0: 00192 _controlfp_s(&control_word,_RC_NEAR,_MCW_RC); 00193 break; 00194 case 1: 00195 _controlfp_s(&control_word,_RC_UP,_MCW_RC); 00196 break; 00197 case 2: 00198 _controlfp_s(&control_word,_RC_CHOP,_MCW_RC); 00199 break; 00200 default: 00201 _controlfp_s(&control_word,_RC_NEAR,_MCW_RC); 00202 } 00203 } 00204 00205 #endif 00206 00207 #if ROUND_C96_SAVE+ROUND_C96_QUICK 00208 #include <fenv96.h> 00209 #define fegetround fegetround96 00210 #define fesetround fesetround96 00211 #endif 00212 00213 00214 namespace cxsc { 00215 00216 inline real addup(const real &x, const real &y) 00217 { 00218 #if ROUND_ASM 00219 double ret; 00220 ret = ra_addu(x.w,y.w); 00221 return real(ret); 00222 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00223 int mode; 00224 mode=fegetround(); 00225 fesetround(FE_UPWARD); 00226 double ret; 00227 ret = x.w + y.w; 00228 fesetround(mode); 00229 return real(ret); 00230 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00231 fesetround(FE_UPWARD); 00232 return( x.w + y.w); 00233 #elif _MSC_VER 00234 unsigned int cw; 00235 unsigned int mask = 0xFFFFF; 00236 _controlfp_s(&cw,0,0); 00237 setround(1); 00238 double ret; 00239 ret = x.w + y.w; 00240 _controlfp_s(&cw,cw,mask); 00241 return real(ret); 00242 #else 00243 a_real ret; 00244 ret = r_addu(_a_real(x.w), _a_real(y.w)); 00245 return _real(ret); 00246 #endif 00247 } 00248 00249 inline real adddown(const real &x, const real &y) 00250 { 00251 #if ROUND_ASM 00252 double ret; 00253 ret = ra_addd(x.w,y.w); 00254 return real(ret); 00255 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00256 int mode; 00257 mode=fegetround(); 00258 fesetround(FE_DOWNWARD); 00259 double ret; 00260 ret = x.w + y.w; 00261 fesetround(mode); 00262 return real(ret); 00263 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00264 fesetround(FE_DOWNWARD); 00265 return( x.w + y.w); 00266 #elif _MSC_VER 00267 unsigned int cw; 00268 unsigned int mask = 0xFFFFF; 00269 _controlfp_s(&cw,0,0); 00270 setround(-1); 00271 double ret; 00272 ret = x.w + y.w; 00273 _controlfp_s(&cw,cw,mask); 00274 return real(ret); 00275 #else 00276 a_real ret; 00277 ret = r_addd(_a_real(x.w), _a_real(y.w)); 00278 return _real(ret); 00279 #endif 00280 } 00281 00282 inline real subup(const real &x, const real &y) 00283 { 00284 #if ROUND_ASM 00285 double ret; 00286 ret = ra_subu(x.w,y.w); 00287 return real(ret); 00288 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00289 int mode; 00290 mode=fegetround(); 00291 fesetround(FE_UPWARD); 00292 double ret; 00293 ret = x.w - y.w; 00294 fesetround(mode); 00295 return real(ret); 00296 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00297 fesetround(FE_UPWARD); 00298 return( x.w - y.w); 00299 #elif _MSC_VER 00300 unsigned int cw; 00301 unsigned int mask = 0xFFFFF; 00302 _controlfp_s(&cw,0,0); 00303 setround(1); 00304 double ret; 00305 ret = x.w - y.w; 00306 _controlfp_s(&cw,cw,mask); 00307 return real(ret); 00308 #else 00309 a_real ret; 00310 ret = r_subu(_a_real(x.w), _a_real(y.w)); 00311 return _real(ret); 00312 #endif 00313 } 00314 00315 inline real subdown(const real &x, const real &y) 00316 { 00317 #if ROUND_ASM 00318 double ret; 00319 ret = ra_subd(x.w,y.w); 00320 return real(ret); 00321 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00322 int mode; 00323 mode=fegetround(); 00324 fesetround(FE_DOWNWARD); 00325 double ret; 00326 ret = x.w - y.w; 00327 fesetround(mode); 00328 return real(ret); 00329 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00330 fesetround(FE_DOWNWARD); 00331 return( x.w - y.w); 00332 #elif _MSC_VER 00333 unsigned int cw; 00334 unsigned int mask = 0xFFFFF; 00335 _controlfp_s(&cw,0,0); 00336 setround(-1); 00337 double ret; 00338 ret = x.w - y.w; 00339 _controlfp_s(&cw,cw,mask); 00340 return real(ret); 00341 #else 00342 a_real ret; 00343 ret = r_subd(_a_real(x.w), _a_real(y.w)); 00344 return _real(ret); 00345 #endif 00346 } 00347 00348 inline real multup(const real &x, const real &y) 00349 { 00350 #if ROUND_ASM 00351 double ret; 00352 ret = ra_mulu(x.w,y.w); 00353 return real(ret); 00354 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00355 int mode; 00356 mode=fegetround(); 00357 fesetround(FE_UPWARD); 00358 double ret; 00359 ret = x.w * y.w; 00360 fesetround(mode); 00361 return real(ret); 00362 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00363 fesetround(FE_UPWARD); 00364 return( x.w * y.w); 00365 #elif _MSC_VER 00366 unsigned int cw; 00367 unsigned int mask = 0xFFFFF; 00368 _controlfp_s(&cw,0,0); 00369 setround(1); 00370 double ret; 00371 ret = x.w * y.w; 00372 _controlfp_s(&cw,cw,mask); 00373 return real(ret); 00374 #else 00375 a_real ret; 00376 ret = r_mulu(_a_real(x.w), _a_real(y.w)); 00377 return _real(ret); 00378 #endif 00379 } 00380 00381 inline real multdown(const real &x, const real &y) 00382 { 00383 #if ROUND_ASM 00384 double ret; 00385 ret = ra_muld(x.w,y.w); 00386 return real(ret); 00387 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00388 int mode; 00389 mode=fegetround(); 00390 fesetround(FE_DOWNWARD); 00391 double ret; 00392 ret = x.w * y.w; 00393 fesetround(mode); 00394 return real(ret); 00395 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00396 fesetround(FE_DOWNWARD); 00397 return( x.w * y.w); 00398 #elif _MSC_VER 00399 unsigned int cw; 00400 unsigned int mask = 0xFFFFF; 00401 _controlfp_s(&cw,0,0); 00402 setround(-1); 00403 double ret; 00404 ret = x.w * y.w; 00405 _controlfp_s(&cw,cw,mask); 00406 return real(ret); 00407 #else 00408 a_real ret; 00409 ret = r_muld(_a_real(x.w), _a_real(y.w)); 00410 return _real(ret); 00411 #endif 00412 } 00413 00414 inline real divup(const real &x, const real &y) 00415 { 00416 #if ROUND_ASM 00417 double ret; 00418 ret = ra_divu(x.w,y.w); 00419 return real(ret); 00420 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00421 int mode; 00422 mode=fegetround(); 00423 fesetround(FE_UPWARD); 00424 double ret; 00425 ret = x.w / y.w; 00426 fesetround(mode); 00427 return real(ret); 00428 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00429 fesetround(FE_UPWARD); 00430 return( x.w / y.w); 00431 #elif _MSC_VER 00432 unsigned int cw; 00433 unsigned int mask = 0xFFFFF; 00434 _controlfp_s(&cw,0,0); 00435 setround(1); 00436 double ret; 00437 ret = x.w / y.w; 00438 _controlfp_s(&cw,cw,mask); 00439 return real(ret); 00440 #else 00441 a_real ret; 00442 ret = r_divu(_a_real(x.w), _a_real(y.w)); 00443 return _real(ret); 00444 #endif 00445 } 00446 00447 inline real divdown(const real &x, const real &y) 00448 { 00449 #if ROUND_ASM 00450 double ret; 00451 ret = ra_divd(x.w,y.w); 00452 return real(ret); 00453 #elif ROUND_C99_SAVE+ROUND_C96_SAVE 00454 int mode; 00455 mode=fegetround(); 00456 fesetround(FE_DOWNWARD); 00457 double ret; 00458 ret = x.w / y.w; 00459 fesetround(mode); 00460 return real(ret); 00461 #elif ROUND_C99_QUICK+ROUND_C96_QUICK 00462 fesetround(FE_DOWNWARD); 00463 return( x.w / y.w); 00464 #elif _MSC_VER 00465 unsigned int cw; 00466 unsigned int mask = 0xFFFFF; 00467 _controlfp_s(&cw,0,0); 00468 setround(-1); 00469 double ret; 00470 ret = x.w / y.w; 00471 _controlfp_s(&cw,cw,mask); 00472 return real(ret); 00473 #else 00474 a_real ret; 00475 ret = r_divd(_a_real(x.w), _a_real(y.w)); 00476 return _real(ret); 00477 #endif 00478 } 00479 00480 //---------------------------------------------------------------------------- 00481 // IsInfinity - prueft ob die uebergebene IEEE 64-Bit Gleitkommazahl 00482 // 'Unendlich' darstellt, dies ist der Fall, falls: 00483 // alle Bits des Exponenten 1 sind 00484 // alle Bits der Mantisse 0 sind 00485 // 00486 inline bool IsInfinity (const real &a) 00487 { 00488 if ((((a_btyp*)&a)[HIGHREAL] & 0x7FF00000L) != 0x7FF00000L) 00489 return false; 00490 if ((((a_btyp*)&a)[HIGHREAL] & 0x000FFFFFL) != 0L) 00491 return false; 00492 if (((a_btyp*)&a)[LOWREAL] != 0L) 00493 return false; 00494 return true; // a ist 'unendlich' 00495 } 00496 00497 //---------------------------------------------------------------------------- 00498 // IsNan - prueft ob die uebergebene IEEE 64-Bit Gleitkommazahl 00499 // eine ungueltige Zahl (Not a number) darstellt, 00500 // dies ist der Fall, falls: 00501 // alle Bits des Exponenten 1 sind 00502 // und nicht alle Bits der Mantisse 0 sind 00503 // 00504 inline bool IsQuietNaN (const real &a) 00505 { 00506 if ((((a_btyp*)&a)[HIGHREAL] & 0x7FF00000L) != 0x7FF00000L) 00507 return false; 00508 00509 if ((((a_btyp*)&a)[HIGHREAL] & 0x000FFFFFL) == 0L) 00510 { 00511 if (((a_btyp*)&a)[LOWREAL] == 0L) 00512 return false; 00513 } 00514 return true; 00515 } 00516 00517 //---------------------------------------------------------------------------- 00518 // IsSignalingNaN - prueft ob die uebergebene IEEE 64-Bit Gleitkommazahl 00519 // undefiniert ist (z.B. Berechnung der Wurzel einer negativen 00520 // Zahl), dies ist der Fall, falls 00521 // das Vorzeichenbit 1 ist 00522 // alle Bits des Exponenten 1 sind 00523 // das hoechste Bit der Mantisse 1 ist 00524 // und alle anderen Bits der Mantisse 0 sind 00525 // 00526 inline bool IsSignalingNaN (const real &a) 00527 { 00528 if ((((a_btyp*)&a)[HIGHREAL] & 0xFFF00000L) != 0xFFF00000L) 00529 return false; 00530 if ((((a_btyp*)&a)[HIGHREAL] & 0x000FFFFFL) != 0x00080000L) 00531 return false; 00532 if (((a_btyp*)&a)[LOWREAL] != 0L) 00533 return false; 00534 return true; 00535 } 00536 00537 } // namespace cxsc 00538 00539 #endif // _CXSC_REAL_HPP_INCLUDED
1.7.5.1