ppALIGN API documentation
src/xdouble.hpp
00001 /****************************************************************************** 00002 Copyright 2009 Stefan Wolfsheimer & Gregory Nuel. 00003 00004 This file is part of ppALIGN 00005 00006 ppALIGN is free software; you can redistribute it and/or modify 00007 it under the terms of the GNU General Public License as published by 00008 the Free Software Foundation; either version 2 of the License, or 00009 (at your option) any later version. 00010 00011 ppALIGN is distributed in the hope that it will be useful, 00012 but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 GNU General Public License for more details. 00015 00016 You should have received a copy of the GNU General Public License 00017 along with ppALIGN; if not, write to the Free Software 00018 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 00019 *******************************************************************************/ 00020 00021 00022 #ifndef _XDOUBLE_HPP_ 00023 #define _XDOUBLE_HPP_ 00024 #include <iostream> 00025 #include <math.h> 00026 #include <stdlib.h> 00027 00028 00031 #define NTL_SNS 00032 #define NTL_BITS_PER_LONG (32) 00033 #define NTL_BITS_PER_INT (32) 00034 #define NTL_BITS_PER_SIZE_T (32) 00035 #define NTL_MAX_INT (2147483647) 00036 #define NTL_MIN_INT (-NTL_MAX_INT - 1) 00037 00038 00039 00040 00046 class xdouble 00047 { 00048 private: 00049 double x; 00050 long e; 00051 00052 xdouble(double xx, long ee) : x(xx), e(ee) { } // internal use only 00053 inline void from_double(double x); 00054 //double to_double(const double & x); 00055 static inline void Error(const char * str); 00056 public: 00057 00058 00061 xdouble() : x(0), e(0) {} 00062 00066 xdouble(const xdouble & a) : x(a.x),e(a.e) { } 00067 00071 xdouble(const double & x) { from_double(x); } 00072 00076 xdouble(const int & x) { from_double(x); } 00077 00081 xdouble(const long & x) { from_double(x); } 00082 00086 xdouble(const unsigned int & x) { from_double(x); } 00087 00091 xdouble(const unsigned long & x) { from_double(x); } 00092 00093 00096 // @{ 00097 00101 inline xdouble& operator=(double a); 00102 00106 inline xdouble& operator=(int a); 00107 00111 inline xdouble& operator=(long a); 00112 00116 inline xdouble& operator=(unsigned int & a); 00117 00121 inline xdouble& operator=(unsigned long & a); 00122 00123 00127 inline xdouble& operator+=(const xdouble& b){ *this = *this + b; return *this; } 00128 00132 inline xdouble& operator-=(const xdouble& b){ *this = *this - b; return *this; } 00133 00137 inline xdouble& operator*=(const xdouble& b){ *this = *this * b; return *this; } 00138 00142 inline xdouble& operator/=(const xdouble& b){ *this = *this / b; return *this; } 00143 // @} 00144 00145 00149 inline double to_double() const; 00150 00151 00152 00154 00156 // @{ 00159 friend inline std::ostream & operator<<(std::ostream & ost, const xdouble & v) { ost << v.to_double(); return ost; } 00160 00163 friend inline std::istream & operator>>(std::istream & ist, xdouble & v) { double temp; ist >> temp; v=temp; return ist; } 00164 // @} 00165 00167 00169 // @{ 00170 00174 friend inline xdouble operator-(const xdouble & a); 00175 // @} 00176 00178 00180 // @{ 00181 00184 // @{ 00188 friend inline xdouble operator+(const xdouble & a,const xdouble & b); 00189 00191 friend inline xdouble operator+(const double & a,const xdouble & b) { return xdouble(a) + b; } 00193 friend inline xdouble operator+(const int & a,const xdouble & b) { return xdouble(a) + b; } 00195 friend inline xdouble operator+(const unsigned int & a,const xdouble & b) { return xdouble(a) + b; } 00197 friend inline xdouble operator+(const long & a,const xdouble & b) { return xdouble(a) + b; } 00199 friend inline xdouble operator+(const unsigned long & a,const xdouble & b) { return xdouble(a) + b; } 00201 friend inline xdouble operator+(const xdouble & a,const double & b) { return a + xdouble(b); } 00203 friend inline xdouble operator+(const xdouble & a,const int & b) { return a + xdouble(b); } 00205 friend inline xdouble operator+(const xdouble & a,const unsigned int & b) { return a + xdouble(b); } 00207 friend inline xdouble operator+(const xdouble & a,const long & b) { return a + xdouble(b); } 00209 friend inline xdouble operator+(const xdouble & a,const unsigned long & b) { return a + xdouble(b); } 00210 // @} 00211 00214 // @{ 00218 friend inline xdouble operator-(const xdouble & a,const xdouble & b); 00219 00221 friend inline xdouble operator-(const double & a,const xdouble & b) { return xdouble(a) - b; } 00223 friend inline xdouble operator-(const int & a,const xdouble & b) { return xdouble(a) - b; } 00225 friend inline xdouble operator-(const unsigned int & a,const xdouble & b) { return xdouble(a) - b; } 00227 friend inline xdouble operator-(const long & a,const xdouble & b) { return xdouble(a) - b; } 00229 friend inline xdouble operator-(const unsigned long & a,const xdouble & b) { return xdouble(a) - b; } 00231 friend inline xdouble operator-(const xdouble & a,const double & b) { return a - xdouble(b); } 00233 friend inline xdouble operator-(const xdouble & a,const int & b) { return a - xdouble(b); } 00235 friend inline xdouble operator-(const xdouble & a,const unsigned int & b) { return a - xdouble(b); } 00237 friend inline xdouble operator-(const xdouble & a,const long & b) { return a - xdouble(b); } 00239 friend inline xdouble operator-(const xdouble & a,const unsigned long & b) { return a - xdouble(b); } 00241 00250 friend inline xdouble operator*(const xdouble & a,const xdouble & b); 00251 00253 friend inline xdouble operator*(const double & a,const xdouble & b) { return xdouble(a) * b; } 00255 friend inline xdouble operator*(const int & a,const xdouble & b) { return xdouble(a) * b; } 00257 friend inline xdouble operator*(const unsigned int & a,const xdouble & b) { return xdouble(a) * b; } 00259 friend inline xdouble operator*(const long & a,const xdouble & b) { return xdouble(a) * b; } 00261 friend inline xdouble operator*(const unsigned long & a,const xdouble & b) { return xdouble(a) * b; } 00263 friend inline xdouble operator*(const xdouble & a,const double & b) { return a * xdouble(b); } 00265 friend inline xdouble operator*(const xdouble & a,const int & b) { return a * xdouble(b); } 00267 friend inline xdouble operator*(const xdouble & a,const unsigned int & b) { return a * xdouble(b); } 00269 friend inline xdouble operator*(const xdouble & a,const long & b) { return a * xdouble(b); } 00271 friend inline xdouble operator*(const xdouble & a,const unsigned long & b) { return a * xdouble(b); } 00273 00274 00281 friend inline xdouble operator/(const xdouble & a,const xdouble & b); 00282 00284 friend inline xdouble operator/(const double & a,const xdouble & b) { return xdouble(a) / b; } 00286 friend inline xdouble operator/(const int & a,const xdouble & b) { return xdouble(a) / b; } 00288 friend inline xdouble operator/(const unsigned int & a,const xdouble & b) { return xdouble(a) / b; } 00290 friend inline xdouble operator/(const long & a,const xdouble & b) { return xdouble(a) / b; } 00292 friend inline xdouble operator/(const unsigned long & a,const xdouble & b) { return xdouble(a) / b; } 00294 friend inline xdouble operator/(const xdouble & a,const double & b) { return a / xdouble(b); } 00296 friend inline xdouble operator/(const xdouble & a,const int & b) { return a / xdouble(b); } 00298 friend inline xdouble operator/(const xdouble & a,const unsigned int & b) { return a / xdouble(b); } 00300 friend inline xdouble operator/(const xdouble & a,const long & b) { return a / xdouble(b); } 00302 friend inline xdouble operator/(const xdouble & a,const unsigned long & b) { return a / xdouble(b); } 00303 // @} 00304 00305 // @} 00306 00307 00310 // @{ 00312 inline xdouble & operator++() { *this+=1.0; return *this; } 00313 00315 inline xdouble & operator--() { *this-=1.0; return *this; } 00316 00318 inline xdouble operator++(int) { xdouble a(*this); *this=*this + 1.0; return a; } 00319 00321 inline xdouble operator--(int) { xdouble a(*this); *this=*this - 1.0; return a; } 00322 // @} 00323 00326 // @{ 00327 friend inline long xcompare(const xdouble& a, const xdouble& b); 00328 friend inline long xcompare(const double& a, const xdouble& b) { return xcompare(xdouble(a),b); } 00329 friend inline long xcompare(const xdouble& a, const double& b) { return xcompare(a,xdouble(b)); } 00330 00331 friend inline long operator==(const xdouble & a,const xdouble & b) { return xcompare(a,b)==0; } 00332 friend inline long operator!=(const xdouble & a,const xdouble & b) { return xcompare(a,b)!=0; } 00333 friend inline long operator<=(const xdouble & a,const xdouble & b) { return xcompare(a,b)<=0; } 00334 friend inline long operator>=(const xdouble & a,const xdouble & b) { return xcompare(a,b)>=0; } 00335 friend inline long operator<(const xdouble & a,const xdouble & b) { return xcompare(a,b)<0; } 00336 friend inline long operator>(const xdouble & a,const xdouble & b) { return xcompare(a,b)>0; } 00337 00338 friend inline long operator==(const double & a,const xdouble & b) { return xcompare(xdouble(a),b)==0; } 00339 friend inline long operator!=(const double & a,const xdouble & b) { return xcompare(xdouble(a),b)!=0; } 00340 friend inline long operator<=(const double & a,const xdouble & b) { return xcompare(xdouble(a),b)<=0; } 00341 friend inline long operator>=(const double & a,const xdouble & b) { return xcompare(xdouble(a),b)>=0; } 00342 friend inline long operator<(const double & a,const xdouble & b) { return xcompare(xdouble(a),b)<0; } 00343 friend inline long operator>(const double & a,const xdouble & b) { return xcompare(xdouble(a),b)>0; } 00344 00345 friend inline long operator==(const int & a,const xdouble & b) { return xcompare(xdouble(a),b)==0; } 00346 friend inline long operator!=(const int & a,const xdouble & b) { return xcompare(xdouble(a),b)!=0; } 00347 friend inline long operator<=(const int & a,const xdouble & b) { return xcompare(xdouble(a),b)<=0; } 00348 friend inline long operator>=(const int & a,const xdouble & b) { return xcompare(xdouble(a),b)>=0; } 00349 friend inline long operator<(const int & a,const xdouble & b) { return xcompare(xdouble(a),b)<0; } 00350 friend inline long operator>(const int & a,const xdouble & b) { return xcompare(xdouble(a),b)>0; } 00351 00352 friend inline long operator==(const unsigned int & a,const xdouble & b) { return xcompare(xdouble(a),b)==0; } 00353 friend inline long operator!=(const unsigned int & a,const xdouble & b) { return xcompare(xdouble(a),b)!=0; } 00354 friend inline long operator<=(const unsigned int & a,const xdouble & b) { return xcompare(xdouble(a),b)<=0; } 00355 friend inline long operator>=(const unsigned int & a,const xdouble & b) { return xcompare(xdouble(a),b)>=0; } 00356 friend inline long operator<(const unsigned int & a,const xdouble & b) { return xcompare(xdouble(a),b)<0; } 00357 friend inline long operator>(const unsigned int & a,const xdouble & b) { return xcompare(xdouble(a),b)>0; } 00358 00359 friend inline long operator==(const long & a,const xdouble & b) { return xcompare(xdouble(a),b)==0; } 00360 friend inline long operator!=(const long & a,const xdouble & b) { return xcompare(xdouble(a),b)!=0; } 00361 friend inline long operator<=(const long & a,const xdouble & b) { return xcompare(xdouble(a),b)<=0; } 00362 friend inline long operator>=(const long & a,const xdouble & b) { return xcompare(xdouble(a),b)>=0; } 00363 friend inline long operator<(const long & a,const xdouble & b) { return xcompare(xdouble(a),b)<0; } 00364 friend inline long operator>(const long & a,const xdouble & b) { return xcompare(xdouble(a),b)>0; } 00365 00366 friend inline long operator==(const unsigned long & a,const xdouble & b) { return xcompare(xdouble(a),b)==0; } 00367 friend inline long operator!=(const unsigned long & a,const xdouble & b) { return xcompare(xdouble(a),b)!=0; } 00368 friend inline long operator<=(const unsigned long & a,const xdouble & b) { return xcompare(xdouble(a),b)<=0; } 00369 friend inline long operator>=(const unsigned long & a,const xdouble & b) { return xcompare(xdouble(a),b)>=0; } 00370 friend inline long operator<(const unsigned long & a,const xdouble & b) { return xcompare(xdouble(a),b)<0; } 00371 friend inline long operator>(const unsigned long & a,const xdouble & b) { return xcompare(xdouble(a),b)>0; } 00372 00373 friend inline long operator==(const xdouble & a,const double & b) { return xcompare(a,xdouble(b))==0; } 00374 friend inline long operator!=(const xdouble & a,const double & b) { return xcompare(a,xdouble(b))!=0; } 00375 friend inline long operator<=(const xdouble & a,const double & b) { return xcompare(a,xdouble(b))<=0; } 00376 friend inline long operator>=(const xdouble & a,const double & b) { return xcompare(a,xdouble(b))>=0; } 00377 friend inline long operator<(const xdouble & a,const double & b) { return xcompare(a,xdouble(b))<0; } 00378 friend inline long operator>(const xdouble & a,const double & b) { return xcompare(a,xdouble(b))>0; } 00379 00380 friend inline long operator==(const xdouble & a,const int & b) { return xcompare(a,xdouble(b))==0; } 00381 friend inline long operator!=(const xdouble & a,const int & b) { return xcompare(a,xdouble(b))!=0; } 00382 friend inline long operator<=(const xdouble & a,const int & b) { return xcompare(a,xdouble(b))<=0; } 00383 friend inline long operator>=(const xdouble & a,const int & b) { return xcompare(a,xdouble(b))>=0; } 00384 friend inline long operator<(const xdouble & a,const int & b) { return xcompare(a,xdouble(b))<0; } 00385 friend inline long operator>(const xdouble & a,const int & b) { return xcompare(a,xdouble(b))>0; } 00386 00387 friend inline long operator==(const xdouble & a,const unsigned int & b) { return xcompare(a,xdouble(b))==0; } 00388 friend inline long operator!=(const xdouble & a,const unsigned int & b) { return xcompare(a,xdouble(b))!=0; } 00389 friend inline long operator<=(const xdouble & a,const unsigned int & b) { return xcompare(a,xdouble(b))<=0; } 00390 friend inline long operator>=(const xdouble & a,const unsigned int & b) { return xcompare(a,xdouble(b))>=0; } 00391 friend inline long operator<(const xdouble & a,const unsigned int & b) { return xcompare(a,xdouble(b))<0; } 00392 friend inline long operator>(const xdouble & a,const unsigned int & b) { return xcompare(a,xdouble(b))>0; } 00393 00394 friend inline long operator==(const xdouble & a,const long & b) { return xcompare(a,xdouble(b))==0; } 00395 friend inline long operator!=(const xdouble & a,const long & b) { return xcompare(a,xdouble(b))!=0; } 00396 friend inline long operator<=(const xdouble & a,const long & b) { return xcompare(a,xdouble(b))<=0; } 00397 friend inline long operator>=(const xdouble & a,const long & b) { return xcompare(a,xdouble(b))>=0; } 00398 friend inline long operator<(const xdouble & a,const long & b) { return xcompare(a,xdouble(b))<0; } 00399 friend inline long operator>(const xdouble & a,const long & b) { return xcompare(a,xdouble(b))>0; } 00400 00401 friend inline long operator==(const xdouble & a,const unsigned long & b) { return xcompare(a,xdouble(b))==0; } 00402 friend inline long operator!=(const xdouble & a,const unsigned long & b) { return xcompare(a,xdouble(b))!=0; } 00403 friend inline long operator<=(const xdouble & a,const unsigned long & b) { return xcompare(a,xdouble(b))<=0; } 00404 friend inline long operator>=(const xdouble & a,const unsigned long & b) { return xcompare(a,xdouble(b))>=0; } 00405 friend inline long operator<(const xdouble & a,const unsigned long & b) { return xcompare(a,xdouble(b))<0; } 00406 friend inline long operator>(const xdouble & a,const unsigned long & b) { return xcompare(a,xdouble(b))>0; } 00408 00415 friend inline xdouble sqrt(const xdouble&); 00416 friend inline long sign(const xdouble&); 00417 friend inline xdouble trunc(const xdouble&); 00418 friend inline xdouble floor(const xdouble&); 00419 friend inline xdouble ceil(const xdouble&); 00420 friend inline xdouble fabs(const xdouble&); 00421 friend inline void power2(xdouble& z, long e); 00422 friend inline xdouble power(const xdouble& z, long e); 00423 friend inline xdouble power2_xdouble(long e); 00424 friend inline double log(const xdouble& a); 00425 friend inline xdouble xexp(double a); 00427 00431 00438 friend inline void MulAdd(xdouble& z, 00439 const xdouble& a, 00440 const xdouble& b, 00441 const xdouble& c); 00442 00449 friend inline void MulSub(xdouble& z, 00450 const xdouble& a, 00451 const xdouble& b, 00452 const xdouble& c) ; 00454 00456 00459 inline void normalize(); 00460 00464 inline double mantissa() const { return x; } 00465 00469 inline long exponent() const { return e; } 00471 00472 00473 private: 00486 static inline long IsFinite(const double *p); 00487 00494 static inline void ForceToMem(double *p) {} 00495 00509 static inline double _ntl_ldexp(double x, long e); 00510 }; 00511 00512 00513 00515 // Implementation 00517 #if (NTL_DOUBLE_PRECISION > NTL_BITS_PER_LONG) 00518 00519 #define NTL_XD_HBOUND (NTL_FDOUBLE_PRECISION*32.0) 00520 #define NTL_XD_HBOUND_LOG (NTL_DOUBLE_PRECISION+4) 00521 00522 #else 00523 00524 #define NTL_XD_HBOUND (double(1L << (NTL_BITS_PER_LONG - 2))*64.0) 00525 #define NTL_XD_HBOUND_LOG (NTL_BITS_PER_LONG+4) 00526 00527 #endif 00528 00529 #define NTL_XD_HBOUND_INV (double(1)/NTL_XD_HBOUND) 00530 00531 #define NTL_XD_BOUND (NTL_XD_HBOUND*NTL_XD_HBOUND) 00532 #define NTL_XD_BOUND_INV (double(1)/NTL_XD_BOUND) 00533 00534 #define NTL_OVFBND (1L << (NTL_BITS_PER_LONG-4)) 00535 #define NTL_OVERFLOW(n, a, b) \ 00536 (((b) >= NTL_OVFBND) || (((long) (n)) > 0 && (((a) >= NTL_OVFBND) || \ 00537 (((long) (n)) >= (NTL_OVFBND-((long)(b))+((long)(a))-1)/((long)(a)))))) 00538 00539 00540 /* 00541 * An IEEE double x is finite if and only if x - x == 0. 00542 * The function _ntl_IsFinite implements this logic; however, 00543 * it does not completely trust that an optimizing compiler 00544 * really implements this correctly, and so it goes out of its way to 00545 * confuse the compiler. For a good compiler that respects IEEE floating 00546 * point arithmetic, this is not be necessary, but it is better 00547 * to be a bit paranoid. 00548 * 00549 * Like the routine _ntl_ForceToMem below, this routine has the 00550 * side effect of forcing its argument into memory. 00551 */ 00552 long xdouble::IsFinite(const double *p) 00553 { 00554 static double _ntl_IsFinite__local; 00555 double *_ntl_IsFinite__ptr1 = &_ntl_IsFinite__local; 00556 double *_ntl_IsFinite__ptr2 = &_ntl_IsFinite__local; 00557 double *_ntl_IsFinite__ptr3 = &_ntl_IsFinite__local; 00558 double *_ntl_IsFinite__ptr4 = &_ntl_IsFinite__local; 00559 *_ntl_IsFinite__ptr1 = *p; 00560 *_ntl_IsFinite__ptr3 = (*_ntl_IsFinite__ptr2 - *p); 00561 if (*_ntl_IsFinite__ptr4 != 0.0) return 0; 00562 return 1; 00563 } 00564 00565 00566 void xdouble::Error(const char * str) 00567 { 00568 std::cerr << str << std::endl; 00569 exit(8); 00570 } 00571 00572 double xdouble::_ntl_ldexp(double x, long e) 00573 { 00574 static const double _ntl_ldexp_zero = 0.0; 00575 if (e > NTL_MAX_INT) 00576 return x/_ntl_ldexp_zero; 00577 else if (e < NTL_MIN_INT) 00578 return x*_ntl_ldexp_zero; 00579 else 00580 return ldexp(x, ((int) e)); 00581 } 00582 00583 00585 void xdouble::from_double(double a) 00586 { 00587 00588 if (a == 0 || a == 1 || (a > 0 && a >= NTL_XD_HBOUND_INV && a <= NTL_XD_HBOUND) 00589 || (a < 0 && a <= -NTL_XD_HBOUND_INV && a >= -NTL_XD_HBOUND)) 00590 { 00591 x = a; 00592 e = 0; 00593 } 00594 if (!xdouble::IsFinite(&a)) 00595 { 00596 xdouble::Error("double to xdouble conversion: non finite value"); 00597 } 00598 x=a; 00599 e=0; 00600 normalize(); 00601 } 00602 00603 void xdouble::normalize() 00604 { 00605 if (x == 0) 00606 e = 0; 00607 else if (x > 0) { 00608 while (x < NTL_XD_HBOUND_INV) { x *= NTL_XD_BOUND; e--; } 00609 while (x > NTL_XD_HBOUND) { x *= NTL_XD_BOUND_INV; e++; } 00610 } 00611 else { 00612 while (x > -NTL_XD_HBOUND_INV) { x *= NTL_XD_BOUND; e--; } 00613 while (x < -NTL_XD_HBOUND) { x *= NTL_XD_BOUND_INV; e++; } 00614 } 00615 00616 if (e >= NTL_OVFBND) 00617 Error("xdouble: overflow"); 00618 00619 if (e <= -NTL_OVFBND) 00620 Error("xdouble: underflow"); 00621 } 00622 00623 double xdouble::to_double() const 00624 { 00625 double tx; 00626 long te; 00627 00628 tx = x; 00629 te = e; 00630 00631 while (te > 0) { tx *= NTL_XD_BOUND; te--; } 00632 while (te < 0) { tx *= NTL_XD_BOUND_INV; te++; } 00633 00634 return tx; 00635 } 00636 00637 00638 xdouble& xdouble::operator=(double a) { from_double(a); return *this; } 00639 xdouble& xdouble::operator=(int a) { from_double(a); return *this;} 00640 xdouble& xdouble::operator=(long a) { from_double(a); return *this; } 00641 xdouble& xdouble::operator=(unsigned int & a) { from_double(a); return *this; } 00642 xdouble& xdouble::operator=(unsigned long & a) { from_double(a); return *this; } 00643 00644 inline xdouble operator+(const xdouble& a, const xdouble& b) 00645 { 00646 xdouble z; 00647 00648 if (a.x == 0) 00649 return b; 00650 00651 if (b.x == 0) 00652 return a; 00653 00654 00655 if (a.e == b.e) { 00656 z.x = a.x + b.x; 00657 z.e = a.e; 00658 z.normalize(); 00659 return z; 00660 } 00661 else if (a.e > b.e) { 00662 if (a.e > b.e+1) 00663 return a; 00664 00665 z.x = a.x + b.x*NTL_XD_BOUND_INV; 00666 z.e = a.e; 00667 z.normalize(); 00668 return z; 00669 } 00670 else { 00671 if (b.e > a.e+1) 00672 return b; 00673 00674 z.x = a.x*NTL_XD_BOUND_INV + b.x; 00675 z.e = b.e; 00676 z.normalize(); 00677 return z; 00678 } 00679 } 00680 00681 00682 inline xdouble operator-(const xdouble& a, const xdouble& b) 00683 { 00684 xdouble z; 00685 00686 if (a.x == 0) 00687 return -b; 00688 00689 if (b.x == 0) 00690 return a; 00691 00692 if (a.e == b.e) { 00693 z.x = a.x - b.x; 00694 z.e = a.e; 00695 z.normalize(); 00696 return z; 00697 } 00698 else if (a.e > b.e) { 00699 if (a.e > b.e+1) 00700 return a; 00701 00702 z.x = a.x - b.x*NTL_XD_BOUND_INV; 00703 z.e = a.e; 00704 z.normalize(); 00705 return z; 00706 } 00707 else { 00708 if (b.e > a.e+1) 00709 return -b; 00710 00711 z.x = a.x*NTL_XD_BOUND_INV - b.x; 00712 z.e = b.e; 00713 z.normalize(); 00714 return z; 00715 } 00716 } 00717 00718 inline xdouble operator-(const xdouble& a) 00719 { 00720 xdouble z; 00721 z.x = -a.x; 00722 z.e = a.e; 00723 return z; 00724 } 00725 00726 inline xdouble operator*(const xdouble& a, const xdouble& b) 00727 { 00728 xdouble z; 00729 00730 z.e = a.e + b.e; 00731 z.x = a.x * b.x; 00732 z.normalize(); 00733 return z; 00734 } 00735 00736 inline xdouble operator/(const xdouble& a, const xdouble& b) 00737 { 00738 xdouble z; 00739 00740 if (b.x == 0) xdouble::Error("xdouble division by 0"); 00741 00742 z.e = a.e - b.e; 00743 z.x = a.x / b.x; 00744 z.normalize(); 00745 return z; 00746 } 00747 00749 inline xdouble sqrt(const xdouble& a) 00750 { 00751 if (a == 0) 00752 return 0; 00753 00754 if (a < 0) 00755 xdouble::Error("xdouble: sqrt of negative number"); 00756 00757 xdouble t; 00758 00759 if (a.e & 1) { 00760 t.e = (a.e - 1)/2; 00761 t.x = sqrt(a.x * NTL_XD_BOUND); 00762 } 00763 else { 00764 t.e = a.e/2; 00765 t.x = sqrt(a.x); 00766 } 00767 00768 t.normalize(); 00769 00770 return t; 00771 } 00772 00773 00774 inline long xcompare(const xdouble& a, const xdouble& b) 00775 { 00776 xdouble z = a - b; 00777 00778 if (z.x < 0) 00779 return -1; 00780 else if (z.x == 0) 00781 return 0; 00782 else 00783 return 1; 00784 } 00785 00786 inline long sign(const xdouble& z) 00787 { 00788 if (z.x < 0) 00789 return -1; 00790 else if (z.x == 0) 00791 return 0; 00792 else 00793 return 1; 00794 } 00795 00796 inline xdouble trunc(const xdouble& a) 00797 { 00798 if (a.x >= 0) 00799 return floor(a); 00800 else 00801 return ceil(a); 00802 } 00803 00804 inline xdouble floor(const xdouble& aa) 00805 { 00806 xdouble z; 00807 00808 xdouble a = aa; 00809 xdouble::ForceToMem(&a.x); 00810 00811 if (a.e == 0) { 00812 z.x = floor(a.x); 00813 z.e = 0; 00814 z.normalize(); 00815 return z; 00816 } 00817 else if (a.e > 0) { 00818 return a; 00819 } 00820 else { 00821 if (a.x < 0) 00822 return -1; 00823 else 00824 return 0; 00825 } 00826 } 00827 00828 00829 inline xdouble ceil(const xdouble& aa) 00830 { 00831 xdouble z; 00832 00833 xdouble a = aa; 00834 xdouble::ForceToMem(&a.x); 00835 00836 if (a.e == 0) { 00837 z.x = ceil(a.x); 00838 z.e = 0; 00839 z.normalize(); 00840 return z; 00841 } 00842 else if (a.e > 0) { 00843 return a; 00844 } 00845 else { 00846 if (a.x < 0) 00847 return 0; 00848 else 00849 return 1; 00850 } 00851 } 00852 00853 00854 inline xdouble fabs(const xdouble& a) 00855 { 00856 xdouble z; 00857 00858 z.e = a.e; 00859 z.x = fabs(a.x); 00860 return z; 00861 } 00862 00863 inline xdouble power(const xdouble& z, long e) 00864 { 00865 xdouble x(z); 00866 power2(x,e); 00867 return x; 00868 } 00869 00870 inline void power2(xdouble& z, long e) 00871 { 00872 long hb = NTL_XD_HBOUND_LOG; 00873 long b = 2*hb; 00874 00875 long q, r; 00876 00877 q = e/b; 00878 r = e%b; 00879 00880 while (r >= hb) { 00881 r -= b; 00882 q++; 00883 } 00884 00885 while (r < -hb) { 00886 r += b; 00887 q--; 00888 } 00889 00890 if (q >= NTL_OVFBND) 00891 xdouble::Error("xdouble: overflow"); 00892 00893 if (q <= -NTL_OVFBND) 00894 xdouble::Error("xdouble: underflow"); 00895 00896 double x = xdouble::_ntl_ldexp(1.0, r); 00897 00898 z.x = x; 00899 z.e = q; 00900 } 00901 00902 inline void MulAdd(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c) 00903 // z = a + b*c 00904 { 00905 double x; 00906 long e; 00907 00908 e = b.e + c.e; 00909 x = b.x * c.x; 00910 00911 if (x == 0) { 00912 z = a; 00913 return; 00914 } 00915 00916 if (a.x == 0) { 00917 z.e = e; 00918 z.x = x; 00919 z.normalize(); 00920 return; 00921 } 00922 00923 00924 if (a.e == e) { 00925 z.x = a.x + x; 00926 z.e = e; 00927 z.normalize(); 00928 return; 00929 } 00930 else if (a.e > e) { 00931 if (a.e > e+1) { 00932 z = a; 00933 return; 00934 } 00935 00936 z.x = a.x + x*NTL_XD_BOUND_INV; 00937 z.e = a.e; 00938 z.normalize(); 00939 return; 00940 } 00941 else { 00942 if (e > a.e+1) { 00943 z.x = x; 00944 z.e = e; 00945 z.normalize(); 00946 return; 00947 } 00948 00949 z.x = a.x*NTL_XD_BOUND_INV + x; 00950 z.e = e; 00951 z.normalize(); 00952 return; 00953 } 00954 } 00955 00956 inline void MulSub(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c) 00957 // z = a - b*c 00958 { 00959 double x; 00960 long e; 00961 00962 e = b.e + c.e; 00963 x = b.x * c.x; 00964 00965 if (x == 0) { 00966 z = a; 00967 return; 00968 } 00969 00970 if (a.x == 0) { 00971 z.e = e; 00972 z.x = -x; 00973 z.normalize(); 00974 return; 00975 } 00976 00977 00978 if (a.e == e) { 00979 z.x = a.x - x; 00980 z.e = e; 00981 z.normalize(); 00982 return; 00983 } 00984 else if (a.e > e) { 00985 if (a.e > e+1) { 00986 z = a; 00987 return; 00988 } 00989 00990 z.x = a.x - x*NTL_XD_BOUND_INV; 00991 z.e = a.e; 00992 z.normalize(); 00993 return; 00994 } 00995 else { 00996 if (e > a.e+1) { 00997 z.x = -x; 00998 z.e = e; 00999 z.normalize(); 01000 return; 01001 } 01002 01003 z.x = a.x*NTL_XD_BOUND_INV - x; 01004 z.e = e; 01005 z.normalize(); 01006 return; 01007 } 01008 } 01009 01010 inline double log(const xdouble& a) 01011 { 01012 static double LogBound = log(NTL_XD_BOUND); 01013 if (a.x <= 0) { 01014 xdouble::Error("log(xdouble): argument must be positive"); 01015 } 01016 01017 return log(a.x) + a.e*LogBound; 01018 } 01019 01020 inline xdouble xexp(double x) 01021 { 01022 const double LogBound = log(NTL_XD_BOUND); 01023 01024 double y = x/LogBound; 01025 double iy = floor(y+0.5); 01026 01027 if (iy >= NTL_OVFBND) 01028 xdouble::Error("xdouble: overflow"); 01029 01030 if (iy <= -NTL_OVFBND) 01031 xdouble::Error("xdouble: underflow"); 01032 01033 01034 double fy = y - iy; 01035 01036 xdouble res; 01037 res.e = long(iy); 01038 res.x = exp(fy*LogBound); 01039 res.normalize(); 01040 return res; 01041 } 01042 01043 01044 inline xdouble power2_xdouble(long e) 01045 { xdouble z; power2(z, e); return z; } 01046 01047 01048 inline void MulAdd(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c); 01049 inline xdouble MulAdd(const xdouble& a, const xdouble& b, 01050 const xdouble& c) 01051 { xdouble z; MulAdd(z, a, b, c); return z; } 01052 01053 01054 inline void MulSub(xdouble& z, const xdouble& a, const xdouble& b, const xdouble& c); 01055 inline xdouble MulSub(const xdouble& a, const xdouble& b, 01056 const xdouble& c) 01057 { xdouble z; MulSub(z, a, b, c); return z; } 01058 01059 01060 #endif