Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <blaze/math/DenseVector.h>
       7             : #include <blaze/math/constraints/SIMDPack.h>
       8             : #include <blaze/math/simd/BasicTypes.h>
       9             : #include <blaze/system/Inline.h>
      10             : #include <blaze/system/Vectorization.h>
      11             : 
      12             : #include "Utilities/Math.hpp"
      13             : 
      14           0 : namespace blaze {
      15             : 
      16             : namespace detail {
      17             : template <typename T>
      18             : BLAZE_ALWAYS_INLINE SIMDdouble integer_pow_impl(const SIMDf64<T>& x,
      19             :                                                 const int e) {
      20             :   ASSERT(e >= 0, "Negative powers are not implemented");
      21             :   int ecount = e;
      22             :   int bitcount = 1;
      23             :   while (ecount >>= 1) {
      24             :     ++bitcount;
      25             :   }
      26             :   SIMDdouble result = blaze::set(1.0);
      27             :   while (bitcount) {
      28             :     result *= result;
      29             :     if ((e >> --bitcount) & 0x1) {
      30             :       result *= x;
      31             :     }
      32             :   }
      33             :   return result;
      34             : }
      35             : }  // namespace detail
      36             : // This vectorized implementation of the integer pow function is necessary
      37             : // because blaze does not offer its own version of a integer pow function.
      38             : template <typename T>
      39           0 : BLAZE_ALWAYS_INLINE SIMDdouble integer_pow(const SIMDf64<T>& b, const int e) {
      40             :   switch (e) {
      41             :     case 0:
      42             :       return blaze::set(1.0);
      43             :     case 1:
      44             :       return b;
      45             :     case 2:
      46             :       return b * b;
      47             :     case 3:
      48             :       return b * b * b;
      49             :     case 4: {
      50             :       const SIMDdouble b2 = b * b;
      51             :       return b2 * b2;
      52             :     }
      53             :     case 5: {
      54             :       const SIMDdouble b2 = b * b;
      55             :       return b2 * b2 * b;
      56             :     }
      57             :     case 6: {
      58             :       const SIMDdouble b2 = b * b;
      59             :       return b2 * b2 * b2;
      60             :     }
      61             :     case 7: {
      62             :       const SIMDdouble b2 = b * b;
      63             :       return b2 * b2 * b2 * b;
      64             :     }
      65             :     case 8: {
      66             :       const SIMDdouble b2 = b * b;
      67             :       const SIMDdouble b4 = b2 * b2;
      68             :       return b4 * b4;
      69             :     }
      70             :     default:
      71             :       return detail::integer_pow_impl(b, e);
      72             :   }
      73             : }
      74             : 
      75           0 : struct IntegerPow {
      76           0 :   int exponent;
      77           0 :   explicit inline IntegerPow(const int e) : exponent(e) {}
      78             : 
      79           0 :   BLAZE_ALWAYS_INLINE double operator()(const double a) const {
      80             :     return ::integer_pow(a, exponent);
      81             :   }
      82             : 
      83             :   template <typename T>
      84           0 :   BLAZE_ALWAYS_INLINE decltype(auto) load(const T& a) const {
      85             :     BLAZE_CONSTRAINT_MUST_BE_SIMD_PACK(T);
      86             :     return integer_pow(a, exponent);
      87             :   }
      88             : };
      89             : }  // namespace blaze
      90             : 
      91             : template <typename VT, bool TF>
      92           0 : BLAZE_ALWAYS_INLINE decltype(auto) integer_pow(
      93             :     const blaze::DenseVector<VT, TF>& vec, const int e) {
      94             :   return map(*vec, blaze::IntegerPow{e});
      95             : }
      96             : 
      97             : template <typename VT, bool TF>
      98           0 : BLAZE_ALWAYS_INLINE decltype(auto) IntegerPow(
      99             :     const blaze::DenseVector<VT, TF>& vec, const int e) {
     100             :   return map(*vec, blaze::IntegerPow{e});
     101             : }