include/qd/dd_inline.h

0001 /*
0002  * include/dd_inline.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-2001
0009  *
0010  * Contains small functions (suitable for inlining) in the double-double
0011  * arithmetic package.
0012  */
0013 #ifndef _QD_DD_INLINE_H
0014 #define _QD_DD_INLINE_H
0015
0016 #include <cmath>
0017 #include <qd/inline.h>
0018
0019 #ifndef QD_INLINE
0020 #define inline
0021 #endif
0022
0023
0024 /*********** Additions ************/
0025 /* double-double = double + double */
0026 inline dd_real dd_real::add(double a, double b) {
0027   double s, e;
0028   s = qd::two_sum(a, b, e);
0029   return dd_real(s, e);
0030 }
0031
0032 /* double-double + double */
0033 inline dd_real operator+(const dd_real &a, double b) {
0034   double s1, s2;
0035   s1 = qd::two_sum(a.x[0], b, s2);
0036   s2 += a.x[1];
0037   s1 = qd::quick_two_sum(s1, s2, s2);
0038   return dd_real(s1, s2);
0039 }
0040
0041 /* double-double + double-double */
0042 inline dd_real dd_real::ieee_add(const dd_real &a, const dd_real &b) {
0043   /* This one satisfies IEEE style error bound,
0044      due to K. Briggs and W. Kahan.                   */
0045   double s1, s2, t1, t2;
0046
0047   s1 = qd::two_sum(a.x[0], b.x[0], s2);
0048   t1 = qd::two_sum(a.x[1], b.x[1], t2);
0049   s2 += t1;
0050   s1 = qd::quick_two_sum(s1, s2, s2);
0051   s2 += t2;
0052   s1 = qd::quick_two_sum(s1, s2, s2);
0053   return dd_real(s1, s2);
0054 }
0055
0056 inline dd_real dd_real::sloppy_add(const dd_real &a, const dd_real &b) {
0057   /* This is the less accurate version ... obeys Cray-style
0058      error bound. */
0059   double s, e;
0060
0061   s = qd::two_sum(a.x[0], b.x[0], e);
0062   e += (a.x[1] + b.x[1]);
0063   s = qd::quick_two_sum(s, e, e);
0064   return dd_real(s, e);
0065 }
0066
0067 inline dd_real operator+(const dd_real &a, const dd_real &b) {
0068 #ifndef QD_IEEE_ADD
0069   return dd_real::sloppy_add(a, b);
0070 #else
0071   return dd_real::ieee_add(a, b);
0072 #endif
0073 }
0074
0075 /* double + double-double */
0076 inline dd_real operator+(double a, const dd_real &b) {
0077   return (b + a);
0078 }
0079
0080
0081 /*********** Self-Additions ************/
0082 /* double-double += double */
0083 inline dd_real &dd_real::operator+=(double a) {
0084   double s1, s2;
0085   s1 = qd::two_sum(x[0], a, s2);
0086   s2 += x[1];
0087   x[0] = qd::quick_two_sum(s1, s2, x[1]);
0088   return *this;
0089 }
0090
0091 /* double-double += double-double */
0092 inline dd_real &dd_real::operator+=(const dd_real &a) {
0093 #ifndef QD_IEEE_ADD
0094   double s, e;
0095   s = qd::two_sum(x[0], a.x[0], e);
0096   e += x[1];
0097   e += a.x[1];
0098   x[0] = qd::quick_two_sum(s, e, x[1]);
0099   return *this;
0100 #else
0101   double s1, s2, t1, t2;
0102   s1 = qd::two_sum(x[0], a.x[0], s2);
0103   t1 = qd::two_sum(x[1], a.x[1], t2);
0104   s2 += t1;
0105   s1 = qd::quick_two_sum(s1, s2, s2);
0106   s2 += t2;
0107   x[0] = qd::quick_two_sum(s1, s2, x[1]);
0108   return *this;
0109 #endif
0110 }
0111
0112 /*********** Subtractions ************/
0113 /* double-double = double - double */
0114 inline dd_real dd_real::sub(double a, double b) {
0115   double s, e;
0116   s = qd::two_diff(a, b, e);
0117   return dd_real(s, e);
0118 }
0119
0120 /* double-double - double */
0121 inline dd_real operator-(const dd_real &a, double b) {
0122   double s1, s2;
0123   s1 = qd::two_diff(a.x[0], b, s2);
0124   s2 += a.x[1];
0125   s1 = qd::quick_two_sum(s1, s2, s2);
0126   return dd_real(s1, s2);
0127 }
0128
0129 /* double-double - double-double */
0130 inline dd_real operator-(const dd_real &a, const dd_real &b) {
0131 #ifndef QD_IEEE_ADD
0132   double s, e;
0133   s = qd::two_diff(a.x[0], b.x[0], e);
0134   e += a.x[1];
0135   e -= b.x[1];
0136   s = qd::quick_two_sum(s, e, e);
0137   return dd_real(s, e);
0138 #else
0139   double s1, s2, t1, t2;
0140   s1 = qd::two_diff(a.x[0], b.x[0], s2);
0141   t1 = qd::two_diff(a.x[1], b.x[1], t2);
0142   s2 += t1;
0143   s1 = qd::quick_two_sum(s1, s2, s2);
0144   s2 += t2;
0145   s1 = qd::quick_two_sum(s1, s2, s2);
0146   return dd_real(s1, s2);
0147 #endif
0148 }
0149
0150 /* double - double-double */
0151 inline dd_real operator-(double a, const dd_real &b) {
0152   double s1, s2;
0153   s1 = qd::two_diff(a, b.x[0], s2);
0154   s2 -= b.x[1];
0155   s1 = qd::quick_two_sum(s1, s2, s2);
0156   return dd_real(s1, s2);
0157 }
0158
0159 /*********** Self-Subtractions ************/
0160 /* double-double -= double */
0161 inline dd_real &dd_real::operator-=(double a) {
0162   double s1, s2;
0163   s1 = qd::two_diff(x[0], a, s2);
0164   s2 += x[1];
0165   x[0] = qd::quick_two_sum(s1, s2, x[1]);
0166   return *this;
0167 }
0168
0169 /* double-double -= double-double */
0170 inline dd_real &dd_real::operator-=(const dd_real &a) {
0171 #ifndef QD_IEEE_ADD
0172   double s, e;
0173   s = qd::two_diff(x[0], a.x[0], e);
0174   e += x[1];
0175   e -= a.x[1];
0176   x[0] = qd::quick_two_sum(s, e, x[1]);
0177   return *this;
0178 #else
0179   double s1, s2, t1, t2;
0180   s1 = qd::two_diff(x[0], a.x[0], s2);
0181   t1 = qd::two_diff(x[1], a.x[1], t2);
0182   s2 += t1;
0183   s1 = qd::quick_two_sum(s1, s2, s2);
0184   s2 += t2;
0185   x[0] = qd::quick_two_sum(s1, s2, x[1]);
0186   return *this;
0187 #endif
0188 }
0189
0190 /*********** Unary Minus ***********/
0191 inline dd_real dd_real::operator-() const {
0192   return dd_real(-x[0], -x[1]);
0193 }
0194
0195 /*********** Multiplications ************/
0196 /* double-double = double * double */
0197 inline dd_real dd_real::mul(double a, double b) {
0198   double p, e;
0199   p = qd::two_prod(a, b, e);
0200   return dd_real(p, e);
0201 }
0202
0203 /* double-double * (2.0 ^ exp) */
0204 inline dd_real ldexp(const dd_real &a, int exp) {
0205   return dd_real(std::ldexp(a.x[0], exp), std::ldexp(a.x[1], exp));
0206 }
0207
0208 /* double-double * double,  where double is a power of 2. */
0209 inline dd_real mul_pwr2(const dd_real &a, double b) {
0210   return dd_real(a.x[0] * b, a.x[1] * b);
0211 }
0212
0213 /* double-double * double */
0214 inline dd_real operator*(const dd_real &a, double b) {
0215   double p1, p2;
0216
0217   p1 = qd::two_prod(a.x[0], b, p2);
0218   p2 += (a.x[1] * b);
0219   p1 = qd::quick_two_sum(p1, p2, p2);
0220   return dd_real(p1, p2);
0221 }
0222
0223 /* double-double * double-double */
0224 inline dd_real operator*(const dd_real &a, const dd_real &b) {
0225   double p1, p2;
0226
0227   p1 = qd::two_prod(a.x[0], b.x[0], p2);
0228   p2 += (a.x[0] * b.x[1] + a.x[1] * b.x[0]);
0229   p1 = qd::quick_two_sum(p1, p2, p2);
0230   return dd_real(p1, p2);
0231 }
0232
0233 /* double * double-double */
0234 inline dd_real operator*(double a, const dd_real &b) {
0235   return (b * a);
0236 }
0237
0238 /*********** Self-Multiplications ************/
0239 /* double-double *= double */
0240 inline dd_real &dd_real::operator*=(double a) {
0241   double p1, p2;
0242   p1 = qd::two_prod(x[0], a, p2);
0243   p2 += x[1] * a;
0244   x[0] = qd::quick_two_sum(p1, p2, x[1]);
0245   return *this;
0246 }
0247
0248 /* double-double *= double-double */
0249 inline dd_real &dd_real::operator*=(const dd_real &a) {
0250   double p1, p2;
0251   p1 = qd::two_prod(x[0], a.x[0], p2);
0252   p2 += a.x[1] * x[0];
0253   p2 += a.x[0] * x[1];
0254   x[0] = qd::quick_two_sum(p1, p2, x[1]);
0255   return *this;
0256 }
0257
0258 /*********** Divisions ************/
0259 inline dd_real dd_real::div(double a, double b) {
0260   double q1, q2;
0261   double p1, p2;
0262   double s, e;
0263
0264   q1 = a / b;
0265
0266   /* Compute  a - q1 * b */
0267   p1 = qd::two_prod(q1, b, p2);
0268   s = qd::two_diff(a, p1, e);
0269   e -= p2;
0270
0271   /* get next approximation */
0272   q2 = (s + e) / b;
0273
0274   s = qd::quick_two_sum(q1, q2, e);
0275
0276   return dd_real(s, e);
0277 }
0278
0279 /* double-double / double */
0280 inline dd_real operator/(const dd_real &a, double b) {
0281
0282   double q1, q2;
0283   double p1, p2;
0284   double s, e;
0285   dd_real r;
0286
0287   q1 = a.x[0] / b;   /* approximate quotient. */
0288
0289   /* Compute  this - q1 * d */
0290   p1 = qd::two_prod(q1, b, p2);
0291   s = qd::two_diff(a.x[0], p1, e);
0292   e += a.x[1];
0293   e -= p2;
0294
0295   /* get next approximation. */
0296   q2 = (s + e) / b;
0297
0298   /* renormalize */
0299   r.x[0] = qd::quick_two_sum(q1, q2, r.x[1]);
0300
0301   return r;
0302 }
0303
0304 inline dd_real dd_real::sloppy_div(const dd_real &a, const dd_real &b) {
0305   double s1, s2;
0306   double q1, q2;
0307   dd_real r;
0308
0309   q1 = a.x[0] / b.x[0];  /* approximate quotient */
0310
0311   /* compute  this - q1 * dd */
0312   r = b * q1;
0313   s1 = qd::two_diff(a.x[0], r.x[0], s2);
0314   s2 -= r.x[1];
0315   s2 += a.x[1];
0316
0317   /* get next approximation */
0318   q2 = (s1 + s2) / b.x[0];
0319
0320   /* renormalize */
0321   r.x[0] = qd::quick_two_sum(q1, q2, r.x[1]);
0322   return r;
0323 }
0324
0325 inline dd_real dd_real::accurate_div(const dd_real &a, const dd_real &b) {
0326   double q1, q2, q3;
0327   dd_real r;
0328
0329   q1 = a.x[0] / b.x[0];  /* approximate quotient */
0330
0331   r = a - q1 * b;
0332
0333   q2 = r.x[0] / b.x[0];
0334   r -= (q2 * b);
0335
0336   q3 = r.x[0] / b.x[0];
0337
0338   q1 = qd::quick_two_sum(q1, q2, q2);
0339   r = dd_real(q1, q2) + q3;
0340   return r;
0341 }
0342
0343 /* double-double / double-double */
0344 inline dd_real operator/(const dd_real &a, const dd_real &b) {
0345 #ifdef QD_SLOPPY_DIV
0346   return dd_real::sloppy_div(a, b);
0347 #else
0348   return dd_real::accurate_div(a, b);
0349 #endif
0350 }
0351
0352 /* double / double-double */
0353 inline dd_real operator/(double a, const dd_real &b) {
0354   return dd_real(a) / b;
0355 }
0356
0357 inline dd_real inv(const dd_real &a) {
0358   return 1.0 / a;
0359 }
0360
0361 /*********** Self-Divisions ************/
0362 /* double-double /= double */
0363 inline dd_real &dd_real::operator/=(double a) {
0364   *this = *this / a;
0365   return *this;
0366 }
0367
0368 /* double-double /= double-double */
0369 inline dd_real &dd_real::operator/=(const dd_real &a) {
0370   *this = *this / a;
0371   return *this;
0372 }
0373
0374 /********** Remainder **********/
0375 inline dd_real drem(const dd_real &a, const dd_real &b) {
0376   dd_real n = nint(a / b);
0377   return (a - n * b);
0378 }
0379
0380 inline dd_real divrem(const dd_real &a, const dd_real &b, dd_real &r) {
0381   dd_real n = nint(a / b);
0382   r = a - n * b;
0383   return n;
0384 }
0385
0386 /*********** Squaring **********/
0387 inline dd_real sqr(const dd_real &a) {
0388   double p1, p2;
0389   double s1, s2;
0390   p1 = qd::two_sqr(a.x[0], p2);
0391   p2 += 2.0 * a.x[0] * a.x[1];
0392   p2 += a.x[1] * a.x[1];
0393   s1 = qd::quick_two_sum(p1, p2, s2);
0394   return dd_real(s1, s2);
0395 }
0396
0397 inline dd_real dd_real::sqr(double a) {
0398   double p1, p2;
0399   p1 = qd::two_sqr(a, p2);
0400   return dd_real(p1, p2);
0401 }
0402
0403
0404 /********** Exponentiation **********/
0405 inline dd_real dd_real::operator^(int n) {
0406   return npwr(*this, n);
0407 }
0408
0409
0410 /*********** Assignments ************/
0411 /* double-double = double */
0412 inline dd_real &dd_real::operator=(double a) {
0413   x[0] = a;
0414   x[1] = 0.0;
0415   return *this;
0416 }
0417
0418 /*********** Equality Comparisons ************/
0419 /* double-double == double */
0420 inline bool operator==(const dd_real &a, double b) {
0421   return (a.x[0] == b && a.x[1] == 0.0);
0422 }
0423
0424 /* double-double == double-double */
0425 inline bool operator==(const dd_real &a, const dd_real &b) {
0426   return (a.x[0] == b.x[0] && a.x[1] == b.x[1]);
0427 }
0428
0429 /* double == double-double */
0430 inline bool operator==(double a, const dd_real &b) {
0431   return (a == b.x[0] && b.x[1] == 0.0);
0432 }
0433
0434 /*********** Greater-Than Comparisons ************/
0435 /* double-double > double */
0436 inline bool operator>(const dd_real &a, double b) {
0437   return (a.x[0] > b || (a.x[0] == b && a.x[1] > 0.0));
0438 }
0439
0440 /* double-double > double-double */
0441 inline bool operator>(const dd_real &a, const dd_real &b) {
0442   return (a.x[0] > b.x[0] || (a.x[0] == b.x[0] && a.x[1] > b.x[1]));
0443 }
0444
0445 /* double > double-double */
0446 inline bool operator>(double a, const dd_real &b) {
0447   return (a > b.x[0] || (a == b.x[0] && b.x[1] < 0.0));
0448 }
0449
0450 /*********** Less-Than Comparisons ************/
0451 /* double-double < double */
0452 inline bool operator<(const dd_real &a, double b) {
0453   return (a.x[0] < b || (a.x[0] == b && a.x[1] < 0.0));
0454 }
0455
0456 /* double-double < double-double */
0457 inline bool operator<(const dd_real &a, const dd_real &b) {
0458   return (a.x[0] < b.x[0] || (a.x[0] == b.x[0] && a.x[1] < b.x[1]));
0459 }
0460
0461 /* double < double-double */
0462 inline bool operator<(double a, const dd_real &b) {
0463   return (a < b.x[0] || (a == b.x[0] && b.x[1] > 0.0));
0464 }
0465
0466 /*********** Greater-Than-Or-Equal-To Comparisons ************/
0467 /* double-double >= double */
0468 inline bool operator>=(const dd_real &a, double b) {
0469   return (a.x[0] > b || (a.x[0] == b && a.x[1] >= 0.0));
0470 }
0471
0472 /* double-double >= double-double */
0473 inline bool operator>=(const dd_real &a, const dd_real &b) {
0474   return (a.x[0] > b.x[0] || (a.x[0] == b.x[0] && a.x[1] >= b.x[1]));
0475 }
0476
0477 /* double >= double-double */
0478 inline bool operator>=(double a, const dd_real &b) {
0479   return (b <= a);
0480 }
0481
0482 /*********** Less-Than-Or-Equal-To Comparisons ************/
0483 /* double-double <= double */
0484 inline bool operator<=(const dd_real &a, double b) {
0485   return (a.x[0] < b || (a.x[0] == b && a.x[1] <= 0.0));
0486 }
0487
0488 /* double-double <= double-double */
0489 inline bool operator<=(const dd_real &a, const dd_real &b) {
0490   return (a.x[0] < b.x[0] || (a.x[0] == b.x[0] && a.x[1] <= b.x[1]));
0491 }
0492
0493 /* double <= double-double */
0494 inline bool operator<=(double a, const dd_real &b) {
0495   return (b >= a);
0496 }
0497
0498 /*********** Not-Equal-To Comparisons ************/
0499 /* double-double != double */
0500 inline bool operator!=(const dd_real &a, double b) {
0501   return (a.x[0] != b || a.x[1] != 0.0);
0502 }
0503
0504 /* double-double != double-double */
0505 inline bool operator!=(const dd_real &a, const dd_real &b) {
0506   return (a.x[0] != b.x[0] || a.x[1] != b.x[1]);
0507 }
0508
0509 /* double != double-double */
0510 inline bool operator!=(double a, const dd_real &b) {
0511   return (a != b.x[0] || b.x[1] != 0.0);
0512 }
0513
0514 /*********** Micellaneous ************/
0515 /*  this == 0 */
0516 inline bool dd_real::is_zero() const {
0517   return (x[0] == 0.0);
0518 }
0519
0520 /*  this == 1 */
0521 inline bool dd_real::is_one() const {
0522   return (x[0] == 1.0 && x[1] == 0.0);
0523 }
0524
0525 /*  this > 0 */
0526 inline bool dd_real::is_positive() const {
0527   return (x[0] > 0.0);
0528 }
0529
0530 /* this < 0 */
0531 inline bool dd_real::is_negative() const {
0532   return (x[0] < 0.0);
0533 }
0534
0535 /* Absolute value */
0536 inline dd_real abs(const dd_real &a) {
0537   return (a.x[0] < 0.0) ? -a : a;
0538 }
0539
0540 inline dd_real fabs(const dd_real &a) {
0541   return abs(a);
0542 }
0543
0544 /* Round to Nearest integer */
0545 inline dd_real nint(const dd_real &a) {
0546   double hi = qd::nint(a.x[0]);
0547   double lo;
0548
0549   if (hi == a.x[0]) {
0550     /* High word is an integer already.  Round the low word.*/
0551     lo = qd::nint(a.x[1]);
0552
0553     /* Renormalize. This is needed if x[0] = some integer, x[1] = 1/2.*/
0554     hi = qd::quick_two_sum(hi, lo, lo);
0555   } else {
0556     /* High word is not an integer. */
0557     lo = 0.0;
0558     if (std::abs(hi-a.x[0]) == 0.5 && a.x[1] < 0.0) {
0559       /* There is a tie in the high word, consult the low word
0560          to break the tie. */
0561       hi -= 1.0;      /* NOTE: This does not cause INEXACT. */
0562     }
0563   }
0564
0565   return dd_real(hi, lo);
0566 }
0567
0568 inline dd_real floor(const dd_real &a) {
0569   double hi = std::floor(a.x[0]);
0570   double lo = 0.0;
0571
0572   if (hi == a.x[0]) {
0573     /* High word is integer already.  Round the low word. */
0574     lo = std::floor(a.x[1]);
0575     hi = qd::quick_two_sum(hi, lo, lo);
0576   }
0577
0578   return dd_real(hi, lo);
0579 }
0580
0581 inline dd_real ceil(const dd_real &a) {
0582   double hi = std::ceil(a.x[0]);
0583   double lo = 0.0;
0584
0585   if (hi == a.x[0]) {
0586     /* High word is integer already.  Round the low word. */
0587     lo = std::ceil(a.x[1]);
0588     hi = qd::quick_two_sum(hi, lo, lo);
0589   }
0590
0591   return dd_real(hi, lo);
0592 }
0593
0594 inline dd_real aint(const dd_real &a) {
0595   return (a.x[0] >= 0.0) ? floor(a) : ceil(a);
0596 }
0597
0598 /* Cast to double. */
0599 inline double to_double(const dd_real &a) {
0600   return a.x[0];
0601 }
0602
0603 /* Cast to int. */
0604 inline int to_int(const dd_real &a) {
0605   return static_cast<int>(a.x[0]);
0606 }
0607
0608 /* Random number generator */
0609 inline dd_real dd_real::rand() {
0610   return ddrand();
0611 }
0612
0613 #endif /* _QD_DD_INLINE_H */