dune-common  2.9.0
quadmath.hh
Go to the documentation of this file.
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 (C) 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 
23 namespace 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 
411 namespace 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
#define DUNE_BINARY_OP(OP)
Definition: debugalign.hh:248
#define DUNE_UNARY_FUNC(name)
#define DUNE_ASSIGN_OP(OP)
Definition: debugalign.hh:207
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:43
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition: bigunsignedint.hh:278
bigunsignedint< k > operator+(const bigunsignedint< k > &x, std::uintmax_t y)
Definition: bigunsignedint.hh:535
bigunsignedint< k > operator-(const bigunsignedint< k > &x, std::uintmax_t y)
Definition: bigunsignedint.hh:542
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
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
T max_value(const AlignedNumber< T, align > &val)
Definition: debugalign.hh:481
T min_value(const AlignedNumber< T, align > &val)
Definition: debugalign.hh:487
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:447
auto max(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:425
Default exception class for range errors.
Definition: exceptions.hh:254