// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@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

#include <limits>
#include "packetmath_test_shared.h"
#include "../Eigen/SpecialFunctions"

using internal::unpacket_traits;

#if EIGEN_ARCH_ARM
// Note: 32-bit arm always flushes subnormals to zero.
#define MAYBE_FLUSH(op)                                                      \
  [](Scalar x) {                                                             \
    Scalar y = static_cast<Scalar>(op(x));                                   \
    if (Eigen::numext::abs(y) < (std::numeric_limits<decltype(y)>::min)()) { \
      y = y * decltype(y)(0); /* Preserve sign. */                           \
    }                                                                        \
    return y;                                                                \
  }
#else
#define MAYBE_FLUSH(op) op
#endif

template <typename Scalar, typename Packet>
void packetmath_real() {
  using std::abs;
  typedef internal::packet_traits<Scalar> PacketTraits;
  const int PacketSize = internal::unpacket_traits<Packet>::size;

  const int size = PacketSize * 4;
  EIGEN_ALIGN_TO_BOUNDARY(unpacket_traits<Packet>::alignment) Scalar data1[PacketSize * 4] = {};
  EIGEN_ALIGN_TO_BOUNDARY(unpacket_traits<Packet>::alignment) Scalar data2[PacketSize * 4] = {};
  EIGEN_ALIGN_TO_BOUNDARY(unpacket_traits<Packet>::alignment) Scalar ref[PacketSize * 4] = {};

  {
    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
    test::packet_helper<internal::packet_traits<Scalar>::HasLGamma, Packet> h;
    h.store(data2, internal::plgamma(h.load(data1)));
    VERIFY((numext::isnan)(data2[0]));
  }
  if (internal::packet_traits<Scalar>::HasErf) {
    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
    data1[1] = std::numeric_limits<Scalar>::infinity();
    data1[2] = -std::numeric_limits<Scalar>::infinity();
    data1[3] = (std::numeric_limits<Scalar>::max)();
    if (size >= 8) {
      data1[4] = -(std::numeric_limits<Scalar>::max)();
      if (sizeof(Scalar) >= 8) {
        data1[5] = Scalar(1e200);
        data1[6] = Scalar(-1e200);
      } else {
        data1[5] = Scalar(1e30f);
        data1[6] = Scalar(-1e30f);
      }
      data1[7] = std::numeric_limits<Scalar>::denorm_min();
    }
    test::packet_helper<internal::packet_traits<Scalar>::HasErf, Packet> h;
    for (int i = 0; i < size; i += PacketSize) {
      h.store(data2 + i, internal::perf(h.load(data1 + i)));
    }
    VERIFY((numext::isnan)(data2[0]));
    VERIFY_IS_EQUAL(data2[1], Scalar(1));
    VERIFY_IS_EQUAL(data2[2], Scalar(-1));
    VERIFY_IS_EQUAL(data2[3], Scalar(1));
    if (size >= 8) {
      VERIFY_IS_EQUAL(data2[4], Scalar(-1));
      VERIFY_IS_EQUAL(data2[5], Scalar(1));
      VERIFY_IS_EQUAL(data2[6], Scalar(-1));
      if (data1[7] > Scalar(0)) {
        VERIFY((data2[7] > Scalar(0)));
      }
    }
  }
  if (internal::packet_traits<Scalar>::HasErfc) {
    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
    data1[1] = std::numeric_limits<Scalar>::infinity();
    data1[2] = -std::numeric_limits<Scalar>::infinity();
    data1[3] = (std::numeric_limits<Scalar>::max)();
    if (size >= 8) {
      data1[4] = -(std::numeric_limits<Scalar>::max)();
      if (sizeof(Scalar) >= 8) {
        data1[5] = Scalar(1e200);
        data1[6] = Scalar(-1e200);
      } else {
        data1[5] = Scalar(1e30f);
        data1[6] = Scalar(-1e30f);
      }
    }
    test::packet_helper<internal::packet_traits<Scalar>::HasErfc, Packet> h;
    for (int i = 0; i < size; i += PacketSize) {
      h.store(data2 + i, internal::perfc(h.load(data1 + i)));
    }
    VERIFY((numext::isnan)(data2[0]));
    VERIFY_IS_EQUAL(data2[1], Scalar(0));
    VERIFY_IS_EQUAL(data2[2], Scalar(2));
    VERIFY_IS_EQUAL(data2[3], Scalar(0));
    if (size >= 8) {
      VERIFY_IS_EQUAL(data2[4], Scalar(2));
      VERIFY_IS_EQUAL(data2[5], Scalar(0));
      VERIFY_IS_EQUAL(data2[6], Scalar(2));
    }
  }
  {
    for (int i = 0; i < size; ++i) {
      data1[i] = internal::random<Scalar>(Scalar(0), Scalar(1));
    }
    CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasNdtri, numext::ndtri, internal::pndtri);
  }

  // For bessel_i*e and bessel_j*, the valid range is negative reals.
  {
    const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10 - 1, 6);
    for (int i = 0; i < size; ++i) {
      data1[i] = internal::random<Scalar>(Scalar(-1), Scalar(1)) *
                 Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-max_exponent), Scalar(max_exponent))));
      data2[i] = internal::random<Scalar>(Scalar(-1), Scalar(1)) *
                 Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-max_exponent), Scalar(max_exponent))));
    }

    CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e);
    CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e);
    CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0);
    CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1);
  }

  // Use a smaller data range for the bessel_i* as these can become very large.
  // Following #1693, we also restrict this range further to avoid inf's due to
  // differences in pexp and exp.
  for (int i = 0; i < size; ++i) {
    data1[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
               Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1), Scalar(2))));
    data2[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
               Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1), Scalar(2))));
  }
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0);
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1);

  // Boundary values for which a naive i0(x) = exp(|x|) * i0e(x) overflows even though the result is
  // finite.  CHECK_CWISE1_IF cannot detect that on its own: its reference is the scalar path through
  // the same generic_i0, so both sides would be +inf and compare equal.  Check finiteness explicitly.
  if (PacketTraits::HasBessel &&
      (internal::is_same<Scalar, float>::value || internal::is_same<Scalar, double>::value)) {
    Scalar boundary = internal::is_same<Scalar, float>::value ? Scalar(90) : Scalar(713);
    for (int i = 0; i < size; ++i) {
      data1[i] = (i % 2 == 0) ? boundary : -boundary;
    }
    CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0);
    for (int i = 0; i < PacketSize; ++i) VERIFY((numext::isfinite)(data2[i]));
    CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1);
    for (int i = 0; i < PacketSize; ++i) VERIFY((numext::isfinite)(data2[i]));
  }

  // y_i, and k_i are valid for x > 0.
  {
    const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10 - 1, 5);
    for (int i = 0; i < size; ++i) {
      data1[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
                 Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-2), Scalar(max_exponent))));
      data2[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
                 Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-2), Scalar(max_exponent))));
    }
  }

  // TODO(srvasude): Re-enable this test once properly investigated why the
  // scalar and vector paths differ.
  // CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y0, internal::pbessel_y0);
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y1, internal::pbessel_y1);
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0e, internal::pbessel_k0e);
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1e, internal::pbessel_k1e);

  // Following #1693, we restrict the range for exp to avoid zeroing out too
  // fast.
  for (int i = 0; i < size; ++i) {
    data1[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
               Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1), Scalar(2))));
    data2[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
               Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1), Scalar(2))));
  }
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0, internal::pbessel_k0);
  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1, internal::pbessel_k1);

  for (int i = 0; i < size; ++i) {
    data1[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
               Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-1), Scalar(2))));
    data2[i] = internal::random<Scalar>(Scalar(0.01), Scalar(1)) *
               Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-1), Scalar(2))));
  }

  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
  // FIXME(rmlarsen): This test occasionally fails due to difference in tiny subnormal results
  // near the underflow boundary. I am not sure which version is correct.
  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, MAYBE_FLUSH(std::erfc), internal::perfc);
}

namespace Eigen {
namespace test {

template <typename Scalar, typename PacketType, bool IsComplex, bool IsInteger>
struct runall {
  static void run() { packetmath_real<Scalar, PacketType>(); }
};

}  // namespace test
}  // namespace Eigen

EIGEN_DECLARE_TEST(special_packetmath) {
  g_first_pass = true;
  for (int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_1(test::runner<float>::run());
    CALL_SUBTEST_2(test::runner<double>::run());
    CALL_SUBTEST_3(test::runner<Eigen::half>::run());
    CALL_SUBTEST_4(test::runner<Eigen::bfloat16>::run());
    g_first_pass = false;
  }
}