Fix random again
diff --git a/Eigen/Core b/Eigen/Core index f9d9974..ed7d353 100644 --- a/Eigen/Core +++ b/Eigen/Core
@@ -178,6 +178,7 @@ #include "src/Core/NumTraits.h" #include "src/Core/MathFunctions.h" +#include "src/Core/RandomImpl.h" #include "src/Core/GenericPacketMath.h" #include "src/Core/MathFunctionsImpl.h" #include "src/Core/arch/Default/ConjHelper.h"
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 3f28068..29e00ff 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h
@@ -604,7 +604,6 @@ struct count_bits_impl { static_assert(std::is_integral<BitsType>::value && std::is_unsigned<BitsType>::value, "BitsType must be an unsigned integer"); - static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) { int n = CHAR_BIT * sizeof(BitsType); int shift = n / 2; @@ -655,7 +654,8 @@ #if EIGEN_COMP_GNUC || EIGEN_COMP_CLANG template <typename BitsType> -struct count_bits_impl<BitsType, std::enable_if_t<sizeof(BitsType) <= sizeof(unsigned int)>> { +struct count_bits_impl< + BitsType, std::enable_if_t<std::is_integral<BitsType>::value && sizeof(BitsType) <= sizeof(unsigned int)>> { static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT); static_assert(std::is_integral<BitsType>::value, "BitsType must be a built-in integer"); static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) { @@ -669,8 +669,9 @@ }; template <typename BitsType> -struct count_bits_impl< - BitsType, std::enable_if_t<sizeof(unsigned int) < sizeof(BitsType) && sizeof(BitsType) <= sizeof(unsigned long)>> { +struct count_bits_impl<BitsType, + std::enable_if_t<std::is_integral<BitsType>::value && sizeof(unsigned int) < sizeof(BitsType) && + sizeof(BitsType) <= sizeof(unsigned long)>> { static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT); static_assert(std::is_integral<BitsType>::value, "BitsType must be a built-in integer"); static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) { @@ -684,8 +685,9 @@ }; template <typename BitsType> -struct count_bits_impl<BitsType, std::enable_if_t<sizeof(unsigned long) < sizeof(BitsType) && - sizeof(BitsType) <= sizeof(unsigned long long)>> { +struct count_bits_impl<BitsType, + std::enable_if_t<std::is_integral<BitsType>::value && sizeof(unsigned long) < sizeof(BitsType) && + sizeof(BitsType) <= sizeof(unsigned long long)>> { static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT); static_assert(std::is_integral<BitsType>::value, "BitsType must be a built-in integer"); static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) { @@ -701,7 +703,8 @@ #elif EIGEN_COMP_MSVC template <typename BitsType> -struct count_bits_impl<BitsType, std::enable_if_t<sizeof(BitsType) <= sizeof(unsigned long)>> { +struct count_bits_impl< + BitsType, std::enable_if_t<std::is_integral<BitsType>::value && sizeof(BitsType) <= sizeof(unsigned long)>> { static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT); static_assert(std::is_integral<BitsType>::value, "BitsType must be a built-in integer"); static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) { @@ -720,8 +723,9 @@ #ifdef _WIN64 template <typename BitsType> -struct count_bits_impl< - BitsType, std::enable_if_t<sizeof(unsigned long) < sizeof(BitsType) && sizeof(BitsType) <= sizeof(__int64)>> { +struct count_bits_impl<BitsType, + std::enable_if_t<std::is_integral<BitsType>::value && sizeof(unsigned long) < sizeof(BitsType) && + sizeof(BitsType) <= sizeof(__int64)>> { static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT); static_assert(std::is_integral<BitsType>::value, "BitsType must be a built-in integer"); static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) { @@ -742,192 +746,27 @@ #endif // EIGEN_COMP_GNUC || EIGEN_COMP_CLANG template <typename BitsType> -int log2_ceil(BitsType x) { - int n = CHAR_BIT * sizeof(BitsType) - clz(x); - bool powerOfTwo = (x & (x - 1)) == 0; - return x == 0 ? 0 : powerOfTwo ? n - 1 : n; +struct log_2_impl { + static constexpr int kTotalBits = sizeof(BitsType) * CHAR_BIT; + static EIGEN_DEVICE_FUNC inline int run_ceil(const BitsType& x) { + const int n = kTotalBits - clz(x); + bool power_of_two = (x & (x - 1)) == 0; + return x == 0 ? 0 : power_of_two ? (n - 1) : n; + } + static EIGEN_DEVICE_FUNC inline int run_floor(const BitsType& x) { + const int n = kTotalBits - clz(x); + return x == 0 ? 0 : n - 1; + } +}; + +template <typename BitsType> +int log2_ceil(const BitsType& x) { + return log_2_impl<BitsType>::run_ceil(x); } template <typename BitsType> -int log2_floor(BitsType x) { - int n = CHAR_BIT * sizeof(BitsType) - clz(x); - return x == 0 ? 0 : n - 1; -} - -/**************************************************************************** - * Implementation of random * - ****************************************************************************/ - -// return a Scalar filled with numRandomBits beginning from the least significant bit -template <typename Scalar> -Scalar getRandomBits(int numRandomBits) { - using BitsType = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type; - enum : int { - StdRandBits = meta_floor_log2<(unsigned int)(RAND_MAX) + 1>::value, - ScalarBits = sizeof(Scalar) * CHAR_BIT - }; - eigen_assert((numRandomBits >= 0) && (numRandomBits <= ScalarBits)); - const BitsType mask = BitsType(-1) >> ((ScalarBits - numRandomBits) & (ScalarBits - 1)); - BitsType randomBits = BitsType(0); - for (int shift = 0; shift < numRandomBits; shift += StdRandBits) { - int r = std::rand(); - randomBits |= static_cast<BitsType>(r) << shift; - } - // clear the excess bits - randomBits &= mask; - return numext::bit_cast<Scalar, BitsType>(randomBits); -} - -template <typename Scalar, bool IsComplex, bool IsInteger> -struct random_default_impl {}; - -template <typename Scalar> -struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {}; - -template <typename Scalar> -struct random_retval { - typedef Scalar type; -}; - -template <typename Scalar> -inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y); -template <typename Scalar> -inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(); - -template <typename Scalar> -struct random_default_impl<Scalar, false, false> { - using BitsType = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type; - static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y, int numRandomBits) { - Scalar half_x = Scalar(0.5) * x; - Scalar half_y = Scalar(0.5) * y; - Scalar result = (half_x + half_y) + (half_y - half_x) * run(numRandomBits); - // result is in the half-open interval [x, y) -- provided that x < y - return result; - } - static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { - const int mantissa_bits = NumTraits<Scalar>::digits() - 1; - return run(x, y, mantissa_bits); - } - static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { - const int mantissa_bits = NumTraits<Scalar>::digits() - 1; - eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissa_bits); - BitsType randomBits = getRandomBits<BitsType>(numRandomBits); - // if fewer than MantissaBits is requested, shift them to the left - randomBits <<= (mantissa_bits - numRandomBits); - // randomBits is in the half-open interval [2,4) - randomBits |= numext::bit_cast<BitsType>(Scalar(2)); - // result is in the half-open interval [-1,1) - Scalar result = numext::bit_cast<Scalar>(randomBits) - Scalar(3); - return result; - } - static EIGEN_DEVICE_FUNC inline Scalar run() { - const int mantissa_bits = NumTraits<Scalar>::digits() - 1; - return run(mantissa_bits); - } -}; - -// TODO: fix this for PPC -template <bool Specialize = sizeof(long double) == 2 * sizeof(uint64_t) && !EIGEN_ARCH_PPC> -struct random_longdouble_impl { - enum : int { - Size = sizeof(long double), - MantissaBits = NumTraits<long double>::digits() - 1, - LowBits = MantissaBits > 64 ? 64 : MantissaBits, - HighBits = MantissaBits > 64 ? MantissaBits - 64 : 0 - }; - static EIGEN_DEVICE_FUNC inline long double run() { - EIGEN_USING_STD(memcpy) - uint64_t randomBits[2]; - long double result = 2.0L; - memcpy(&randomBits, &result, Size); - randomBits[0] |= getRandomBits<uint64_t>(LowBits); - randomBits[1] |= getRandomBits<uint64_t>(HighBits); - memcpy(&result, &randomBits, Size); - result -= 3.0L; - return result; - } -}; - -// GPUs treat long double as double. -#ifndef EIGEN_GPU_COMPILE_PHASE -template <> -struct random_longdouble_impl<false> { - using Impl = random_impl<double>; - static EIGEN_DEVICE_FUNC inline long double run() { return static_cast<long double>(Impl::run()); } -}; - -template <> -struct random_impl<long double> { - static EIGEN_DEVICE_FUNC inline long double run(const long double& x, const long double& y) { - long double half_x = 0.5L * x; - long double half_y = 0.5L * y; - long double result = (half_x + half_y) + (half_y - half_x) * run(); - return result; - } - static EIGEN_DEVICE_FUNC inline long double run() { return random_longdouble_impl<>::run(); } -}; -#endif - -template <typename Scalar> -struct random_default_impl<Scalar, false, true> { - using BitsType = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type; - enum : int { ScalarBits = sizeof(Scalar) * CHAR_BIT }; - static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { - if (y <= x) return x; - const BitsType range = static_cast<BitsType>(y) - static_cast<BitsType>(x) + 1; - // handle edge case where [x,y] spans the entire range of Scalar - if (range == 0) return getRandomBits<Scalar>(ScalarBits); - // calculate the number of random bits needed to fill range - const int numRandomBits = log2_ceil(range); - BitsType randomBits; - do { - randomBits = getRandomBits<BitsType>(numRandomBits); - // if the random draw is outside [0, range), try again (rejection sampling) - // in the worst-case scenario, the probability of rejection is: 1/2 - 1/2^numRandomBits < 50% - } while (randomBits >= range); - // Avoid overflow in the case where `x` is negative and there is a large range so - // `randomBits` would also be negative if cast to `Scalar` first. - Scalar result = static_cast<Scalar>(static_cast<BitsType>(x) + randomBits); - return result; - } - - static EIGEN_DEVICE_FUNC inline Scalar run() { -#ifdef EIGEN_MAKING_DOCS - return run(Scalar(NumTraits<Scalar>::IsSigned ? -10 : 0), Scalar(10)); -#else - return getRandomBits<Scalar>(ScalarBits); -#endif - } -}; - -template <> -struct random_impl<bool> { - static EIGEN_DEVICE_FUNC inline bool run(const bool& x, const bool& y) { - if (y <= x) return x; - return run(); - } - static EIGEN_DEVICE_FUNC inline bool run() { return getRandomBits<int>(1) ? true : false; } -}; - -template <typename Scalar> -struct random_default_impl<Scalar, true, false> { - static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { - return Scalar(random(x.real(), y.real()), random(x.imag(), y.imag())); - } - static EIGEN_DEVICE_FUNC inline Scalar run() { - typedef typename NumTraits<Scalar>::Real RealScalar; - return Scalar(random<RealScalar>(), random<RealScalar>()); - } -}; - -template <typename Scalar> -inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y) { - return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y); -} - -template <typename Scalar> -inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random() { - return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(); +int log2_floor(const BitsType& x) { + return log_2_impl<BitsType>::run_floor(x); } // Implementation of is* functions
diff --git a/Eigen/src/Core/RandomImpl.h b/Eigen/src/Core/RandomImpl.h new file mode 100644 index 0000000..d5cd335 --- /dev/null +++ b/Eigen/src/Core/RandomImpl.h
@@ -0,0 +1,253 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2024 Charles Schlosser <cs.schlosser@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/. + +#ifndef EIGEN_RANDOM_IMPL_H +#define EIGEN_RANDOM_IMPL_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +namespace internal { + +/**************************************************************************** + * Implementation of random * + ****************************************************************************/ + +template <typename Scalar, bool IsComplex, bool IsInteger> +struct random_default_impl {}; + +template <typename Scalar> +struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {}; + +template <typename Scalar> +struct random_retval { + typedef Scalar type; +}; + +template <typename Scalar> +inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y) { + return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y); +} + +template <typename Scalar> +inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random() { + return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(); +} + +// TODO: replace or provide alternatives to this, e.g. std::random_device +struct eigen_random_device { + using ReturnType = int; + static constexpr int Entropy = meta_floor_log2<(unsigned int)(RAND_MAX) + 1>::value; + static constexpr ReturnType Highest = RAND_MAX; + static EIGEN_DEVICE_FUNC inline ReturnType run() { return std::rand(); }; +}; + +// Fill a built-in unsigned integer with numRandomBits beginning with the least significant bit +template <typename Scalar> +struct random_bits_impl { + EIGEN_STATIC_ASSERT(std::is_unsigned<Scalar>::value, SCALAR MUST BE A BUILT - IN UNSIGNED INTEGER) + using RandomDevice = eigen_random_device; + using RandomReturnType = typename RandomDevice::ReturnType; + static constexpr int kEntropy = RandomDevice::Entropy; + static constexpr int kTotalBits = sizeof(Scalar) * CHAR_BIT; + // return a Scalar filled with numRandomBits beginning from the least significant bit + static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { + eigen_assert((numRandomBits >= 0) && (numRandomBits <= kTotalBits)); + const Scalar mask = Scalar(-1) >> ((kTotalBits - numRandomBits) & (kTotalBits - 1)); + Scalar randomBits = 0; + for (int shift = 0; shift < numRandomBits; shift += kEntropy) { + RandomReturnType r = RandomDevice::run(); + randomBits |= static_cast<Scalar>(r) << shift; + } + // clear the excess bits + randomBits &= mask; + return randomBits; + } +}; + +template <typename BitsType> +EIGEN_DEVICE_FUNC inline BitsType getRandomBits(int numRandomBits) { + return random_bits_impl<BitsType>::run(numRandomBits); +} + +// random implementation for a built-in floating point type +template <typename Scalar, bool BuiltIn = std::is_floating_point<Scalar>::value> +struct random_float_impl { + using BitsType = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type; + static constexpr EIGEN_DEVICE_FUNC inline int mantissaBits() { + const int digits = NumTraits<Scalar>::digits(); + return digits - 1; + } + static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { + eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissaBits()); + BitsType randomBits = getRandomBits<BitsType>(numRandomBits); + // if fewer than MantissaBits is requested, shift them to the left + randomBits <<= (mantissaBits() - numRandomBits); + // randomBits is in the half-open interval [2,4) + randomBits |= numext::bit_cast<BitsType>(Scalar(2)); + // result is in the half-open interval [-1,1) + Scalar result = numext::bit_cast<Scalar>(randomBits) - Scalar(3); + return result; + } +}; +// random implementation for a custom floating point type +// uses double as the implementation with a mantissa with a size equal to either the target scalar's mantissa or that of +// double, whichever is smaller +template <typename Scalar> +struct random_float_impl<Scalar, false> { + static EIGEN_DEVICE_FUNC inline int mantissaBits() { + const int digits = NumTraits<Scalar>::digits(); + constexpr int kDoubleDigits = NumTraits<double>::digits(); + return numext::mini(digits, kDoubleDigits) - 1; + } + static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { + eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissaBits()); + Scalar result = static_cast<Scalar>(random_float_impl<double>::run(numRandomBits)); + return result; + } +}; + +// random implementation for long double +// this specialization is not compatible with double-double scalars +template <bool Specialize = (sizeof(long double) == 2 * sizeof(uint64_t)) && + ((std::numeric_limits<long double>::digits != (2 * std::numeric_limits<double>::digits)))> +struct random_longdouble_impl { + static constexpr int Size = sizeof(long double); + static constexpr EIGEN_DEVICE_FUNC inline int mantissaBits() { return NumTraits<long double>::digits() - 1; } + static EIGEN_DEVICE_FUNC inline long double run(int numRandomBits) { + eigen_assert(numRandomBits >= 0 && numRandomBits <= mantissaBits()); + EIGEN_USING_STD(memcpy); + int numLowBits = numext::mini(numRandomBits, 64); + int numHighBits = numext::maxi(numRandomBits - 64, 0); + uint64_t randomBits[2]; + long double result = 2.0L; + memcpy(&randomBits, &result, Size); + randomBits[0] |= getRandomBits<uint64_t>(numLowBits); + randomBits[1] |= getRandomBits<uint64_t>(numHighBits); + memcpy(&result, &randomBits, Size); + result -= 3.0L; + return result; + } +}; +template <> +struct random_longdouble_impl<false> { + static constexpr EIGEN_DEVICE_FUNC inline int mantissaBits() { return NumTraits<long double>::digits() - 1; } + static EIGEN_DEVICE_FUNC inline long double run(int numRandomBits) { + return static_cast<long double>(random_float_impl<double>::run(numRandomBits)); + } +}; +template <> +struct random_float_impl<long double> : random_longdouble_impl<> {}; + +template <typename Scalar> +struct random_default_impl<Scalar, false, false> { + using Impl = random_float_impl<Scalar>; + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y, int numRandomBits) { + Scalar half_x = Scalar(0.5) * x; + Scalar half_y = Scalar(0.5) * y; + Scalar result = (half_x + half_y) + (half_y - half_x) * run(numRandomBits); + // result is in the half-open interval [x, y) -- provided that x < y + return result; + } + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { + return run(x, y, Impl::mantissaBits()); + } + static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { return Impl::run(numRandomBits); } + static EIGEN_DEVICE_FUNC inline Scalar run() { return run(Impl::mantissaBits()); } +}; + +template <typename Scalar, bool IsSigned = NumTraits<Scalar>::IsSigned, bool BuiltIn = std::is_integral<Scalar>::value> +struct random_int_impl; + +// random implementation for a built-in unsigned integer type +template <typename Scalar> +struct random_int_impl<Scalar, false, true> { + static constexpr int kTotalBits = sizeof(Scalar) * CHAR_BIT; + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { + if (y <= x) return x; + Scalar range = y - x; + // handle edge case where [x,y] spans the entire range of Scalar + if (range == NumTraits<Scalar>::highest()) return run(); + Scalar count = range + 1; + // calculate the number of random bits needed to fill range + int numRandomBits = log2_ceil(count); + Scalar randomBits; + do { + randomBits = getRandomBits<Scalar>(numRandomBits); + // if the random draw is outside [0, range), try again (rejection sampling) + // in the worst-case scenario, the probability of rejection is: 1/2 - 1/2^numRandomBits < 50% + } while (randomBits >= count); + Scalar result = x + randomBits; + return result; + } + static EIGEN_DEVICE_FUNC inline Scalar run() { return getRandomBits<Scalar>(kTotalBits); } +}; + +// random implementation for a built-in signed integer type +template <typename Scalar> +struct random_int_impl<Scalar, true, true> { + static constexpr int kTotalBits = sizeof(Scalar) * CHAR_BIT; + using BitsType = typename make_unsigned<Scalar>::type; + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { + if (y <= x) return x; + // Avoid overflow by representing `range` as an unsigned type + BitsType range = static_cast<BitsType>(y) - static_cast<BitsType>(x); + BitsType randomBits = random_int_impl<BitsType>::run(0, range); + // Avoid overflow in the case where `x` is negative and there is a large range so + // `randomBits` would also be negative if cast to `Scalar` first. + Scalar result = static_cast<Scalar>(static_cast<BitsType>(x) + randomBits); + return result; + } + static EIGEN_DEVICE_FUNC inline Scalar run() { return static_cast<Scalar>(getRandomBits<BitsType>(kTotalBits)); } +}; + +// todo: custom integers +template <typename Scalar, bool IsSigned> +struct random_int_impl<Scalar, IsSigned, false> { + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar&, const Scalar&) { return run(); } + static EIGEN_DEVICE_FUNC inline Scalar run() { + eigen_assert(std::false_type::value && "RANDOM FOR CUSTOM INTEGERS NOT YET SUPPORTED"); + return Scalar(0); + } +}; + +template <typename Scalar> +struct random_default_impl<Scalar, false, true> : random_int_impl<Scalar> {}; + +template <> +struct random_impl<bool> { + static EIGEN_DEVICE_FUNC inline bool run(const bool& x, const bool& y) { + if (y <= x) return x; + return run(); + } + static EIGEN_DEVICE_FUNC inline bool run() { return getRandomBits<unsigned>(1) ? true : false; } +}; + +template <typename Scalar> +struct random_default_impl<Scalar, true, false> { + typedef typename NumTraits<Scalar>::Real RealScalar; + using Impl = random_impl<RealScalar>; + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y, int numRandomBits) { + return Scalar(Impl::run(x.real(), y.real(), numRandomBits), Impl::run(x.imag(), y.imag(), numRandomBits)); + } + static EIGEN_DEVICE_FUNC inline Scalar run(const Scalar& x, const Scalar& y) { + return Scalar(Impl::run(x.real(), y.real()), Impl::run(x.imag(), y.imag())); + } + static EIGEN_DEVICE_FUNC inline Scalar run(int numRandomBits) { + return Scalar(Impl::run(numRandomBits), Impl::run(numRandomBits)); + } + static EIGEN_DEVICE_FUNC inline Scalar run() { return Scalar(Impl::run(), Impl::run()); } +}; + +} // namespace internal +} // namespace Eigen + +#endif // EIGEN_RANDOM_IMPL_H
diff --git a/test/AnnoyingScalar.h b/test/AnnoyingScalar.h index 637fdbf..00a20c7 100644 --- a/test/AnnoyingScalar.h +++ b/test/AnnoyingScalar.h
@@ -184,19 +184,6 @@ return *x.v; } -template <> -struct random_impl<AnnoyingScalar> { - using Impl = random_impl<float>; - static EIGEN_DEVICE_FUNC inline AnnoyingScalar run(const AnnoyingScalar& x, const AnnoyingScalar& y) { - float result = Impl::run(*x.v, *y.v); - return AnnoyingScalar(result); - } - static EIGEN_DEVICE_FUNC inline AnnoyingScalar run() { - float result = Impl::run(); - return AnnoyingScalar(result); - } -}; - } // namespace internal } // namespace Eigen
diff --git a/test/MovableScalar.h b/test/MovableScalar.h index 56a873e..c8bf546 100644 --- a/test/MovableScalar.h +++ b/test/MovableScalar.h
@@ -26,24 +26,10 @@ operator Scalar() const { return this->size() > 0 ? this->back() : Scalar(); } }; -template <> -struct NumTraits<MovableScalar<float>> : GenericNumTraits<float> {}; - -namespace internal { -template <typename T> -struct random_impl<MovableScalar<T>> { - using MoveableT = MovableScalar<T>; - using Impl = random_impl<T>; - static EIGEN_DEVICE_FUNC inline MoveableT run(const MoveableT& x, const MoveableT& y) { - T result = Impl::run(x, y); - return MoveableT(result); - } - static EIGEN_DEVICE_FUNC inline MoveableT run() { - T result = Impl::run(); - return MoveableT(result); - } +template <typename Scalar> +struct NumTraits<MovableScalar<Scalar>> : GenericNumTraits<Scalar> { + enum { RequireInitialization = 1 }; }; -} // namespace internal } // namespace Eigen
diff --git a/test/SafeScalar.h b/test/SafeScalar.h index 4f4da56..33a54c5 100644 --- a/test/SafeScalar.h +++ b/test/SafeScalar.h
@@ -4,43 +4,30 @@ class SafeScalar { public: SafeScalar() : initialized_(false) {} - SafeScalar(const SafeScalar& other) { *this = other; } - SafeScalar& operator=(const SafeScalar& other) { - val_ = T(other); - initialized_ = true; - return *this; - } - SafeScalar(T val) : val_(val), initialized_(true) {} - SafeScalar& operator=(T val) { - val_ = val; - initialized_ = true; - } + SafeScalar(const T& val) : val_(val), initialized_(true) {} + + template <typename Source> + explicit SafeScalar(const Source& val) : SafeScalar(T(val)) {} operator T() const { VERIFY(initialized_ && "Uninitialized access."); return val_; } + template <typename Target> + explicit operator Target() const { + return Target(this->operator T()); + } + private: T val_; bool initialized_; }; namespace Eigen { -namespace internal { template <typename T> -struct random_impl<SafeScalar<T>> { - using SafeT = SafeScalar<T>; - using Impl = random_impl<T>; - static EIGEN_DEVICE_FUNC inline SafeT run(const SafeT& x, const SafeT& y) { - T result = Impl::run(x, y); - return SafeT(result); - } - static EIGEN_DEVICE_FUNC inline SafeT run() { - T result = Impl::run(); - return SafeT(result); - } +struct NumTraits<SafeScalar<T>> : GenericNumTraits<T> { + enum { RequireInitialization = 1 }; }; -} // namespace internal -} // namespace Eigen +} // namespace Eigen \ No newline at end of file
diff --git a/test/rand.cpp b/test/rand.cpp index 6a7c316..4131f38 100644 --- a/test/rand.cpp +++ b/test/rand.cpp
@@ -9,6 +9,10 @@ #include <cstdlib> #include "main.h" +#include "SafeScalar.h" + +// SafeScalar<T> is used to simulate custom Scalar types, which use a more generalized approach to generate random +// numbers // For GCC-6, if this function is inlined then there seems to be an optimization // bug that triggers a failure. This failure goes away if you access `r` in @@ -25,15 +29,28 @@ template <typename Scalar> void check_all_in_range(Scalar x, Scalar y) { - Array<int, 1, Dynamic> mask(y - x + 1); - mask.fill(0); - int64_t n = (y - x + 1) * 32; - for (int64_t k = 0; k < n; ++k) { - mask(check_in_range(x, y) - x)++; - } + constexpr int repeats = 32; + uint64_t count = static_cast<uint64_t>(y) - static_cast<uint64_t>(x) + 1; + ArrayX<bool> mask(count); + // ensure that `count` does not overflow the return type of `mask.size()` + VERIFY(count == static_cast<uint64_t>(mask.size())); + mask.setConstant(false); + for (uint64_t k = 0; k < count; k++) + for (int repeat = 0; repeat < repeats; repeat++) { + Scalar r = check_in_range(x, y); + Index i = static_cast<Index>(r) - static_cast<Index>(x); + mask(i) = true; + } for (Index i = 0; i < mask.size(); ++i) - if (mask(i) == 0) std::cout << "WARNING: value " << x + i << " not reached." << std::endl; - VERIFY((mask > 0).all()); + if (mask(i) == false) std::cout << "WARNING: value " << x + i << " not reached." << std::endl; + VERIFY(mask.cwiseEqual(true).all()); +} + +template <typename Scalar> +void check_all_in_range() { + const Scalar x = NumTraits<Scalar>::lowest(); + const Scalar y = NumTraits<Scalar>::highest(); + check_all_in_range(x, y); } template <typename Scalar, typename EnableIf = void> @@ -66,10 +83,16 @@ double bin_width_; }; +// helper class to avoid extending std:: namespace +template <typename T> +struct get_range_type : internal::make_unsigned<T> {}; +template <typename T> +struct get_range_type<SafeScalar<T>> : internal::make_unsigned<T> {}; + template <typename Scalar> class HistogramHelper<Scalar, std::enable_if_t<Eigen::NumTraits<Scalar>::IsInteger>> { public: - using RangeType = typename Eigen::internal::make_unsigned<Scalar>::type; + using RangeType = typename get_range_type<Scalar>::type; HistogramHelper(int nbins) : HistogramHelper(Eigen::NumTraits<Scalar>::lowest(), Eigen::NumTraits<Scalar>::highest(), nbins) {} HistogramHelper(Scalar lower, Scalar upper, int nbins) @@ -109,38 +132,59 @@ template <typename Scalar> void check_histogram(Scalar x, Scalar y, int bins) { + constexpr int repeats = 10000; + double count = double(bins) * double(repeats); Eigen::VectorXd hist = Eigen::VectorXd::Zero(bins); HistogramHelper<Scalar> hist_helper(x, y, bins); - int64_t n = static_cast<int64_t>(bins) * 10000; // Approx 10000 per bin. - for (int64_t k = 0; k < n; ++k) { - Scalar r = check_in_range(x, y); - int bin = hist_helper.bin(r); - hist(bin)++; - } - // Normalize bins by probability. + for (int k = 0; k < bins; k++) + for (int repeat = 0; repeat < repeats; repeat++) { + Scalar r = check_in_range(x, y); + int bin = hist_helper.bin(r); + hist(bin)++; + } + // Normalize bins by probability. + hist /= count; for (int i = 0; i < bins; ++i) { - hist(i) = hist(i) / n / hist_helper.uniform_bin_probability(i); + hist(i) = hist(i) / hist_helper.uniform_bin_probability(i); } VERIFY(((hist.array() - 1.0).abs() < 0.05).all()); } template <typename Scalar> void check_histogram(int bins) { + constexpr int repeats = 10000; + double count = double(bins) * double(repeats); Eigen::VectorXd hist = Eigen::VectorXd::Zero(bins); HistogramHelper<Scalar> hist_helper(bins); - int64_t n = static_cast<int64_t>(bins) * 10000; // Approx 10000 per bin. - for (int64_t k = 0; k < n; ++k) { - Scalar r = Eigen::internal::random<Scalar>(); - int bin = hist_helper.bin(r); - hist(bin)++; - } - // Normalize bins by probability. + for (int k = 0; k < bins; k++) + for (int repeat = 0; repeat < repeats; repeat++) { + Scalar r = Eigen::internal::random<Scalar>(); + int bin = hist_helper.bin(r); + hist(bin)++; + } + // Normalize bins by probability. + hist /= count; for (int i = 0; i < bins; ++i) { - hist(i) = hist(i) / n / hist_helper.uniform_bin_probability(i); + hist(i) = hist(i) / hist_helper.uniform_bin_probability(i); } VERIFY(((hist.array() - 1.0).abs() < 0.05).all()); } +template <> +void check_histogram<bool>(int) { + constexpr int bins = 2; + constexpr int repeats = 10000; + double count = double(bins) * double(repeats); + double true_count = 0.0; + for (int k = 0; k < bins; k++) + for (int repeat = 0; repeat < repeats; repeat++) { + bool r = Eigen::internal::random<bool>(); + if (r) true_count += 1.0; + } + double p = true_count / count; + VERIFY(numext::abs(p - 0.5) < 0.05); +} + EIGEN_DECLARE_TEST(rand) { int64_t int64_ref = NumTraits<int64_t>::highest() / 10; // the minimum guarantees that these conversions are safe @@ -191,14 +235,16 @@ CALL_SUBTEST_7(check_all_in_range<int8_t>(-11 - int8t_offset, -11)); CALL_SUBTEST_7(check_all_in_range<int8_t>(-126, -126 + int8t_offset)); CALL_SUBTEST_7(check_all_in_range<int8_t>(126 - int8t_offset, 126)); - CALL_SUBTEST_7(check_all_in_range<int8_t>(-126, 126)); + CALL_SUBTEST_7(check_all_in_range<int8_t>()); + CALL_SUBTEST_7(check_all_in_range<uint8_t>()); CALL_SUBTEST_8(check_all_in_range<int16_t>(11, 11)); CALL_SUBTEST_8(check_all_in_range<int16_t>(11, 11 + int16t_offset)); CALL_SUBTEST_8(check_all_in_range<int16_t>(-5, 5)); CALL_SUBTEST_8(check_all_in_range<int16_t>(-11 - int16t_offset, -11)); CALL_SUBTEST_8(check_all_in_range<int16_t>(-24345, -24345 + int16t_offset)); - CALL_SUBTEST_8(check_all_in_range<int16_t>(24345, 24345 + int16t_offset)); + CALL_SUBTEST_8(check_all_in_range<int16_t>()); + CALL_SUBTEST_8(check_all_in_range<uint16_t>()); CALL_SUBTEST_9(check_all_in_range<int32_t>(11, 11)); CALL_SUBTEST_9(check_all_in_range<int32_t>(11, 11 + g_repeat)); @@ -223,6 +269,7 @@ CALL_SUBTEST_11(check_histogram<int32_t>(-RAND_MAX + 10, -int64_t(RAND_MAX) + 10 + bins * (2 * int64_t(RAND_MAX) / bins) - 1, bins)); + CALL_SUBTEST_12(check_histogram<bool>(/*bins=*/2)); CALL_SUBTEST_12(check_histogram<uint8_t>(/*bins=*/16)); CALL_SUBTEST_12(check_histogram<uint16_t>(/*bins=*/1024)); CALL_SUBTEST_12(check_histogram<uint32_t>(/*bins=*/1024)); @@ -238,10 +285,16 @@ CALL_SUBTEST_14(check_histogram<long double>(-10.0L, 10.0L, /*bins=*/1024)); CALL_SUBTEST_14(check_histogram<half>(half(-10.0f), half(10.0f), /*bins=*/512)); CALL_SUBTEST_14(check_histogram<bfloat16>(bfloat16(-10.0f), bfloat16(10.0f), /*bins=*/64)); + CALL_SUBTEST_14(check_histogram<SafeScalar<float>>(-10.0f, 10.0f, /*bins=*/1024)); + CALL_SUBTEST_14(check_histogram<SafeScalar<half>>(half(-10.0f), half(10.0f), /*bins=*/512)); + CALL_SUBTEST_14(check_histogram<SafeScalar<bfloat16>>(bfloat16(-10.0f), bfloat16(10.0f), /*bins=*/64)); CALL_SUBTEST_15(check_histogram<float>(/*bins=*/1024)); CALL_SUBTEST_15(check_histogram<double>(/*bins=*/1024)); CALL_SUBTEST_15(check_histogram<long double>(/*bins=*/1024)); CALL_SUBTEST_15(check_histogram<half>(/*bins=*/512)); CALL_SUBTEST_15(check_histogram<bfloat16>(/*bins=*/64)); + CALL_SUBTEST_15(check_histogram<SafeScalar<float>>(/*bins=*/1024)); + CALL_SUBTEST_15(check_histogram<SafeScalar<half>>(/*bins=*/512)); + CALL_SUBTEST_15(check_histogram<SafeScalar<bfloat16>>(/*bins=*/64)); }