MATH-1293

New inner class "CombinatoricsUtils.FactorialLog" that allows to precompute a configurable number of values.

Thanks to Aleksei Dievskii for the original feature request and implementation.
The proposed code was further modified to enhance usage flexibility and performance and to include additional unit tests.

src/changes/changes.xml

@@ -54,6 +54,10 @@ If the output is not quite correct, check for invisible trailing spaces!
     </release>
 
     <release version="4.0" date="XXXX-XX-XX" description="">
+      <action dev="erans" type="add" issue="MATH-1295" due-to="Aleksei Dievskii">
+        "FactorialLog": Cache-based computation of the "log factorial" function (implemented
+        as an inner class of "CombinatoricsUtils" in package "o.a.c.m.util").
+      </action>
       <action dev="erans" type="fix" issue="MATH-1295" due-to="Luke Lindsay"> <!-- backported to 3.6 -->
         Increased default value for number of allowed evaluations in
         "o.a.c.m.optim.univariate.BracketFinder".

src/main/java/org/apache/commons/math4/util/CombinatoricsUtils.java

@@ -23,6 +23,7 @@ import org.apache.commons.math4.exception.MathArithmeticException;
 import org.apache.commons.math4.exception.NotPositiveException;
 import org.apache.commons.math4.exception.NumberIsTooLargeException;
 import org.apache.commons.math4.exception.util.LocalizedFormats;
+import org.apache.commons.math4.special.Gamma;
 
 /**
  * Combinatorial utilities.
@@ -307,6 +308,15 @@ public final class CombinatoricsUtils {
         return FastMath.floor(FastMath.exp(CombinatoricsUtils.factorialLog(n)) + 0.5);
     }
 
+    /**
+     * Default implementation of {@link #factorialLog(int)} method:
+     * <ul>
+     *  <li>No pre-computation</li>
+     *  <li>No cache allocation</li>
+     * </ul>
+     */
+    private static final FactorialLog FACTORIAL_LOG_NO_CACHE = FactorialLog.create();
+
     /**
      * Compute the natural logarithm of the factorial of {@code n}.
      *
@@ -315,18 +325,7 @@ public final class CombinatoricsUtils {
      * @throws NotPositiveException if {@code n < 0}.
      */
     public static double factorialLog(final int n) throws NotPositiveException {
-        if (n < 0) {
+        return FACTORIAL_LOG_NO_CACHE.value(n);
-            throw new NotPositiveException(LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
-                                           n);
-        }
-        if (n < 21) {
-            return FastMath.log(FACTORIALS[n]);
-        }
-        double logSum = 0;
-        for (int i = 2; i <= n; i++) {
-            logSum += FastMath.log(i);
-        }
-        return logSum;
     }
 
     /**
@@ -459,4 +458,95 @@ public final class CombinatoricsUtils {
             throw new NotPositiveException(LocalizedFormats.BINOMIAL_NEGATIVE_PARAMETER, n);
         }
     }
+
+    /**
+     * Class for computing the natural logarithm of the factorial of {@code n}.
+     * It allows to allocate a cache of precomputed values.
+     * In case of cache miss, computation is preformed by a call to
+     * {@link Gamma#logGamma(double)}.
+     */
+    public static final class FactorialLog {
+        /**
+         * Precomputed values of the function:
+         * {@code LOG_FACTORIALS[i] = log(i!)}.
+         */
+        private final double[] LOG_FACTORIALS;
+
+        /**
+         * Creates an instance, reusing the already computed values if available.
+         *
+         * @param numValues Number of values of the function to compute.
+         * @param cache Existing cache.
+         * @throw NotPositiveException if {@code n < 0}.
+         */
+        private FactorialLog(int numValues,
+                             double[] cache) {
+            if (numValues < 0) {
+                throw new NotPositiveException(numValues);
+            }
+
+            LOG_FACTORIALS = new double[numValues];
+
+            final int beginCopy = 2;
+            final int endCopy = cache == null || cache.length <= beginCopy ?
+                beginCopy : cache.length <= numValues ?
+                cache.length : numValues;
+
+            // Copy available values.
+            for (int i = beginCopy; i < endCopy; i++) {
+                LOG_FACTORIALS[i] = cache[i];
+            }
+
+            // Precompute.
+            for (int i = endCopy; i < numValues; i++) {
+                LOG_FACTORIALS[i] = LOG_FACTORIALS[i - 1] + FastMath.log(i);
+            }
+        }
+
+        /**
+         * Creates an instance with no precomputed values.
+         */
+        public static FactorialLog create() {
+            return new FactorialLog(0, null);
+        }
+
+        /**
+         * Creates an instance with the specified cache size.
+         *
+         * @param cacheSize Number of precomputed values of the function.
+         * @return an new instance where {@code cacheSize} values have been
+         * precomputed.
+         * @throw NotPositiveException if {@code n < 0}.
+         */
+        public FactorialLog withCache(final int cacheSize) {
+            return new FactorialLog(cacheSize, LOG_FACTORIALS);
+        }
+
+        /**
+         * Computes {@code log(n!)}.
+         *
+         * @param n Argument.
+         * @return {@code log(n!)}.
+         * @throws NotPositiveException if {@code n < 0}.
+         */
+        public double value(final int n) {
+            if (n < 0) {
+                throw new NotPositiveException(LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
+                                               n);
+            }
+
+            // Use cache of precomputed values.
+            if (n < LOG_FACTORIALS.length) {
+                return LOG_FACTORIALS[n];
+            }
+
+            // Use cache of precomputed factorial values.
+            if (n < FACTORIALS.length) {
+                return FastMath.log(FACTORIALS[n]);
+            }
+
+            // Delegate.
+            return Gamma.logGamma(n + 1);
+        }
+    }
 }

