support factorials fac>20! with long double.

evaleev · Feb 23, 2024 · 4708c2a · 4708c2a
1 parent cbe5144
commit 4708c2a
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 8 deletions.
diff --git a/include/libint2/shell.h b/include/libint2/shell.h
@@ -60,6 +60,12 @@ static constexpr std::array<int64_t, 21> fac = {{1LL,
                                                  6402373705728000LL,
                                                  121645100408832000LL,
                                                  2432902008176640000LL}};
+
+/// fac_ld[k] = k! if 21<=k<=25
+static constexpr std::array<long double, 5> fac_ld = {
+    {5.109094217170944e+19, 1.12400072777760768e+21, 2.585201673888497664e+22,
+     6.2044840173323943936e+23, 1.5511210043330985984e+25}};
+
 /// df_Kminus1[k] = (k-1)!!
 static constexpr std::array<int64_t, 31> df_Kminus1 = {{1LL,
                                                         1LL,

diff --git a/include/libint2/solidharmonics.h b/include/libint2/solidharmonics.h
@@ -83,8 +83,8 @@ class SolidHarmonicsCoefficients {
   static const SolidHarmonicsCoefficients& instance(unsigned int l) {
     static std::vector<SolidHarmonicsCoefficients> shg_coefs(
         SolidHarmonicsCoefficients::CtorHelperIter(0),
-        SolidHarmonicsCoefficients::CtorHelperIter(11));
-    assert(l <= 10);  // see coeff() for explanation of the upper limit on l
+        SolidHarmonicsCoefficients::CtorHelperIter(13));
+    assert(l <= 12);  // see coeff() for explanation of the upper limit on l
     return shg_coefs[l];
   }
 
@@ -111,6 +111,7 @@ class SolidHarmonicsCoefficients {
     using libint2::math::bc;
     using libint2::math::df_Kminus1;
     using libint2::math::fac;
+    using libint2::math::fac_ld;
 
     auto abs_m = std::abs(m);
     if ((lx + ly - abs_m) % 2) return 0.0;
@@ -127,12 +128,21 @@ class SolidHarmonicsCoefficients {
     auto i = abs_m - lx;
     if (comp != parity(abs(i))) return 0.0;
 
-    assert(l <= 10);  // libint2::math::fac[] is only defined up to 20
-    Real pfac =
-        sqrt(((Real(fac[2 * lx]) * Real(fac[2 * ly]) * Real(fac[2 * lz])) /
-              fac[2 * l]) *
-             ((Real(fac[l - abs_m])) / (fac[l])) * (Real(1) / fac[l + abs_m]) *
-             (Real(1) / (fac[lx] * fac[ly] * fac[lz])));
+    assert(l <= 12);  // libint2::math::fac[] is only defined up to 20
+    Real pfac;
+    if (l <= 10)
+      pfac = sqrt(((Real(fac[2 * lx]) * Real(fac[2 * ly]) * Real(fac[2 * lz])) /
+                   fac[2 * l]) *
+                  ((Real(fac[l - abs_m])) / (fac[l])) *
+                  (Real(1) / fac[l + abs_m]) *
+                  (Real(1) / (fac[lx] * fac[ly] * fac[lz])));
+    else
+      pfac = sqrt(((Real(fac[2 * lx]) * Real(fac[2 * ly]) * Real(fac[2 * lz])) /
+                   Real(fac_ld[2 * l - 21])) *
+                  ((Real(fac[l - abs_m])) / (fac[l])) *
+                  (Real(1) / fac[l + abs_m]) *
+                  (Real(1) / (fac[lx] * fac[ly] * fac[lz])));
+
     /*  pfac = sqrt(fac[l-abs_m]/(fac[l]*fac[l]*fac[l+abs_m]));*/
     pfac /= (1L << l);
     if (m < 0)