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 #include <numpy/utils.h>
29
30 #ifdef __cplusplus
31 extern "C" {
32 #endif
33
34 #define _DD_SPLITTER 134217729.0 // = 2^27 + 1
35 #define _DD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996
36
37 /*
38 ************************************************************************
39 The basic routines taking double arguments, returning 1 (or 2) doubles
40 ************************************************************************
41 */
42
43 /* Computes fl(a+b) and err(a+b). Assumes |a| >= |b|. */
44 static NPY_INLINE double
quick_two_sum(double a,double b,double * err)45 quick_two_sum(double a, double b, double *err)
46 {
47 volatile double s = a + b;
48 volatile double c = s - a;
49 *err = b - c;
50 return s;
51 }
52
53 /* Computes fl(a-b) and err(a-b). Assumes |a| >= |b| */
54 static NPY_INLINE double
quick_two_diff(double a,double b,double * err)55 quick_two_diff(double a, double b, double *err)
56 {
57 volatile double s = a - b;
58 volatile double c = a - s;
59 *err = c - b;
60 return s;
61 }
62
63 /* Computes fl(a+b) and err(a+b). */
64 static NPY_INLINE double
two_sum(double a,double b,double * err)65 two_sum(double a, double b, double *err)
66 {
67 volatile double s = a + b;
68 volatile double c = s - a;
69 volatile double d = b - c;
70 volatile double e = s - c;
71 *err = (a - e) + d;
72 return s;
73 }
74
75 /* Computes fl(a-b) and err(a-b). */
76 static NPY_INLINE double
two_diff(double a,double b,double * err)77 two_diff(double a, double b, double *err)
78 {
79 volatile double s = a - b;
80 volatile double c = s - a;
81 volatile double d = b + c;
82 volatile double e = s - c;
83 *err = (a - e) - d;
84 return s;
85 }
86
87 /* Computes high word and lo word of a */
88 static NPY_INLINE void
two_split(double a,double * hi,double * lo)89 two_split(double a, double *hi, double *lo)
90 {
91 volatile double temp, tempma;
92 if (a > _DD_SPLIT_THRESH || a < -_DD_SPLIT_THRESH) {
93 a *= 3.7252902984619140625e-09; // 2^-28
94 temp = _DD_SPLITTER * a;
95 tempma = temp - a;
96 *hi = temp - tempma;
97 *lo = a - *hi;
98 *hi *= 268435456.0; // 2^28
99 *lo *= 268435456.0; // 2^28
100 }
101 else {
102 temp = _DD_SPLITTER * a;
103 tempma = temp - a;
104 *hi = temp - tempma;
105 *lo = a - *hi;
106 }
107 }
108
109 /* Computes fl(a*b) and err(a*b). */
110 static NPY_INLINE double
two_prod(double a,double b,double * err)111 two_prod(double a, double b, double *err)
112 {
113 #ifdef DD_FMS
114 volatile double p = a * b;
115 *err = DD_FMS(a, b, p);
116 return p;
117 #else
118 double a_hi, a_lo, b_hi, b_lo;
119 double p = a * b;
120 volatile double c, d;
121 two_split(a, &a_hi, &a_lo);
122 two_split(b, &b_hi, &b_lo);
123 c = a_hi * b_hi - p;
124 d = c + a_hi * b_lo + a_lo * b_hi;
125 *err = d + a_lo * b_lo;
126 return p;
127 #endif /* DD_FMA */
128 }
129
130 /* Computes fl(a*a) and err(a*a). Faster than the above method. */
131 static NPY_INLINE double
two_sqr(double a,double * err)132 two_sqr(double a, double *err)
133 {
134 #ifdef DD_FMS
135 volatile double p = a * a;
136 *err = DD_FMS(a, a, p);
137 return p;
138 #else
139 double hi, lo;
140 volatile double c;
141 double q = a * a;
142 two_split(a, &hi, &lo);
143 c = hi * hi - q;
144 *err = (c + 2.0 * hi * lo) + lo * lo;
145 return q;
146 #endif /* DD_FMS */
147 }
148
149 static NPY_INLINE double
two_div(double a,double b,double * err)150 two_div(double a, double b, double *err)
151 {
152 volatile double q1, q2;
153 double p1, p2;
154 double s, e;
155
156 q1 = a / b;
157
158 /* Compute a - q1 * b */
159 p1 = two_prod(q1, b, &p2);
160 s = two_diff(a, p1, &e);
161 e -= p2;
162
163 /* get next approximation */
164 q2 = (s + e) / b;
165
166 return quick_two_sum(q1, q2, err);
167 }
168
169 /* Computes the nearest integer to d. */
170 static NPY_INLINE double
two_nint(double d)171 two_nint(double d)
172 {
173 if (d == floor(d)) {
174 return d;
175 }
176 return floor(d + 0.5);
177 }
178
179 /* Computes the truncated integer. */
180 static NPY_INLINE double
two_aint(double d)181 two_aint(double d)
182 {
183 return (d >= 0.0 ? floor(d) : ceil(d));
184 }
185
186
187 /* Compare a and b */
188 static NPY_INLINE int
two_comp(const double a,const double b)189 two_comp(const double a, const double b)
190 {
191 /* Works for non-NAN inputs */
192 return (a < b ? -1 : (a > b ? 1 : 0));
193 }
194
195
196 #ifdef __cplusplus
197 }
198 #endif
199
200 #endif /* _DD_IDEFS_H_ */
201