src/test/java/org/apache/commons/math4/util/FactorialLogTest.java

@@ -0,0 +1,111 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math4.util;
+
+import org.apache.commons.math4.exception.NotPositiveException;
+import org.apache.commons.math4.special.Gamma;
+
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test cases for the {@link CombinatoricsUtils.FactorialLog} class.
+ */
+public class FactorialLogTest {
+
+    @Test(expected=NotPositiveException.class)
+    public void testPrecondition1() {
+        CombinatoricsUtils.FactorialLog.create().withCache(-1);
+    }
+
+    @Test(expected=NotPositiveException.class)
+    public void testNonPositiveArgument() {
+        final CombinatoricsUtils.FactorialLog f = CombinatoricsUtils.FactorialLog.create();
+        f.value(-1);
+    }
+
+    @Test
+    public void testDelegation() {
+        final CombinatoricsUtils.FactorialLog f = CombinatoricsUtils.FactorialLog.create();
+
+        // Starting at 21 because for smaller arguments, there is no delegation to the
+        // "Gamma" class.
+        for (int i = 21; i < 10000; i++) {
+            final double expected = Gamma.logGamma(i + 1);
+            Assert.assertEquals(i + "! ",
+                                expected, f.value(i), 0d);
+        }
+    }
+
+    @Test
+    public void testCompareDirectWithoutCache() {
+        // This test shows that delegating to the "Gamma" class will also lead to a
+        // less accurate result.
+
+        final int max = 100;
+        final CombinatoricsUtils.FactorialLog f = CombinatoricsUtils.FactorialLog.create();
+
+        for (int i = 0; i < max; i++) {
+            final double expected = factorialLog(i);
+            Assert.assertEquals(i + "! ",
+                                expected, f.value(i), 2 * Math.ulp(expected));
+        }
+    }
+
+    @Test
+    public void testCompareDirectWithCache() {
+        final int max = 1000;
+        final CombinatoricsUtils.FactorialLog f = CombinatoricsUtils.FactorialLog.create().withCache(max);
+
+        for (int i = 0; i < max; i++) {
+            final double expected = factorialLog(i);
+            Assert.assertEquals(i + "! ",
+                                expected, f.value(i), 0d);
+        }
+    }
+
+    @Test
+    public void testCacheIncrease() {
+        final int max = 100;
+        final CombinatoricsUtils.FactorialLog f1 = CombinatoricsUtils.FactorialLog.create().withCache(max);
+        final CombinatoricsUtils.FactorialLog f2 = f1.withCache(2 * max);
+
+        final int val = max + max / 2;
+        final double expected = factorialLog(val);
+        Assert.assertEquals(expected, f2.value(val), 0d);
+    }
+
+    @Test
+    public void testCacheDecrease() {
+        final int max = 100;
+        final CombinatoricsUtils.FactorialLog f1 = CombinatoricsUtils.FactorialLog.create().withCache(max);
+        final CombinatoricsUtils.FactorialLog f2 = f1.withCache(max / 2);
+
+        final int val = max / 4;
+        final double expected = factorialLog(val);
+        Assert.assertEquals(expected, f2.value(val), 0d);
+    }
+
+    // Direct implementation.
+    private double factorialLog(final int n) {
+        double logSum = 0;
+        for (int i = 2; i <= n; i++) {
+            logSum += FastMath.log(i);
+        }
+        return logSum;
+    }
+}