quadmath.hh 17.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_QUADMATH_HH
#define DUNE_QUADMATH_HH

#if HAVE_QUADMATH
#include <quadmath.h>

#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstdlib> // abs
#include <istream>
#include <ostream>
#include <type_traits>
#include <utility>

#include <dune/common/typetraits.hh>

namespace Dune
{
  namespace Impl
  {
    // forward declaration
    class Float128;

  } // end namespace Impl

  using Impl::Float128;

  // The purpose of this namespace is to move the `<cmath>` function overloads
  // out of namespace `Dune`, see AlignedNumber in debugalign.hh.
  namespace Impl
  {
    using float128_t = __float128;

    /// Wrapper for quad-precision type __float128
    class Float128
    {
      float128_t value_ = 0.0q;

    public:
      constexpr Float128() = default;
      constexpr Float128(const float128_t& value) noexcept
        : value_(value)
      {}

      // constructor from any floating-point or integer type
      template <class T,
        std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>
      constexpr Float128(const T& value) noexcept
        : value_(value)
      {}

      // constructor from pointer to null-terminated byte string
      Float128(const char* str) noexcept
        : value_(strtoflt128(str, NULL))
      {}

      // accessors
      constexpr operator float128_t() const noexcept { return value_; }

      constexpr float128_t const& value() const noexcept { return value_; }
      constexpr float128_t&       value() noexcept       { return value_; }

      // I/O
      template<class CharT, class Traits>
      friend std::basic_istream<CharT, Traits>&
      operator>>(std::basic_istream<CharT, Traits>& in, Float128& x)
      {
        std::string buf;
        buf.reserve(128);
        in >> buf;
        x.value() = strtoflt128(buf.c_str(), NULL);
        return in;
      }

      template<class CharT, class Traits>
      friend std::basic_ostream<CharT, Traits>&
      operator<<(std::basic_ostream<CharT, Traits>& out, const Float128& x)
      {
        const std::size_t bufSize = 128;
        CharT buf[128];

        std::string format = "%." + std::to_string(out.precision()) + "Q" +
                              ((out.flags() | std::ios_base::scientific) ? "e" : "f");
        const int numChars = quadmath_snprintf(buf, bufSize, format.c_str(), x.value());
        if (std::size_t(numChars) >= bufSize) {
          DUNE_THROW(Dune::RangeError, "Failed to print Float128 value: buffer overflow");
        }
        out << buf;
        return out;
      }

      // Increment, decrement
      constexpr Float128& operator++() noexcept { ++value_; return *this; }
      constexpr Float128& operator--() noexcept { --value_; return *this; }

      constexpr Float128 operator++(int) noexcept { Float128 tmp{*this}; ++value_; return tmp; }
      constexpr Float128 operator--(int) noexcept { Float128 tmp{*this}; --value_; return tmp; }

      // unary operators
      constexpr Float128 operator+() const noexcept { return Float128{+value_}; }
      constexpr Float128 operator-() const noexcept { return Float128{-value_}; }

      // assignment operators
#define DUNE_ASSIGN_OP(OP)                                              \
      constexpr Float128& operator OP(const Float128& u) noexcept       \
      {                                                                 \
        value_ OP float128_t(u);                                        \
        return *this;                                                   \
      }                                                                 \
      static_assert(true, "Require semicolon to unconfuse editors")

      DUNE_ASSIGN_OP(+=);
      DUNE_ASSIGN_OP(-=);

      DUNE_ASSIGN_OP(*=);
      DUNE_ASSIGN_OP(/=);

#undef DUNE_ASSIGN_OP

    }; // end class Float128

    // binary operators:
    // For symmetry provide overloads with arithmetic types
    // in the first or second argument.
#define DUNE_BINARY_OP(OP)                                              \
    constexpr Float128 operator OP(const Float128& t,                   \
                                   const Float128& u) noexcept          \
    {                                                                   \
      return Float128{float128_t(t) OP float128_t(u)};                  \
    }                                                                   \
    template <class T,                                                  \
      std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>          \
    constexpr Float128 operator OP(const T& t,                          \
                                   const Float128& u) noexcept          \
    {                                                                   \
      return Float128{float128_t(t) OP float128_t(u)};                  \
    }                                                                   \
    template <class U,                                                  \
      std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>          \
    constexpr Float128 operator OP(const Float128& t,                   \
                                   const U& u) noexcept                 \
    {                                                                   \
      return Float128{float128_t(t) OP float128_t(u)};                  \
    }                                                                   \
    static_assert(true, "Require semicolon to unconfuse editors")

