MAINT: Increase logfactorial lookup table size

zoj613 · Mar 26, 2021 · d3ef3d3 · d3ef3d3
1 parent 2af0edb
commit d3ef3d3
Showing 1 changed file with 30 additions and 7 deletions.
diff --git a/src/pgm_common.h b/src/pgm_common.h
@@ -125,9 +125,8 @@ pgm_erfc(double x)
 NPY_INLINE double
 pgm_lgamma(double z)
 {
-    /* lookup table for integer values of log-gamma function where 1<=z<=126
-     * courtesy of NumPy devs. */
-    static const double logfactorial[126] = {
+    /* lookup table for integer values of log-gamma function where 1<=z<=200 */
+    static const double logfactorial[200] = {
         0.000000000000000, 0.0000000000000000, 0.69314718055994529,
         1.791759469228055, 3.1780538303479458, 4.7874917427820458,
         6.5792512120101012, 8.5251613610654147, 10.604602902745251,
@@ -169,12 +168,36 @@ pgm_lgamma(double z)
         429.21439186665157, 433.95932399501481, 438.71291418612117,
         443.47508812091894, 448.24577274538461, 453.02489623849613,
         457.81238798127816, 462.60817852687489, 467.4121995716082,
-        472.22438392698058, 477.04466549258564, 481.87297922988796
+        472.22438392698058, 477.04466549258564, 481.87297922988796,
+        486.70926113683936, 491.55344822329801, 496.40547848721764,
+        501.26529089157924, 506.13282534203483, 511.00802266523596,
+        515.89082458782241, 520.78117371604412, 525.67901351599517,
+        530.58428829443358, 535.49694318016952, 540.41692410599762,
+        545.34417779115483, 550.27865172428551, 555.22029414689484,
+        560.16905403727310, 565.12488109487424, 570.08772572513419,
+        575.05753902471020, 580.03427276713080, 585.01787938883899,
+        590.00831197561786, 595.00552424938201, 600.00947055532743,
+        605.02010584942377, 610.03738568623862, 615.06126620708494,
+        620.09170412847732, 625.12865673089095, 630.17208184781020,
+        635.22193785505965, 640.27818366040810, 645.34077869343503,
+        650.40968289565524, 655.48485671088906, 660.56626107587351,
+        665.65385741110595, 670.74760761191271, 675.84747403973688,
+        680.95341951363753, 686.06540730199413, 691.18340111441080,
+        696.30736509381404, 701.43726380873704, 706.57306224578736,
+        711.71472580228999, 716.86222027910355, 722.01551187360133,
+        727.17456717281584, 732.33935314673920, 737.50983714177733,
+        742.68598687435122, 747.86777042464337, 753.05515623048404,
+        758.24811308137441, 763.44661011264009, 768.65061679971711,
+        773.86010295255835, 779.07503871016729, 784.29539453524569,
+        789.52114120895885, 794.75224982581346, 799.98869178864345,
+        805.23043880370301, 810.47746287586358, 815.72973630391016,
+        820.98723167593789, 826.24992186484292, 831.51778002390620,
+        836.79077958246978, 842.06889424170038, 847.35209797043842,
+        852.64036500113298, 857.93366982585735,
     };
 
-    size_t zz = (size_t)z;
-    if (z < 127 && z == zz) {
-        return logfactorial[zz - 1];
+    if (z < 201 && z == (size_t)z) {
+        return logfactorial[(size_t)z - 1];
     }
     else if (z > 12) {
         static const double a1 = 0.08333333333333333;  // 1 / 12