Eigen/src/Core/arch/Default/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) {
  typedef typename unpacket_traits<Packet>::integer_packet PacketI;

  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 = pset1frombits<Packet>(static_cast<Eigen::numext::uint32_t>(0xff800000u));
  const Packet cst_pos_inf = pset1frombits<Packet>(static_cast<Eigen::numext::uint32_t>(0x7f800000u));
  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) {
    typedef typename unpacket_traits<Packet>::integer_packet PacketI;
    // 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 = pset1frombits<Packet>(static_cast<uint64_t>(0xfff0000000000000ull));
  const Packet cst_pos_inf = pset1frombits<Packet>(static_cast<uint64_t>(0x7ff0000000000000ull));

  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);
}

/** \internal \returns log(1 + x) for single precision float.
    Computes log(1+x) using plog_core_float for the core range reduction
    and polynomial evaluation. The rounding error from forming u = fl(1+x)
    is recovered as dx = x - (u - 1), and folded in as a first-order
    correction dx/u after the polynomial evaluation.
 */
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 = pset1frombits<Packet>(static_cast<Eigen::numext::uint32_t>(0xff800000u));
  const Packet cst_pos_inf = pset1frombits<Packet>(static_cast<Eigen::numext::uint32_t>(0x7f800000u));

  // 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;
}

/** \internal \returns log(1 + x) for double precision.
    Computes log(1+x) using plog_core_double for the core range reduction and
    polynomial evaluation.  The rounding error from forming u = fl(1+x) is
    recovered as dx = x - (u - 1) and folded in as a first-order correction
    dx/u after the polynomial evaluation.
 */
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 = pset1frombits<Packet>(static_cast<uint64_t>(0xfff0000000000000ull));
  const Packet cst_pos_inf = pset1frombits<Packet>(static_cast<uint64_t>(0x7ff0000000000000ull));

  // 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;
}

/** \internal \returns log(1 + x) computed using W. Kahan's formula.
    See: http://www.plunk.org/~hatch/rightway.php
    This is the generic fallback for types without a specialized implementation.
 */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_log1p(const Packet& x) {
  typedef typename unpacket_traits<Packet>::type ScalarType;
  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);
}

/** \internal \returns exp(x)-1 computed using W. Kahan's formula.
    See: http://www.plunk.org/~hatch/rightway.php
 */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_expm1(const Packet& x) {
  typedef typename unpacket_traits<Packet>::type ScalarType;
  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));
}

// 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) {
  typedef typename unpacket_traits<Packet>::integer_packet PacketI;

  const Packet cst_one = pset1<Packet>(1.0f);
  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);
  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);

  // 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);

  // Evaluate the 6th order polynomial approximation to exp(r)
  // with r in the interval [-ln(2)/2;ln(2)/2].
  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);
  Packet y = pmadd(r, p_odd, p_even);
  y = pmadd(r2, y, p_low);

  // Construct 2^m by directly manipulating the exponent bits.
  // After clamping, m is in [-150, 128], so biased exponent m+127 is in [-23, 255].
  // We only need the lower clamp to 0 (the upper bound 255 is exact).
  const PacketI cst_bias = pset1<PacketI>(127);
  PacketI mi = pcast<Packet, PacketI>(m);
  mi = pmax(padd(mi, cst_bias), pzero(mi));
  const Packet pow2m = preinterpret<Packet>(plogical_shift_left<23>(mi));
  y = pmul(y, pow2m);

  if (!IsFinite) {
    // Handle NaN: exp(nan) = nan. Use pmax to propagate NaN from input.
    y = pmax(y, _x);
  }
  return y;
}

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp_double(const Packet _x) {
  const Packet cst_zero = pset1<Packet>(0.0);
  const Packet cst_1 = pset1<Packet>(1.0);
  const Packet cst_2 = pset1<Packet>(2.0);

  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_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 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);

  Packet x2 = pmul(x, x);

  // 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, x);

  // 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).
  x = pdiv(px, psub(qx, px));
  x = pmadd(cst_2, x, cst_1);

  // 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 pselect(zero_mask, cst_zero, pmax(pldexp(x, fx), _x));
}

// This function computes exp2(x) = exp(ln(2) * x).
// To improve accuracy, the product ln(2)*x is computed using the twoprod
// algorithm, such that ln(2) * x = p_hi + p_lo holds exactly. Then exp2(x) is
// computed as exp2(x) = exp(p_hi) * exp(p_lo) ~= exp(p_hi) * (1 + p_lo). This
// correction step this reduces the maximum absolute error as follows:
//
// type   | max error (simple product) | max error (twoprod) |
// -----------------------------------------------------------
// float  |       35 ulps              |       4 ulps        |
// double |      363 ulps              |     110 ulps        |
//
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_exp2(const Packet& _x) {
  typedef typename unpacket_traits<Packet>::type Scalar;
  constexpr int max_exponent = std::numeric_limits<Scalar>::max_exponent;
  constexpr int digits = std::numeric_limits<Scalar>::digits;
  constexpr Scalar max_cap = Scalar(max_exponent + 1);
  constexpr Scalar min_cap = -Scalar(max_exponent + digits - 1);
  Packet x = pmax(pmin(_x, pset1<Packet>(max_cap)), pset1<Packet>(min_cap));
  Packet p_hi, p_lo;
  twoprod(pset1<Packet>(Scalar(EIGEN_LN2)), x, p_hi, p_lo);
  Packet exp2_hi = pexp(p_hi);
  Packet exp2_lo = padd(pset1<Packet>(Scalar(1)), p_lo);
  return pmul(exp2_hi, exp2_lo);
}

/** \internal \returns log10(x) for single precision float.
    Computed as log(x) * log10(e).
    Simply multiplying by a single float constant loses accuracy because
    float(log10(e)) has rounding error. We use a hi+lo split instead:
    log10(x) ~= log(x) * hi + log(x) * lo, computed via fma. */
template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog10_float(const Packet& x) {
  typedef typename unpacket_traits<Packet>::type Scalar;
  // 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 = pset1<Packet>(NumTraits<Scalar>::infinity());
  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);
}

/** \internal \returns log10(x) for double precision float.
    Computed as log(x) * log10(e). For double, single-constant rounding error
    is ~1e-17 (sub-ULP), so the simple form is precise enough. */
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>
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>::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 abs_a = pabs(a);
    const Packet sign_mask = pandnot(a, abs_a);
    const Packet nonzero_mask = pcmp_lt(cst_zero, abs_a);

    return pselect(nonzero_mask, por(sign_mask, cst_one), abs_a);
  }
};

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) {
    typedef typename unpacket_traits<Packet>::type Scalar;
    typedef typename Scalar::value_type RealScalar;
    typedef typename unpacket_traits<Packet>::as_real RealPacket;

    // 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 = pset1<Packet>(static_cast<Scalar>(-0.0));
  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