Automatic Differentiation
 
Loading...
Searching...
No Matches
double_d.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_OPENCL_DOUBLE_D_HPP
2#define STAN_MATH_OPENCL_DOUBLE_D_HPP
3
5
6#include <string>
7#include <cmath>
8#include <cfloat>
9#include <limits>
10
14#define STRINGIFY2(A) \
15#A; \
16 A
17namespace stan {
18namespace math {
19namespace internal {
20
25typedef struct double_d {
26 double high;
27 double low;
28
30 double_d(double a) {
31 high = a;
32 low = 0;
33 }
34 double_d(double h, double l) {
35 high = h;
36 low = l;
37 }
38} double_d;
39
40using std::copysign;
41using std::fabs;
42using std::isfinite;
43using std::isinf;
44using std::isnan;
45using std::ldexp;
46
47inline double_d mul_d_d(double a, double b) {
48 double_d res;
49 // for some reason the std::fma only gives accurate result if compiling for
50 // instruction set that includes a fma instruction
51#if defined(__FMA__)
52 res.high = a * b;
53 if (isfinite(res.high)) {
54 res.low = fma(a, b, -res.high);
55 } else {
56 res.low = 0;
57 }
58#else
59 const double C = (2 << 26) + 1;
60 double p = a * C;
61 if (!isfinite(p)) {
62 res.high = p;
63 res.low = 0;
64 return res;
65 }
66 double hx = a - p + p;
67 double tx = a - hx;
68 p = b * C;
69 if (!isfinite(p)) {
70 res.high = p;
71 res.low = 0;
72 return res;
73 }
74 double hy = b - p + p;
75 double ty = b - hy;
76 p = hx * hy;
77 double q = hx * ty + tx * hy;
78 res.high = p + q;
79 if (isfinite(res.high)) {
80 res.low = p - res.high + q + tx * ty;
81 } else {
82 res.low = 0;
83 }
84#endif
85 return res;
86}
87
88// \cond
89static constexpr const char* double_d_src = STRINGIFY(
90 // \endcond
91 typedef struct {
92 double high;
93 double low;
94 } double_d;
95
96 inline double_d mul_d_d(double a, double b) {
97 double_d res;
98 res.high = a * b;
99 if (isfinite(res.high)) {
100 res.low = fma(a, b, -res.high);
101 } else {
102 res.low = 0;
103 }
104 return res;
105 }
106 // \cond
107 )
109 // \endcond
110
112 double_d res;
113 res.high = -a.high;
114 res.low = -a.low;
115 return res;
116 }
117
119 double_d res;
120 double high = a.high + b.high;
121 if (isfinite(high)) {
122 double low;
123 if (fabs(a.high) > fabs(b.high)) {
124 low = a.high - high + b.high;
125 } else {
126 low = b.high - high + a.high;
127 }
128 low += a.low + b.low;
129 res.high = high + low;
130 res.low = high - res.high + low;
131 } else {
132 res.high = high;
133 res.low = 0;
134 }
135 return res;
136 }
137
138 inline double_d add_dd_d(double_d a, double b) {
139 double_d res;
140 double high = a.high + b;
141 if (isfinite(high)) {
142 double low;
143 if (fabs(a.high) > fabs(b)) {
144 low = a.high - high + b;
145 } else {
146 low = b - high + a.high;
147 }
148 low += a.low;
149 res.high = high + low;
150 res.low = high - res.high + low;
151 } else {
152 res.high = high;
153 res.low = 0;
154 }
155 return res;
156 }
157
159 return add_dd_dd(a, neg(b));
160 }
161
162 inline double_d sub_dd_d(double_d a, double b) {
163 return add_dd_d(a, -b);
164 }
165
166 inline double_d sub_d_dd(double a, double_d b) {
167 return add_dd_d(neg(b), a);
168 }
169
171 double_d high = mul_d_d(a.high, b.high);
172 double_d res;
173 if (isfinite(high.high)) {
174 high.low += a.high * b.low + a.low * b.high;
175 res.high = high.high + high.low;
176 res.low = high.high - res.high + high.low;
177 } else {
178 high.low = 0;
179 return high;
180 }
181 return res;
182 }
183
184 inline double_d mul_dd_d(double_d a, double b) {
185 double_d high = mul_d_d(a.high, b);
186 double_d res;
187 if (isfinite(high.high)) {
188 high.low += a.low * b;
189 res.high = high.high + high.low;
190 res.low = high.high - res.high + high.low;
191 return res;
192 } else {
193 high.low = 0;
194 return high;
195 }
196 }
197
199 double high = a.high / b.high;
200 double_d res;
201 if (isfinite(high)) {
202 double_d u = mul_d_d(high, b.high);
203 if (isfinite(u.high)) {
204 double low
205 = (a.high - u.high - u.low + a.low - high * b.low) / b.high;
206 res.high = high + low;
207 res.low = high - res.high + low;
208 return res;
209 }
210 }
211 res.high = high;
212 res.low = 0;
213 return res;
214 }
215
216 inline double_d div_dd_d(double_d a, double b) {
217 double high = a.high / b;
218 double_d res;
219 if (isfinite(high)) {
220 double_d u = mul_d_d(high, b);
221 double low = (a.high - u.high - u.low + a.low) / b;
222 res.high = high + low;
223 res.low = high - res.high + low;
224 } else {
225 res.high = high;
226 res.low = 0;
227 }
228 return res;
229 }
230
231 inline double_d div_d_dd(double a, double_d b) {
232 double high = a / b.high;
233 double_d res;
234 if (isfinite(high)) {
235 double_d u = mul_d_d(high, b.high);
236 double low = (a - u.high - u.low - high * b.low) / b.high;
237 res.high = high + low;
238 res.low = high - res.high + low;
239 } else {
240 res.high = high;
241 res.low = 0;
242 }
243 return res;
244 }
245
246 inline double copysign_d_dd(double a, double_d b) {
247 return copysign(a, b.high);
248 }
249
250 inline double_d abs_dd(double_d a) { return a.high > 0 ? a : neg(a); }
251
252 inline bool isnan_dd(double_d a) { return isnan(a.high); }
253
254 inline bool isinf_dd(double_d a) { return isinf(a.high); }
255
256 inline bool lt_dd_dd(double_d a, double_d b) {
257 return a.high < b.high || (a.high == b.high && a.low < b.low);
258 }
259
260 inline bool lt_d_dd(double a, double_d b) {
261 return lt_dd_dd((double_d){a, 0}, b);
262 }
263
264 inline bool lt_dd_d(double_d a, double b) {
265 return lt_dd_dd(a, (double_d){b, 0});
266 }
267
268 inline bool gt_dd_dd(double_d a, double_d b) {
269 return a.high > b.high || (a.high == b.high && a.low > b.low);
270 }
271
272 inline bool gt_d_dd(double a, double_d b) {
273 return gt_dd_dd((double_d){a, 0}, b);
274 }
275
276 inline bool gt_dd_d(double_d a, double b) {
277 return gt_dd_dd(a, (double_d){b, 0});
278 }
279
280 inline bool le_dd_dd(double_d a, double_d b) {
281 return !gt_dd_dd(a, b);
282 } inline bool le_d_dd(double a, double_d b) {
283 return !gt_d_dd(a, b);
284 } inline bool le_dd_d(double_d a, double b) { return !gt_dd_d(a, b); }
285
286 inline bool ge_dd_dd(double_d a, double_d b) {
287 return !lt_dd_dd(a, b);
288 } inline bool ge_d_dd(double a, double_d b) {
289 return !lt_d_dd(a, b);
290 } inline bool ge_dd_d(double_d a, double b) { return !lt_dd_d(a, b); }
291 // \cond
292 ); // NOLINT(whitespace/parens)
293// \endcond
294
295inline double_d operator+(double_d a, double b) { return add_dd_d(a, b); }
296inline double_d operator+(double a, double_d b) { return add_dd_d(b, a); }
297inline double_d operator+(double_d a, double_d b) { return add_dd_dd(a, b); }
298
299inline double_d operator-(double a, double_d b) { return sub_d_dd(a, b); }
300inline double_d operator-(double_d a, double b) { return sub_dd_d(a, b); }
301inline double_d operator-(double_d a, double_d b) { return sub_dd_dd(a, b); }
302
303inline double_d operator*(double_d a, double b) { return mul_dd_d(a, b); }
304inline double_d operator*(double a, double_d b) { return mul_dd_d(b, a); }
305inline double_d operator*(double_d a, double_d b) { return mul_dd_dd(a, b); }
306
307inline double_d operator/(double_d a, double b) { return div_dd_d(a, b); }
308inline double_d operator/(double a, double_d b) { return div_d_dd(a, b); }
309inline double_d operator/(double_d a, double_d b) { return div_dd_dd(a, b); }
310
311inline double_d operator-(double_d a) { return neg(a); }
312
313inline bool operator<(double_d a, double_d b) { return lt_dd_dd(a, b); }
314inline bool operator<(double a, double_d b) { return lt_d_dd(a, b); }
315inline bool operator<(double_d a, double b) { return lt_dd_d(a, b); }
316inline bool operator>(double_d a, double_d b) { return gt_dd_dd(a, b); }
317inline bool operator>(double a, double_d b) { return gt_d_dd(a, b); }
318inline bool operator>(double_d a, double b) { return gt_dd_d(a, b); }
319inline bool operator<=(double_d a, double_d b) { return le_dd_dd(a, b); }
320inline bool operator<=(double a, double_d b) { return le_d_dd(a, b); }
321inline bool operator<=(double_d a, double b) { return le_dd_d(a, b); }
322inline bool operator>=(double_d a, double_d b) { return ge_dd_dd(a, b); }
323inline bool operator>=(double a, double_d b) { return ge_d_dd(a, b); }
324inline bool operator>=(double_d a, double b) { return ge_dd_d(a, b); }
325
326inline double_d fabs(double_d a) { return abs_dd(a); }
327inline bool isnan(double_d a) { return isnan_dd(a); }
328inline bool isinf(double_d a) { return isinf_dd(a); }
329inline double copysign(double a, double_d b) { return copysign_d_dd(a, b); }
330
331} // namespace internal
332} // namespace math
333} // namespace stan
334
335#endif
#define STRINGIFY2(A)
Alternative stringify, that also exposes the code for use in C++ host code.
Definition double_d.hpp:14
isfinite_< as_operation_cl_t< T > > isfinite(T &&a)
double_d operator*(double_d a, double b)
Definition double_d.hpp:303
bool le_dd_d(double_d a, double b)
Definition double_d.hpp:284
double_d operator/(double_d a, double b)
Definition double_d.hpp:307
double_d div_dd_dd(double_d a, double_d b)
Definition double_d.hpp:198
double_d mul_dd_d(double_d a, double b)
Definition double_d.hpp:184
double_d sub_d_dd(double a, double_d b)
Definition double_d.hpp:166
double_d mul_d_d(double a, double b)
Definition double_d.hpp:47
double_d abs_dd(double_d a)
Definition double_d.hpp:250
bool isinf(double_d a)
Definition double_d.hpp:328
double_d add_dd_dd(double_d a, double_d b)
Definition double_d.hpp:118
bool lt_d_dd(double a, double_d b)
Definition double_d.hpp:260
bool operator>(double_d a, double_d b)
Definition double_d.hpp:316
double_d sub_dd_dd(double_d a, double_d b)
Definition double_d.hpp:158
double_d operator-(double a, double_d b)
Definition double_d.hpp:299
double_d operator+(double_d a, double b)
Definition double_d.hpp:295
double_d sub_dd_d(double_d a, double b)
Definition double_d.hpp:162
bool le_d_dd(double a, double_d b)
Definition double_d.hpp:282
bool operator<=(double_d a, double_d b)
Definition double_d.hpp:319
double copysign(double a, double_d b)
Definition double_d.hpp:329
double_d div_dd_d(double_d a, double b)
Definition double_d.hpp:216
bool gt_dd_d(double_d a, double b)
Definition double_d.hpp:276
bool operator<(double_d a, double_d b)
Definition double_d.hpp:313
double_d div_d_dd(double a, double_d b)
Definition double_d.hpp:231
bool isinf_dd(double_d a)
Definition double_d.hpp:254
double copysign_d_dd(double a, double_d b)
Definition double_d.hpp:246
bool isnan_dd(double_d a)
Definition double_d.hpp:252
double_d neg(double_d a)
Definition double_d.hpp:111
double_d mul_dd_dd(double_d a, double_d b)
Definition double_d.hpp:170
bool lt_dd_dd(double_d a, double_d b)
Definition double_d.hpp:256
bool ge_d_dd(double a, double_d b)
Definition double_d.hpp:288
bool isnan(double_d a)
Definition double_d.hpp:327
bool gt_d_dd(double a, double_d b)
Definition double_d.hpp:272
bool ge_dd_d(double_d a, double b)
Definition double_d.hpp:290
bool ge_dd_dd(double_d a, double_d b)
Definition double_d.hpp:286
bool le_dd_dd(double_d a, double_d b)
Definition double_d.hpp:280
double_d add_dd_d(double_d a, double b)
Definition double_d.hpp:138
bool lt_dd_d(double_d a, double b)
Definition double_d.hpp:264
bool operator>=(double_d a, double_d b)
Definition double_d.hpp:322
bool gt_dd_dd(double_d a, double_d b)
Definition double_d.hpp:268
fvar< return_type_t< T1, T2, T3 > > fma(const fvar< T1 > &x1, const fvar< T2 > &x2, const fvar< T3 > &x3)
The fused multiply-add operation (C99).
Definition fma.hpp:60
fvar< T > fabs(const fvar< T > &x)
Definition fabs.hpp:15
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
bool isnan(const stan::math::var &a)
Checks if the given number is NaN.
Definition std_isnan.hpp:18
bool isinf(const stan::math::var &a)
Return 1 if the specified argument is positive infinity or negative infinity and 0 otherwise.
Definition std_isinf.hpp:16
#define STRINGIFY(...)
Definition stringify.hpp:9
double low
Definition double_d.hpp:27
double_d(double h, double l)
Definition double_d.hpp:34
double_d(double a)
Definition double_d.hpp:30
double high
Definition double_d.hpp:26
Double double - a 128 bit floating point number defined as an exact sum of 2 doubles.
Definition double_d.hpp:25