include/qd/inline.h

0001 /*
0002  * include/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  * This file contains the basic functions used both by double-double
0011  * and quad-double package.  These are declared as inline functions as
0012  * they are the smallest building blocks of the double-double and
0013  * quad-double arithmetic.
0014  */
0015 #ifndef _QD_INLINE_H
0016 #define _QD_INLINE_H
0017
0018 #define _QD_SPLITTER 134217729.0               // = 2^27 + 1
0019 #define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996
0020
0021 #ifdef QD_VACPP_BUILTINS_H
0022 /* For VisualAge C++ __fmadd */
0023 #include <builtins.h>
0024 #endif
0025
0026 #include <cmath>
0027 #include <limits>
0028
0029 namespace qd {
0030
0031 static const double _d_nan = std::numeric_limits<double>::quiet_NaN();
0032 static const double _d_inf = std::numeric_limits<double>::infinity();
0033
0034 /*********** Basic Functions ************/
0035 /* Computes fl(a+b) and err(a+b).  Assumes |a| >= |b|. */
0036 inline double quick_two_sum(double a, double b, double &err) {
0037   double s = a + b;
0038   err = b - (s - a);
0039   return s;
0040 }
0041
0042 /* Computes fl(a-b) and err(a-b).  Assumes |a| >= |b| */
0043 inline double quick_two_diff(double a, double b, double &err) {
0044   double s = a - b;
0045   err = (a - s) - b;
0046   return s;
0047 }
0048
0049 /* Computes fl(a+b) and err(a+b).  */
0050 inline double two_sum(double a, double b, double &err) {
0051   double s = a + b;
0052   double bb = s - a;
0053   err = (a - (s - bb)) + (b - bb);
0054   return s;
0055 }
0056
0057 /* Computes fl(a-b) and err(a-b).  */
0058 inline double two_diff(double a, double b, double &err) {
0059   double s = a - b;
0060   double bb = s - a;
0061   err = (a - (s - bb)) - (b + bb);
0062   return s;
0063 }
0064
0065 #ifndef QD_FMS
0066 /* Computes high word and lo word of a */
0067 inline void split(double a, double &hi, double &lo) {
0068   double temp;
0069   if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {
0070     a *= 3.7252902984619140625e-09;  // 2^-28
0071     temp = _QD_SPLITTER * a;
0072     hi = temp - (temp - a);
0073     lo = a - hi;
0074     hi *= 268435456.0;          // 2^28
0075     lo *= 268435456.0;          // 2^28
0076   } else {
0077     temp = _QD_SPLITTER * a;
0078     hi = temp - (temp - a);
0079     lo = a - hi;
0080   }
0081 }
0082 #endif
0083
0084 /* Computes fl(a*b) and err(a*b). */
0085 inline double two_prod(double a, double b, double &err) {
0086 #ifdef QD_FMS
0087   double p = a * b;
0088   err = QD_FMS(a, b, p);
0089   return p;
0090 #else
0091   double a_hi, a_lo, b_hi, b_lo;
0092   double p = a * b;
0093   split(a, a_hi, a_lo);
0094   split(b, b_hi, b_lo);
0095   err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
0096   return p;
0097 #endif
0098 }
0099
0100 /* Computes fl(a*a) and err(a*a).  Faster than the above method. */
0101 inline double two_sqr(double a, double &err) {
0102 #ifdef QD_FMS
0103   double p = a * a;
0104   err = QD_FMS(a, a, p);
0105   return p;
0106 #else
0107   double hi, lo;
0108   double q = a * a;
0109   split(a, hi, lo);
0110   err = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;
0111   return q;
0112 #endif
0113 }
0114
0115 /* Computes the nearest integer to d. */
0116 inline double nint(double d) {
0117   if (d == std::floor(d))
0118     return d;
0119   return std::floor(d + 0.5);
0120 }
0121
0122 /* Computes the truncated integer. */
0123 inline double aint(double d) {
0124   return (d >= 0.0) ? std::floor(d) : std::ceil(d);
0125 }
0126
0127 /* These are provided to give consistent
0128    interface for double with double-double and quad-double. */
0129 inline void sincosh(double t, double &sinh_t, double &cosh_t) {
0130   sinh_t = std::sinh(t);
0131   cosh_t = std::cosh(t);
0132 }
0133
0134 inline double sqr(double t) {
0135   return t * t;
0136 }
0137
0138 inline double to_double(double a) { return a; }
0139 inline int    to_int(double a) { return static_cast<int>(a); }
0140
0141 }
0142
0143 #endif /* _QD_INLINE_H */