GenericPacketMathFunctions.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2007 Julien Pommier
// Copyright (C) 2014 Pedro Gonnet (pedro.gonnet@gmail.com)
// Copyright (C) 2009-2019 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2018-2025 Rasmus Munk Larsen <rmlarsen@gmail.com>
//
// 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 exp and log functions of this file initially come from
 * Julien Pommier's sse math library: http://gruntthepeon.free.fr/ssemath/
 */
 
#ifndef EIGEN_ARCH_GENERIC_PACKET_MATH_FUNCTIONS_H
#define EIGEN_ARCH_GENERIC_PACKET_MATH_FUNCTIONS_H
 
// IWYU pragma: private
#include "../../InternalHeaderCheck.h"
#include "GenericPacketMathPolynomials.h"
#include "GenericPacketMathFrexpLdexp.h"
#include "GenericPacketMathDoubleWord.h"
 
namespace Eigen {
namespace internal {
 
//----------------------------------------------------------------------
// Exponential and Logarithmic Functions
//----------------------------------------------------------------------
 
// Core range reduction and polynomial evaluation for float logarithm.
//
// Given a positive float value v (may be denormal), decomposes it as
// v = 2^e * (1+f) with f in [sqrt(0.5)-1, sqrt(2)-1], then evaluates
// log(1+f) ≈ f + f^2 * P(f) using a degree-7 minimax polynomial.
//
// Returns the approximation of log(v_mantissa) in log_mantissa and the
// integer exponent in e. The caller combines these as appropriate
// (e.g. e*ln2 + log_mantissa for natural log, or log_mantissa*log2e + e
// for log2).
//
// Range reduction uses integer bit manipulation (musl-inspired) instead of the
// heavier pfrexp_generic, saving ~12 ops. The minimax polynomial was found via
// Sollya's fpminimax, giving faithfully-rounded results (max 1 ULP for log).
template <typename Packet>
EIGEN_STRONG_INLINE void plog_core_float(const Packet v, Packet& log_mantissa, Packet& e) {
  using PacketI = typename unpacket_traits<Packet>::integer_packet;
 
  const PacketI cst_min_normal = pset1<PacketI>(0x00800000);
  const PacketI cst_mant_mask = pset1<PacketI>(0x007fffff);
  // Adding this offset to the integer representation biases the exponent so
  // that values near 1 (0x3f800000) map to exponent 0, and values below
  // sqrt(0.5) get folded into the previous exponent.  The magic constant is
  // 0x3f800000 - 0x3f3504f3 = 0x004afb0d, where 0x3f3504f3 ≈ sqrt(0.5).
  const PacketI cst_sqrt_half_offset = pset1<PacketI>(0x004afb0d);
  const PacketI cst_exp_bias = pset1<PacketI>(0x7f);         // 127
  const PacketI cst_half_mant = pset1<PacketI>(0x3f3504f3);  // sqrt(0.5)
 
  // Normalize denormals by multiplying by 2^23.
  PacketI vi = preinterpret<PacketI>(v);
  PacketI is_denormal = pcmp_lt(vi, cst_min_normal);
  Packet v_normalized = pmul(v, pset1<Packet>(8388608.0f));  // 2^23
  vi = pselect(is_denormal, preinterpret<PacketI>(v_normalized), vi);
  // Denormal exponent adjustment: subtract 23 from exponent.
  PacketI denorm_adj = pand(is_denormal, pset1<PacketI>(23));
 
  // Combined range reduction: bias integer representation so that exponent
  // extraction automatically shifts mantissa to [sqrt(0.5), sqrt(2)).
  PacketI vi_biased = padd(vi, cst_sqrt_half_offset);
  // Extract exponent as integer, subtract bias and denormal adjustment.
  PacketI e_int = psub(psub(plogical_shift_right<23>(vi_biased), cst_exp_bias), denorm_adj);
  e = pcast<PacketI, Packet>(e_int);
  // Reconstruct mantissa in [sqrt(0.5), sqrt(2)). The integer addition of the
  // masked mantissa with 0x3f3504f3 (sqrt(0.5)) naturally produces carry into
  // the exponent field, yielding values in [sqrt(0.5), 1) or [1, sqrt(2)).
  // Then subtract 1 to center on 0 → f in [sqrt(0.5)-1, sqrt(2)-1].
  Packet f = psub(preinterpret<Packet>(padd(pand(vi_biased, cst_mant_mask), cst_half_mant)), pset1<Packet>(1.0f));
 
  // Minimax degree-7 polynomial for g(f) = (log(1+f) - f) / f^2 on
  // [sqrt(0.5)-1, sqrt(2)-1], so log(1+f) ≈ f + f^2 * P(f).
  // Generated by Sollya: fpminimax(g, 7, [|single...|], [lo;hi])
  // Mathematical approximation error: max |log(1+f) - (f + f^2*P(f))| < 2.04e-8.
  // Coefficients stored in reverse order for ppolevl (highest degree first).
  constexpr float coeffs[] = {
      8.8758550584316254e-02f,   //  c7 (x^7)
      -1.4199858903884888e-01f,  //  c6 (x^6)
      1.4824025332927704e-01f,   //  c5 (x^5)
      -1.6583317518234253e-01f,  //  c4 (x^4)
      1.9972395896911621e-01f,   //  c3 (x^3)
      -2.5001299381256104e-01f,  //  c2 (x^2)
      3.3333668112754822e-01f,   //  c1 (x^1)
      -4.9999997019767761e-01f,  //  c0 (x^0)
  };
 
  // Evaluate P(f) via Horner's method, then log(1+f) ≈ f + f^2 * P(f).
  Packet f2 = pmul(f, f);
  Packet p = ppolevl<Packet, 7>::run(f, coeffs);
  log_mantissa = pmadd(p, f2, f);
}
 
// Natural or base-2 logarithm for float packets.
//
// Computes log(x) as e*C + log(m), where x = 2^e * m with m in [sqrt(1/2), sqrt(2))
// and C = ln(2) for natural log, C = 1 for log2.
template <typename Packet, bool base2>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_impl_float(const Packet _x) {
  Packet log_mantissa, e;
  plog_core_float(_x, log_mantissa, e);
 
  // Add the logarithm of the exponent back to the result.
  Packet x;
  if (base2) {
    const Packet cst_log2e = pset1<Packet>(static_cast<float>(EIGEN_LOG2E));
    x = pmadd(log_mantissa, cst_log2e, e);
  } else {
    const Packet cst_ln2 = pset1<Packet>(static_cast<float>(EIGEN_LN2));
    x = pmadd(e, cst_ln2, log_mantissa);
  }
 
  // Filter out invalid inputs:
  //  - negative arg → NAN
  //  - 0 → -INF
  //  - +INF → +INF
  const Packet cst_minus_inf = por(psignmask<Packet>(), pinf<Packet>());
  const Packet cst_pos_inf = pinf<Packet>();
  Packet invalid_mask = pcmp_lt_or_nan(_x, pzero(_x));
  Packet iszero_mask = pcmp_eq(_x, pzero(_x));
  Packet pos_inf_mask = pcmp_eq(_x, cst_pos_inf);
  return pselect(iszero_mask, cst_minus_inf, por(pselect(pos_inf_mask, cst_pos_inf, x), invalid_mask));
}
 
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_float(const Packet _x) {
  return plog_impl_float<Packet, /* base2 */ false>(_x);
}
 
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog2_float(const Packet _x) {
  return plog_impl_float<Packet, /* base2 */ true>(_x);
}
 
// -----------------------------------------------------------------------
// Double logarithm: shared polynomial + two range-reduction backends
// -----------------------------------------------------------------------
 
// Cephes rational-polynomial approximation of log(1+f) for
// f in [sqrt(0.5)-1, sqrt(2)-1].
// Evaluates x - 0.5*x^2 + x^3 * P(x)/Q(x) where P and Q are degree-5.
// See: http://www.netlib.org/cephes/
template <typename Packet>
EIGEN_STRONG_INLINE Packet plog_mantissa_double(const Packet x) {
  const Packet cst_cephes_log_p0 = pset1<Packet>(1.01875663804580931796E-4);
  const Packet cst_cephes_log_p1 = pset1<Packet>(4.97494994976747001425E-1);
  const Packet cst_cephes_log_p2 = pset1<Packet>(4.70579119878881725854E0);
  const Packet cst_cephes_log_p3 = pset1<Packet>(1.44989225341610930846E1);
  const Packet cst_cephes_log_p4 = pset1<Packet>(1.79368678507819816313E1);
  const Packet cst_cephes_log_p5 = pset1<Packet>(7.70838733755885391666E0);
  // Q0 = 1.0; pmadd(1, x, q1) simplifies to padd(x, q1).
  const Packet cst_cephes_log_q1 = pset1<Packet>(1.12873587189167450590E1);
  const Packet cst_cephes_log_q2 = pset1<Packet>(4.52279145837532221105E1);
  const Packet cst_cephes_log_q3 = pset1<Packet>(8.29875266912776603211E1);
  const Packet cst_cephes_log_q4 = pset1<Packet>(7.11544750618563894466E1);
  const Packet cst_cephes_log_q5 = pset1<Packet>(2.31251620126765340583E1);
 
  Packet x2 = pmul(x, x);
  Packet x3 = pmul(x2, x);
 
  // Evaluate P and Q simultaneously for better ILP.
  Packet y, y1, y_;
  y = pmadd(cst_cephes_log_p0, x, cst_cephes_log_p1);
  y1 = pmadd(cst_cephes_log_p3, x, cst_cephes_log_p4);
  y = pmadd(y, x, cst_cephes_log_p2);
  y1 = pmadd(y1, x, cst_cephes_log_p5);
  y_ = pmadd(y, x3, y1);
 
  y = padd(x, cst_cephes_log_q1);
  y1 = pmadd(cst_cephes_log_q3, x, cst_cephes_log_q4);
  y = pmadd(y, x, cst_cephes_log_q2);
  y1 = pmadd(y1, x, cst_cephes_log_q5);
  y = pmadd(y, x3, y1);
 
  y_ = pmul(y_, x3);
  y = pdiv(y_, y);
  y = pnmadd(pset1<Packet>(0.5), x2, y);
  return padd(x, y);
}
 
// Detect whether unpacket_traits<Packet>::integer_packet is defined.
template <typename Packet, typename = void>
struct packet_has_integer_packet : std::false_type {};
template <typename Packet>
struct packet_has_integer_packet<Packet, void_t<typename unpacket_traits<Packet>::integer_packet>> : std::true_type {};
 
// Dispatch struct for double-precision range reduction.
// Primary template: pfrexp-based fallback (used when integer_packet is absent).
template <typename Packet, bool UseIntegerPacket>
struct plog_range_reduce_double {
  EIGEN_STRONG_INLINE static void run(const Packet v, Packet& f, Packet& e) {
    const Packet one = pset1<Packet>(1.0);
    const Packet cst_cephes_SQRTHF = pset1<Packet>(0.70710678118654752440E0);
    // pfrexp: f in [0.5, 1), e = unbiased exponent as double.
    f = pfrexp(v, e);
    // Shift [0.5,1) -> [sqrt(0.5)-1, sqrt(2)-1] with exponent correction:
    //   if f < sqrt(0.5): f = f + f - 1, e -= 1   (giving f in [0, sqrt(2)-1))
    //   else:             f = f - 1                (giving f in [sqrt(0.5)-1, 0))
    Packet mask = pcmp_lt(f, cst_cephes_SQRTHF);
    Packet tmp = pand(f, mask);
    f = psub(f, one);
    e = psub(e, pand(one, mask));
    f = padd(f, tmp);
  }
};
 
// Specialisation: fast integer-bit-manipulation path (musl-inspired).
// Requires unpacket_traits<Packet>::integer_packet to be a 64-bit integer packet.
template <typename Packet>
struct plog_range_reduce_double<Packet, true> {
  EIGEN_STRONG_INLINE static void run(const Packet v, Packet& f, Packet& e) {
    using PacketI = typename unpacket_traits<Packet>::integer_packet;
    // 2^-1022: smallest positive normal double.
    const PacketI cst_min_normal = pset1<PacketI>(static_cast<int64_t>(0x0010000000000000LL));
    // Lower 52-bit mask (IEEE mantissa field).
    const PacketI cst_mant_mask = pset1<PacketI>(static_cast<int64_t>(0x000FFFFFFFFFFFFFLL));
    // Offset = 1.0_bits - sqrt(0.5)_bits.  Adding this to the integer
    // representation shifts the exponent field so that the [sqrt(0.5), sqrt(2))
    // half-octave boundary falls on an exact biased-exponent boundary, letting
    // us extract e with a single right shift.  The constant is:
    //   0x3FF0000000000000 - 0x3FE6A09E667F3BCD = 0x00095F619980C433
    const PacketI cst_sqrt_half_offset =
        pset1<PacketI>(static_cast<int64_t>(0x3FF0000000000000LL - 0x3FE6A09E667F3BCDLL));
    // IEEE double exponent bias (1023).
    const PacketI cst_exp_bias = pset1<PacketI>(static_cast<int64_t>(1023));
    // sqrt(0.5) IEEE bits — used to reconstruct f from biased mantissa.
    const PacketI cst_half_mant = pset1<PacketI>(static_cast<int64_t>(0x3FE6A09E667F3BCDLL));
 
    // Reinterpret v as a 64-bit integer vector.
    PacketI vi = preinterpret<PacketI>(v);
 
    // Normalise denormals: multiply by 2^52 and correct the exponent by -52.
    PacketI is_denormal = pcmp_lt(vi, cst_min_normal);
    // 2^52 via bit pattern: biased exponent = 52 + 1023 = 0x433, mantissa = 0.
    Packet v_norm = pmul(v, pset1frombits<Packet>(static_cast<uint64_t>(int64_t(52 + 0x3ff) << 52)));
    vi = pselect(is_denormal, preinterpret<PacketI>(v_norm), vi);
    PacketI denorm_adj = pand(is_denormal, pset1<PacketI>(static_cast<int64_t>(52)));
 
    // Bias the integer representation so the exponent field directly encodes
    // the half-octave index.
    PacketI vi_biased = padd(vi, cst_sqrt_half_offset);
    // Extract unbiased exponent: shift out mantissa bits, subtract IEEE bias
    // and denormal adjustment.
    PacketI e_int = psub(psub(plogical_shift_right<52>(vi_biased), cst_exp_bias), denorm_adj);
    // Convert integer exponent to floating-point.
    e = pcast<PacketI, Packet>(e_int);
 
    // Reconstruct mantissa in [sqrt(0.5), sqrt(2)) via integer arithmetic.
    // The integer addition of the masked mantissa bits and the sqrt(0.5) bit
    // pattern carries into the exponent field, yielding a value in that range.
    // Then subtract 1 to centre on 0: f in [sqrt(0.5)-1, sqrt(2)-1].
    f = psub(preinterpret<Packet>(padd(pand(vi_biased, cst_mant_mask), cst_half_mant)), pset1<Packet>(1.0));
  }
};
 
// Core range reduction and polynomial for double logarithm.
// Input:  v > 0 (zero / negative / inf / nan are handled by the caller).
// Output: log_mantissa ≈ log(mantissa of v in [sqrt(0.5), sqrt(2))),
//         e            = unbiased exponent of v as a double.
// Selects the fast integer path when integer_packet is available, otherwise
// falls back to pfrexp.
template <typename Packet>
EIGEN_STRONG_INLINE void plog_core_double(const Packet v, Packet& log_mantissa, Packet& e) {
  Packet f;
  plog_range_reduce_double<Packet, packet_has_integer_packet<Packet>::value>::run(v, f, e);
  log_mantissa = plog_mantissa_double(f);
}
 
/* Returns the base e (2.718...) or base 2 logarithm of x.
 * The argument is separated into its exponent and fractional parts.
 * The logarithm of the fraction in the interval [sqrt(1/2), sqrt(2)],
 * is approximated by
 *
 *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
 *
 * for more detail see: http://www.netlib.org/cephes/
 */
template <typename Packet, bool base2>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_impl_double(const Packet _x) {
  const Packet cst_minus_inf = por(psignmask<Packet>(), pinf<Packet>());
  const Packet cst_pos_inf = pinf<Packet>();
 
  Packet log_mantissa, e;
  plog_core_double(_x, log_mantissa, e);
 
  // Combine: log(x) = e * ln2 + log(mantissa), or log2(x) = log(mantissa)*log2e + e.
  Packet x;
  if (base2) {
    const Packet cst_log2e = pset1<Packet>(static_cast<double>(EIGEN_LOG2E));
    x = pmadd(log_mantissa, cst_log2e, e);
  } else {
    const Packet cst_ln2 = pset1<Packet>(static_cast<double>(EIGEN_LN2));
    x = pmadd(e, cst_ln2, log_mantissa);
  }
 
  Packet invalid_mask = pcmp_lt_or_nan(_x, pzero(_x));
  Packet iszero_mask = pcmp_eq(_x, pzero(_x));
  Packet pos_inf_mask = pcmp_eq(_x, cst_pos_inf);
  // Filter out invalid inputs:
  //  - negative arg → NAN
  //  - 0            → -INF
  //  - +INF         → +INF
  return pselect(iszero_mask, cst_minus_inf, por(pselect(pos_inf_mask, cst_pos_inf, x), invalid_mask));
}
 
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog_double(const Packet _x) {
  return plog_impl_double<Packet, /* base2 */ false>(_x);
}
 
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog2_double(const Packet _x) {
  return plog_impl_double<Packet, /* base2 */ true>(_x);
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_log1p_float(const Packet& x) {
  const Packet one = pset1<Packet>(1.0f);
  const Packet cst_minus_inf = por(psignmask<Packet>(), pinf<Packet>());
  const Packet cst_pos_inf = pinf<Packet>();
 
  // u = 1 + x, with rounding. Recover the lost low bits: dx = x - (u - 1).
  Packet u = padd(one, x);
  Packet dx = psub(x, psub(u, one));
 
  // For |x| tiny enough that u rounds to 1, return x directly.
  Packet small_mask = pcmp_eq(u, one);
  // For u = +inf (x very large), return +inf.
  Packet inf_mask = pcmp_eq(u, cst_pos_inf);
 
  // Core range reduction and polynomial on u.
  Packet log_u, e;
  plog_core_float(u, log_u, e);
 
  // result = e * ln2 + log(u) + dx/u.
  // The dx/u term corrects for the rounding error in u = fl(1+x).
  const Packet cst_ln2 = pset1<Packet>(static_cast<float>(EIGEN_LN2));
  Packet result = pmadd(e, cst_ln2, padd(log_u, pdiv(dx, u)));
 
  // Handle special cases.
  Packet neg_mask = pcmp_lt(u, pzero(u));
  Packet zero_mask = pcmp_eq(x, pset1<Packet>(-1.0f));
  result = pselect(small_mask, x, result);
  result = pselect(inf_mask, cst_pos_inf, result);
  result = pselect(zero_mask, cst_minus_inf, result);
  result = por(neg_mask, result);  // NaN for x < -1
  return result;
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_log1p_double(const Packet& x) {
  const Packet one = pset1<Packet>(1.0);
  const Packet cst_minus_inf = por(psignmask<Packet>(), pinf<Packet>());
  const Packet cst_pos_inf = pinf<Packet>();
 
  // u = 1 + x, with rounding.  Recover the lost low bits: dx = x - (u - 1).
  Packet u = padd(one, x);
  Packet dx = psub(x, psub(u, one));
 
  // For |x| tiny enough that u rounds to 1, return x directly.
  Packet small_mask = pcmp_eq(u, one);
  // For u = +inf (x very large), return +inf.
  Packet inf_mask = pcmp_eq(u, cst_pos_inf);
 
  // Core range reduction and polynomial on u.
  Packet log_u, e;
  plog_core_double(u, log_u, e);
 
  // result = e * ln2 + log(u) + dx/u.
  // The dx/u term corrects for the rounding error in u = fl(1+x).
  const Packet cst_ln2 = pset1<Packet>(static_cast<double>(EIGEN_LN2));
  Packet result = pmadd(e, cst_ln2, padd(log_u, pdiv(dx, u)));
 
  // Handle special cases.
  Packet neg_mask = pcmp_lt(u, pzero(u));
  Packet zero_mask = pcmp_eq(x, pset1<Packet>(-1.0));
  result = pselect(small_mask, x, result);
  result = pselect(inf_mask, cst_pos_inf, result);
  result = pselect(zero_mask, cst_minus_inf, result);
  result = por(neg_mask, result);  // NaN for x < -1
  return result;
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_log1p(const Packet& x) {
  using ScalarType = typename unpacket_traits<Packet>::type;
  const Packet one = pset1<Packet>(ScalarType(1));
  Packet xp1 = padd(x, one);
  Packet small_mask = pcmp_eq(xp1, one);
  Packet log1 = plog(xp1);
  Packet inf_mask = pcmp_eq(xp1, log1);
  Packet log_large = pmul(x, pdiv(log1, psub(xp1, one)));
  return pselect(por(small_mask, inf_mask), x, log_large);
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_expm1(const Packet& x) {
  using ScalarType = typename unpacket_traits<Packet>::type;
  const Packet one = pset1<Packet>(ScalarType(1));
  const Packet neg_one = pset1<Packet>(ScalarType(-1));
  Packet u = pexp(x);
  Packet one_mask = pcmp_eq(u, one);
  Packet u_minus_one = psub(u, one);
  Packet neg_one_mask = pcmp_eq(u_minus_one, neg_one);
  Packet logu = plog(u);
  // The following comparison is to catch the case where
  // exp(x) = +inf. It is written in this way to avoid having
  // to form the constant +inf, which depends on the packet
  // type.
  Packet pos_inf_mask = pcmp_eq(logu, u);
  Packet expm1 = pmul(u_minus_one, pdiv(x, logu));
  expm1 = pselect(pos_inf_mask, u, expm1);
  return pselect(one_mask, x, pselect(neg_one_mask, neg_one, expm1));
}
 
// exp(r) for r in [-ln(2)/2, ln(2)/2], by a 6th order minimax polynomial.
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pexp_float_reduced(const Packet& r) {
  const Packet cst_one = pset1<Packet>(1.0f);
  const Packet cst_p2 = pset1<Packet>(0.49999988079071044921875f);
  const Packet cst_p3 = pset1<Packet>(0.16666518151760101318359375f);
  const Packet cst_p4 = pset1<Packet>(4.166965186595916748046875e-2f);
  const Packet cst_p5 = pset1<Packet>(8.36894474923610687255859375e-3f);
  const Packet cst_p6 = pset1<Packet>(1.37449637986719608306884765625e-3f);
 
  const Packet r2 = pmul(r, r);
  Packet p_even = pmadd(r2, cst_p6, cst_p4);
  const Packet p_odd = pmadd(r2, cst_p5, cst_p3);
  p_even = pmadd(r2, p_even, cst_p2);
  const Packet p_low = padd(r, cst_one);
  const Packet y = pmadd(r, p_odd, p_even);
  return pmadd(r2, y, p_low);
}
 
// Exponential function. Works by writing "x = m*log(2) + r" where
// "m = rint(x/log(2))" and "r" is the remainder. The result is then
// "exp(x) = 2^m*exp(r)" where exp(r) is in the range [-1,1).
// exp(r) is computed using a 6th order minimax polynomial approximation.
template <typename Packet, bool IsFinite>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp_float(const Packet _x) {
  using PacketI = typename unpacket_traits<Packet>::integer_packet;
 
  const Packet cst_exp_hi = pset1<Packet>(88.723f);
  const Packet cst_exp_lo = pset1<Packet>(-104.f);
 
  const Packet cst_cephes_LOG2EF = pset1<Packet>(1.44269504088896341f);
 
  // Clamp x to prevent overflow/underflow.
  Packet x = pmin(pmax(_x, cst_exp_lo), cst_exp_hi);
 
  // Express exp(x) as exp(m*ln(2) + r), start by extracting
  // m = rint(x/ln(2)).
  Packet m = print(pmul(x, cst_cephes_LOG2EF));
 
  // Get r = x - m*ln(2). m*ln(2) is subtracted out in two parts,
  // m*C1+m*C2 = m*ln(2), to avoid accumulating truncation errors.
  const Packet cst_cephes_exp_C1 = pset1<Packet>(-0.693359375f);
  const Packet cst_cephes_exp_C2 = pset1<Packet>(2.12194440e-4f);
  Packet r = pmadd(m, cst_cephes_exp_C1, x);
  r = pmadd(m, cst_cephes_exp_C2, r);
 
  Packet y = pexp_float_reduced(r);
 
  // Construct the result y * 2^m via a 2-way exponent split. Writing
  //   2^m = 2^floor(m/2) * 2^(m - floor(m/2))
  // keeps each constructed power-of-two within the normal float range (since
  // |m| <= 150 after the input clamp implies each half lies in [-75, 65]), so
  // the IEEE-correct rounding of normal*normal handles the subnormal output
  // range without special-casing.
  //
  // To minimize integer ops we fold the bias into the split:
  //   biased_sum = m + 254                      (always nonnegative)
  //   biased_hi  = biased_sum >> 1  =  m/2 + 127
  //   biased_lo  = biased_sum - biased_hi = (m - m/2) + 127
  // so each <<23 directly yields the float bit pattern for 2^(m/2) and
  // 2^(m - m/2).
  const PacketI cst_double_bias = pset1<PacketI>(254);
  const PacketI mi = pcast<Packet, PacketI>(m);
  const PacketI biased_sum = padd(mi, cst_double_bias);
  const PacketI biased_hi = plogical_shift_right<1>(biased_sum);  // m/2 + 127
  const PacketI biased_lo = psub(biased_sum, biased_hi);          // (m - m/2) + 127
  const Packet pow2_hi = preinterpret<Packet>(plogical_shift_left<23>(biased_hi));
  const Packet pow2_lo = preinterpret<Packet>(plogical_shift_left<23>(biased_lo));
  y = pmul(pmul(y, pow2_hi), pow2_lo);
 
  EIGEN_IF_CONSTEXPR (!IsFinite) {
    // Handle NaN: exp(nan) = nan. Use pmax to propagate NaN from input.
    y = pmax(y, _x);
  }
  return y;
}
 
// exp(g) for g in [-ln(2)/2, ln(2)/2], by the rational interpolant exp(g) = 1 + 2 g P(g^2) / (Q(g^2) - g P(g^2)).
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet pexp_double_reduced(const Packet& g) {
  const Packet cst_1 = pset1<Packet>(1.0);
  const Packet cst_2 = pset1<Packet>(2.0);
  const Packet cst_cephes_exp_p0 = pset1<Packet>(1.26177193074810590878e-4);
  const Packet cst_cephes_exp_p1 = pset1<Packet>(3.02994407707441961300e-2);
  const Packet cst_cephes_exp_p2 = pset1<Packet>(9.99999999999999999910e-1);
  const Packet cst_cephes_exp_q0 = pset1<Packet>(3.00198505138664455042e-6);
  const Packet cst_cephes_exp_q1 = pset1<Packet>(2.52448340349684104192e-3);
  const Packet cst_cephes_exp_q2 = pset1<Packet>(2.27265548208155028766e-1);
  const Packet cst_cephes_exp_q3 = pset1<Packet>(2.00000000000000000009e0);
 
  const Packet x2 = pmul(g, g);
 
  // Evaluate the numerator polynomial of the rational interpolant.
  Packet px = cst_cephes_exp_p0;
  px = pmadd(px, x2, cst_cephes_exp_p1);
  px = pmadd(px, x2, cst_cephes_exp_p2);
  px = pmul(px, g);
 
  // Evaluate the denominator polynomial of the rational interpolant.
  Packet qx = cst_cephes_exp_q0;
  qx = pmadd(qx, x2, cst_cephes_exp_q1);
  qx = pmadd(qx, x2, cst_cephes_exp_q2);
  qx = pmadd(qx, x2, cst_cephes_exp_q3);
 
  // exp(g) = 1 + 2*px/(qx - px).
  return pmadd(cst_2, pdiv(px, psub(qx, px)), cst_1);
}
 
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp_double(const Packet _x) {
  const Packet cst_exp_hi = pset1<Packet>(709.784);
  const Packet cst_exp_lo = pset1<Packet>(-745.519);
  const Packet cst_pldexp_threshold = pset1<Packet>(708.0);
  const Packet cst_cephes_LOG2EF = pset1<Packet>(1.4426950408889634073599);
  const Packet cst_cephes_exp_C1 = pset1<Packet>(0.693145751953125);
  const Packet cst_cephes_exp_C2 = pset1<Packet>(1.42860682030941723212e-6);
 
  // Clamp x.
  Packet zero_mask = pcmp_lt(_x, cst_exp_lo);
  Packet x = pmin(_x, cst_exp_hi);
 
  // Express exp(x) as exp(g + n*log(2)).
  // n = rint(x / ln(2)).
  Packet fx = print(pmul(x, cst_cephes_LOG2EF));
 
  // Get the remainder modulo log(2), i.e. the "g" described above. Subtract
  // n*log(2) out in two steps, i.e. n*C1 + n*C2, C1+C2=log2 to get the last
  // digits right.
  x = pnmadd(fx, cst_cephes_exp_C1, x);
  x = pnmadd(fx, cst_cephes_exp_C2, x);
 
  x = pexp_double_reduced(x);
 
  // Construct the result 2^n * exp(g) = e * x. The max is used to catch
  // non-finite values in the input.
  const Packet fast_pldexp_unsafe = pcmp_lt(cst_pldexp_threshold, pabs(_x));
  if (!predux_any(fast_pldexp_unsafe)) {
    // For |x| <= 708, we know the result is not zero or inf, and we can safely use
    // the fast version of pldexp.
    return pmax(pldexp_fast(x, fx), _x);
  }
  return pandnot(pmax(pldexp(x, fx), _x), zero_mask);
}
 
template <typename Packet, typename Scalar = typename unpacket_traits<Packet>::type>
struct pexp_reduced_impl;
template <typename Packet>
struct pexp_reduced_impl<Packet, float> {
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& r) { return pexp_float_reduced(r); }
};
template <typename Packet>
struct pexp_reduced_impl<Packet, double> {
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet& r) { return pexp_double_reduced(r); }
};
 
// 2^f for |f| <= 1/2: f ln(2) = r_hi + r_lo with ln(2) = ln2_hi + ln2_lo, and exp(r) = exp(r_hi) (1 + r_lo) +
// O(r_lo^2).
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_exp2_reduced(const Packet& f) {
  using Scalar = typename unpacket_traits<Packet>::type;
  // ln(2) - ln2_hi, from (double(ln(2)) - ln2_hi) exact in double plus ln(2) - double(ln(2)).
  const Scalar ln2_hi = Scalar(EIGEN_LN2);
  const Scalar ln2_lo = Scalar((static_cast<double>(EIGEN_LN2) - double(ln2_hi)) + 2.3190468138462996e-17);
  Packet r_hi, r_lo;
  twoprod(pset1<Packet>(ln2_hi), f, r_hi, r_lo);
  r_lo = pmadd(pset1<Packet>(ln2_lo), f, r_lo);
  const Packet e = pexp_reduced_impl<Packet>::run(r_hi);
  return pmadd(e, r_lo, e);
}
 
// exp2(x) = 2^n 2^f with n = rint(x) and f = x - n, which is exact with |f| <= 1/2. So exp2(n) = 2^n exactly, and
// neither the reduction nor the scaling adds an error that grows with |x|.
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_exp2(const Packet& _x) {
  using Scalar = typename unpacket_traits<Packet>::type;
  constexpr int max_exponent = std::numeric_limits<Scalar>::max_exponent;
  constexpr int digits = std::numeric_limits<Scalar>::digits;
  // Beyond these, the result is inf or rounds to 0; within them, n / 2 and n - n / 2 are normal exponents.
  constexpr Scalar max_cap = Scalar(max_exponent + 1);
  constexpr Scalar min_cap = -Scalar(max_exponent + digits - 1);
  const Packet x = pmax(pmin(_x, pset1<Packet>(max_cap)), pset1<Packet>(min_cap));
  const Packet n = print(x);
  const Packet y = generic_exp2_reduced(psub(x, n));
  // One pldexp_fast while 2^n is normal. Near the ends of the range, y 2^n = (y 2^h) 2^(n - h) with h = rint(n / 2):
  // both powers are normal and the first product is exact, so only the last one rounds, also into the subnormal range.
  // The barrier keeps -ffast-math from forming 2^h 2^(n - h), which overflows for n = max_exponent.
  if (!predux_any(pcmp_lt(pset1<Packet>(Scalar(max_exponent - 2)), pabs(n)))) return pldexp_fast(y, n);
  const Packet h = print(pmul(n, pset1<Packet>(Scalar(0.5))));
  Packet out = pldexp_fast(y, h);
  EIGEN_OPTIMIZATION_BARRIER(out)
  return pldexp_fast(out, psub(n, h));
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_float(const Packet& x) {
  // log10(e) in higher precision, split into hi+lo so log_x*(hi+lo) is reconstructed via FMA.
  // hi = round-to-nearest-float of log10(e); lo = float(log10(e) - hi).
  const Packet cst_log10e_hi = pset1<Packet>(0.4342944920063018f);
  const Packet cst_log10e_lo = pset1<Packet>(-1.0103049952192578e-08f);
  const Packet cst_inf = pinf<Packet>();
  const Packet cst_zero = pzero(x);
 
  const Packet log_x = plog(x);
  const Packet finite_mask = pcmp_lt(pabs(log_x), cst_inf);
  const Packet finite_log_x = pselect(finite_mask, log_x, cst_zero);
  const Packet split_log10_x = pmadd(finite_log_x, cst_log10e_hi, pmul(finite_log_x, cst_log10e_lo));
  return pselect(finite_mask, split_log10_x, log_x);
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_double(const Packet& x) {
  const Packet cst_log10e = pset1<Packet>(0.43429448190325182);
  return pmul(plog(x), cst_log10e);
}
 
}  // end namespace internal
}  // end namespace Eigen
 
// Include the split-out sections. Order matters: Pow depends on exp/log and FrexpLdexp,
// Trig depends on exp (for ptanh_float), Complex depends on Trig (for psincos_selector).
#include "GenericPacketMathPow.h"
#include "GenericPacketMathTrig.h"
#include "GenericPacketMathComplex.h"
 
namespace Eigen {
namespace internal {
 
//----------------------------------------------------------------------
// Sign Function
//----------------------------------------------------------------------
 
template <typename Packet, typename Scalar = typename unpacket_traits<Packet>::type>
struct psign_uses_binary_encoding
    : bool_constant<!is_scalar<Packet>::value && packet_has_integer_packet<Packet>::value &&
                    (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value)> {};
 
template <typename Packet>
struct psign_impl<Packet, std::enable_if_t<psign_uses_binary_encoding<Packet>::value>> {
  static EIGEN_DEVICE_FUNC inline Packet run(const Packet& a) {
    using Scalar = typename unpacket_traits<Packet>::type;
    using PacketI = typename unpacket_traits<Packet>::integer_packet;
    using IntScalar = typename unpacket_traits<PacketI>::type;
    using Binary = binary_floating_point_traits<Scalar>;
 
    // Magnitude is a nonnegative integer. Comparing it with 1 avoids Clang
    // folding an integer == 0 test into a DAZ-sensitive FP comparison.
    const PacketI magnitude = pand(preinterpret<PacketI>(a), pset1<PacketI>(IntScalar(~Binary::kSignBit)));
    const Packet is_zero = preinterpret<Packet>(pcmp_lt(magnitude, pset1<PacketI>(IntScalar(1))));
    const Packet is_nan = preinterpret<Packet>(pcmp_lt(pset1<PacketI>(IntScalar(Binary::kExponentMask)), magnitude));
    const Packet signed_one = por(pandnot(a, pabs(a)), pset1<Packet>(Scalar(1)));
    return pselect(is_nan, a, pandnot(signed_one, is_zero));
  }
};
 
template <typename Packet>
struct psign_impl<Packet, std::enable_if_t<!psign_uses_binary_encoding<Packet>::value && !is_scalar<Packet>::value &&
                                           !NumTraits<typename unpacket_traits<Packet>::type>::IsComplex &&
                                           !NumTraits<typename unpacket_traits<Packet>::type>::IsInteger>> {
  static EIGEN_DEVICE_FUNC inline Packet run(const Packet& a) {
    using Scalar = typename unpacket_traits<Packet>::type;
    const Packet cst_one = pset1<Packet>(Scalar(1));
 
    const Packet abs_a = pabs(a);
    const Packet sign_mask = pandnot(a, abs_a);
    const Packet zero_mask = pcmp_eq(abs_a, pzero(abs_a));
 
    return pselect(pisnan(a), a, pselect(zero_mask, abs_a, por(sign_mask, cst_one)));
  }
};
 
template <typename Packet>
struct psign_impl<Packet, std::enable_if_t<!is_scalar<Packet>::value &&
                                           !NumTraits<typename unpacket_traits<Packet>::type>::IsComplex &&
                                           NumTraits<typename unpacket_traits<Packet>::type>::IsSigned &&
                                           NumTraits<typename unpacket_traits<Packet>::type>::IsInteger>> {
  static EIGEN_DEVICE_FUNC inline Packet run(const Packet& a) {
    using Scalar = typename unpacket_traits<Packet>::type;
    const Packet cst_one = pset1<Packet>(Scalar(1));
    const Packet cst_minus_one = pset1<Packet>(Scalar(-1));
    const Packet cst_zero = pzero(a);
 
    const Packet positive_mask = pcmp_lt(cst_zero, a);
    const Packet positive = pand(positive_mask, cst_one);
    const Packet negative_mask = pcmp_lt(a, cst_zero);
    const Packet negative = pand(negative_mask, cst_minus_one);
 
    return por(positive, negative);
  }
};
 
template <typename Packet>
struct psign_impl<Packet, std::enable_if_t<!is_scalar<Packet>::value &&
                                           !NumTraits<typename unpacket_traits<Packet>::type>::IsComplex &&
                                           !NumTraits<typename unpacket_traits<Packet>::type>::IsSigned &&
                                           NumTraits<typename unpacket_traits<Packet>::type>::IsInteger>> {
  static EIGEN_DEVICE_FUNC inline Packet run(const Packet& a) {
    using Scalar = typename unpacket_traits<Packet>::type;
    const Packet cst_one = pset1<Packet>(Scalar(1));
    const Packet cst_zero = pzero(a);
 
    const Packet zero_mask = pcmp_eq(cst_zero, a);
    return pandnot(cst_one, zero_mask);
  }
};
 
// \internal \returns the sign of a complex number z, defined as z / abs(z).
template <typename Packet>
struct psign_impl<Packet, std::enable_if_t<!is_scalar<Packet>::value &&
                                           NumTraits<typename unpacket_traits<Packet>::type>::IsComplex &&
                                           unpacket_traits<Packet>::vectorizable>> {
  static EIGEN_DEVICE_FUNC inline Packet run(const Packet& a) {
    using Scalar = typename unpacket_traits<Packet>::type;
    using RealScalar = typename Scalar::value_type;
    using RealPacket = typename unpacket_traits<Packet>::as_real;
 
    // Step 1. Compute (for each element z = x + i*y in a)
    //     l = abs(z) = sqrt(x^2 + y^2).
    // To avoid over- and underflow, we use the stable formula for each hypotenuse
    //    l = (zmin == 0 ? zmax : zmax * sqrt(1 + (zmin/zmax)**2)),
    // where zmax = max(|x|, |y|), zmin = min(|x|, |y|),
    RealPacket a_abs = pabs(a.v);
    RealPacket a_abs_flip = pcplxflip(Packet(a_abs)).v;
    RealPacket a_max = pmax(a_abs, a_abs_flip);
    RealPacket a_min = pmin(a_abs, a_abs_flip);
    RealPacket a_min_zero_mask = pcmp_eq(a_min, pzero(a_min));
    RealPacket a_max_zero_mask = pcmp_eq(a_max, pzero(a_max));
    RealPacket r = pdiv(a_min, a_max);
    const RealPacket cst_one = pset1<RealPacket>(RealScalar(1));
    RealPacket l = pmul(a_max, psqrt(padd(cst_one, pmul(r, r))));  // [l0, l0, l1, l1]
    // Set l to a_max if a_min is zero, since the roundtrip sqrt(a_max^2) may be
    // lossy.
    l = pselect(a_min_zero_mask, a_max, l);
    // Step 2 compute a / abs(a).
    RealPacket sign_as_real = pandnot(pdiv(a.v, l), a_max_zero_mask);
    Packet sign;
    sign.v = sign_as_real;
    return sign;
  }
};
 
//----------------------------------------------------------------------
// Rounding Functions
//----------------------------------------------------------------------
 
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_rint(const Packet& a) {
  using Scalar = typename unpacket_traits<Packet>::type;
  using IntType = typename numext::get_integer_by_size<sizeof(Scalar)>::signed_type;
  // Adds and subtracts signum(a) * 2^kMantissaBits to force rounding.
  const IntType kLimit = IntType(1) << (NumTraits<Scalar>::digits() - 1);
  const Packet cst_limit = pset1<Packet>(static_cast<Scalar>(kLimit));
  Packet abs_a = pabs(a);
  Packet sign_a = pandnot(a, abs_a);
  Packet rint_a = padd(abs_a, cst_limit);
  // Don't compile-away addition and subtraction.
  EIGEN_OPTIMIZATION_BARRIER(rint_a);
  rint_a = psub(rint_a, cst_limit);
  rint_a = por(rint_a, sign_a);
  // If greater than limit (or NaN), simply return a.
  Packet mask = pcmp_lt(abs_a, cst_limit);
  Packet result = pselect(mask, rint_a, a);
  return result;
}
 
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_floor(const Packet& a) {
  using Scalar = typename unpacket_traits<Packet>::type;
  const Packet cst_1 = pset1<Packet>(Scalar(1));
  Packet rint_a = generic_rint(a);
  // if a < rint(a), then rint(a) == ceil(a)
  Packet mask = pcmp_lt(a, rint_a);
  Packet offset = pand(cst_1, mask);
  Packet result = psub(rint_a, offset);
  return result;
}
 
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_ceil(const Packet& a) {
  using Scalar = typename unpacket_traits<Packet>::type;
  const Packet cst_1 = pset1<Packet>(Scalar(1));
  const Packet sign_mask = psignmask<Packet>();
  Packet rint_a = generic_rint(a);
  // if rint(a) < a, then rint(a) == floor(a)
  Packet mask = pcmp_lt(rint_a, a);
  Packet offset = pand(cst_1, mask);
  Packet result = padd(rint_a, offset);
  // Signed zero must remain signed (e.g. ceil(-0.02) == -0).
  result = por(result, pand(sign_mask, a));
  return result;
}
 
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_trunc(const Packet& a) {
  Packet abs_a = pabs(a);
  Packet sign_a = pandnot(a, abs_a);
  Packet floor_abs_a = generic_floor(abs_a);
  Packet result = por(floor_abs_a, sign_a);
  return result;
}
 
template <typename Packet>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_round(const Packet& a) {
  using Scalar = typename unpacket_traits<Packet>::type;
  const Packet cst_half = pset1<Packet>(Scalar(0.5));
  const Packet cst_1 = pset1<Packet>(Scalar(1));
  Packet abs_a = pabs(a);
  Packet sign_a = pandnot(a, abs_a);
  Packet floor_abs_a = generic_floor(abs_a);
  Packet diff = psub(abs_a, floor_abs_a);
  Packet mask = pcmp_le(cst_half, diff);
  Packet offset = pand(cst_1, mask);
  Packet result = padd(floor_abs_a, offset);
  result = por(result, sign_a);
  return result;
}
 
template <typename Packet>
struct nearest_integer_packetop_impl<Packet, /*IsScalar*/ false, /*IsInteger*/ false> {
  using Scalar = typename unpacket_traits<Packet>::type;
  static_assert(packet_traits<Scalar>::HasRound, "Generic nearest integer functions are disabled for this type.");
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_floor(const Packet& x) { return generic_floor(x); }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_ceil(const Packet& x) { return generic_ceil(x); }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_rint(const Packet& x) { return generic_rint(x); }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_round(const Packet& x) { return generic_round(x); }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_trunc(const Packet& x) { return generic_trunc(x); }
};
 
template <typename Packet>
struct nearest_integer_packetop_impl<Packet, /*IsScalar*/ false, /*IsInteger*/ true> {
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_floor(const Packet& x) { return x; }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_ceil(const Packet& x) { return x; }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_rint(const Packet& x) { return x; }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_round(const Packet& x) { return x; }
  static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run_trunc(const Packet& x) { return x; }
};
 
}  // end namespace internal
}  // end namespace Eigen
 
#endif  // EIGEN_ARCH_GENERIC_PACKET_MATH_FUNCTIONS_H