dd_real_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_REAL_IDEFS_H_
22 #define _DD_REAL_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 #include "dd_idefs.h"
33 
34 /*
35  ************************************************************************
36  Now for the double2 routines
37  ************************************************************************
38 */
39 
40 static inline double
41 dd_hi(const double2 a)
42 {
43  return a.x[0];
44 }
45 
46 static inline double
47 dd_lo(const double2 a)
48 {
49  return a.x[1];
50 }
51 
52 static inline int
54 {
55  return isfinite(a.x[0]);
56 }
57 
58 static inline int
60 {
61  return isinf(a.x[0]);
62 }
63 
64 static inline int
66 {
67  return (a.x[0] == 0.0);
68 }
69 
70 static inline int
72 {
73  return (a.x[0] == 1.0 && a.x[1] == 0.0);
74 }
75 
76 static inline int
78 {
79  return (a.x[0] > 0.0);
80 }
81 
82 static inline int
84 {
85  return (a.x[0] < 0.0);
86 }
87 
88 /* Cast to double. */
89 static inline double
91 {
92  return a.x[0];
93 }
94 
95 /* Cast to int. */
96 static inline int
98 {
99  return DD_STATIC_CAST(int, a.x[0]);
100 }
101 
102 /*********** Equality and Other Comparisons ************/
103 static inline int
104 dd_comp(const double2 a, const double2 b)
105 {
106  int cmp = two_comp(a.x[0], b.x[0]);
107  if (cmp == 0) {
108  cmp = two_comp(a.x[1], b.x[1]);
109  }
110  return cmp;
111 }
112 
113 static inline int
114 dd_comp_dd_d(const double2 a, double b)
115 {
116  int cmp = two_comp(a.x[0], b);
117  if (cmp == 0) {
118  cmp = two_comp(a.x[1], 0);
119  }
120  return cmp;
121 }
122 
123 static inline int
124 dd_comp_d_dd(double a, const double2 b)
125 {
126  int cmp = two_comp(a, b.x[0]);
127  if (cmp == 0) {
128  cmp = two_comp(0.0, b.x[1]);
129  }
130  return cmp;
131 }
132 
133 
134 /*********** Creation ************/
135 static inline double2
136 dd_create(double hi, double lo)
137 {
138  double2 ret = {{hi, lo}};
139  return ret;
140 }
141 
142 static inline double2
143 dd_zero(void)
144 {
145  return DD_C_ZERO;
146 }
147 
148 static inline double2
149 dd_create_d(double hi)
150 {
151  double2 ret = {{hi, 0.0}};
152  return ret;
153 }
154 
155 static inline double2
156 dd_create_i(int hi)
157 {
158  double2 ret = {{DD_STATIC_CAST(double, hi), 0.0}};
159  return ret;
160 }
161 
162 static inline double2
163 dd_create_dp(const double *d)
164 {
165  double2 ret = {{d[0], d[1]}};
166  return ret;
167 }
168 
169 
170 /*********** Unary Minus ***********/
171 static inline double2
173 {
174  double2 ret = {{-a.x[0], -a.x[1]}};
175  return ret;
176 }
177 
178 /*********** Rounding ************/
179 /* Round to Nearest integer */
180 static inline double2
182 {
183  double hi = two_nint(a.x[0]);
184  double lo;
185 
186  if (hi == a.x[0]) {
187  /* High word is an integer already. Round the low word.*/
188  lo = two_nint(a.x[1]);
189 
190  /* Renormalize. This is needed if x[0] = some integer, x[1] = 1/2.*/
191  hi = quick_two_sum(hi, lo, &lo);
192  }
193  else {
194  /* High word is not an integer. */
195  lo = 0.0;
196  if (fabs(hi - a.x[0]) == 0.5 && a.x[1] < 0.0) {
197  /* There is a tie in the high word, consult the low word
198  to break the tie. */
199  hi -= 1.0; /* NOTE: This does not cause INEXACT. */
200  }
201  }
202 
203  return dd_create(hi, lo);
204 }
205 
206 static inline double2
208 {
209  double hi = floor(a.x[0]);
210  double lo = 0.0;
211 
212  if (hi == a.x[0]) {
213  /* High word is integer already. Round the low word. */
214  lo = floor(a.x[1]);
215  hi = quick_two_sum(hi, lo, &lo);
216  }
217 
218  return dd_create(hi, lo);
219 }
220 
221 static inline double2
223 {
224  double hi = ceil(a.x[0]);
225  double lo = 0.0;
226 
227  if (hi == a.x[0]) {
228  /* High word is integer already. Round the low word. */
229  lo = ceil(a.x[1]);
230  hi = quick_two_sum(hi, lo, &lo);
231  }
232 
233  return dd_create(hi, lo);
234 }
235 
236 static inline double2
238 {
239  return (a.x[0] >= 0.0) ? dd_floor(a) : dd_ceil(a);
240 }
241 
242 /* Absolute value */
243 static inline double2
245 {
246  return (a.x[0] < 0.0 ? dd_neg(a) : a);
247 }
248 
249 static inline double2
251 {
252  return dd_abs(a);
253 }
254 
255 
256 /*********** Normalizing ***********/
257 /* double-double * (2.0 ^ expt) */
258 static inline double2
259 dd_ldexp(const double2 a, int expt)
260 {
261  return dd_create(ldexp(a.x[0], expt), ldexp(a.x[1], expt));
262 }
263 
264 static inline double2
265 dd_frexp(const double2 a, int *expt)
266 {
267 // r"""return b and l s.t. 0.5<=|b|<1 and 2^l == a
268 // 0.5<=|b[0]|<1.0 or |b[0]| == 1.0 and b[0]*b[1]<0
269 // """
270  int exponent;
271  double man = frexp(a.x[0], &exponent);
272  double b1 = ldexp(a.x[1], -exponent);
273  if (fabs(man) == 0.5 && man * b1 < 0)
274  {
275  man *=2;
276  b1 *= 2;
277  exponent -= 1;
278  }
279  *expt = exponent;
280  return dd_create(man, b1);
281 }
282 
283 
284 /*********** Additions ************/
285 static inline double2
286 dd_add_d_d(double a, double b)
287 {
288  double s, e;
289  s = two_sum(a, b, &e);
290  return dd_create(s, e);
291 }
292 
293 static inline double2
294 dd_add_dd_d(const double2 a, double b)
295 {
296  double s1, s2;
297  s1 = two_sum(a.x[0], b, &s2);
298  s2 += a.x[1];
299  s1 = quick_two_sum(s1, s2, &s2);
300  return dd_create(s1, s2);
301 }
302 
303 static inline double2
304 dd_add_d_dd(double a, const double2 b)
305 {
306  double s1, s2;
307  s1 = two_sum(a, b.x[0], &s2);
308  s2 += b.x[1];
309  s1 = quick_two_sum(s1, s2, &s2);
310  return dd_create(s1, s2);
311 }
312 
313 static inline double2
315 {
316  /* This one satisfies IEEE style error bound,
317  due to K. Briggs and W. Kahan. */
318  double s1, s2, t1, t2;
319 
320  s1 = two_sum(a.x[0], b.x[0], &s2);
321  t1 = two_sum(a.x[1], b.x[1], &t2);
322  s2 += t1;
323  s1 = quick_two_sum(s1, s2, &s2);
324  s2 += t2;
325  s1 = quick_two_sum(s1, s2, &s2);
326  return dd_create(s1, s2);
327 }
328 
329 static inline double2
331 {
332  /* This is the less accurate version ... obeys Cray-style
333  error bound. */
334  double s, e;
335 
336  s = two_sum(a.x[0], b.x[0], &e);
337  e += (a.x[1] + b.x[1]);
338  s = quick_two_sum(s, e, &e);
339  return dd_create(s, e);
340 }
341 
342 static inline double2
343 dd_add(const double2 a, const double2 b)
344 {
345  /* Always require IEEE-style error bounds */
346  return dd_ieee_add(a, b);
347 }
348 
349 /*********** Subtractions ************/
350 /* double-double = double - double */
351 static inline double2
352 dd_sub_d_d(double a, double b)
353 {
354  double s, e;
355  s = two_diff(a, b, &e);
356  return dd_create(s, e);
357 }
358 
359 static inline double2
360 dd_sub(const double2 a, const double2 b)
361 {
362  return dd_ieee_add(a, dd_neg(b));
363 }
364 
365 static inline double2
366 dd_sub_dd_d(const double2 a, double b)
367 {
368  double s1, s2;
369  s1 = two_sum(a.x[0], -b, &s2);
370  s2 += a.x[1];
371  s1 = quick_two_sum(s1, s2, &s2);
372  return dd_create(s1, s2);
373 }
374 
375 static inline double2
376 dd_sub_d_dd(double a, const double2 b)
377 {
378  double s1, s2;
379  s1 = two_sum(a, -b.x[0], &s2);
380  s2 -= b.x[1];
381  s1 = quick_two_sum(s1, s2, &s2);
382  return dd_create(s1, s2);
383 }
384 
385 
386 /*********** Multiplications ************/
387 /* double-double = double * double */
388 static inline double2
389 dd_mul_d_d(double a, double b)
390 {
391  double p, e;
392  p = two_prod(a, b, &e);
393  return dd_create(p, e);
394 }
395 
396 /* double-double * double, where double is a power of 2. */
397 static inline double2
398 dd_mul_pwr2(const double2 a, double b)
399 {
400  return dd_create(a.x[0] * b, a.x[1] * b);
401 }
402 
403 static inline double2
404 dd_mul(const double2 a, const double2 b)
405 {
406  double p1, p2;
407  p1 = two_prod(a.x[0], b.x[0], &p2);
408  p2 += (a.x[0] * b.x[1] + a.x[1] * b.x[0]);
409  p1 = quick_two_sum(p1, p2, &p2);
410  return dd_create(p1, p2);
411 }
412 
413 static inline double2
414 dd_mul_dd_d(const double2 a, double b)
415 {
416  double p1, p2, e1, e2;
417  p1 = two_prod(a.x[0], b, &e1);
418  p2 = two_prod(a.x[1], b, &e2);
419  p1 = quick_two_sum(p1, e2 + p2 + e1, &e1);
420  return dd_create(p1, e1);
421 }
422 
423 static inline double2
424 dd_mul_d_dd(double a, const double2 b)
425 {
426  double p1, p2, e1, e2;
427  p1 = two_prod(a, b.x[0], &e1);
428  p2 = two_prod(a, b.x[1], &e2);
429  p1 = quick_two_sum(p1, e2 + p2 + e1, &e1);
430  return dd_create(p1, e1);
431 }
432 
433 
434 /*********** Divisions ************/
435 static inline double2
437 {
438  double s1, s2;
439  double q1, q2;
440  double2 r;
441 
442  q1 = a.x[0] / b.x[0]; /* approximate quotient */
443 
444  /* compute this - q1 * dd */
445  r = dd_sub(a, dd_mul_dd_d(b, q1));
446  s1 = two_diff(a.x[0], r.x[0], &s2);
447  s2 -= r.x[1];
448  s2 += a.x[1];
449 
450  /* get next approximation */
451  q2 = (s1 + s2) / b.x[0];
452 
453  /* renormalize */
454  r.x[0] = quick_two_sum(q1, q2, &r.x[1]);
455  return r;
456 }
457 
458 static inline double2
460 {
461  double q1, q2, q3;
462  double2 r;
463 
464  q1 = a.x[0] / b.x[0]; /* approximate quotient */
465 
466  r = dd_sub(a, dd_mul_dd_d(b, q1));
467 
468  q2 = r.x[0] / b.x[0];
469  r = dd_sub(r, dd_mul_dd_d(b, q2));
470 
471  q3 = r.x[0] / b.x[0];
472 
473  q1 = quick_two_sum(q1, q2, &q2);
474  r = dd_add_dd_d(dd_create(q1, q2), q3);
475  return r;
476 }
477 
478 static inline double2
479 dd_div(const double2 a, const double2 b)
480 {
481  return dd_accurate_div(a, b);
482 }
483 
484 static inline double2
485 dd_div_d_d(double a, double b)
486 {
488 }
489 
490 static inline double2
491 dd_div_dd_d(const double2 a, double b)
492 {
493  return dd_accurate_div(a, dd_create_d(b));
494 }
495 
496 static inline double2
497 dd_div_d_dd(double a, const double2 b)
498 {
499  return dd_accurate_div(dd_create_d(a), b);
500 }
501 
502 static inline double2
504 {
505  return dd_div(DD_C_ONE, a);
506 }
507 
508 
509 /********** Remainder **********/
510 static inline double2
511 dd_drem(const double2 a, const double2 b)
512 {
513  double2 n = dd_nint(dd_div(a, b));
514  return dd_sub(a, dd_mul(n, b));
515 }
516 
517 static inline double2
518 dd_divrem(const double2 a, const double2 b, double2 *r)
519 {
520  double2 n = dd_nint(dd_div(a, b));
521  *r = dd_sub(a, dd_mul(n, b));
522  return n;
523 }
524 
525 static inline double2
526 dd_fmod(const double2 a, const double2 b)
527 {
528  double2 n = dd_aint(dd_div(a, b));
529  return dd_sub(a, dd_mul(b, n));
530 }
531 
532 /*********** Squaring **********/
533 static inline double2
535 {
536  double p1, p2;
537  double s1, s2;
538  p1 = two_sqr(a.x[0], &p2);
539  p2 += 2.0 * a.x[0] * a.x[1];
540  p2 += a.x[1] * a.x[1];
541  s1 = quick_two_sum(p1, p2, &s2);
542  return dd_create(s1, s2);
543 }
544 
545 static inline double2
546 dd_sqr_d(double a)
547 {
548  double p1, p2;
549  p1 = two_sqr(a, &p2);
550  return dd_create(p1, p2);
551 }
552 
553 #ifdef __cplusplus
554 }
555 #endif
556 
557 #endif /* _DD_REAL_IDEFS_H_ */
dd_floor
static double2 dd_floor(const double2 a)
Definition: dd_real_idefs.h:207
dd_mul_pwr2
static double2 dd_mul_pwr2(const double2 a, double b)
Definition: dd_real_idefs.h:398
dd_to_int
static int dd_to_int(const double2 a)
Definition: dd_real_idefs.h:97
DD_C_ONE
const double2 DD_C_ONE
Definition: dd_real.c:47
dd_add_d_dd
static double2 dd_add_d_dd(double a, const double2 b)
Definition: dd_real_idefs.h:304
two_comp
static int two_comp(const double a, const double b)
Definition: dd_idefs.h:187
simple_graph::b1
Vector2 b1(2, -1)
dd_sloppy_add
static double2 dd_sloppy_add(const double2 a, const double2 b)
Definition: dd_real_idefs.h:330
double2
Definition: dd_real.h:74
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.)
dd_zero
static double2 dd_zero(void)
Definition: dd_real_idefs.h:143
d
static const double d[K][N]
Definition: igam.h:11
dd_divrem
static double2 dd_divrem(const double2 a, const double2 b, double2 *r)
Definition: dd_real_idefs.h:518
dd_comp_d_dd
static int dd_comp_d_dd(double a, const double2 b)
Definition: dd_real_idefs.h:124
dd_mul
static double2 dd_mul(const double2 a, const double2 b)
Definition: dd_real_idefs.h:404
b
Scalar * b
Definition: benchVecAdd.cpp:17
dd_neg
static double2 dd_neg(const double2 a)
Definition: dd_real_idefs.h:172
dd_div_d_d
static double2 dd_div_d_d(double a, double b)
Definition: dd_real_idefs.h:485
dd_drem
static double2 dd_drem(const double2 a, const double2 b)
Definition: dd_real_idefs.h:511
ret
DenseIndex ret
Definition: level1_cplx_impl.h:44
dd_idefs.h
dd_comp
static int dd_comp(const double2 a, const double2 b)
Definition: dd_real_idefs.h:104
dd_fabs
static double2 dd_fabs(const double2 a)
Definition: dd_real_idefs.h:250
dd_fmod
static double2 dd_fmod(const double2 a, const double2 b)
Definition: dd_real_idefs.h:526
dd_div
static double2 dd_div(const double2 a, const double2 b)
Definition: dd_real_idefs.h:479
dd_inv
static double2 dd_inv(const double2 a)
Definition: dd_real_idefs.h:503
dd_mul_d_dd
static double2 dd_mul_d_dd(double a, const double2 b)
Definition: dd_real_idefs.h:424
DD_C_ZERO
const double2 DD_C_ZERO
Definition: dd_real.c:46
dd_add_d_d
static double2 dd_add_d_d(double a, double b)
Definition: dd_real_idefs.h:286
DD_STATIC_CAST
#define DD_STATIC_CAST(T, X)
Definition: dd_real.h:69
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
n
int n
Definition: BiCGSTAB_simple.cpp:1
simple::p2
static Point3 p2
Definition: testInitializePose3.cpp:51
double2::x
double x[2]
Definition: dd_real.h:76
dd_ceil
static double2 dd_ceil(const double2 a)
Definition: dd_real_idefs.h:222
dd_accurate_div
static double2 dd_accurate_div(const double2 a, const double2 b)
Definition: dd_real_idefs.h:459
dd_ldexp
static double2 dd_ldexp(const double2 a, int expt)
Definition: dd_real_idefs.h:259
dd_div_dd_d
static double2 dd_div_dd_d(const double2 a, double b)
Definition: dd_real_idefs.h:491
dd_sqr_d
static double2 dd_sqr_d(double a)
Definition: dd_real_idefs.h:546
dd_to_double
static double dd_to_double(const double2 a)
Definition: dd_real_idefs.h:90
isfinite
#define isfinite(X)
Definition: main.h:95
dd_frexp
static double2 dd_frexp(const double2 a, int *expt)
Definition: dd_real_idefs.h:265
dd_hi
static double dd_hi(const double2 a)
Definition: dd_real_idefs.h:41
dd_create_dp
static double2 dd_create_dp(const double *d)
Definition: dd_real_idefs.h:163
dd_isfinite
static int dd_isfinite(const double2 a)
Definition: dd_real_idefs.h:53
dd_div_d_dd
static double2 dd_div_d_dd(double a, const double2 b)
Definition: dd_real_idefs.h:497
dd_is_zero
static int dd_is_zero(const double2 a)
Definition: dd_real_idefs.h:65
dd_sqr
static double2 dd_sqr(const double2 a)
Definition: dd_real_idefs.h:534
dd_nint
static double2 dd_nint(const double2 a)
Definition: dd_real_idefs.h:181
dd_add
static double2 dd_add(const double2 a, const double2 b)
Definition: dd_real_idefs.h:343
two_sum
static double two_sum(double a, double b, double *err)
Definition: dd_idefs.h:63
dd_isinf
static int dd_isinf(const double2 a)
Definition: dd_real_idefs.h:59
dd_lo
static double dd_lo(const double2 a)
Definition: dd_real_idefs.h:47
dd_abs
static double2 dd_abs(const double2 a)
Definition: dd_real_idefs.h:244
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
dd_sloppy_div
static double2 dd_sloppy_div(const double2 a, const double2 b)
Definition: dd_real_idefs.h:436
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
two_prod
static double two_prod(double a, double b, double *err)
Definition: dd_idefs.h:109
dd_sub
static double2 dd_sub(const double2 a, const double2 b)
Definition: dd_real_idefs.h:360
dd_is_positive
static int dd_is_positive(const double2 a)
Definition: dd_real_idefs.h:77
dd_is_one
static int dd_is_one(const double2 a)
Definition: dd_real_idefs.h:71
dd_mul_d_d
static double2 dd_mul_d_d(double a, double b)
Definition: dd_real_idefs.h:389
p
float * p
Definition: Tutorial_Map_using.cpp:9
dd_sub_dd_d
static double2 dd_sub_dd_d(const double2 a, double b)
Definition: dd_real_idefs.h:366
dd_sub_d_dd
static double2 dd_sub_d_dd(double a, const double2 b)
Definition: dd_real_idefs.h:376
dd_create_d
static double2 dd_create_d(double hi)
Definition: dd_real_idefs.h:149
dd_mul_dd_d
static double2 dd_mul_dd_d(const double2 a, double b)
Definition: dd_real_idefs.h:414
ceil
const EIGEN_DEVICE_FUNC CeilReturnType ceil() const
Definition: ArrayCwiseUnaryOps.h:495
dd_create_i
static double2 dd_create_i(int hi)
Definition: dd_real_idefs.h:156
dd_ieee_add
static double2 dd_ieee_add(const double2 a, const double2 b)
Definition: dd_real_idefs.h:314
two_diff
static double two_diff(double a, double b, double *err)
Definition: dd_idefs.h:75
dd_aint
static double2 dd_aint(const double2 a)
Definition: dd_real_idefs.h:237
dd_create
static double2 dd_create(double hi, double lo)
Definition: dd_real_idefs.h:136
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
dd_add_dd_d
static double2 dd_add_dd_d(const double2 a, double b)
Definition: dd_real_idefs.h:294
dd_sub_d_d
static double2 dd_sub_d_d(double a, double b)
Definition: dd_real_idefs.h:352
two_sqr
static double two_sqr(double a, double *err)
Definition: dd_idefs.h:130
dd_comp_dd_d
static int dd_comp_dd_d(const double2 a, double b)
Definition: dd_real_idefs.h:114
dd_is_negative
static int dd_is_negative(const double2 a)
Definition: dd_real_idefs.h:83
isinf
#define isinf(X)
Definition: main.h:94


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