    DUNE_BINARY_OP(+);
    DUNE_BINARY_OP(-);
    DUNE_BINARY_OP(*);
    DUNE_BINARY_OP(/);

#undef DUNE_BINARY_OP

    // logical operators:
    // For symmetry provide overloads with arithmetic types
    // in the first or second argument.
#define DUNE_BINARY_BOOL_OP(OP)                                         \
    constexpr bool operator OP(const Float128& t,                       \
                               const Float128& u) noexcept              \
    {                                                                   \
      return float128_t(t) OP float128_t(u);                            \
    }                                                                   \
    template <class T,                                                  \
      std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>          \
    constexpr bool operator OP(const T& t,                              \
                               const Float128& u) noexcept              \
    {                                                                   \
      return float128_t(t) OP float128_t(u);                            \
    }                                                                   \
    template <class U,                                                  \
      std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>          \
    constexpr bool operator OP(const Float128& t,                       \
                               const U& u) noexcept                     \
    {                                                                   \
      return float128_t(t) OP float128_t(u);                            \
    }                                                                   \
    static_assert(true, "Require semicolon to unconfuse editors")

    DUNE_BINARY_BOOL_OP(==);
    DUNE_BINARY_BOOL_OP(!=);
    DUNE_BINARY_BOOL_OP(<);
    DUNE_BINARY_BOOL_OP(>);
    DUNE_BINARY_BOOL_OP(<=);
    DUNE_BINARY_BOOL_OP(>=);

#undef DUNE_BINARY_BOOL_OP

    // Overloads for the cmath functions

    // function with name `name` redirects to quadmath function `func`
#define DUNE_UNARY_FUNC(name,func)                                         \
    inline Float128 name(const Float128& u) noexcept                       \
    {                                                                      \
      return Float128{func (float128_t(u))};                               \
    }                                                                      \
    static_assert(true, "Require semicolon to unconfuse editors")

    // like DUNE_UNARY_FUNC but with cutom return type
#define DUNE_CUSTOM_UNARY_FUNC(type,name,func)                             \
    inline type name(const Float128& u) noexcept                           \
    {                                                                      \
      return (type)(func (float128_t(u)));                                 \
    }                                                                      \
    static_assert(true, "Require semicolon to unconfuse editors")

    // redirects to quadmath function with two arguments
#define DUNE_BINARY_FUNC(name,func)                                        \
    inline Float128 name(const Float128& t,                                \
                         const Float128& u) noexcept                       \
    {                                                                      \
      return Float128{func (float128_t(t), float128_t(u))};                \
    }                                                                      \
    static_assert(true, "Require semicolon to unconfuse editors")

    DUNE_UNARY_FUNC(abs, fabsq);
    DUNE_UNARY_FUNC(acos, acosq);
    DUNE_UNARY_FUNC(acosh, acoshq);
    DUNE_UNARY_FUNC(asin, asinq);
    DUNE_UNARY_FUNC(asinh, asinhq);
    DUNE_UNARY_FUNC(atan, atanq);
    DUNE_UNARY_FUNC(atanh, atanhq);
    DUNE_UNARY_FUNC(cbrt, cbrtq);
    DUNE_UNARY_FUNC(ceil, ceilq);
    DUNE_UNARY_FUNC(cos, cosq);
    DUNE_UNARY_FUNC(cosh, coshq);
    DUNE_UNARY_FUNC(erf, erfq);
    DUNE_UNARY_FUNC(erfc, erfcq);
    DUNE_UNARY_FUNC(exp, expq);
    DUNE_UNARY_FUNC(expm1, expm1q);
    DUNE_UNARY_FUNC(fabs, fabsq);
    DUNE_UNARY_FUNC(floor, floorq);
    DUNE_CUSTOM_UNARY_FUNC(int, ilogb, ilogbq);
    DUNE_UNARY_FUNC(lgamma, lgammaq);
    DUNE_CUSTOM_UNARY_FUNC(long long int, llrint, llrintq);
    DUNE_CUSTOM_UNARY_FUNC(long long int, llround, llroundq);
    DUNE_UNARY_FUNC(log, logq);
    DUNE_UNARY_FUNC(log10, log10q);
    DUNE_UNARY_FUNC(log1p, log1pq);
    DUNE_UNARY_FUNC(log2, log2q);
    // DUNE_UNARY_FUNC(logb, logbq); // not available in gcc5
    DUNE_CUSTOM_UNARY_FUNC(long int, lrint, lrintq);
    DUNE_CUSTOM_UNARY_FUNC(long int, lround, lroundq);
    DUNE_UNARY_FUNC(nearbyint, nearbyintq);
    DUNE_BINARY_FUNC(nextafter, nextafterq);
    DUNE_BINARY_FUNC(pow, powq); // overload for integer argument see below
    DUNE_UNARY_FUNC(rint, rintq);
    DUNE_UNARY_FUNC(round, roundq);
    DUNE_UNARY_FUNC(sin, sinq);
    DUNE_UNARY_FUNC(sinh, sinhq);
    DUNE_UNARY_FUNC(sqrt, sqrtq);
    DUNE_UNARY_FUNC(tan, tanq);
    DUNE_UNARY_FUNC(tanh, tanhq);
    DUNE_UNARY_FUNC(tgamma, tgammaq);
    DUNE_UNARY_FUNC(trunc, truncq);

