include/qd/qd_real.h

0001 /*
0002  * include/qd_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  * Quad-double precision (>= 212-bit significand) floating point arithmetic
0011  * package, written in ANSI C++, taking full advantage of operator overloading.
0012  * Uses similar techniques as that of David Bailey's double-double package
0013  * and that of Jonathan Shewchuk's adaptive precision floating point
0014  * arithmetic package.  See
0015  *
0016  *   http://www.nersc.gov/~dhbailey/mpdist/mpdist.html
0017  *   http://www.cs.cmu.edu/~quake/robust.html
0018  *
0019  * for more details.
0020  *
0021  * Yozo Hida
0022  */
0023 #ifndef _QD_QD_REAL_H
0024 #define _QD_QD_REAL_H
0025
0026 #include <iostream>
0027 #include <string>
0028 #include <limits>
0029 #include <qd/qd_config.h>
0030 #include <qd/dd_real.h>
0031
0032 struct QD_API qd_real {
0033   double x[4];    /* The Components. */
0034
0035   /* Eliminates any zeros in the middle component(s). */
0036   void zero_elim();
0037   void zero_elim(double &e);
0038
0039   void renorm();
0040   void renorm(double &e);
0041
0042   void quick_accum(double d, double &e);
0043   void quick_prod_accum(double a, double b, double &e);
0044
0045   qd_real(double x0, double x1, double x2, double x3);
0046   explicit qd_real(const double *xx);
0047
0048   static const qd_real _2pi;
0049   static const qd_real _pi;
0050   static const qd_real _3pi4;
0051   static const qd_real _pi2;
0052   static const qd_real _pi4;
0053   static const qd_real _e;
0054   static const qd_real _log2;
0055   static const qd_real _log10;
0056   static const qd_real _nan;
0057   static const qd_real _inf;
0058
0059   static const double _eps;
0060   static const double _min_normalized;
0061   static const qd_real _max;
0062   static const qd_real _safe_max;
0063   static const int _ndigits;
0064
0065   qd_real();
0066   qd_real(const char *s);
0067   qd_real(const dd_real &dd);
0068   qd_real(double d);
0069   qd_real(int i);
0070
0071   double operator[](int i) const;
0072   double &operator[](int i);
0073
0074   static void error(const char *msg);
0075
0076   bool isnan() const;
0077   bool isfinite() const { return QD_ISFINITE(x[0]); }
0078   bool isinf() const { return QD_ISINF(x[0]); }
0079
0080   static qd_real ieee_add(const qd_real &a, const qd_real &b);
0081   static qd_real sloppy_add(const qd_real &a, const qd_real &b);
0082
0083   qd_real &operator+=(double a);
0084   qd_real &operator+=(const dd_real &a);
0085   qd_real &operator+=(const qd_real &a);
0086
0087   qd_real &operator-=(double a);
0088   qd_real &operator-=(const dd_real &a);
0089   qd_real &operator-=(const qd_real &a);
0090
0091   static qd_real sloppy_mul(const qd_real &a, const qd_real &b);
0092   static qd_real accurate_mul(const qd_real &a, const qd_real &b);
0093
0094   qd_real &operator*=(double a);
0095   qd_real &operator*=(const dd_real &a);
0096   qd_real &operator*=(const qd_real &a);
0097
0098   static qd_real sloppy_div(const qd_real &a, const dd_real &b);
0099   static qd_real accurate_div(const qd_real &a, const dd_real &b);
0100   static qd_real sloppy_div(const qd_real &a, const qd_real &b);
0101   static qd_real accurate_div(const qd_real &a, const qd_real &b);
0102
0103   qd_real &operator/=(double a);
0104   qd_real &operator/=(const dd_real &a);
0105   qd_real &operator/=(const qd_real &a);
0106
0107   qd_real operator^(int n) const;
0108
0109   qd_real operator-() const;
0110
0111   qd_real &operator=(double a);
0112   qd_real &operator=(const dd_real &a);
0113   qd_real &operator=(const char *s);
0114
0115   bool is_zero() const;
0116   bool is_one() const;
0117   bool is_positive() const;
0118   bool is_negative() const;
0119
0120   static qd_real rand(void);
0121
0122   void to_digits(char *s, int &expn, int precision = _ndigits) const;
0123   void write(char *s, int len, int precision = _ndigits,
0124       bool showpos = false, bool uppercase = false) const;
0125   std::string to_string(int precision = _ndigits, int width = 0,
0126       std::ios_base::fmtflags fmt = static_cast<std::ios_base::fmtflags>(0),
0127       bool showpos = false, bool uppercase = false, char fill = ' ') const;
0128   static int read(const char *s, qd_real &a);
0129
0130   /* Debugging methods */
0131   void dump(const std::string &name = "", std::ostream &os = std::cerr) const;
0132   void dump_bits(const std::string &name = "",
0133                  std::ostream &os = std::cerr) const;
0134
0135   static qd_real debug_rand();
0136
0137 };
0138
0139 namespace std {
0140   template <>
0141   class numeric_limits<qd_real> : public numeric_limits<double> {
0142   public:
0143     inline static double epsilon() { return qd_real::_eps; }
0144     inline static double min() { return qd_real::_min_normalized; }
0145     inline static qd_real max() { return qd_real::_max; }
0146     inline static qd_real safe_max() { return qd_real::_safe_max; }
0147     static const int digits = 209;
0148     static const int digits10 = 62;
0149   };
0150 }
0151
0152 QD_API qd_real polyeval(const qd_real *c, int n, const qd_real &x);
0153 QD_API qd_real polyroot(const qd_real *c, int n,
0154     const qd_real &x0, int max_iter = 64, double thresh = 0.0);
0155
0156 QD_API qd_real qdrand(void);
0157 QD_API qd_real sqrt(const qd_real &a);
0158
0159 QD_API inline bool isnan(const qd_real &a) { return a.isnan(); }
0160 QD_API inline bool isfinite(const qd_real &a) { return a.isfinite(); }
0161 QD_API inline bool isinf(const qd_real &a) { return a.isinf(); }
0162
0163 /* Computes  qd * d  where d is known to be a power of 2.
0164    This can be done component wise.                      */
0165 QD_API qd_real mul_pwr2(const qd_real &qd, double d);
0166
0167 QD_API qd_real operator+(const qd_real &a, const qd_real &b);
0168 QD_API qd_real operator+(const dd_real &a, const qd_real &b);
0169 QD_API qd_real operator+(const qd_real &a, const dd_real &b);
0170 QD_API qd_real operator+(const qd_real &a, double b);
0171 QD_API qd_real operator+(double a, const qd_real &b);
0172
0173 QD_API qd_real operator-(const qd_real &a, const qd_real &b);
0174 QD_API qd_real operator-(const dd_real &a, const qd_real &b);
0175 QD_API qd_real operator-(const qd_real &a, const dd_real &b);
0176 QD_API qd_real operator-(const qd_real &a, double b);
0177 QD_API qd_real operator-(double a, const qd_real &b);
0178
0179 QD_API qd_real operator*(const qd_real &a, const qd_real &b);
0180 QD_API qd_real operator*(const dd_real &a, const qd_real &b);
0181 QD_API qd_real operator*(const qd_real &a, const dd_real &b);
0182 QD_API qd_real operator*(const qd_real &a, double b);
0183 QD_API qd_real operator*(double a, const qd_real &b);
0184
0185 QD_API qd_real operator/(const qd_real &a, const qd_real &b);
0186 QD_API qd_real operator/(const dd_real &a, const qd_real &b);
0187 QD_API qd_real operator/(const qd_real &a, const dd_real &b);
0188 QD_API qd_real operator/(const qd_real &a, double b);
0189 QD_API qd_real operator/(double a, const qd_real &b);
0190
0191 QD_API qd_real sqr(const qd_real &a);
0192 QD_API qd_real sqrt(const qd_real &a);
0193 QD_API qd_real fsqrt(const qd_real &a, int &flag);
0194 QD_API qd_real pow(const qd_real &a, int n);
0195 QD_API qd_real pow(const qd_real &a, const qd_real &b);
0196 QD_API qd_real npwr(const qd_real &a, int n);
0197
0198 QD_API qd_real nroot(const qd_real &a, int n);
0199
0200 QD_API qd_real rem(const qd_real &a, const qd_real &b);
0201 QD_API qd_real drem(const qd_real &a, const qd_real &b);
0202 QD_API qd_real divrem(const qd_real &a, const qd_real &b, qd_real &r);
0203
0204 dd_real to_dd_real(const qd_real &a);
0205 double  to_double(const qd_real &a);
0206 int     to_int(const qd_real &a);
0207
0208 QD_API bool operator==(const qd_real &a, const qd_real &b);
0209 QD_API bool operator==(const qd_real &a, const dd_real &b);
0210 QD_API bool operator==(const dd_real &a, const qd_real &b);
0211 QD_API bool operator==(double a, const qd_real &b);
0212 QD_API bool operator==(const qd_real &a, double b);
0213
0214 QD_API bool operator<(const qd_real &a, const qd_real &b);
0215 QD_API bool operator<(const qd_real &a, const dd_real &b);
0216 QD_API bool operator<(const dd_real &a, const qd_real &b);
0217 QD_API bool operator<(double a, const qd_real &b);
0218 QD_API bool operator<(const qd_real &a, double b);
0219
0220 QD_API bool operator>(const qd_real &a, const qd_real &b);
0221 QD_API bool operator>(const qd_real &a, const dd_real &b);
0222 QD_API bool operator>(const dd_real &a, const qd_real &b);
0223 QD_API bool operator>(double a, const qd_real &b);
0224 QD_API bool operator>(const qd_real &a, double b);
0225
0226 QD_API bool operator<=(const qd_real &a, const qd_real &b);
0227 QD_API bool operator<=(const qd_real &a, const dd_real &b);
0228 QD_API bool operator<=(const dd_real &a, const qd_real &b);
0229 QD_API bool operator<=(double a, const qd_real &b);
0230 QD_API bool operator<=(const qd_real &a, double b);
0231
0232 QD_API bool operator>=(const qd_real &a, const qd_real &b);
0233 QD_API bool operator>=(const qd_real &a, const dd_real &b);
0234 QD_API bool operator>=(const dd_real &a, const qd_real &b);
0235 QD_API bool operator>=(double a, const qd_real &b);
0236 QD_API bool operator>=(const qd_real &a, double b);
0237
0238 QD_API bool operator!=(const qd_real &a, const qd_real &b);
0239 QD_API bool operator!=(const qd_real &a, const dd_real &b);
0240 QD_API bool operator!=(const dd_real &a, const qd_real &b);
0241 QD_API bool operator!=(double a, const qd_real &b);
0242 QD_API bool operator!=(const qd_real &a, double b);
0243
0244 QD_API qd_real fabs(const qd_real &a);
0245 QD_API qd_real abs(const qd_real &a);    /* same as fabs */
0246
0247 QD_API qd_real ldexp(const qd_real &a, int n);
0248
0249 QD_API qd_real nint(const qd_real &a);
0250 QD_API qd_real quick_nint(const qd_real &a);
0251 QD_API qd_real floor(const qd_real &a);
0252 QD_API qd_real ceil(const qd_real &a);
0253 QD_API qd_real aint(const qd_real &a);
0254
0255 QD_API qd_real sin(const qd_real &a);
0256 QD_API qd_real cos(const qd_real &a);
0257 QD_API qd_real tan(const qd_real &a);
0258 QD_API void sincos(const qd_real &a, qd_real &s, qd_real &c);
0259
0260 QD_API qd_real asin(const qd_real &a);
0261 QD_API qd_real acos(const qd_real &a);
0262 QD_API qd_real atan(const qd_real &a);
0263 QD_API qd_real atan2(const qd_real &y, const qd_real &x);
0264
0265 QD_API qd_real exp(const qd_real &a);
0266 QD_API qd_real log(const qd_real &a);
0267 QD_API qd_real log10(const qd_real &a);
0268
0269 QD_API qd_real sinh(const qd_real &a);
0270 QD_API qd_real cosh(const qd_real &a);
0271 QD_API qd_real tanh(const qd_real &a);
0272 QD_API void sincosh(const qd_real &a, qd_real &sin_qd, qd_real &cos_qd);
0273
0274 QD_API qd_real asinh(const qd_real &a);
0275 QD_API qd_real acosh(const qd_real &a);
0276 QD_API qd_real atanh(const qd_real &a);
0277
0278 QD_API qd_real qdrand(void);
0279
0280 QD_API qd_real max(const qd_real &a, const qd_real &b);
0281 QD_API qd_real max(const qd_real &a, const qd_real &b, const qd_real &c);
0282 QD_API qd_real min(const qd_real &a, const qd_real &b);
0283 QD_API qd_real min(const qd_real &a, const qd_real &b, const qd_real &c);
0284
0285 QD_API qd_real fmod(const qd_real &a, const qd_real &b);
0286
0287 QD_API std::ostream &operator<<(std::ostream &s, const qd_real &a);
0288 QD_API std::istream &operator>>(std::istream &s, qd_real &a);
0289 #ifdef QD_INLINE
0290 #include <qd/qd_inline.h>
0291 #endif
0292
0293 #endif /* _QD_QD_REAL_H */
0294