Remove FactorialList

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 {

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

Tests/UnitTests/Core/Detector/PrecomputedTest.h

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