    DUNE_CUSTOM_UNARY_FUNC(bool, isfinite, finiteq);
    DUNE_CUSTOM_UNARY_FUNC(bool, isinf, isinfq);
    DUNE_CUSTOM_UNARY_FUNC(bool, isnan, isnanq);
    DUNE_CUSTOM_UNARY_FUNC(bool, signbit, signbitq);

#undef DUNE_UNARY_FUNC
#undef DUNE_CUSTOM_UNARY_FUNC
#undef DUNE_BINARY_FUNC

    // like DUNE_BINARY_FUNC but provide overloads with arithmetic
    // types in the first or second argument.
#define DUNE_BINARY_ARITHMETIC_FUNC(name,func)                          \
    inline Float128 name(const Float128& t,                             \
                         const Float128& u) noexcept                    \
    {                                                                   \
      return Float128{func (float128_t(t), float128_t(u))};             \
    }                                                                   \
    template <class T,                                                  \
      std::enable_if_t<std::is_arithmetic<T>::value, int> = 0>          \
    inline Float128 name(const T& t,                                    \
                         const Float128& u) noexcept                    \
    {                                                                   \
      return Float128{func (float128_t(t), float128_t(u))};             \
    }                                                                   \
    template <class U,                                                  \
      std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>          \
    inline Float128 name(const Float128& t,                             \
                         const U& u) noexcept                           \
    {                                                                   \
      return Float128{func (float128_t(t), float128_t(u))};             \
    }                                                                   \
    static_assert(true, "Require semicolon to unconfuse editors")

    DUNE_BINARY_ARITHMETIC_FUNC(atan2,atan2q);
    DUNE_BINARY_ARITHMETIC_FUNC(copysign,copysignq);
    DUNE_BINARY_ARITHMETIC_FUNC(fdim,fdimq);
    DUNE_BINARY_ARITHMETIC_FUNC(fmax,fmaxq);
    DUNE_BINARY_ARITHMETIC_FUNC(fmin,fminq);
    DUNE_BINARY_ARITHMETIC_FUNC(fmod,fmodq);
    DUNE_BINARY_ARITHMETIC_FUNC(hypot,hypotq);
    DUNE_BINARY_ARITHMETIC_FUNC(remainder,remainderq);

#undef DUNE_BINARY_ARITHMETIC_FUNC

    // some more cmath functions with special signature

    inline Float128 fma(const Float128& t, const Float128& u, const Float128& v)
    {
      return Float128{fmaq(float128_t(t),float128_t(u),float128_t(v))};
    }

    inline Float128 frexp(const Float128& u, int* p)
    {
      return Float128{frexpq(float128_t(u), p)};
    }

    inline Float128 ldexp(const Float128& u, int p)
    {
      return Float128{ldexpq(float128_t(u), p)};
    }

    inline Float128 remquo(const Float128& t, const Float128& u, int* quo)
    {
      return Float128{remquoq(float128_t(t), float128_t(u), quo)};
    }

    inline Float128 scalbln(const Float128& u, long int e)
    {
      return Float128{scalblnq(float128_t(u), e)};
    }

    inline Float128 scalbn(const Float128& u, int e)
    {
      return Float128{scalbnq(float128_t(u), e)};
    }

