dune-common  2.8.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 #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 
21 namespace 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 
409 namespace 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
#define DUNE_BINARY_OP(OP)
Definition: debugalign.hh:235
#define DUNE_UNARY_FUNC(name)
#define DUNE_ASSIGN_OP(OP)
Definition: debugalign.hh:194
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
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition: bigunsignedint.hh:273
bigunsignedint< k > operator+(const bigunsignedint< k > &x, std::uintmax_t y)
Definition: bigunsignedint.hh:530
bigunsignedint< k > operator-(const bigunsignedint< k > &x, std::uintmax_t y)
Definition: bigunsignedint.hh:537
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
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:11
T max_value(const AlignedNumber< T, align > &val)
Definition: debugalign.hh:468
T min_value(const AlignedNumber< T, align > &val)
Definition: debugalign.hh:474
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:434
auto max(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:412
Default exception class for range errors.
Definition: exceptions.hh:252