Line data    Source code

       1             : // SPDX-License-Identifier: BSD-3-Clause
       2             : // Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
       3             : /// @file
       4             : /// @author Simon Heybrock
       5             : #pragma once
       6             : 
       7             : #include <algorithm>
       8             : #include <cmath>
       9             : #include <limits>
      10             : #include <type_traits>
      11             : 
      12             : #include "scipp/common/index.h"
      13             : 
      14             : namespace scipp::numeric {
      15             : 
      16       64102 : template <class Range> bool islinspace(const Range &range) {
      17       64102 :   if (scipp::size(range) < 2)
      18           0 :     return false;
      19       64102 :   if (range.back() <= range.front())
      20           2 :     return false;
      21             : 
      22             :   using T = typename Range::value_type;
      23             :   if constexpr (std::is_floating_point_v<T>) {
      24       53754 :     const T delta = (range.back() - range.front()) / (scipp::size(range) - 1);
      25       53754 :     constexpr int32_t ulp = 4;
      26       53754 :     const T epsilon = std::numeric_limits<T>::epsilon() *
      27       53754 :                       (std::abs(range.front()) + std::abs(range.back())) * ulp;
      28       53754 :     return std::adjacent_find(range.begin(), range.end(),
      29     1591446 :                               [epsilon, delta](const auto &a, const auto &b) {
      30      530482 :                                 return std::abs(std::abs(b - a) - delta) >
      31      530482 :                                        epsilon;
      32       53754 :                               }) == range.end();
      33             :   } else {
      34       10346 :     const auto delta = range[1] - range[0];
      35       10346 :     return std::adjacent_find(range.begin(), range.end(),
      36       78121 :                               [delta](const auto &a, const auto &b) {
      37       39357 :                                 return std::abs(b - a) != delta;
      38       10346 :                               }) == range.end();
      39             :   }
      40             : }
      41             : 
      42        4831 : template <class Range> bool isarange(const Range &range) {
      43             :   static_assert(std::is_integral_v<typename Range::value_type>);
      44        4831 :   if (scipp::size(range) < 2)
      45         639 :     return true;
      46        4192 :   if (range.back() <= range.front())
      47           8 :     return false;
      48        4184 :   return std::adjacent_find(range.begin(), range.end(),
      49     3268844 :                             [](const auto &a, const auto &b) {
      50     3268844 :                               return std::abs(b - a) != 1;
      51        4184 :                             }) == range.end();
      52             : }
      53             : 
      54             : // Division like Python's __truediv__
      55      656426 : template <class T, class U> auto true_divide(const T &a, const U &b) {
      56             :   if constexpr (std::is_integral_v<T> && std::is_integral_v<U>)
      57          91 :     return static_cast<double>(a) / static_cast<double>(b);
      58             :   else
      59      656335 :     return a / b;
      60             : }
      61             : 
      62             : // Division like Python's __floordiv__
      63             : template <class T, class U>
      64    48203103 : std::common_type_t<T, U> floor_divide(const T &a, const U &b) {
      65             :   using std::floor;
      66             :   if constexpr (std::is_integral_v<T> && std::is_integral_v<U>)
      67    48203065 :     return b == 0 ? 0 : floor(static_cast<double>(a) / static_cast<double>(b));
      68             :   else
      69          38 :     return floor(a / b);
      70             : }
      71             : 
      72             : // Remainder like Python's __mod__, complementary to floor_divide.
      73    24099583 : template <class T, class U> auto remainder(const T &a, const U &b) {
      74             :   if constexpr (std::is_floating_point_v<T> || std::is_floating_point_v<U>) {
      75          16 :     return b == 0 ? NAN : a - floor_divide(a, b) * b;
      76             :   } else {
      77    24099567 :     return b == 0 ? 0 : a - floor_divide(a, b) * b;
      78             :   }
      79             : }
      80             : 
      81             : namespace {
      82             : template <class B, class E>
      83         660 : constexpr auto integer_pow_pos_exponent(const B &base,
      84             :                                         const E exponent) noexcept {
      85             :   static_assert(std::is_integral_v<std::decay_t<E>>);
      86             : 
      87         660 :   if (exponent == 0)
      88         166 :     return static_cast<B>(1);
      89         494 :   if (exponent == 1)
      90         245 :     return base;
      91             : 
      92         249 :   const auto aux = integer_pow_pos_exponent(base, exponent / 2);
      93         249 :   if (exponent % 2 == 0)
      94         217 :     return aux * aux;
      95          32 :   return base * aux * aux;
      96             : }
      97             : } // namespace
      98             : 
      99             : template <class B, class E>
     100         669 : constexpr auto pow(const B base, const E exponent) noexcept {
     101             :   if constexpr (std::is_integral_v<std::decay_t<E>>) {
     102         411 :     if (exponent >= 0) {
     103         409 :       return integer_pow_pos_exponent(base, exponent);
     104             :     } else {
     105             :       return static_cast<std::decay_t<B>>(1) /
     106           2 :              integer_pow_pos_exponent(base, -exponent);
     107             :     }
     108             :   } else {
     109             :     using std::pow;
     110         258 :     return pow(base, exponent);
     111             :   }
     112             : }
     113             : 
     114   198404288 : template <class T> bool isnan([[maybe_unused]] const T &x) {
     115             :   // This is not completely failsafe, e.g., std::complex<float> would not be
     116             :   // caught, but it should help to avoid the most common issues.
     117             :   static_assert(sizeof(T) <= 8, "Likely missing specialization of isnan");
     118             :   if constexpr (std::is_floating_point_v<std::decay_t<T>>)
     119    22910560 :     return std::isnan(x);
     120             :   else
     121   175493728 :     return false;
     122             : }
     123             : 
     124         356 : template <class T> bool isinf([[maybe_unused]] const T &x) {
     125             :   static_assert(sizeof(T) <= 8, "Likely missing specialization of isinf");
     126             :   if constexpr (std::is_floating_point_v<std::decay_t<T>>)
     127         356 :     return std::isinf(x);
     128             :   else
     129           0 :     return false;
     130             : }
     131             : 
     132         856 : template <class T> bool isfinite([[maybe_unused]] const T &x) {
     133             :   static_assert(sizeof(T) <= 8, "Likely missing specialization of isfinite");
     134             :   if constexpr (std::is_floating_point_v<std::decay_t<T>>)
     135         724 :     return std::isfinite(x);
     136             :   else
     137         132 :     return true;
     138             : }
     139             : 
     140          12 : template <class T> bool signbit([[maybe_unused]] T x) {
     141             :   if constexpr (std::is_floating_point_v<std::decay_t<T>>)
     142          12 :     return std::signbit(x);
     143             :   else
     144           0 :     return std::signbit(double(x));
     145             : }
     146             : 
     147             : } // namespace scipp::numeric