include/qd/dd_real.h

0001 /*
0002  * include/dd_real.h
0003  *
0004  * This work was supported by the Director, Office of Science, Division
0005  * of Mathematical, Information, and Computational Sciences of the
0006  * U.S. Department of Energy under contract number DE-AC03-76SF00098.
0007  *
0008  * Copyright (c) 2000-2007
0009  *
0010  * Double-double precision (>= 106-bit significand) floating point
0011  * arithmetic package based on David Bailey's Fortran-90 double-double
0012  * package, with some changes. See
0013  *
0014  *   http://www.nersc.gov/~dhbailey/mpdist/mpdist.html
0015  *
0016  * for the original Fortran-90 version.
0017  *
0018  * Overall structure is similar to that of Keith Brigg's C++ double-double
0019  * package.  See
0020  *
0021  *   http://www-epidem.plansci.cam.ac.uk/~kbriggs/doubledouble.html
0022  *
0023  * for more details.  In particular, the fix for x86 computers is borrowed
0024  * from his code.
0025  *
0026  * Yozo Hida
0027  */
0028
0029 #ifndef _QD_DD_REAL_H
0030 #define _QD_DD_REAL_H
0031
0032 #include <cmath>
0033 #include <iostream>
0034 #include <string>
0035 #include <limits>
0036 #include <qd/qd_config.h>
0037 #include <qd/fpu.h>
0038
0039 // Some compilers define isnan, isfinite, and isinf as macros, even for
0040 // C++ codes, which cause havoc when overloading these functions.  We undef
0041 // them here.
0042 #ifdef isnan
0043 #undef isnan
0044 #endif
0045
0046 #ifdef isfinite
0047 #undef isfinite
0048 #endif
0049
0050 #ifdef isinf
0051 #undef isinf
0052 #endif
0053
0054 #ifdef max
0055 #undef max
0056 #endif
0057
0058 #ifdef min
0059 #undef min
0060 #endif
0061
0062 struct QD_API dd_real {
0063   double x[2];
0064
0065   dd_real(double hi, double lo) { x[0] = hi; x[1] = lo; }
0066   dd_real() {x[0] = 0.0; x[1] = 0.0; }
0067   dd_real(double h) { x[0] = h; x[1] = 0.0; }
0068   dd_real(int h) {
0069     x[0] = (static_cast<double>(h));
0070     x[1] = 0.0;
0071   }
0072
0073   dd_real (const char *s);
0074   explicit dd_real (const double *d) {
0075     x[0] = d[0]; x[1] = d[1];
0076   }
0077
0078   static void error(const char *msg);
0079
0080   double _hi() const { return x[0]; }
0081   double _lo() const { return x[1]; }
0082
0083   static const dd_real _2pi;
0084   static const dd_real _pi;
0085   static const dd_real _3pi4;
0086   static const dd_real _pi2;
0087   static const dd_real _pi4;
0088   static const dd_real _e;
0089   static const dd_real _log2;
0090   static const dd_real _log10;
0091   static const dd_real _nan;
0092   static const dd_real _inf;
0093
0094   static const double _eps;
0095   static const double _min_normalized;
0096   static const dd_real _max;
0097   static const dd_real _safe_max;
0098   static const int _ndigits;
0099
0100   bool isnan() const { return QD_ISNAN(x[0]) || QD_ISNAN(x[1]); }
0101   bool isfinite() const { return QD_ISFINITE(x[0]); }
0102   bool isinf() const { return QD_ISINF(x[0]); }
0103
0104   static dd_real add(double a, double b);
0105   static dd_real ieee_add(const dd_real &a, const dd_real &b);
0106   static dd_real sloppy_add(const dd_real &a, const dd_real &b);
0107
0108   dd_real &operator+=(double a);
0109   dd_real &operator+=(const dd_real &a);
0110
0111   static dd_real sub(double a, double b);
0112
0113   dd_real &operator-=(double a);
0114   dd_real &operator-=(const dd_real &a);
0115
0116   dd_real operator-() const;
0117
0118   static dd_real mul(double a, double b);
0119
0120   dd_real &operator*=(double a);
0121   dd_real &operator*=(const dd_real &a);
0122
0123   static dd_real div(double a, double b);
0124   static dd_real sloppy_div(const dd_real &a, const dd_real &b);
0125   static dd_real accurate_div(const dd_real &a, const dd_real &b);
0126
0127   dd_real &operator/=(double a);
0128   dd_real &operator/=(const dd_real &a);
0129
0130   dd_real &operator=(double a);
0131   dd_real &operator=(const char *s);
0132
0133   dd_real operator^(int n);
0134   static dd_real sqr(double d);
0135
0136   static dd_real sqrt(double a);
0137
0138   bool is_zero() const;
0139   bool is_one() const;
0140   bool is_positive() const;
0141   bool is_negative() const;
0142
0143   static dd_real rand(void);
0144
0145   void to_digits(char *s, int &expn, int precision = _ndigits) const;
0146   void write(char *s, int len, int precision = _ndigits,
0147       bool showpos = false, bool uppercase = false) const;
0148   std::string to_string(int precision = _ndigits, int width = 0,
0149       std::ios_base::fmtflags fmt = static_cast<std::ios_base::fmtflags>(0),
0150       bool showpos = false, bool uppercase = false, char fill = ' ') const;
0151   int read(const char *s, dd_real &a);
0152
0153   /* Debugging Methods */
0154   void dump(const std::string &name = "", std::ostream &os = std::cerr) const;
0155   void dump_bits(const std::string &name = "",
0156                  std::ostream &os = std::cerr) const;
0157
0158   static dd_real debug_rand();
0159 };
0160
0161
0162 namespace std {
0163   template <>
0164   class numeric_limits<dd_real> : public numeric_limits<double> {
0165   public:
0166     inline static double epsilon() { return dd_real::_eps; }
0167     inline static dd_real max() { return dd_real::_max; }
0168     inline static dd_real safe_max() { return dd_real::_safe_max; }
0169     inline static double min() { return dd_real::_min_normalized; }
0170     static const int digits = 104;
0171     static const int digits10 = 31;
0172   };
0173 }
0174
0175 QD_API dd_real ddrand(void);
0176 QD_API dd_real sqrt(const dd_real &a);
0177
0178 QD_API dd_real polyeval(const dd_real *c, int n, const dd_real &x);
0179 QD_API dd_real polyroot(const dd_real *c, int n,
0180     const dd_real &x0, int max_iter = 32, double thresh = 0.0);
0181
0182 QD_API inline bool isnan(const dd_real &a) { return a.isnan(); }
0183 QD_API inline bool isfinite(const dd_real &a) { return a.isfinite(); }
0184 QD_API inline bool isinf(const dd_real &a) { return a.isinf(); }
0185
0186 /* Computes  dd * d  where d is known to be a power of 2. */
0187 QD_API dd_real mul_pwr2(const dd_real &dd, double d);
0188
0189 QD_API dd_real operator+(const dd_real &a, double b);
0190 QD_API dd_real operator+(double a, const dd_real &b);
0191 QD_API dd_real operator+(const dd_real &a, const dd_real &b);
0192
0193 QD_API dd_real operator-(const dd_real &a, double b);
0194 QD_API dd_real operator-(double a, const dd_real &b);
0195 QD_API dd_real operator-(const dd_real &a, const dd_real &b);
0196
0197 QD_API dd_real operator*(const dd_real &a, double b);
0198 QD_API dd_real operator*(double a, const dd_real &b);
0199 QD_API dd_real operator*(const dd_real &a, const dd_real &b);
0200
0201 QD_API dd_real operator/(const dd_real &a, double b);
0202 QD_API dd_real operator/(double a, const dd_real &b);
0203 QD_API dd_real operator/(const dd_real &a, const dd_real &b);
0204
0205 QD_API dd_real inv(const dd_real &a);
0206
0207 QD_API dd_real rem(const dd_real &a, const dd_real &b);
0208 QD_API dd_real drem(const dd_real &a, const dd_real &b);
0209 QD_API dd_real divrem(const dd_real &a, const dd_real &b, dd_real &r);
0210
0211 QD_API dd_real pow(const dd_real &a, int n);
0212 QD_API dd_real pow(const dd_real &a, const dd_real &b);
0213 QD_API dd_real npwr(const dd_real &a, int n);
0214 QD_API dd_real sqr(const dd_real &a);
0215
0216 QD_API dd_real sqrt(const dd_real &a);
0217 QD_API dd_real nroot(const dd_real &a, int n);
0218
0219 QD_API bool operator==(const dd_real &a, double b);
0220 QD_API bool operator==(double a, const dd_real &b);
0221 QD_API bool operator==(const dd_real &a, const dd_real &b);
0222
0223 QD_API bool operator<=(const dd_real &a, double b);
0224 QD_API bool operator<=(double a, const dd_real &b);
0225 QD_API bool operator<=(const dd_real &a, const dd_real &b);
0226
0227 QD_API bool operator>=(const dd_real &a, double b);
0228 QD_API bool operator>=(double a, const dd_real &b);
0229 QD_API bool operator>=(const dd_real &a, const dd_real &b);
0230
0231 QD_API bool operator<(const dd_real &a, double b);
0232 QD_API bool operator<(double a, const dd_real &b);
0233 QD_API bool operator<(const dd_real &a, const dd_real &b);
0234
0235 QD_API bool operator>(const dd_real &a, double b);
0236 QD_API bool operator>(double a, const dd_real &b);
0237 QD_API bool operator>(const dd_real &a, const dd_real &b);
0238
0239 QD_API bool operator!=(const dd_real &a, double b);
0240 QD_API bool operator!=(double a, const dd_real &b);
0241 QD_API bool operator!=(const dd_real &a, const dd_real &b);
0242
0243 QD_API dd_real nint(const dd_real &a);
0244 QD_API dd_real floor(const dd_real &a);
0245 QD_API dd_real ceil(const dd_real &a);
0246 QD_API dd_real aint(const dd_real &a);
0247
0248 QD_API dd_real ddrand(void);
0249
0250 double to_double(const dd_real &a);
0251 int    to_int(const dd_real &a);
0252
0253 QD_API dd_real exp(const dd_real &a);
0254 QD_API dd_real ldexp(const dd_real &a, int exp);
0255 QD_API dd_real log(const dd_real &a);
0256 QD_API dd_real log10(const dd_real &a);
0257
0258 QD_API dd_real sin(const dd_real &a);
0259 QD_API dd_real cos(const dd_real &a);
0260 QD_API dd_real tan(const dd_real &a);
0261 QD_API void sincos(const dd_real &a, dd_real &sin_a, dd_real &cos_a);
0262
0263 QD_API dd_real asin(const dd_real &a);
0264 QD_API dd_real acos(const dd_real &a);
0265 QD_API dd_real atan(const dd_real &a);
0266 QD_API dd_real atan2(const dd_real &y, const dd_real &x);
0267
0268 QD_API dd_real sinh(const dd_real &a);
0269 QD_API dd_real cosh(const dd_real &a);
0270 QD_API dd_real tanh(const dd_real &a);
0271 QD_API void sincosh(const dd_real &a,
0272                       dd_real &sinh_a, dd_real &cosh_a);
0273
0274 QD_API dd_real asinh(const dd_real &a);
0275 QD_API dd_real acosh(const dd_real &a);
0276 QD_API dd_real atanh(const dd_real &a);
0277
0278 QD_API dd_real fabs(const dd_real &a);
0279 QD_API dd_real abs(const dd_real &a);   /* same as fabs */
0280
0281 QD_API dd_real fmod(const dd_real &a, const dd_real &b);
0282
0283 QD_API std::ostream& operator<<(std::ostream &s, const dd_real &a);
0284 QD_API std::istream& operator>>(std::istream &s, dd_real &a);
0285 #ifdef QD_INLINE
0286 #include <qd/dd_inline.h>
0287 #endif
0288
0289 #endif /* _QD_DD_REAL_H */
0290