text: /***************************************************************************
  * Copyright (c) 2016, Johan Mabille and Sylvain Corlay                     *
  *                                                                          *
  * Distributed under the terms of the BSD 3-Clause License.                 *
  *                                                                          *
  * The full license is in the file LICENSE, distributed with this software. *
  ****************************************************************************/
  
  #ifndef XSIMD_POWER_HPP
  #define XSIMD_POWER_HPP
  
  #include "xsimd_basic_math.hpp"
  #include "xsimd_exponential.hpp"
  #include "xsimd_fp_manipulation.hpp"
  #include "xsimd_fp_sign.hpp"
  #include "xsimd_horner.hpp"
  #include "xsimd_logarithm.hpp"
  #include "xsimd_numerical_constant.hpp"
  
  namespace xsimd
  {
  
      /**
       * Computes the value of the batch \c x raised to the power
       * \c y.
       * @param x batch of floating point values.
       * @param y batch of floating point values.
       * @return \c x raised to the power \c y.
       */
      template <class T, std::size_t N>
      batch<T, N> pow(const batch<T, N>& x, const batch<T, N>& y);
  
      /**
       * Computes the cubic root of the batch \c x.
       * @param x batch of floating point values.
       * @return the cubic root of \c x.
       */
      template <class T, std::size_t N>
      batch<T, N> cbrt(const batch<T, N>& x);
  
      /**
       * Computes the square root of the sum of the squares of the batches
       * \c x, and \c y.
       * @param x batch of floating point values.
       * @param y batch of floating point values.
       * @return the square root of the sum of the squares of \c x and \c y.
       */
      template <class T, std::size_t N>
      batch<T, N> hypot(const batch<T, N>& x, const batch<T, N>& y);
  
      /**********************
       * pow implementation *
       **********************/
  
      /* origin: boost/simd/arch/common/simd/function/pow.hpp*/
      /*
       * ====================================================
       * copyright 2016 NumScale SAS
       *
       * Distributed under the Boost Software License, Version 1.0.
       * (See copy at http://boost.org/LICENSE_1_0.txt)
       * ====================================================
       */
  
      namespace detail
      {
          template <class B>
          struct pow_kernel
          {
              static inline B compute(const B& x, const B& y)
              {
                  using b_type = B;
                  auto negx = x < b_type(0.);
                  b_type z = exp(y * log(abs(x)));
                  z = select(is_odd(y) && negx, -z, z);
                  auto invalid = negx && !(is_flint(y) || isinf(y));
                  return select(invalid, nan<b_type>(), z);
              }
          };
  
        template <class T0, class T1>
        inline T0
        ipow(const T0& t0, const T1& t1)
        {
            static_assert(std::is_integral<T1>::value, "second argument must be an integer");
            T0 a = t0;
            T1 b = t1;
            bool const recip = b < 0;
            T0 r{static_cast<T0>(1)};
            while (1)
            {
                if (b & 1)
                {
                    r *= a;
                }
                b /= 2;
                if (b == 0)
                {
                    break;
                }
                a *= a;
            }
            return recip ? 1 / r : r;
        }
      }
  
      template <class T, std::size_t N>
      inline batch<T, N> pow(const batch<T, N>& x, const batch<T, N>& y)
      {
          return detail::pow_kernel<batch<T, N>>::compute(x, y);
      }
  
      template <class T0, std::size_t N, class T1>
      inline typename std::enable_if<std::is_integral<T1>::value, batch<T0, N>>::type
      pow(const batch<T0, N>& t0, const T1& t1)
      {
          return detail::ipow(t0, t1);
      }
  
      /***********************
       * cbrt implementation *
       ***********************/
  
      namespace detail
      {
          /* origin: boost/simd/arch/common/simd/function/cbrt.hpp */
          /*
           * ====================================================
           * copyright 2016 NumScale SAS
           *
           * Distributed under the Boost Software License, Version 1.0.
           * (See copy at http://boost.org/LICENSE_1_0.txt)
           * ====================================================
           */
  
          template <class B, class T = typename B::value_type>
          struct cbrt_kernel;
  
          template <class B>
          struct cbrt_kernel<B, float>
          {
              static inline B compute(const B& a)
              {
                  B z = abs(a);
  #ifndef XSIMD_NO_DENORMALS
                  auto denormal = z < smallestposval<B>();
                  z = select(denormal, z * twotonmb<B>(), z);
                  B f = select(denormal, twotonmbo3<B>(), B(1.));
  #endif
                  static const B CBRT2 = B(detail::caster32_t(0x3fa14518).f);
                  static const B CBRT4 = B(detail::caster32_t(0x3fcb2ff5).f);
                  static const B CBRT2I = B(detail::caster32_t(0x3f4b2ff5).f);
                  static const B CBRT4I = B(detail::caster32_t(0x3f214518).f);
                  using i_type = as_integer_t<B>;
                  i_type e;
                  B x = frexp(z, e);
                  x = horner<B,
                             0x3ece0609,
                             0x3f91eb77,
                             0xbf745265,
                             0x3f0bf0fe,
                             0xbe09e49a>(x);
                  auto flag = e >= i_type(0);
                  i_type e1 = abs(e);
                  i_type rem = e1;
                  e1 /= i_type(3);
                  rem -= e1 * i_type(3);
                  e = e1 * sign(e);
                  const B cbrt2 = select(bool_cast(flag), CBRT2, CBRT2I);
                  const B cbrt4 = select(bool_cast(flag), CBRT4, CBRT4I);
                  B fact = select(bool_cast(rem == i_type(1)), cbrt2, B(1.));
                  fact = select(bool_cast(rem == i_type(2)), cbrt4, fact);
                  x = ldexp(x * fact, e);
                  x -= (x - z / (x * x)) * B(1.f / 3.f);
  #ifndef XSIMD_NO_DENORMALS
                  x = (x | bitofsign(a)) * f;
  #else
                  x = x | bitofsign(a);
  #endif
  #ifndef XSIMD_NO_INFINITIES
                  return select(a == B(0.) || isinf(a), a, x);
  #else
                  return select(a == B(0.), a, x);
  #endif
              }
          };
  
          template <class B>
          struct cbrt_kernel<B, double>
          {
              static inline B compute(const B& a)
              {
                  B z = abs(a);
  #ifndef XSIMD_NO_DENORMALS
                  auto denormal = z < smallestposval<B>();
                  z = select(denormal, z * twotonmb<B>(), z);
                  B f = select(denormal, twotonmbo3<B>(), B(1.));
  #endif
                  static const B CBRT2 = B(detail::caster64_t(int64_t(0x3ff428a2f98d728b)).f);
                  static const B CBRT4 = B(detail::caster64_t(int64_t(0x3ff965fea53d6e3d)).f);
                  static const B CBRT2I = B(detail::caster64_t(int64_t(0x3fe965fea53d6e3d)).f);
                  static const B CBRT4I = B(detail::caster64_t(int64_t(0x3fe428a2f98d728b)).f);
                  using i_type = as_integer_t<B>;
                  i_type e;
                  B x = frexp(z, e);
                  x = horner<B,
                             0x3fd9c0c12122a4fell,
                             0x3ff23d6ee505873all,
                             0xbfee8a4ca3ba37b8ll,
                             0x3fe17e1fc7e59d58ll,
                             0xbfc13c93386fdff6ll>(x);
                  auto flag = e >= zero<i_type>();
                  i_type e1 = abs(e);
                  i_type rem = e1;
                  e1 /= i_type(3);
                  rem -= e1 * i_type(3);
                  e = e1 * sign(e);
                  const B cbrt2 = select(bool_cast(flag), CBRT2, CBRT2I);
                  const B cbrt4 = select(bool_cast(flag), CBRT4, CBRT4I);
                  B fact = select(bool_cast(rem == i_type(1)), cbrt2, B(1.));
                  fact = select(bool_cast(rem == i_type(2)), cbrt4, fact);
                  x = ldexp(x * fact, e);
                  x -= (x - z / (x * x)) * B(1. / 3.);
                  x -= (x - z / (x * x)) * B(1. / 3.);
  #ifndef XSIMD_NO_DENORMALS
                  x = (x | bitofsign(a)) * f;
  #else
                  x = x | bitofsign(a);
  #endif
  #ifndef XSIMD_NO_INFINITIES
                  return select(a == B(0.) || isinf(a), a, x);
  #else
                  return select(a == B(0.), a, x);
  #endif
              }
          };
      }
  
      template <class T, std::size_t N>
      inline batch<T, N> cbrt(const batch<T, N>& x)
      {
          return detail::cbrt_kernel<batch<T, N>>::compute(x);
      }
  
      /************************
       * hypot implementation *
       ************************/
  
      template <class T, std::size_t N>
      inline batch<T, N> hypot(const batch<T, N>& x, const batch<T, N>& y)
      {
          return sqrt(fma(x, x, y * y));
      }
  }
  
  #endif
  