DUNE PDELab (git)

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