Eigen/src/Core/arch/SSE/MathFunctions.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2007 Julien Pommier
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

/* The sin and cos functions of this file come from
 * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
 */

#ifndef EIGEN_MATH_FUNCTIONS_SSE_H
#define EIGEN_MATH_FUNCTIONS_SSE_H

// IWYU pragma: private
#include "../../InternalHeaderCheck.h"

namespace Eigen {

namespace internal {

EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_FLOAT(Packet4f)
EIGEN_INSTANTIATE_GENERIC_MATH_FUNCS_DOUBLE(Packet2d)

// Notice that for newer processors, it is counterproductive to use Newton
// iteration for square root. In particular, Skylake and Zen2 processors
// have approximately doubled throughput of the _mm_sqrt_ps instruction
// compared to their predecessors.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet4f psqrt<Packet4f>(const Packet4f& x) {
  return _mm_sqrt_ps(x);
}
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet2d psqrt<Packet2d>(const Packet2d& x) {
  return _mm_sqrt_pd(x);
}
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet16b psqrt<Packet16b>(const Packet16b& x) {
  return x;
}

#if EIGEN_FAST_MATH
// Even on Skylake, using Newton iteration is a win for reciprocal square root. rsqrtps reads a subnormal x as zero,
// giving +-inf, so the subnormal lanes, 1 <= bits(|x|) <= 0x007fffff, take 1 / sqrt(x) instead. Adding 0x7f800000
// moves exactly that range above 0x7f800000 and wraps every larger value negative. With SSE this test of x alone is
// faster than the test of x * estimate in prsqrt<Packet8f>, which wins with AVX. Under DAZ the fallback reads a
// subnormal as zero too, and returns the estimate's +-inf.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f prsqrt<Packet4f>(const Packet4f& x) {
  const Packet4f r = generic_rsqrt_newton_step<Packet4f, /*Steps=*/1>::run(x, _mm_rsqrt_ps(x));
  const Packet4i bias = pset1<Packet4i>(0x7f800000);
  const Packet4i abs_bits = pand(Packet4i(_mm_castps_si128(x)), pset1<Packet4i>(0x7fffffff));
  const Packet4f redo = _mm_castsi128_ps(pcmp_lt(bias, padd(abs_bits, bias)));
  return predux_any(redo) ? pselect(redo, pdiv(pset1<Packet4f>(1.0f), psqrt(x)), r) : r;
}

#ifdef EIGEN_VECTORIZE_FMA
// Trying to speed up reciprocal using Newton-Raphson is counterproductive
// unless FMA is available. Without FMA pdiv(pset1<Packet>(Scalar(1),a)) is
// 30% faster.
template <>
EIGEN_STRONG_INLINE Packet4f preciprocal<Packet4f>(const Packet4f& x) {
#ifdef EIGEN_VECTORIZE_AVX
  // generic_reciprocal_newton_step, with its NaN-or-zero test as a single compare: r == 0 or unordered.
  const Packet4f one = pset1<Packet4f>(1.0f);
  const Packet4f r0 = _mm_rcp_ps(x);
  const Packet4f refined = pmadd(r0, pnmadd(x, r0, one), r0);
  const Packet4f redo = _mm_cmp_ps(refined, _mm_setzero_ps(), _CMP_EQ_UQ);
  return predux_any(redo) ? pselect(redo, pdiv(one, x), refined) : refined;
#else
  return generic_reciprocal_newton_step<Packet4f, /*Steps=*/1>::run(x, _mm_rcp_ps(x));
#endif
}
#endif

#endif

}  // end namespace internal

namespace numext {

template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float sqrt(const float& x) {
  return internal::pfirst(internal::Packet4f(_mm_sqrt_ss(_mm_set_ss(x))));
}

template <>
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double sqrt(const double& x) {
#if EIGEN_COMP_GNUC_STRICT
  // This works around a GCC bug generating poor code for _mm_sqrt_pd
  // See https://gitlab.com/libeigen/eigen/commit/8dca9f97e38970
  return internal::pfirst(internal::Packet2d(__builtin_ia32_sqrtsd(_mm_set_sd(x))));
#else
  return internal::pfirst(internal::Packet2d(_mm_sqrt_pd(_mm_set_sd(x))));
#endif
}

}  // namespace numext

}  // end namespace Eigen

#endif  // EIGEN_MATH_FUNCTIONS_SSE_H