| // SPDX-FileCopyrightText: The Eigen Authors |
| // SPDX-License-Identifier: MPL-2.0 |
| |
| #include "main.h" |
| |
| #include <cfenv> |
| #include <cmath> |
| #include <complex> |
| #include <limits> |
| #include <vector> |
| |
| // use_double_word<packet_traits<Scalar>::type> names a vector type such as __m512d as a template argument. |
| #if EIGEN_GNUC_STRICT_AT_LEAST(6, 0, 0) |
| #pragma GCC diagnostic ignored "-Wignored-attributes" |
| #endif |
| |
| // Integer powers of binary32 and binary64 bases, real and complex, are computed by repeated squaring in |
| // double-word arithmetic: correctly rounded up to the rare cases where the exact power lies within about |
| // 7 * n * u^2 of a rounding boundary. These tests check that bound against independent references, on the |
| // vectorized path and on the scalar path, together with the exceptional values and signed zeros. |
| |
| template <typename Scalar> |
| Scalar ulp_of(const Scalar& x) { |
| Scalar a = numext::abs(x); |
| return numext::nextafter(a, std::numeric_limits<Scalar>::infinity()) - a; |
| } |
| |
| // |a - b| in units of ulp(b), with non-finite and zero references matched exactly (including the sign of zero). |
| template <typename Scalar> |
| bool within_ulps(const Scalar& a, const Scalar& b, double ulps) { |
| if ((numext::isnan)(a) || (numext::isnan)(b)) return (numext::isnan)(a) && (numext::isnan)(b); |
| if ((numext::isinf)(b) || b == Scalar(0)) return a == b && bool(numext::signbit(a)) == bool(numext::signbit(b)); |
| if ((numext::isinf)(a)) return false; |
| return double(numext::abs(a - b)) <= ulps * double(ulp_of(b)); |
| } |
| |
| // Log-uniform magnitudes over [2^lo, 2^hi], both signs. |
| template <typename Scalar> |
| std::vector<Scalar> log_uniform_bases(int lo, int hi, int count) { |
| std::vector<Scalar> bases; |
| for (int i = 0; i < count; ++i) { |
| double exponent = internal::random<double>(double(lo), double(hi)); |
| Scalar x = Scalar(std::ldexp(internal::random<double>(1.0, 2.0), int(std::floor(exponent)))); |
| bases.push_back(i % 2 ? -x : x); |
| } |
| return bases; |
| } |
| |
| // y = x.pow(n), returning which of FE_INVALID and FE_OVERFLOW it raised on a target whose packets flush subnormals, |
| // where the results below the smallest normal are rebuilt from integer bits. Elsewhere it reports none: the double |
| // word can raise FE_INVALID for a zero base, which an AVX2 partial packet holds as padding. |
| template <typename ArrayType, typename Exponent> |
| int pow_exceptions(const ArrayType& x, const Exponent& n, ArrayType& y) { |
| using Packet = typename internal::packet_traits<typename ArrayType::Scalar>::type; |
| std::fenv_t saved_environment; |
| const bool environment_saved = |
| internal::unary_pow::flushes_subnormals<Packet>::value && std::feholdexcept(&saved_environment) == 0; |
| y = x.pow(n); |
| const int exceptions = environment_saved ? std::fetestexcept(FE_INVALID | FE_OVERFLOW) : 0; |
| if (environment_saved) std::fesetenv(&saved_environment); |
| return exceptions; |
| } |
| |
| // Reference: pow in the next wider format, so its own rounding is negligible for float; for double, the |
| // platform pow, which is within 1 ulp on every supported libm. |
| template <typename Scalar> |
| Scalar reference_pow(const Scalar& x, long n) { |
| using Wide = std::conditional_t<std::is_same<Scalar, float>::value, double, Scalar>; |
| return Scalar(std::pow(Wide(x), Wide(n))); |
| } |
| |
| template <typename Scalar, typename Exponent> |
| void check_real_pow(const std::vector<Scalar>& bases, const std::vector<long>& exponents, double ulps) { |
| Index size = 2 * internal::packet_traits<Scalar>::size + 1; |
| ArrayX<Scalar> x(size), y(size); |
| for (long n : exponents) { |
| Exponent exponent = Exponent(n); |
| for (const Scalar& base : bases) { |
| x.setConstant(base); |
| y = x.pow(exponent); |
| Scalar reference = reference_pow(base, n); |
| for (Index i = 0; i < size; ++i) { |
| bool ok = within_ulps(y(i), reference, ulps); |
| if (!ok) |
| std::cout << "pow(" << base << ", " << n << ") = " << y(i) << " != " << reference << " (" << ulps << " ulp)" |
| << std::endl; |
| VERIFY(ok); |
| } |
| } |
| } |
| } |
| |
| template <typename Scalar> |
| void real_pow_test() { |
| // Exponents on both sides of the squaring cutoffs: the plain-squaring range [-3, 7] of other scalar types, the |
| // fallback to generic_pow above 2^12 (float) or 2^20 (double), and integer exponents beyond the significand, |
| // which generic_pow could not represent. |
| long cutoff = long(internal::unary_pow::max_squaring_exponent<Scalar>()); |
| const std::vector<long> exponents = {2, 3, 4, 5, 7, 8, 9, 15, 16, 17, 31, 100, 127, 1000, |
| -2, -3, -4, -5, -7, -8, -9, -15, -100, -127, -1000, 4096, cutoff, -cutoff}; |
| const std::vector<long> huge_exponents = {cutoff + 1, |
| -cutoff - 1, |
| 1 << 24, |
| -(1 << 24), |
| (1 << 24) + 1, |
| -((1 << 24) + 1), |
| (std::numeric_limits<int>::max)(), |
| (std::numeric_limits<int>::min)()}; |
| constexpr int min_exponent = std::numeric_limits<Scalar>::min_exponent; |
| constexpr int max_exponent = std::numeric_limits<Scalar>::max_exponent; |
| |
| // Correctly rounded within 2 ulps only holds for double-word repeated squaring, whose fallback -- plain |
| // squaring -- is just "a few ulps" and exercised elsewhere. Only the MSA and HVX packets lack the path, so |
| // require it everywhere else rather than skipping on the trait that selects it: losing it has to fail here. |
| using Packet = typename internal::packet_traits<Scalar>::type; |
| constexpr bool has_double_word = internal::unary_pow::use_double_word<Packet>::value; |
| #if !defined(EIGEN_VECTORIZE_MSA) && !defined(EIGEN_VECTORIZE_HVX) |
| VERIFY(has_double_word); |
| #endif |
| if (has_double_word) { |
| for (long n : exponents) { |
| int limit = numext::mini((max_exponent - 2) / int(numext::abs(n)), max_exponent - 2); |
| const std::vector<Scalar> bases = log_uniform_bases<Scalar>(-limit, limit, 64); |
| check_real_pow<Scalar, int>(bases, {n}, 2.0); |
| check_real_pow<Scalar, Scalar>(bases, {n}, 2.0); |
| } |
| } |
| // Results near the underflow and overflow thresholds: the final scaling rounds once more into the subnormal |
| // range, and the reference rounds too. |
| { |
| const std::vector<Scalar> bases = log_uniform_bases<Scalar>(min_exponent - 4, 4, 32); |
| check_real_pow<Scalar, int>(bases, {2, 3, 8, -2, -3, -8}, 3.0); |
| check_real_pow<Scalar, Scalar>(bases, {2, 3, 8, -2, -3, -8}, 3.0); |
| const std::vector<Scalar> tiny = log_uniform_bases<Scalar>(min_exponent - 20, min_exponent / 8, 32); |
| check_real_pow<Scalar, int>(tiny, {8, 9, 100, -8, -9, -100}, 3.0); |
| check_real_pow<Scalar, Scalar>(tiny, {8, 9, 100, -8, -9, -100}, 3.0); |
| } |
| // x^-1 and x^2 round once, as 1/x and x * x do, also where their results are subnormal. |
| { |
| const Index size = 2 * internal::packet_traits<Scalar>::size + 1; |
| for (const Scalar& base : log_uniform_bases<Scalar>(max_exponent - 3, max_exponent, 64)) { |
| ArrayX<Scalar> y = ArrayX<Scalar>::Constant(size, base).pow(-1); |
| for (Index k = 0; k < size; ++k) VERIFY(within_ulps(y(k), Scalar(1) / base, 0.0)); |
| } |
| for (const Scalar& base : log_uniform_bases<Scalar>(min_exponent / 2 - 13, min_exponent / 2, 64)) { |
| ArrayX<Scalar> y = ArrayX<Scalar>::Constant(size, base).pow(2); |
| for (Index k = 0; k < size; ++k) VERIFY(within_ulps(y(k), base * base, 0.0)); |
| } |
| } |
| // A subnormal power beside normal ones, 2^-n and 1^n in one packet, raises neither FE_INVALID nor FE_OVERFLOW. |
| { |
| const int n = 1 - min_exponent + std::numeric_limits<Scalar>::digits / 2; |
| const Scalar tiny = Scalar(std::ldexp(1.0, -n)); |
| ArrayX<Scalar> x = ArrayX<Scalar>::Ones(2 * internal::packet_traits<Scalar>::size + 1), y; |
| x(0) = Scalar(0.5); |
| VERIFY_IS_EQUAL(pow_exceptions(x, n, y), 0); |
| for (Index k = 0; k < x.size(); ++k) VERIFY(within_ulps(y(k), k == 0 ? tiny : Scalar(1), 0.0)); |
| VERIFY_IS_EQUAL(pow_exceptions(x, Scalar(n), y), 0); |
| for (Index k = 0; k < x.size(); ++k) VERIFY(within_ulps(y(k), k == 0 ? tiny : Scalar(1), 0.0)); |
| } |
| // Exponents beyond the squaring cutoff, on bases within a few ulps of one, where the power stays finite for |
| // every exponent up to INT_MAX. The reference pow then carries the full rounding error of its own exponent |
| // conversion, which is why the double reference is allowed a wider margin. |
| { |
| std::vector<Scalar> bases; |
| for (int i = -16; i <= 16; ++i) bases.push_back(Scalar(1) + Scalar(i) * NumTraits<Scalar>::epsilon()); |
| check_real_pow<Scalar, int>(bases, huge_exponents, 3.0); |
| for (long n : huge_exponents) |
| if (n == long(Scalar(n))) check_real_pow<Scalar, Scalar>(bases, {n}, 3.0); |
| } |
| // Exceptional values and signed zeros, exactly as std::pow. |
| { |
| Scalar inf = std::numeric_limits<Scalar>::infinity(); |
| Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); |
| Scalar subnormal = std::numeric_limits<Scalar>::denorm_min() * Scalar(3); |
| const std::vector<Scalar> specials = {Scalar(0), |
| -Scalar(0), |
| inf, |
| -inf, |
| nan, |
| subnormal, |
| -subnormal, |
| (std::numeric_limits<Scalar>::max)(), |
| (std::numeric_limits<Scalar>::min)(), |
| Scalar(1), |
| Scalar(-1)}; |
| const std::vector<long> small = {1, 2, 3, 4, 5, 8, 9, 100, -1, -2, -3, -4, -5, -8, -9, -100}; |
| check_real_pow<Scalar, int>(specials, small, 1.0); |
| check_real_pow<Scalar, Scalar>(specials, small, 1.0); |
| } |
| } |
| |
| // Complex reference: componentwise double-word repeated multiplication, no scaling, so the bases are kept where |
| // nothing under- or overflows. |
| template <typename Real> |
| struct double_word { |
| Real hi, lo; |
| static double_word product(const double_word& a, const double_word& b) { |
| Real p = a.hi * b.hi; |
| Real e = numext::fma(a.hi, b.hi, -p) + (a.hi * b.lo + a.lo * b.hi); |
| Real s = p + e; |
| return {s, e - (s - p)}; |
| } |
| static double_word sum(const double_word& a, const double_word& b) { |
| Real s = a.hi + b.hi; |
| Real bb = s - a.hi; |
| Real e = ((a.hi - (s - bb)) + (b.hi - bb)) + a.lo + b.lo; |
| Real t = s + e; |
| return {t, e - (t - s)}; |
| } |
| static double_word negate(const double_word& a) { return {-a.hi, -a.lo}; } |
| }; |
| |
| template <typename Real> |
| std::complex<Real> reference_complex_pow(const std::complex<Real>& z, long n) { |
| using DW = double_word<Real>; |
| DW re{Real(1), Real(0)}, im{Real(0), Real(0)}; |
| DW zr{numext::real(z), Real(0)}, zi{numext::imag(z), Real(0)}; |
| for (long k = 0; k < numext::abs(n); ++k) { |
| DW new_re = DW::sum(DW::product(re, zr), DW::negate(DW::product(im, zi))); |
| im = DW::sum(DW::product(re, zi), DW::product(im, zr)); |
| re = new_re; |
| } |
| if (n < 0) { |
| // 1/(re + im i) = (re - im i)/(re^2 + im^2) in double-word arithmetic, then rounded. |
| using LongDouble = long double; |
| DW norm = DW::sum(DW::product(re, re), DW::product(im, im)); |
| LongDouble inv = 1.0L / (LongDouble(norm.hi) + LongDouble(norm.lo)); |
| return std::complex<Real>(Real((LongDouble(re.hi) + LongDouble(re.lo)) * inv), |
| Real(-(LongDouble(im.hi) + LongDouble(im.lo)) * inv)); |
| } |
| return std::complex<Real>(Real(re.hi + re.lo), Real(im.hi + im.lo)); |
| } |
| |
| // The reciprocal above rounds in long double, which is a reference only where that is wider than Real: not where |
| // long double is binary64 (MSVC, Apple arm64), for Real = double. |
| template <typename Real> |
| bool has_complex_reference(long n) { |
| return n >= 0 || sizeof(long double) > sizeof(Real); |
| } |
| |
| // The error of a complex result is measured against the magnitude of the reference: a component that nearly |
| // cancels is not required to be accurate on its own. |
| template <typename Real> |
| bool complex_within_ulps(const std::complex<Real>& a, const std::complex<Real>& b, double ulps) { |
| double magnitude = double(numext::abs(b)); |
| double error = double(numext::abs(a - b)); |
| return error <= ulps * double(ulp_of(Real(magnitude))); |
| } |
| |
| template <typename Real, typename Exponent> |
| void check_complex_pow(const std::vector<std::complex<Real>>& bases, const std::vector<long>& exponents, double ulps) { |
| using Complex = std::complex<Real>; |
| Index size = 2 * internal::packet_traits<Complex>::size + 1; |
| ArrayX<Complex> z(size), y(size); |
| for (long n : exponents) { |
| if (!has_complex_reference<Real>(n)) continue; |
| Exponent exponent = Exponent(n); |
| for (const Complex& base : bases) { |
| z.setConstant(base); |
| y = z.pow(exponent); |
| Complex reference = reference_complex_pow(base, n); |
| for (Index i = 0; i < size; ++i) { |
| bool ok = complex_within_ulps(y(i), reference, ulps); |
| if (!ok) |
| std::cout << "pow(" << base << ", " << n << ") = " << y(i) << " != " << reference << " (" << ulps << " ulp)" |
| << std::endl; |
| VERIFY(ok); |
| } |
| } |
| } |
| } |
| |
| template <typename Real> |
| void complex_pow_test() { |
| using Complex = std::complex<Real>; |
| Index size = 2 * internal::packet_traits<Complex>::size + 1; |
| |
| // Small integer inputs have exactly representable powers. |
| { |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, Complex(3, 4)); |
| VERIFY((z.pow(3) == Complex(-117, 44)).all()); |
| VERIFY((z.pow(Real(3)) == Complex(-117, 44)).all()); |
| VERIFY((z.pow(4) == Complex(-527, -336)).all()); |
| z.setConstant(Complex(1, 1)); |
| VERIFY((z.pow(2) == Complex(0, 2)).all()); |
| VERIFY((z.pow(Real(2)) == Complex(0, 2)).all()); |
| VERIFY((z.pow(4) == Complex(-4, 0)).all()); |
| VERIFY((z.pow(-2) == Complex(0, Real(-0.5))).all()); |
| VERIFY((z.pow(-4) == Complex(Real(-0.25), 0)).all()); |
| VERIFY((z.pow(0) == Complex(1, 0)).all()); |
| VERIFY((z.pow(1) == Complex(1, 1)).all()); |
| // (1+i)^(8k+7) = 2^(4k+3) (1 - i): exact, with a magnitude of sqrt(2) times a finite component. The unscaled |
| // loop must not admit it once the exponent budget is exhausted, or its intermediates overflow. |
| constexpr int k = (std::numeric_limits<Real>::max_exponent - 8) / 4; |
| Real big = Real(std::ldexp(1.0, 4 * k + 3)); |
| VERIFY((z.pow(8 * k + 7) == Complex(big, -big)).all()); |
| VERIFY((z.pow(Real(8 * k + 7)) == Complex(big, -big)).all()); |
| } |
| // Components so far apart that the smaller one underflows when both take the larger one's scale, while the |
| // power's imaginary part, 2ab, is normal: each component is checked on its own, as the magnitude-relative check |
| // would not notice the imaginary part missing. |
| { |
| Real big = Real(std::ldexp(1.0, std::numeric_limits<Real>::max_exponent / 2 - 2)); |
| Real tiny = Real(std::ldexp(1.0, std::numeric_limits<Real>::min_exponent + std::numeric_limits<Real>::digits)); |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, Complex(big, tiny)); |
| ArrayX<Complex> y = z.pow(2); |
| for (Index k = 0; k < size; ++k) { |
| VERIFY(within_ulps(numext::real(y(k)), big * big, 1.0)); |
| VERIFY(within_ulps(numext::imag(y(k)), Real(2) * big * tiny, 1.0)); |
| } |
| } |
| // Powers with both components subnormal, exact: 2^(5s) times (1 + i)^5 = -4 (1 + i) and (1 + i)^-5 = -(1 - i) / 8. |
| // Compared coefficient by coefficient, as a vectorized comparison may flush the expected value to zero too. |
| { |
| const int s = (std::numeric_limits<Real>::min_exponent - 5) / 5; |
| const Complex tiny = Real(std::ldexp(1.0, 5 * s + 2)) * Complex(-1, -1); |
| const Complex tinier = Real(std::ldexp(1.0, 5 * s - 3)) * Complex(-1, 1); |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, Real(std::ldexp(1.0, s)) * Complex(1, 1)); |
| ArrayX<Complex> y = z.pow(5), y_real = z.pow(Real(5)); |
| for (Index k = 0; k < size; ++k) VERIFY(y(k) == tiny && y_real(k) == tiny); |
| z.setConstant(Real(std::ldexp(1.0, -s)) * Complex(1, 1)); |
| y = z.pow(-5); |
| y_real = z.pow(Real(-5)); |
| for (Index k = 0; k < size; ++k) VERIFY(y(k) == tinier && y_real(k) == tinier); |
| } |
| // A zero component beside nonzero ones, (1 + i)^2 = 2i and (1 + 2i)^2 = -3 + 4i in one packet, raises neither |
| // FE_INVALID nor FE_OVERFLOW. |
| { |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, Complex(1, 2)), y; |
| z(0) = Complex(1, 1); |
| VERIFY_IS_EQUAL(pow_exceptions(z, 2, y), 0); |
| for (Index k = 0; k < size; ++k) VERIFY(y(k) == (k == 0 ? Complex(0, 2) : Complex(-3, 4))); |
| VERIFY_IS_EQUAL(pow_exceptions(z, Real(2), y), 0); |
| for (Index k = 0; k < size; ++k) VERIFY(y(k) == (k == 0 ? Complex(0, 2) : Complex(-3, 4))); |
| } |
| // Reciprocals of normal bases whose norm exceeds 1 / min, c (1 + i) with |z|^2 = 4.5 * 2^(2e), or whose smaller |
| // component is subnormal: 1/(2^p + 2^-30 i) = 2^-p - 2^(-30 - 2p) i to working precision. |
| { |
| const int e = (std::numeric_limits<Real>::max_exponent - 1) / 2 - 1; |
| const Real c = Real(std::ldexp(1.5, e)); |
| ArrayX<Complex> y = ArrayX<Complex>::Constant(size, Complex(c, c)).pow(-1); |
| for (Index k = 0; k < size; ++k) { |
| VERIFY(within_ulps(numext::real(y(k)), Real(1) / (Real(2) * c), 2.0)); |
| VERIFY(within_ulps(numext::imag(y(k)), Real(-1) / (Real(2) * c), 2.0)); |
| } |
| const int p = (std::numeric_limits<Real>::max_exponent - 1) / 2 - 13; |
| y = ArrayX<Complex>::Constant(size, Complex(Real(std::ldexp(1.0, p)), Real(std::ldexp(1.0, -30)))).pow(-1); |
| for (Index k = 0; k < size; ++k) { |
| VERIFY(within_ulps(numext::real(y(k)), Real(std::ldexp(1.0, -p)), 2.0)); |
| VERIFY(within_ulps(numext::imag(y(k)), Real(-std::ldexp(1.0, -30 - 2 * p)), 2.0)); |
| } |
| } |
| // Separated components near |z| = 1, with a large |n|: (a + bi)^n = a^n (1 + i n b/a) and |
| // (b + ai)^n = i^n a^n (1 - i n b/a) to working precision, both components to double-word accuracy. The smaller |
| // component is subnormal, or normal but small enough that a negative power, which shrinks between |
| // renormalizations, would take its residuals below the normal range. |
| const Complex i_powers[4] = {Complex(1, 0), Complex(0, 1), Complex(-1, 0), Complex(0, -1)}; |
| for (int b_exponent : {std::numeric_limits<Real>::min_exponent - 6, std::numeric_limits<Real>::min_exponent + 4}) { |
| Real a = Real(1) + Real(std::ldexp(1.0, -10)); |
| Real b = Real(std::ldexp(1.0, b_exponent)); |
| for (long n : {100L, 101L, 1000L, 1003L, -99L, -100L, -1000L, -1001L}) { |
| if (!has_complex_reference<Real>(n)) continue; |
| Real power = numext::real(reference_complex_pow(Complex(a), n)); |
| Real term = Real(n) * numext::real(reference_complex_pow(Complex(a), n - 1)) * b; |
| for (bool swapped : {false, true}) { |
| Complex expected = swapped ? i_powers[((n % 4) + 4) % 4] * Complex(power, -term) : Complex(power, term); |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, swapped ? Complex(b, a) : Complex(a, b)); |
| ArrayX<Complex> y = z.pow(int(n)); |
| for (Index k = 0; k < size; ++k) { |
| VERIFY(within_ulps(numext::real(y(k)), numext::real(expected), 2.0)); |
| VERIFY(within_ulps(numext::imag(y(k)), numext::imag(expected), 2.0)); |
| } |
| } |
| } |
| } |
| // (1 + 2^-40 i)^n with n = 2^30 + 1, where n r = 2^-10 makes first order 8 ulp off for float: such an n |
| // renormalizes at every step instead. The reference is the binomial series in double, whose terms |
| // C(n, k) r^k fall by 2^-10 each. |
| EIGEN_IF_CONSTEXPR (std::is_same<Real, float>::value) { |
| const long n = (1L << 30) + 1; |
| const double r = std::ldexp(1.0, -40); |
| double re = 0, im = 0, term = 1; |
| for (int k = 0; k < 12; ++k) { |
| if (k > 0) term *= double(n - k + 1) / k * r; |
| double signed_term = (k / 2) % 2 == 0 ? term : -term; |
| (k % 2 == 0 ? re : im) += signed_term; |
| } |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, Complex(1, Real(r))); |
| ArrayX<Complex> y = z.pow(int(n)); |
| for (Index k = 0; k < size; ++k) { |
| VERIFY(within_ulps(numext::real(y(k)), Real(re), 1.0)); |
| VERIFY(within_ulps(numext::imag(y(k)), Real(im), 1.0)); |
| } |
| } |
| // n = 2^digits + 1 is inexact in Real, and the first-order coefficient n b / a must not round it: the |
| // imaginary part of (1 + 1.5 * 2^e i)^(+-n) is +-(3 * 2^(digits - 1) + 1.5) 2^e, correctly rounded to |
| // +-(3 * 2^(digits - 1) + 2) 2^e, and the real part is 1. |
| { |
| constexpr int digits = std::numeric_limits<Real>::digits; |
| constexpr int e = std::numeric_limits<Real>::min_exponent - 2; |
| const long long n = (1LL << digits) + 1; |
| Real im = Real(std::ldexp(3.0 * std::ldexp(1.0, digits - 1) + 2.0, e)); |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, Complex(1, Real(std::ldexp(1.5, e)))); |
| for (long long sign : {1LL, -1LL}) { |
| ArrayX<Complex> y = z.pow(sign * n); |
| for (Index k = 0; k < size; ++k) VERIFY(y(k) == Complex(1, Real(sign) * im)); |
| } |
| } |
| // An element's power does not depend on its packet neighbours, one of which here forces the scaled loop: the |
| // unscaled loop, whose residuals would underflow in the smaller component (here 8 a b^7, subnormal), requires |
| // both components in range, so this base takes the scaled loop either way. |
| { |
| constexpr int min_exponent = std::numeric_limits<Real>::min_exponent; |
| constexpr int e = min_exponent / 9; |
| for (int i = 0; i < 16; ++i) { |
| Real a = std::ldexp(internal::random<Real>(Real(1), Real(2)), min_exponent - 20 - 7 * e); |
| Real b = std::ldexp(internal::random<Real>(Real(-2), Real(-1)), e); |
| ArrayX<Complex> alone = ArrayX<Complex>::Constant(size, Complex(a, b)), mixed = alone; |
| for (Index k = 1; k < size; k += 2) mixed(k) = Complex((std::numeric_limits<Real>::max)() / 4, 1); |
| ArrayX<Complex> y_alone = alone.pow(8), y_mixed = mixed.pow(8); |
| for (Index k = 0; k < size; k += 2) VERIFY(y_alone(k) == y_mixed(k)); |
| } |
| } |
| // Random bases at magnitudes where the reference cannot under- or overflow; the exponent runs past the |
| // renormalization interval of the implementation. |
| { |
| const std::vector<long> exponents = {2, 3, 4, 5, 7, 8, 9, 15, 16, 17, 31, 40, -2, -3, -4, -5, -8, -9, -17, -40}; |
| std::vector<Complex> bases; |
| for (int i = 0; i < 64; ++i) { |
| Real magnitude = internal::random<Real>(Real(0.5), Real(2)); |
| Real angle = internal::random<Real>(Real(-3.14159), Real(3.14159)); |
| bases.push_back(Complex(magnitude * std::cos(angle), magnitude * std::sin(angle))); |
| } |
| check_complex_pow<Real, int>(bases, exponents, 2.0); |
| check_complex_pow<Real, Real>(bases, exponents, 2.0); |
| } |
| // Magnitudes across the range, against the same base scaled by a power of two, which only moves the exponent. |
| // The scale runs through the range in steps, ends included: at its extremes |z| exceeds what the reciprocal's |
| // norm can represent. |
| { |
| constexpr int max_exponent = std::numeric_limits<Real>::max_exponent; |
| for (long n : {1L, 3L, 8L, 100L, -1L, -2L, -3L, -8L, -100L}) { |
| int extreme = (max_exponent - 4) / int(numext::abs(n)) - 1; // |base|^n within 2^(+-|n|) stays finite |
| for (int i = 0; i <= 32; ++i) { |
| int e = -extreme + (2 * extreme * i) / 32; |
| Real scale = Real(std::ldexp(1.0, e)); |
| Real magnitude = internal::random<Real>(Real(0.5), Real(1.5)); |
| Real angle = internal::random<Real>(Real(-3.14159), Real(3.14159)); |
| Complex base(magnitude * std::cos(angle), magnitude * std::sin(angle)); |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, base), w = ArrayX<Complex>::Constant(size, base * scale); |
| ArrayX<Complex> zn = z.pow(int(n)), wn = w.pow(int(n)); |
| Real power_of_scale = Real(std::pow(double(scale), double(n))); |
| for (Index k = 0; k < size; ++k) { |
| Complex expected = zn(k) * power_of_scale; |
| VERIFY(complex_within_ulps(wn(k), expected, 2.0)); |
| } |
| } |
| } |
| } |
| // Zero and infinite bases are the points 0 and infinity of the extended plane, any value with an infinite |
| // component being the infinity +inf + 0i of std::proj: a zero base gives +0 + 0i for n > 1 and +inf + 0i for n < 0, |
| // an infinite base the reverse, whatever the signs of their zeros; a NaN component without an infinite one gives |
| // NaN + NaN i. std::pow leaves all of these to the implementation. In a packet, in the scalar tail, and next to |
| // ordinary bases in the same packet, for integer and floating exponents; n = 1 returns the base. |
| { |
| Real inf = std::numeric_limits<Real>::infinity(); |
| Real nan = std::numeric_limits<Real>::quiet_NaN(); |
| auto same = [](const Complex& a, const Complex& b) { |
| auto same_real = [](const Real& x, const Real& y) { |
| return ((numext::isnan)(x) && (numext::isnan)(y)) || (x == y && std::signbit(x) == std::signbit(y)); |
| }; |
| return same_real(numext::real(a), numext::real(b)) && same_real(numext::imag(a), numext::imag(b)); |
| }; |
| std::vector<Complex> zeros = {Complex(0, 0), Complex(-Real(0), 0), Complex(0, -Real(0)), |
| Complex(-Real(0), -Real(0))}; |
| std::vector<Complex> infinities = {Complex(inf, 0), Complex(-inf, 2), Complex(0, inf), Complex(-2, -inf), |
| Complex(inf, inf), Complex(-inf, inf), Complex(inf, nan), Complex(nan, -inf)}; |
| std::vector<Complex> nans = {Complex(nan, 1), Complex(1, nan), Complex(nan, 0), Complex(nan, nan)}; |
| for (int n : {2, 3, 4, 5, -1, -2, -3, -4}) { |
| for (int kind = 0; kind < 3; ++kind) { |
| const std::vector<Complex>& bases = kind == 0 ? zeros : kind == 1 ? infinities : nans; |
| Complex expected = kind == 2 ? Complex(nan, nan) : (kind == 0) == (n > 0) ? Complex(0, 0) : Complex(inf, 0); |
| for (const Complex& base : bases) { |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, base); |
| for (Index k = 1; k < size; k += 2) z(k) = Complex(Real(0.75), Real(-1.25)); |
| ArrayX<Complex> y = z.pow(n), y_real = z.pow(Real(n)), ordinary = ArrayX<Complex>::Constant(1, z(1)).pow(n); |
| for (Index k = 0; k < size; k += 2) VERIFY(same(y(k), expected) && same(y_real(k), expected)); |
| for (Index k = 1; k < size; k += 2) VERIFY(same(y(k), ordinary(0))); |
| VERIFY(same(ArrayX<Complex>::Constant(1, base).pow(1)(0), base)); |
| } |
| } |
| } |
| } |
| // A base with one zero component is the limit of a base whose component is tiny and of the same sign, and the |
| // zero component of its power takes the sign of that limit: here of the power of the perturbed base, computed |
| // in double by repeated multiplication, where the tiny component cannot cancel. |
| for (Real a : {Real(2), Real(-2), Real(0.5), Real(-3)}) { |
| for (Real s : {Real(0), -Real(0)}) { |
| for (bool swapped : {false, true}) { |
| for (int n = -12; n <= 12; ++n) { |
| if (n == 0) continue; |
| double tiny = std::copysign(std::ldexp(1.0, -60), double(s)); |
| std::complex<double> w = swapped ? std::complex<double>(tiny, a) : std::complex<double>(a, tiny); |
| if (n < 0) w = std::conj(w) / std::norm(w); |
| std::complex<double> reference = w; |
| for (int k = 1; k < numext::abs(n); ++k) reference *= w; |
| bool real_is_zero = std::abs(reference.real()) < std::abs(reference.imag()); |
| double zero_sign = real_is_zero ? reference.real() : reference.imag(); |
| ArrayX<Complex> z = ArrayX<Complex>::Constant(size, swapped ? Complex(s, a) : Complex(a, s)); |
| ArrayX<Complex> y = z.pow(n); |
| for (Index k = 0; k < size; ++k) { |
| Real zero = real_is_zero ? numext::real(y(k)) : numext::imag(y(k)); |
| VERIFY(zero == Real(0) && std::signbit(zero) == std::signbit(zero_sign)); |
| } |
| } |
| } |
| } |
| } |
| } |
| |
| // half and bfloat16 keep plain repeated squaring, whose error grows like n ulp, for integer exponents and small |
| // floating-point ones, and generic_pow (pow through float) beyond; with F16C their packets take the vectorized |
| // path. Check both against float. |
| template <typename Scalar> |
| void narrow_pow_test() { |
| Index size = 2 * internal::packet_traits<Scalar>::size + 1; |
| ArrayX<Scalar> x(size), y(size); |
| for (long n : {1L, 2L, 3L, 5L, 7L, 8L, 16L, 100L, 4097L, -1L, -2L, -3L, -8L, -100L}) { |
| // Powers within the normal range of the type: |x|^n in 2^(+-(max_exponent - 2)). |
| double limit = double(std::numeric_limits<Scalar>::max_exponent - 2) / double(numext::abs(n)); |
| for (int i = 0; i < 32; ++i) { |
| float base = float(std::exp2(internal::random<double>(-limit, limit))); |
| x.setConstant(Scalar(i % 2 ? -base : base)); |
| Scalar reference = Scalar(std::pow(float(x(0)), float(n))); |
| double ulps = 2.0 + double(numext::abs(n)); |
| y = x.pow(int(n)); |
| for (Index k = 0; k < size; ++k) { |
| if (!within_ulps(y(k), reference, ulps)) |
| std::cout << "pow(" << x(0) << ", " << n << ") = " << y(k) << " != " << reference << std::endl; |
| VERIFY(within_ulps(y(k), reference, ulps)); |
| } |
| if (n == long(Scalar(float(n)))) { |
| y = x.pow(Scalar(float(n))); |
| for (Index k = 0; k < size; ++k) VERIFY(within_ulps(y(k), reference, ulps)); |
| } |
| } |
| } |
| } |
| |
| // complex<long double> keeps plain squaring; its reciprocal must stay range-safe, as the double-word path's |
| // conj(z)/|z|^2 would not be for long double's wider exponent range. |
| void complex_long_double_range_test() { |
| using Complex = std::complex<long double>; |
| constexpr int max_exponent = std::numeric_limits<long double>::max_exponent; |
| if (max_exponent <= std::numeric_limits<double>::max_exponent) return; |
| ArrayX<Complex> z(9), y(9); |
| auto close = [](const Complex& a, const Complex& b, long double ulps) { |
| return std::abs(a - b) <= ulps * std::numeric_limits<long double>::epsilon() * std::abs(b); |
| }; |
| // Built with ldexp rather than written as literals: 1e3000L and the like are out of range, and so a |
| // compile-time error, wherever long double is binary64, which the run-time guard above cannot prevent. |
| for (int e : {(4 * max_exponent) / 5, -(4 * max_exponent) / 5, (2 * max_exponent) / 5, -(2 * max_exponent) / 5}) { |
| long double magnitude = std::ldexp(1.0L, e); |
| z.setConstant(Complex(magnitude, magnitude / 3)); |
| y = z.pow(-1); |
| Complex expected = Complex(1) / z(0); |
| for (Index k = 0; k < 9; ++k) VERIFY(close(y(k), expected, 4.0L)); |
| if (2 * numext::abs(e) >= max_exponent) continue; // the square would leave the range |
| y = z.pow(-2); |
| expected = Complex(1) / (z(0) * z(0)); |
| for (Index k = 0; k < 9; ++k) VERIFY(close(y(k), expected, 8.0L)); |
| } |
| } |
| |
| EIGEN_DECLARE_TEST(unary_pow) { |
| for (int i = 0; i < g_repeat; i++) { |
| CALL_SUBTEST_1(real_pow_test<float>()); |
| CALL_SUBTEST_2(real_pow_test<double>()); |
| CALL_SUBTEST_3(complex_pow_test<float>()); |
| CALL_SUBTEST_4(complex_pow_test<double>()); |
| CALL_SUBTEST_5(narrow_pow_test<half>()); |
| CALL_SUBTEST_6(narrow_pow_test<bfloat16>()); |
| CALL_SUBTEST_7(complex_long_double_range_test()); |
| } |
| } |
