diff --git a/Core/HardParticle/FormFactorPolyhedron.cpp b/Core/HardParticle/FormFactorPolyhedron.cpp index 05476d484cffe34d55c4a6f54f659a40dca41a51..ad14bd70bf306a9cf0386ce1a2a64a919eaa1ecf 100644 --- a/Core/HardParticle/FormFactorPolyhedron.cpp +++ b/Core/HardParticle/FormFactorPolyhedron.cpp @@ -27,6 +27,7 @@ namespace { const complex_t I = {0.,1.}; const double eps = 2e-16; + constexpr auto ReciprocalFactorialArray = Precomputed::GenerateReciprocalFactorialArray<171>(); } double PolyhedralFace::qpa_limit_series = 3e-2; @@ -64,20 +65,20 @@ complex_t PolyhedralEdge::contrib(int M, cvector_t qpa, complex_t qrperp) const if( M&1 ) // M is odd return 0.; else - return Precomputed::ReciprocalFactorialArray[M] * ( pow(u, M)/(M+1.) - pow(v1, M) ); + return ReciprocalFactorialArray[M] * ( pow(u, M)/(M+1.) - pow(v1, M) ); } complex_t ret = 0; // the l=0 term, minus (qperp.R)^M, which cancels under the sum over E*contrib() if ( v1==0. ) { - ret = Precomputed::ReciprocalFactorialArray[M] * pow(v2, M); + ret = ReciprocalFactorialArray[M] * pow(v2, M); } else if ( v2==0. ) { ; // leave ret=0 } else { // binomial expansion for( int mm=1; mm<=M; ++mm ) { complex_t term = - Precomputed::ReciprocalFactorialArray[mm] * - Precomputed::ReciprocalFactorialArray[M-mm] * + ReciprocalFactorialArray[mm] * + ReciprocalFactorialArray[M-mm] * pow(v2, mm) * pow(v1, M-mm); ret += term; #ifdef POLYHEDRAL_DIAGNOSTIC @@ -90,8 +91,8 @@ complex_t PolyhedralEdge::contrib(int M, cvector_t qpa, complex_t qrperp) const return ret; for( int l=1; l<=M/2; ++l ) { complex_t term = - Precomputed::ReciprocalFactorialArray[M-2*l] * - Precomputed::ReciprocalFactorialArray[2*l+1] * + ReciprocalFactorialArray[M-2*l] * + ReciprocalFactorialArray[2*l+1] * pow(u, 2*l) * pow(v, M-2*l); ret += term; #ifdef POLYHEDRAL_DIAGNOSTIC @@ -250,7 +251,7 @@ complex_t PolyhedralFace::ff_n( int n, cvector_t q ) const decompose_q( q, qperp, qpa ); double qpa_mag2 = qpa.mag2(); if ( qpa_mag2==0. ) { - return qn * pow(qperp*m_rperp, n) * m_area / Precomputed::FactorialArray[n]; + return qn * pow(qperp*m_rperp, n) * m_area * ReciprocalFactorialArray[n]; } else if ( sym_S2 ) { return qn * ( ff_n_core( n, qpa, qperp ) + ff_n_core( n, -qpa, qperp ) ) / qpa_mag2; } else { diff --git a/Core/Tools/Precomputed.h b/Core/Tools/Precomputed.h index 78a309bc9dc871e626cce617c842d3b31ec44109..ee498060e8d498bf9eefbdd5ac34650c16ffe2c5 100644 --- a/Core/Tools/Precomputed.h +++ b/Core/Tools/Precomputed.h @@ -22,48 +22,32 @@ #include <utility> #include <vector> -//! Compile-time generated std::arrays of factorials and their reciprocals +//! Compile-time generated std::array of reciprocal factorials namespace Precomputed { template<size_t N> -struct Factorial +struct ReciprocalFactorial { - static constexpr double value = N*Factorial<N - 1>::value; + static constexpr double value = ReciprocalFactorial<N-1>::value/N; }; template<> -struct Factorial<0> +struct ReciprocalFactorial<0> { static constexpr double value = 1.0; }; -template<size_t N> -struct ReciprocalFactorial -{ - static constexpr double value = 1.0/Factorial<N>::value; -}; - template<template<size_t> class F, size_t... I> constexpr std::array<double, sizeof...(I)> GenerateArrayHelper(std::index_sequence<I...>) { return { F<I>::value... }; }; -template<size_t N, typename Indices = std::make_index_sequence<N>> -constexpr std::array<double, N> GenerateFactorialArray() -{ - return GenerateArrayHelper<Factorial>(Indices{}); -}; - template<size_t N, typename Indices = std::make_index_sequence<N>> constexpr std::array<double, N> GenerateReciprocalFactorialArray() { return GenerateArrayHelper<ReciprocalFactorial>(Indices{}); }; - -static constexpr auto FactorialArray = GenerateFactorialArray<171>(); -static constexpr auto ReciprocalFactorialArray = GenerateReciprocalFactorialArray<171>(); - } // namespace Precomputed #endif // PRECOMPUTED_H diff --git a/Tests/UnitTests/Core/Detector/PrecomputedTest.h b/Tests/UnitTests/Core/Detector/PrecomputedTest.h index 3452d71a64673997da9ade1b872a455738a9b57a..ec5be695daab38b28bb8829b0e8855d1295cca61 100644 --- a/Tests/UnitTests/Core/Detector/PrecomputedTest.h +++ b/Tests/UnitTests/Core/Detector/PrecomputedTest.h @@ -1,23 +1,28 @@ #include "Precomputed.h" #include <memory> +namespace { + constexpr auto ReciprocalFactorialArray = Precomputed::GenerateReciprocalFactorialArray<171>(); +} + class PrecomputedTest : public ::testing::Test { public: }; -TEST_F(PrecomputedTest, Factorial) +TEST_F(PrecomputedTest, ReciprocalFactorial) { const double eps = 2.3e-16; // about the machine precision - EXPECT_TRUE(Precomputed::FactorialArray.size()>150); - EXPECT_DOUBLE_EQ(Precomputed::FactorialArray[0], 1.); - EXPECT_DOUBLE_EQ(Precomputed::FactorialArray[1], 1.); - EXPECT_DOUBLE_EQ(Precomputed::FactorialArray[2], 2.); - EXPECT_DOUBLE_EQ(Precomputed::FactorialArray[3], 6.); + EXPECT_TRUE(ReciprocalFactorialArray.size()>150); + EXPECT_DOUBLE_EQ(ReciprocalFactorialArray[0], 1.); + EXPECT_DOUBLE_EQ(ReciprocalFactorialArray[1], 1.); + EXPECT_DOUBLE_EQ(ReciprocalFactorialArray[2], 0.5); + EXPECT_DOUBLE_EQ(ReciprocalFactorialArray[3], 1.0/6); /* the following disabled because tgamma is too unprecise under old versions of glibc (at leat up to 2.12, but less than 2.22) for( size_t k=4; k<precomputed.factorial.size(); ++k ) EXPECT_NEAR(precomputed.factorial[k], tgamma(k+1.), 12*eps*tgamma(k+1.) ); */ - EXPECT_NEAR(Precomputed::FactorialArray[150], 5.71338395644585459e262, 4*eps*Precomputed::FactorialArray[150]); + EXPECT_NEAR(ReciprocalFactorialArray[150], 1.75027620692601519e-263, + 4*eps*ReciprocalFactorialArray[150]); }