1  // Copyright (C) 2020-2023 Free Software Foundation, Inc.
       2  //
       3  // This file is part of the GNU ISO C++ Library.  This library is free
       4  // software; you can redistribute it and/or modify it under the
       5  // terms of the GNU General Public License as published by the
       6  // Free Software Foundation; either version 3, or (at your option)
       7  // any later version.
       8  //
       9  // This library is distributed in the hope that it will be useful,
      10  // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11  // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      12  // GNU General Public License for more details.
      13  //
      14  // You should have received a copy of the GNU General Public License along
      15  // with this library; see the file COPYING3.  If not see
      16  // <http://www.gnu.org/licenses/>.
      17  
      18  #ifndef ULP_H
      19  #define ULP_H
      20  
      21  #include <cmath>
      22  #include <experimental/simd>
      23  #include <type_traits>
      24  #include <cfenv>
      25  
      26  namespace vir {
      27    namespace test {
      28      template <typename T, typename R = typename T::value_type>
      29        R
      30        value_type_impl(int);
      31  
      32      template <typename T>
      33        T
      34        value_type_impl(float);
      35  
      36      template <typename T>
      37        using value_type_t = decltype(value_type_impl<T>(int()));
      38  
      39      template <typename T>
      40        inline T
      41        ulp_distance(const T& val_, const T& ref_)
      42        {
      43          if constexpr (std::is_floating_point_v<value_type_t<T>>)
      44            {
      45              const int fp_exceptions = std::fetestexcept(FE_ALL_EXCEPT);
      46              T val = val_;
      47              T ref = ref_;
      48  
      49              T diff = T();
      50  
      51              using std::abs;
      52              using std::fpclassify;
      53              using std::frexp;
      54              using std::isnan;
      55              using std::isinf;
      56              using std::ldexp;
      57              using std::max;
      58              using std::experimental::where;
      59              using TT = value_type_t<T>;
      60  
      61              where(ref == 0, val) = abs(val);
      62              where(ref == 0, diff) = 1;
      63              where(ref == 0, ref) = std::__norm_min_v<TT>;
      64              where(isinf(ref) && ref == val, ref)
      65                = 0; // where(val_ == ref_) = 0 below will fix it up
      66  
      67              where(val == 0, ref) = abs(ref);
      68              where(val == 0, diff) += 1;
      69              where(val == 0, val) = std::__norm_min_v<TT>;
      70  
      71              using I = decltype(fpclassify(std::declval<T>()));
      72              I exp = {};
      73              frexp(ref, &exp);
      74              // lower bound for exp must be min_exponent to scale the resulting
      75              // difference from a denormal correctly
      76              exp = max(exp, I(std::__min_exponent_v<TT>));
      77              diff += ldexp(abs(ref - val), std::__digits_v<TT> - exp);
      78              where(val_ == ref_ || (isnan(val_) && isnan(ref_)), diff) = T();
      79              std::feclearexcept(FE_ALL_EXCEPT ^ fp_exceptions);
      80              return diff;
      81            }
      82          else
      83            {
      84              if (val_ > ref_)
      85                return val_ - ref_;
      86              else
      87                return ref_ - val_;
      88            }
      89        }
      90  
      91      template <typename T>
      92        inline T
      93        ulp_distance_signed(const T& _val, const T& _ref)
      94        {
      95          using std::copysign;
      96          return copysign(ulp_distance(_val, _ref), _val - _ref);
      97        }
      98    } // namespace test
      99  } // namespace vir
     100  
     101  #endif // ULP_H