diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 28ebacc2e..d4fac5506 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -54,6 +54,10 @@ If the output is not quite correct, check for invisible trailing spaces! + + "FactorialLog": Cache-based computation of the "log factorial" function (implemented + as an inner class of "CombinatoricsUtils" in package "o.a.c.m.util"). + Increased default value for number of allowed evaluations in "o.a.c.m.optim.univariate.BracketFinder". diff --git a/src/main/java/org/apache/commons/math4/util/CombinatoricsUtils.java b/src/main/java/org/apache/commons/math4/util/CombinatoricsUtils.java index 0ca2e7527..019a80f70 100644 --- a/src/main/java/org/apache/commons/math4/util/CombinatoricsUtils.java +++ b/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: + * + */ + 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) { - 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; + return FACTORIAL_LOG_NO_CACHE.value(n); } /** @@ -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); + } + } } diff --git a/src/test/java/org/apache/commons/math4/util/FactorialLogTest.java b/src/test/java/org/apache/commons/math4/util/FactorialLogTest.java new file mode 100644 index 000000000..980a834b2 --- /dev/null +++ b/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; + } +}