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              : 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              : }