    /// \brief Overload of `pow` function for integer exponents.
    // NOTE: This is much faster than a pow(x, Float128(p)) call
    // NOTE: Derived from the boost::math::cstdfloat::detail::pown implementation
    template <class Int,
      std::enable_if_t<std::is_integral<Int>::value, int> = 0>
    inline Float128 pow(const Float128& x, const Int p)
    {
      static const Float128 max_value = FLT128_MAX;
      static const Float128 min_value = FLT128_MIN;
      static const Float128 inf_value = float128_t{1} / float128_t{0};

      const bool isneg = (x < 0);
      const bool isnan = (x != x);
      const bool isinf = (isneg ? bool(-x > max_value) : bool(+x > max_value));

      if (isnan) { return x; }
      if (isinf) { return Float128{nanq("")}; }

      const Float128 abs_x = (isneg ? -x : x);
      if (p < Int(0)) {
        if (abs_x < min_value)
          return (isneg ? -inf_value : +inf_value);
        else
          return Float128(1) / pow(x, Int(-p));
      }

      if (p == Int(0)) { return Float128(1); }
      if (p == Int(1)) { return x; }
      if (abs_x > max_value)
        return (isneg ? -inf_value : +inf_value);

      if (p == Int(2)) { return  (x * x); }
      if (p == Int(3)) { return ((x * x) * x); }
      if (p == Int(4)) { const Float128 x2 = (x * x); return (x2 * x2); }

      Float128 result = ((p % Int(2)) != Int(0)) ? x : Float128(1);
      Float128 xn     = x;  // binary powers of x

      Int p2 = p;
      while (Int(p2 /= 2) != Int(0)) {
        xn *= xn;  // Square xn for each binary power

        const bool has_binary_power = (Int(p2 % Int(2)) != Int(0));
        if (has_binary_power)
          result *= xn;
      }

      return result;
    }


  } // end namespace Impl

  template <>
  struct IsNumber<Impl::Float128>
      : public std::true_type {};

} // end namespace Dune

namespace std
{
#ifndef NO_STD_NUMERIC_LIMITS_SPECIALIZATION
  template <>
  class numeric_limits<Dune::Impl::Float128>
  {
    using Float128 = Dune::Impl::Float128;
    using float128_t = Dune::Impl::float128_t;

  public:
    static constexpr bool is_specialized = true;
    static constexpr Float128 min() noexcept { return FLT128_MIN; }
    static constexpr Float128 max() noexcept { return FLT128_MAX; }
    static constexpr Float128 lowest() noexcept { return -FLT128_MAX; }
    static constexpr int digits = FLT128_MANT_DIG;
    static constexpr int digits10 = 34;
    static constexpr int max_digits10 = 36;
    static constexpr bool is_signed = true;
    static constexpr bool is_integer = false;
    static constexpr bool is_exact = false;
    static constexpr int radix = 2;
    static constexpr Float128 epsilon() noexcept { return FLT128_EPSILON; }
    static constexpr Float128 round_error() noexcept { return float128_t{0.5}; }
    static constexpr int min_exponent = FLT128_MIN_EXP;
    static constexpr int min_exponent10 = FLT128_MIN_10_EXP;
    static constexpr int max_exponent = FLT128_MAX_EXP;
    static constexpr int max_exponent10 = FLT128_MAX_10_EXP;
    static constexpr bool has_infinity = true;
    static constexpr bool has_quiet_NaN = true;
    static constexpr bool has_signaling_NaN = false;
    static constexpr float_denorm_style has_denorm = denorm_present;
    static constexpr bool has_denorm_loss = false;
    static constexpr Float128 infinity() noexcept { return float128_t{1}/float128_t{0}; }
    static Float128 quiet_NaN() noexcept { return nanq(""); }
    static constexpr Float128 signaling_NaN() noexcept { return float128_t{}; }
    static constexpr Float128 denorm_min() noexcept { return FLT128_DENORM_MIN; }
    static constexpr bool is_iec559 = true;
    static constexpr bool is_bounded = false;
    static constexpr bool is_modulo = false;
    static constexpr bool traps = false;
    static constexpr bool tinyness_before = false;
    static constexpr float_round_style round_style = round_to_nearest;
  };
#endif
} // end namespace std

#endif // HAVE_QUADMATH
#endif // DUNE_QUADMATH_HH