Dune Core Modules (2.8.0)

quadmath.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_QUADMATH_HH
4#define DUNE_QUADMATH_HH
5
6#if HAVE_QUADMATH
7#include <quadmath.h>
8
9#include <cmath>
10#include <cstddef>
11#include <cstdint>
12#include <cstdlib> // abs
13#include <istream>
14#include <ostream>
15#include <type_traits>
16#include <utility>
17
20
21namespace Dune
22{
23 namespace Impl
24 {
25 // forward declaration
26 class Float128;
27
28 } // end namespace Impl
29
30 using Impl::Float128;
31
32 // The purpose of this namespace is to move the `<cmath>` function overloads
33 // out of namespace `Dune`, see AlignedNumber in debugalign.hh.
34 namespace Impl
35 {
36 using float128_t = __float128;
37
39 class Float128
40 {
41 float128_t value_ = 0.0q;
42
43 public:
44 constexpr Float128() = default;
45 constexpr Float128(const float128_t& value) noexcept
46 : value_(value)
47 {}
48
49 // constructor from any floating-point or integer type
50 template <class T,
51 std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
52 constexpr Float128(const T& value) noexcept
53 : value_(value)
54 {}
55
56 // constructor from pointer to null-terminated byte string
57 explicit Float128(const char* str) noexcept
58 : value_(strtoflt128(str, NULL))
59 {}
60
61 // accessors
62 constexpr operator float128_t() const noexcept { return value_; }
63
64 constexpr float128_t const& value() const noexcept { return value_; }
65 constexpr float128_t& value() noexcept { return value_; }
66
67 // I/O
68 template<class CharT, class Traits>
69 friend std::basic_istream<CharT, Traits>&
70 operator>>(std::basic_istream<CharT, Traits>& in, Float128& x)
71 {
72 std::string buf;
73 buf.reserve(128);
74 in >> buf;
75 x.value() = strtoflt128(buf.c_str(), NULL);
76 return in;
77 }
78
79 template<class CharT, class Traits>
80 friend std::basic_ostream<CharT, Traits>&
81 operator<<(std::basic_ostream<CharT, Traits>& out, const Float128& x)
82 {
83 const std::size_t bufSize = 128;
84 CharT buf[128];
85
86 std::string format = "%." + std::to_string(out.precision()) + "Q" +
87 ((out.flags() | std::ios_base::scientific) ? "e" : "f");
88 const int numChars = quadmath_snprintf(buf, bufSize, format.c_str(), x.value());
89 if (std::size_t(numChars) >= bufSize) {
90 DUNE_THROW(Dune::RangeError, "Failed to print Float128 value: buffer overflow");
91 }
92 out << buf;
93 return out;
94 }
95
96 // Increment, decrement
97 constexpr Float128& operator++() noexcept { ++value_; return *this; }
98 constexpr Float128& operator--() noexcept { --value_; return *this; }
99
100 constexpr Float128 operator++(int) noexcept { Float128 tmp{*this}; ++value_; return tmp; }
101 constexpr Float128 operator--(int) noexcept { Float128 tmp{*this}; --value_; return tmp; }
102
103 // unary operators
104 constexpr Float128 operator+() const noexcept { return Float128{+value_}; }
105 constexpr Float128 operator-() const noexcept { return Float128{-value_}; }
106
107 // assignment operators
108#define DUNE_ASSIGN_OP(OP) \
109 constexpr Float128& operator OP(const Float128& u) noexcept \
110 { \
111 value_ OP float128_t(u); \
112 return *this; \
113 } \
114 static_assert(true, "Require semicolon to unconfuse editors")
115
116 DUNE_ASSIGN_OP(+=);
117 DUNE_ASSIGN_OP(-=);
118
119 DUNE_ASSIGN_OP(*=);
120 DUNE_ASSIGN_OP(/=);
121
122#undef DUNE_ASSIGN_OP
123
124 }; // end class Float128
125
126 // binary operators:
127 // For symmetry provide overloads with arithmetic types
128 // in the first or second argument.
129#define DUNE_BINARY_OP(OP) \
130 constexpr Float128 operator OP(const Float128& t, \
131 const Float128& u) noexcept \
132 { \
133 return Float128{float128_t(t) OP float128_t(u)}; \
134 } \
135 constexpr Float128 operator OP(const float128_t& t, \
136 const Float128& u) noexcept \
137 { \
138 return Float128{t OP float128_t(u)}; \
139 } \
140 constexpr Float128 operator OP(const Float128& t, \
141 const float128_t& u) noexcept \
142 { \
143 return Float128{float128_t(t) OP u}; \
144 } \
145 template <class T, \
146 std::enable_if_t<std::is_arithmetic<T>::value, int> = 0> \
147 constexpr Float128 operator OP(const T& t, \
148 const Float128& u) noexcept \
149 { \
150 return Float128{float128_t(t) OP float128_t(u)}; \
151 } \
152 template <class U, \
153 std::enable_if_t<std::is_arithmetic<U>::value, int> = 0> \
154 constexpr Float128 operator OP(const Float128& t, \
155 const U& u) noexcept \
156 { \
157 return Float128{float128_t(t) OP float128_t(u)}; \
158 } \
159 static_assert(true, "Require semicolon to unconfuse editors")
160
161 DUNE_BINARY_OP(+);
162 DUNE_BINARY_OP(-);
163 DUNE_BINARY_OP(*);
164 DUNE_BINARY_OP(/);
165
166#undef DUNE_BINARY_OP
167
168 // logical operators:
169 // For symmetry provide overloads with arithmetic types
170 // in the first or second argument.
171#define DUNE_BINARY_BOOL_OP(OP) \
172 constexpr bool operator OP(const Float128& t, \
173 const Float128& u) noexcept \
174 { \
175 return float128_t(t) OP float128_t(u); \
176 } \
177 template <class T, \
178 std::enable_if_t<std::is_arithmetic<T>::value, int> = 0> \
179 constexpr bool operator OP(const T& t, \
180 const Float128& u) noexcept \
181 { \
182 return float128_t(t) OP float128_t(u); \
183 } \
184 template <class U, \
185 std::enable_if_t<std::is_arithmetic<U>::value, int> = 0> \
186 constexpr bool operator OP(const Float128& t, \
187 const U& u) noexcept \
188 { \
189 return float128_t(t) OP float128_t(u); \
190 } \
191 static_assert(true, "Require semicolon to unconfuse editors")
192
193 DUNE_BINARY_BOOL_OP(==);
194 DUNE_BINARY_BOOL_OP(!=);
195 DUNE_BINARY_BOOL_OP(<);
196 DUNE_BINARY_BOOL_OP(>);
197 DUNE_BINARY_BOOL_OP(<=);
198 DUNE_BINARY_BOOL_OP(>=);
199
200#undef DUNE_BINARY_BOOL_OP
201
202 // Overloads for the cmath functions
203
204 // function with name `name` redirects to quadmath function `func`
205#define DUNE_UNARY_FUNC(name,func) \
206 inline Float128 name(const Float128& u) noexcept \
207 { \
208 return Float128{func (float128_t(u))}; \
209 } \
210 static_assert(true, "Require semicolon to unconfuse editors")
211
212 // like DUNE_UNARY_FUNC but with cutom return type
213#define DUNE_CUSTOM_UNARY_FUNC(type,name,func) \
214 inline type name(const Float128& u) noexcept \
215 { \
216 return (type)(func (float128_t(u))); \
217 } \
218 static_assert(true, "Require semicolon to unconfuse editors")
219
220 // redirects to quadmath function with two arguments
221#define DUNE_BINARY_FUNC(name,func) \
222 inline Float128 name(const Float128& t, \
223 const Float128& u) noexcept \
224 { \
225 return Float128{func (float128_t(t), float128_t(u))}; \
226 } \
227 static_assert(true, "Require semicolon to unconfuse editors")
228
229 DUNE_UNARY_FUNC(abs, fabsq);
230 DUNE_UNARY_FUNC(acos, acosq);
231 DUNE_UNARY_FUNC(acosh, acoshq);
232 DUNE_UNARY_FUNC(asin, asinq);
233 DUNE_UNARY_FUNC(asinh, asinhq);
234 DUNE_UNARY_FUNC(atan, atanq);
235 DUNE_UNARY_FUNC(atanh, atanhq);
236 DUNE_UNARY_FUNC(cbrt, cbrtq);
237 DUNE_UNARY_FUNC(ceil, ceilq);
238 DUNE_UNARY_FUNC(cos, cosq);
239 DUNE_UNARY_FUNC(cosh, coshq);
240 DUNE_UNARY_FUNC(erf, erfq);
241 DUNE_UNARY_FUNC(erfc, erfcq);
242 DUNE_UNARY_FUNC(exp, expq);
243 DUNE_UNARY_FUNC(expm1, expm1q);
244 DUNE_UNARY_FUNC(fabs, fabsq);
245 DUNE_UNARY_FUNC(floor, floorq);
246 DUNE_CUSTOM_UNARY_FUNC(int, ilogb, ilogbq);
247 DUNE_UNARY_FUNC(lgamma, lgammaq);
248 DUNE_CUSTOM_UNARY_FUNC(long long int, llrint, llrintq);
249 DUNE_CUSTOM_UNARY_FUNC(long long int, llround, llroundq);
250 DUNE_UNARY_FUNC(log, logq);
251 DUNE_UNARY_FUNC(log10, log10q);
252 DUNE_UNARY_FUNC(log1p, log1pq);
253 DUNE_UNARY_FUNC(log2, log2q);
254 // DUNE_UNARY_FUNC(logb, logbq); // not available in gcc5
255 DUNE_CUSTOM_UNARY_FUNC(long int, lrint, lrintq);
256 DUNE_CUSTOM_UNARY_FUNC(long int, lround, lroundq);
257 DUNE_UNARY_FUNC(nearbyint, nearbyintq);
258 DUNE_BINARY_FUNC(nextafter, nextafterq);
259 DUNE_BINARY_FUNC(pow, powq); // overload for integer argument see below
260 DUNE_UNARY_FUNC(rint, rintq);
261 DUNE_UNARY_FUNC(round, roundq);
262 DUNE_UNARY_FUNC(sin, sinq);
263 DUNE_UNARY_FUNC(sinh, sinhq);
264 DUNE_UNARY_FUNC(sqrt, sqrtq);
265 DUNE_UNARY_FUNC(tan, tanq);
266 DUNE_UNARY_FUNC(tanh, tanhq);
267 DUNE_UNARY_FUNC(tgamma, tgammaq);
268 DUNE_UNARY_FUNC(trunc, truncq);
269
270 DUNE_CUSTOM_UNARY_FUNC(bool, isfinite, finiteq);
271 DUNE_CUSTOM_UNARY_FUNC(bool, isinf, isinfq);
272 DUNE_CUSTOM_UNARY_FUNC(bool, isnan, isnanq);
273 DUNE_CUSTOM_UNARY_FUNC(bool, signbit, signbitq);
274
275#undef DUNE_UNARY_FUNC
276#undef DUNE_CUSTOM_UNARY_FUNC
277#undef DUNE_BINARY_FUNC
278
279 // like DUNE_BINARY_FUNC but provide overloads with arithmetic
280 // types in the first or second argument.
281#define DUNE_BINARY_ARITHMETIC_FUNC(name,func) \
282 inline Float128 name(const Float128& t, \
283 const Float128& u) noexcept \
284 { \
285 return Float128{func (float128_t(t), float128_t(u))}; \
286 } \
287 template <class T, \
288 std::enable_if_t<std::is_arithmetic<T>::value, int> = 0> \
289 inline Float128 name(const T& t, \
290 const Float128& u) noexcept \
291 { \
292 return Float128{func (float128_t(t), float128_t(u))}; \
293 } \
294 template <class U, \
295 std::enable_if_t<std::is_arithmetic<U>::value, int> = 0> \
296 inline Float128 name(const Float128& t, \
297 const U& u) noexcept \
298 { \
299 return Float128{func (float128_t(t), float128_t(u))}; \
300 } \
301 static_assert(true, "Require semicolon to unconfuse editors")
302
303 DUNE_BINARY_ARITHMETIC_FUNC(atan2,atan2q);
304 DUNE_BINARY_ARITHMETIC_FUNC(copysign,copysignq);
305 DUNE_BINARY_ARITHMETIC_FUNC(fdim,fdimq);
306 DUNE_BINARY_ARITHMETIC_FUNC(fmax,fmaxq);
307 DUNE_BINARY_ARITHMETIC_FUNC(fmin,fminq);
308 DUNE_BINARY_ARITHMETIC_FUNC(fmod,fmodq);
309 DUNE_BINARY_ARITHMETIC_FUNC(hypot,hypotq);
310 DUNE_BINARY_ARITHMETIC_FUNC(remainder,remainderq);
311
312#undef DUNE_BINARY_ARITHMETIC_FUNC
313
314 // some more cmath functions with special signature
315
316 inline Float128 fma(const Float128& t, const Float128& u, const Float128& v)
317 {
318 return Float128{fmaq(float128_t(t),float128_t(u),float128_t(v))};
319 }
320
321 inline Float128 frexp(const Float128& u, int* p)
322 {
323 return Float128{frexpq(float128_t(u), p)};
324 }
325
326 inline Float128 ldexp(const Float128& u, int p)
327 {
328 return Float128{ldexpq(float128_t(u), p)};
329 }
330
331 inline Float128 remquo(const Float128& t, const Float128& u, int* quo)
332 {
333 return Float128{remquoq(float128_t(t), float128_t(u), quo)};
334 }
335
336 inline Float128 scalbln(const Float128& u, long int e)
337 {
338 return Float128{scalblnq(float128_t(u), e)};
339 }
340
341 inline Float128 scalbn(const Float128& u, int e)
342 {
343 return Float128{scalbnq(float128_t(u), e)};
344 }
345
347 // NOTE: This is much faster than a pow(x, Float128(p)) call
348 // NOTE: This is a modified version of boost::math::cstdfloat::detail::pown
349 // (adapted to the type Float128) that is part of the Boost 1.65 Math toolkit 2.8.0
350 // and is implemented by Christopher Kormanyos, John Maddock, and Paul A. Bristow,
351 // distributed under the Boost Software License, Version 1.0
352 // (See http://www.boost.org/LICENSE_1_0.txt)
353 template <class Int,
354 std::enable_if_t<std::is_integral<Int>::value, int> = 0>
355 inline Float128 pow(const Float128& x, const Int p)
356 {
357 static const Float128 max_value = FLT128_MAX;
358 static const Float128 min_value = FLT128_MIN;
359 static const Float128 inf_value = float128_t{1} / float128_t{0};
360
361 const bool isneg = (x < 0);
362 const bool isnan = (x != x);
363 const bool isinf = (isneg ? bool(-x > max_value) : bool(+x > max_value));
364
365 if (isnan) { return x; }
366 if (isinf) { return Float128{nanq("")}; }
367
368 const Float128 abs_x = (isneg ? -x : x);
369 if (p < Int(0)) {
370 if (abs_x < min_value)
371 return (isneg ? -inf_value : +inf_value);
372 else
373 return Float128(1) / pow(x, Int(-p));
374 }
375
376 if (p == Int(0)) { return Float128(1); }
377 if (p == Int(1)) { return x; }
378 if (abs_x > max_value)
379 return (isneg ? -inf_value : +inf_value);
380
381 if (p == Int(2)) { return (x * x); }
382 if (p == Int(3)) { return ((x * x) * x); }
383 if (p == Int(4)) { const Float128 x2 = (x * x); return (x2 * x2); }
384
385 Float128 result = ((p % Int(2)) != Int(0)) ? x : Float128(1);
386 Float128 xn = x; // binary powers of x
387
388 Int p2 = p;
389 while (Int(p2 /= 2) != Int(0)) {
390 xn *= xn; // Square xn for each binary power
391
392 const bool has_binary_power = (Int(p2 % Int(2)) != Int(0));
393 if (has_binary_power)
394 result *= xn;
395 }
396
397 return result;
398 }
399
400
401 } // end namespace Impl
402
403 template <>
404 struct IsNumber<Impl::Float128>
405 : public std::true_type {};
406
407} // end namespace Dune
408
409namespace std
410{
411#ifndef NO_STD_NUMERIC_LIMITS_SPECIALIZATION
412 template <>
413 class numeric_limits<Dune::Impl::Float128>
414 {
415 using Float128 = Dune::Impl::Float128;
416 using float128_t = Dune::Impl::float128_t;
417
418 public:
419 static constexpr bool is_specialized = true;
420 static constexpr Float128 min() noexcept { return FLT128_MIN; }
421 static constexpr Float128 max() noexcept { return FLT128_MAX; }
422 static constexpr Float128 lowest() noexcept { return -FLT128_MAX; }
423 static constexpr int digits = FLT128_MANT_DIG;
424 static constexpr int digits10 = 34;
425 static constexpr int max_digits10 = 36;
426 static constexpr bool is_signed = true;
427 static constexpr bool is_integer = false;
428 static constexpr bool is_exact = false;
429 static constexpr int radix = 2;
430 static constexpr Float128 epsilon() noexcept { return FLT128_EPSILON; }
431 static constexpr Float128 round_error() noexcept { return float128_t{0.5}; }
432 static constexpr int min_exponent = FLT128_MIN_EXP;
433 static constexpr int min_exponent10 = FLT128_MIN_10_EXP;
434 static constexpr int max_exponent = FLT128_MAX_EXP;
435 static constexpr int max_exponent10 = FLT128_MAX_10_EXP;
436 static constexpr bool has_infinity = true;
437 static constexpr bool has_quiet_NaN = true;
438 static constexpr bool has_signaling_NaN = false;
439 static constexpr float_denorm_style has_denorm = denorm_present;
440 static constexpr bool has_denorm_loss = false;
441 static constexpr Float128 infinity() noexcept { return float128_t{1}/float128_t{0}; }
442 static Float128 quiet_NaN() noexcept { return nanq(""); }
443 static constexpr Float128 signaling_NaN() noexcept { return float128_t{}; }
444 static constexpr Float128 denorm_min() noexcept { return FLT128_DENORM_MIN; }
445 static constexpr bool is_iec559 = true;
446 static constexpr bool is_bounded = false;
447 static constexpr bool is_modulo = false;
448 static constexpr bool traps = false;
449 static constexpr bool tinyness_before = false;
450 static constexpr float_round_style round_style = round_to_nearest;
451 };
452#endif
453} // end namespace std
454
455#endif // HAVE_QUADMATH
456#endif // DUNE_QUADMATH_HH
Default exception class for range errors.
Definition: exceptions.hh:252
A few common exception classes.
Traits for type conversions and type information.
Stream & operator>>(Stream &stream, std::tuple< Ts... > &t)
Read a std::tuple.
Definition: streamoperators.hh:41
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
I round(const T &val, typename EpsilonType< T >::Type epsilon)
round using epsilon
Definition: float_cmp.cc:309
I trunc(const T &val, typename EpsilonType< T >::Type epsilon)
truncate using epsilon
Definition: float_cmp.cc:405
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Dune namespace.
Definition: alignedallocator.hh:11
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)