dd_idefs.h
Go to the documentation of this file.
1 /*
2  * include/dd_inline.h
3  *
4  * This work was supported by the Director, Office of Science, Division
5  * of Mathematical, Information, and Computational Sciences of the
6  * U.S. Department of Energy under contract numbers DE-AC03-76SF00098 and
7  * DE-AC02-05CH11231.
8  *
9  * Copyright (c) 2003-2009, The Regents of the University of California,
10  * through Lawrence Berkeley National Laboratory (subject to receipt of
11  * any required approvals from U.S. Dept. of Energy) All rights reserved.
12  *
13  * By downloading or using this software you are agreeing to the modified
14  * BSD license "BSD-LBNL-License.doc" (see LICENSE.txt).
15  */
16 /*
17  * Contains small functions (suitable for inlining) in the double-double
18  * arithmetic package.
19  */
20 
21 #ifndef _DD_IDEFS_H_
22 #define _DD_IDEFS_H_ 1
23 
24 #include <float.h>
25 #include <limits.h>
26 #include <math.h>
27 
28 #ifdef __cplusplus
29 extern "C" {
30 #endif
31 
32 #define _DD_SPLITTER 134217729.0 // = 2^27 + 1
33 #define _DD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996
34 
35 /*
36  ************************************************************************
37  The basic routines taking double arguments, returning 1 (or 2) doubles
38  ************************************************************************
39 */
40 
41 /* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */
42 static inline double
43 quick_two_sum(double a, double b, double *err)
44 {
45  volatile double s = a + b;
46  volatile double c = s - a;
47  *err = b - c;
48  return s;
49 }
50 
51 /* Computes fl(a-b) and err(a-b). Assumes |a| >= |b| */
52 static inline double
53 quick_two_diff(double a, double b, double *err)
54 {
55  volatile double s = a - b;
56  volatile double c = a - s;
57  *err = c - b;
58  return s;
59 }
60 
61 /* Computes fl(a+b) and err(a+b). */
62 static inline double
63 two_sum(double a, double b, double *err)
64 {
65  volatile double s = a + b;
66  volatile double c = s - a;
67  volatile double d = b - c;
68  volatile double e = s - c;
69  *err = (a - e) + d;
70  return s;
71 }
72 
73 /* Computes fl(a-b) and err(a-b). */
74 static inline double
75 two_diff(double a, double b, double *err)
76 {
77  volatile double s = a - b;
78  volatile double c = s - a;
79  volatile double d = b + c;
80  volatile double e = s - c;
81  *err = (a - e) - d;
82  return s;
83 }
84 
85 /* Computes high word and lo word of a */
86 static inline void
87 two_split(double a, double *hi, double *lo)
88 {
89  volatile double temp, tempma;
90  if (a > _DD_SPLIT_THRESH || a < -_DD_SPLIT_THRESH) {
91  a *= 3.7252902984619140625e-09; // 2^-28
92  temp = _DD_SPLITTER * a;
93  tempma = temp - a;
94  *hi = temp - tempma;
95  *lo = a - *hi;
96  *hi *= 268435456.0; // 2^28
97  *lo *= 268435456.0; // 2^28
98  }
99  else {
100  temp = _DD_SPLITTER * a;
101  tempma = temp - a;
102  *hi = temp - tempma;
103  *lo = a - *hi;
104  }
105 }
106 
107 /* Computes fl(a*b) and err(a*b). */
108 static inline double
109 two_prod(double a, double b, double *err)
110 {
111 #ifdef DD_FMS
112  volatile double p = a * b;
113  *err = DD_FMS(a, b, p);
114  return p;
115 #else
116  double a_hi, a_lo, b_hi, b_lo;
117  double p = a * b;
118  volatile double c, d;
119  two_split(a, &a_hi, &a_lo);
120  two_split(b, &b_hi, &b_lo);
121  c = a_hi * b_hi - p;
122  d = c + a_hi * b_lo + a_lo * b_hi;
123  *err = d + a_lo * b_lo;
124  return p;
125 #endif /* DD_FMA */
126 }
127 
128 /* Computes fl(a*a) and err(a*a). Faster than the above method. */
129 static inline double
130 two_sqr(double a, double *err)
131 {
132 #ifdef DD_FMS
133  volatile double p = a * a;
134  *err = DD_FMS(a, a, p);
135  return p;
136 #else
137  double hi, lo;
138  volatile double c;
139  double q = a * a;
140  two_split(a, &hi, &lo);
141  c = hi * hi - q;
142  *err = (c + 2.0 * hi * lo) + lo * lo;
143  return q;
144 #endif /* DD_FMS */
145 }
146 
147 static inline double
148 two_div(double a, double b, double *err)
149 {
150  volatile double q1, q2;
151  double p1, p2;
152  double s, e;
153 
154  q1 = a / b;
155 
156  /* Compute a - q1 * b */
157  p1 = two_prod(q1, b, &p2);
158  s = two_diff(a, p1, &e);
159  e -= p2;
160 
161  /* get next approximation */
162  q2 = (s + e) / b;
163 
164  return quick_two_sum(q1, q2, err);
165 }
166 
167 /* Computes the nearest integer to d. */
168 static inline double
169 two_nint(double d)
170 {
171  if (d == floor(d)) {
172  return d;
173  }
174  return floor(d + 0.5);
175 }
176 
177 /* Computes the truncated integer. */
178 static inline double
179 two_aint(double d)
180 {
181  return (d >= 0.0 ? floor(d) : ceil(d));
182 }
183 
184 
185 /* Compare a and b */
186 static inline int
187 two_comp(const double a, const double b)
188 {
189  /* Works for non-NAN inputs */
190  return (a < b ? -1 : (a > b ? 1 : 0));
191 }
192 
193 
194 #ifdef __cplusplus
195 }
196 #endif
197 
198 #endif /* _DD_IDEFS_H_ */
two_comp
static int two_comp(const double a, const double b)
Definition: dd_idefs.h:187
two_nint
static double two_nint(double d)
Definition: dd_idefs.h:169
s
RealScalar s
Definition: level1_cplx_impl.h:126
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
d
static const double d[K][N]
Definition: igam.h:11
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
two_split
static void two_split(double a, double *hi, double *lo)
Definition: dd_idefs.h:87
simple::p2
static Point3 p2
Definition: testInitializePose3.cpp:51
quick_two_diff
static double quick_two_diff(double a, double b, double *err)
Definition: dd_idefs.h:53
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
two_sum
static double two_sum(double a, double b, double *err)
Definition: dd_idefs.h:63
p1
Vector3f p1
Definition: MatrixBase_all.cpp:2
quick_two_sum
static double quick_two_sum(double a, double b, double *err)
Definition: dd_idefs.h:43
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
two_div
static double two_div(double a, double b, double *err)
Definition: dd_idefs.h:148
two_prod
static double two_prod(double a, double b, double *err)
Definition: dd_idefs.h:109
two_aint
static double two_aint(double d)
Definition: dd_idefs.h:179
p
float * p
Definition: Tutorial_Map_using.cpp:9
_DD_SPLITTER
#define _DD_SPLITTER
Definition: dd_idefs.h:32
_DD_SPLIT_THRESH
#define _DD_SPLIT_THRESH
Definition: dd_idefs.h:33
ceil
const EIGEN_DEVICE_FUNC CeilReturnType ceil() const
Definition: ArrayCwiseUnaryOps.h:495
two_diff
static double two_diff(double a, double b, double *err)
Definition: dd_idefs.h:75
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
two_sqr
static double two_sqr(double a, double *err)
Definition: dd_idefs.h:130


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:10