CombinatoricsUtils.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.util;

  24. import java.lang.reflect.Array;
  25. import java.util.Iterator;
  26. import java.util.List;
  27. import java.util.Spliterator;
  28. import java.util.Spliterators;
  29. import java.util.concurrent.atomic.AtomicReference;
  30. import java.util.stream.Stream;
  31. import java.util.stream.StreamSupport;

  33. import org.hipparchus.exception.LocalizedCoreFormats;
  34. import org.hipparchus.exception.MathIllegalArgumentException;
  35. import org.hipparchus.exception.MathRuntimeException;
  36. import org.hipparchus.special.Gamma;

  38. /**
  39.  * Combinatorial utilities.
  40.  */
  41. public final class CombinatoricsUtils {

  43.     /** Maximum index of Bell number that fits into a long.
  44.      * @since 2.2
  45.      */
  46.     public static final int MAX_BELL = 25;

  48.     /** All long-representable factorials */
  49.     static final long[] FACTORIALS = {
  50.                        1L,                  1L,                   2L,
  51.                        6L,                 24L,                 120L,
  52.                      720L,               5040L,               40320L,
  53.                   362880L,            3628800L,            39916800L,
  54.                479001600L,         6227020800L,         87178291200L,
  55.            1307674368000L,     20922789888000L,     355687428096000L,
  56.         6402373705728000L, 121645100408832000L, 2432902008176640000L };

  58.     /** Stirling numbers of the second kind. */
  59.     static final AtomicReference<long[][]> STIRLING_S2 = new AtomicReference<> (null);

  61.     /**
  62.      * Default implementation of {@link #factorialLog(int)} method:
  63.      * <ul>
  64.      *  <li>No pre-computation</li>
  65.      *  <li>No cache allocation</li>
  66.      * </ul>
  67.      */
  68.     private static final FactorialLog FACTORIAL_LOG_NO_CACHE = FactorialLog.create();

  70.     /** Bell numbers.
  71.      * @since 2.2
  72.      */
  73.     private static final AtomicReference<long[]> BELL = new AtomicReference<> (null);

  75.     /** Private constructor (class contains only static methods). */
  76.     private CombinatoricsUtils() {}


  79.     /**
  80.      * Returns an exact representation of the <a
  81.      * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
  82.      * Coefficient</a>, "{@code n choose k}", the number of
  83.      * {@code k}-element subsets that can be selected from an
  84.      * {@code n}-element set.
  85.      * <p><Strong>Preconditions</strong>:</p>
  86.      * <ul>
  87.      * <li> {@code 0 <= k <= n } (otherwise
  88.      * {@code MathIllegalArgumentException} is thrown)</li>
  89.      * <li> The result is small enough to fit into a {@code long}. The
  90.      * largest value of {@code n} for which all coefficients are
  91.      * {@code  < Long.MAX_VALUE} is 66. If the computed value exceeds
  92.      * {@code Long.MAX_VALUE} a {@code MathRuntimeException} is
  93.      * thrown.</li>
  94.      * </ul>
  95.      *
  96.      * @param n the size of the set
  97.      * @param k the size of the subsets to be counted
  98.      * @return {@code n choose k}
  99.      * @throws MathIllegalArgumentException if {@code n < 0}.
  100.      * @throws MathIllegalArgumentException if {@code k > n}.
  101.      * @throws MathRuntimeException if the result is too large to be
  102.      * represented by a long integer.
  103.      */
  104.     public static long binomialCoefficient(final int n, final int k)
  105.         throws MathIllegalArgumentException, MathRuntimeException {
  106.         checkBinomial(n, k);
  107.         if ((n == k) || (k == 0)) {
  108.             return 1;
  109.         }
  110.         if ((k == 1) || (k == n - 1)) {
  111.             return n;
  112.         }
  113.         // Use symmetry for large k
  114.         if (k > n / 2) {
  115.             return binomialCoefficient(n, n - k);
  116.         }

  118.         // We use the formula
  119.         // (n choose k) = n! / (n-k)! / k!
  120.         // (n choose k) == ((n-k+1)*...*n) / (1*...*k)
  121.         // which could be written
  122.         // (n choose k) == (n-1 choose k-1) * n / k
  123.         long result = 1;
  124.         if (n <= 61) {
  125.             // For n <= 61, the naive implementation cannot overflow.
  126.             int i = n - k + 1;
  127.             for (int j = 1; j <= k; j++) {
  128.                 result = result * i / j;
  129.                 i++;
  130.             }
  131.         } else if (n <= 66) {
  132.             // For n > 61 but n <= 66, the result cannot overflow,
  133.             // but we must take care not to overflow intermediate values.
  134.             int i = n - k + 1;
  135.             for (int j = 1; j <= k; j++) {
  136.                 // We know that (result * i) is divisible by j,
  137.                 // but (result * i) may overflow, so we split j:
  138.                 // Filter out the gcd, d, so j/d and i/d are integer.
  139.                 // result is divisible by (j/d) because (j/d)
  140.                 // is relative prime to (i/d) and is a divisor of
  141.                 // result * (i/d).
  142.                 final long d = ArithmeticUtils.gcd(i, j);
  143.                 result = (result / (j / d)) * (i / d);
  144.                 i++;
  145.             }
  146.         } else {
  147.             // For n > 66, a result overflow might occur, so we check
  148.             // the multiplication, taking care to not overflow
  149.             // unnecessary.
  150.             int i = n - k + 1;
  151.             for (int j = 1; j <= k; j++) {
  152.                 final long d = ArithmeticUtils.gcd(i, j);
  153.                 result = ArithmeticUtils.mulAndCheck(result / (j / d), i / d);
  154.                 i++;
  155.             }
  156.         }
  157.         return result;
  158.     }

  160.     /**
  161.      * Returns a {@code double} representation of the <a
  162.      * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
  163.      * Coefficient</a>, "{@code n choose k}", the number of
  164.      * {@code k}-element subsets that can be selected from an
  165.      * {@code n}-element set.
  166.      * <p>* <Strong>Preconditions</strong>:</p>
  167.      * <ul>
  168.      * <li> {@code 0 <= k <= n } (otherwise
  169.      * {@code IllegalArgumentException} is thrown)</li>
  170.      * <li> The result is small enough to fit into a {@code double}. The
  171.      * largest value of {@code n} for which all coefficients are &lt;
  172.      * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
  173.      * Double.POSITIVE_INFINITY is returned</li>
  174.      * </ul>
  175.      *
  176.      * @param n the size of the set
  177.      * @param k the size of the subsets to be counted
  178.      * @return {@code n choose k}
  179.      * @throws MathIllegalArgumentException if {@code n < 0}.
  180.      * @throws MathIllegalArgumentException if {@code k > n}.
  181.      * @throws MathRuntimeException if the result is too large to be
  182.      * represented by a long integer.
  183.      */
  184.     public static double binomialCoefficientDouble(final int n, final int k)
  185.         throws MathIllegalArgumentException, MathRuntimeException {
  186.         checkBinomial(n, k);
  187.         if ((n == k) || (k == 0)) {
  188.             return 1d;
  189.         }
  190.         if ((k == 1) || (k == n - 1)) {
  191.             return n;
  192.         }
  193.         if (k > n/2) {
  194.             return binomialCoefficientDouble(n, n - k);
  195.         }
  196.         if (n < 67) {
  197.             return binomialCoefficient(n,k);
  198.         }

  200.         double result = 1d;
  201.         for (int i = 1; i <= k; i++) {
  202.              result *= ((double) (n - k + i)) / i;
  203.         }

  205.         return FastMath.floor(result + 0.5);
  206.     }

  208.     /**
  209.      * Returns the natural {@code log} of the <a
  210.      * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
  211.      * Coefficient</a>, "{@code n choose k}", the number of
  212.      * {@code k}-element subsets that can be selected from an
  213.      * {@code n}-element set.
  214.      * <p>* <Strong>Preconditions</strong>:</p>
  215.      * <ul>
  216.      * <li> {@code 0 <= k <= n } (otherwise
  217.      * {@code MathIllegalArgumentException} is thrown)</li>
  218.      * </ul>
  219.      *
  220.      * @param n the size of the set
  221.      * @param k the size of the subsets to be counted
  222.      * @return {@code n choose k}
  223.      * @throws MathIllegalArgumentException if {@code n < 0}.
  224.      * @throws MathIllegalArgumentException if {@code k > n}.
  225.      * @throws MathRuntimeException if the result is too large to be
  226.      * represented by a long integer.
  227.      */
  228.     public static double binomialCoefficientLog(final int n, final int k)
  229.         throws MathIllegalArgumentException, MathRuntimeException {
  230.         checkBinomial(n, k);
  231.         if ((n == k) || (k == 0)) {
  232.             return 0;
  233.         }
  234.         if ((k == 1) || (k == n - 1)) {
  235.             return FastMath.log(n);
  236.         }

  238.         /*
  239.          * For values small enough to do exact integer computation,
  240.          * return the log of the exact value
  241.          */
  242.         if (n < 67) {
  243.             return FastMath.log(binomialCoefficient(n,k));
  244.         }

  246.         /*
  247.          * Return the log of binomialCoefficientDouble for values that will not
  248.          * overflow binomialCoefficientDouble
  249.          */
  250.         if (n < 1030) {
  251.             return FastMath.log(binomialCoefficientDouble(n, k));
  252.         }

  254.         if (k > n / 2) {
  255.             return binomialCoefficientLog(n, n - k);
  256.         }

  258.         /*
  259.          * Sum logs for values that could overflow
  260.          */
  261.         double logSum = 0;

  263.         // n!/(n-k)!
  264.         for (int i = n - k + 1; i <= n; i++) {
  265.             logSum += FastMath.log(i);
  266.         }

  268.         // divide by k!
  269.         for (int i = 2; i <= k; i++) {
  270.             logSum -= FastMath.log(i);
  271.         }

  273.         return logSum;
  274.     }

  276.     /**
  277.      * Returns n!. Shorthand for {@code n} <a
  278.      * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
  279.      * product of the numbers {@code 1,...,n}.
  280.      * <p>* <Strong>Preconditions</strong>:</p>
  281.      * <ul>
  282.      * <li> {@code n >= 0} (otherwise
  283.      * {@code MathIllegalArgumentException} is thrown)</li>
  284.      * <li> The result is small enough to fit into a {@code long}. The
  285.      * largest value of {@code n} for which {@code n!} does not exceed
  286.      * Long.MAX_VALUE} is 20. If the computed value exceeds {@code Long.MAX_VALUE}
  287.      * an {@code MathRuntimeException } is thrown.</li>
  288.      * </ul>
  289.      *
  290.      * @param n argument
  291.      * @return {@code n!}
  292.      * @throws MathRuntimeException if the result is too large to be represented
  293.      * by a {@code long}.
  294.      * @throws MathIllegalArgumentException if {@code n < 0}.
  295.      * @throws MathIllegalArgumentException if {@code n > 20}: The factorial value is too
  296.      * large to fit in a {@code long}.
  297.      */
  298.     public static long factorial(final int n) throws MathIllegalArgumentException {
  299.         if (n < 0) {
  300.             throw new MathIllegalArgumentException(LocalizedCoreFormats.FACTORIAL_NEGATIVE_PARAMETER, n);
  301.         }
  302.         if (n > 20) {
  303.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE, n, 20);
  304.         }
  305.         return FACTORIALS[n];
  306.     }

  308.     /**
  309.      * Compute n!, the<a href="http://mathworld.wolfram.com/Factorial.html">
  310.      * factorial</a> of {@code n} (the product of the numbers 1 to n), as a
  311.      * {@code double}.
  312.      * The result should be small enough to fit into a {@code double}: The
  313.      * largest {@code n} for which {@code n!} does not exceed
  314.      * {@code Double.MAX_VALUE} is 170. If the computed value exceeds
  315.      * {@code Double.MAX_VALUE}, {@code Double.POSITIVE_INFINITY} is returned.
  316.      *
  317.      * @param n Argument.
  318.      * @return {@code n!}
  319.      * @throws MathIllegalArgumentException if {@code n < 0}.
  320.      */
  321.     public static double factorialDouble(final int n) throws MathIllegalArgumentException {
  322.         if (n < 0) {
  323.             throw new MathIllegalArgumentException(LocalizedCoreFormats.FACTORIAL_NEGATIVE_PARAMETER,
  324.                                            n);
  325.         }
  326.         if (n < 21) {
  327.             return FACTORIALS[n];
  328.         }
  329.         return FastMath.floor(FastMath.exp(factorialLog(n)) + 0.5);
  330.     }

  332.     /**
  333.      * Compute the natural logarithm of the factorial of {@code n}.
  334.      *
  335.      * @param n Argument.
  336.      * @return {@code log(n!)}
  337.      * @throws MathIllegalArgumentException if {@code n < 0}.
  338.      */
  339.     public static double factorialLog(final int n) throws MathIllegalArgumentException {
  340.         return FACTORIAL_LOG_NO_CACHE.value(n);
  341.     }

  343.     /**
  344.      * Returns the <a
  345.      * href="http://mathworld.wolfram.com/StirlingNumberoftheSecondKind.html">
  346.      * Stirling number of the second kind</a>, "{@code S(n,k)}", the number of
  347.      * ways of partitioning an {@code n}-element set into {@code k} non-empty
  348.      * subsets.
  349.      * <p>
  350.      * The preconditions are {@code 0 <= k <= n } (otherwise
  351.      * {@code MathIllegalArgumentException} is thrown)
  352.      * </p>
  353.      * @param n the size of the set
  354.      * @param k the number of non-empty subsets
  355.      * @return {@code S(n,k)}
  356.      * @throws MathIllegalArgumentException if {@code k < 0}.
  357.      * @throws MathIllegalArgumentException if {@code k > n}.
  358.      * @throws MathRuntimeException if some overflow happens, typically for n exceeding 25 and
  359.      * k between 20 and n-2 (S(n,n-1) is handled specifically and does not overflow)
  360.      */
  361.     public static long stirlingS2(final int n, final int k)
  362.         throws MathIllegalArgumentException, MathRuntimeException {
  363.         if (k < 0) {
  364.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL, k, 0);
  365.         }
  366.         if (k > n) {
  367.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
  368.                                                    k, n);
  369.         }

  371.         long[][] stirlingS2 = STIRLING_S2.get();

  373.         if (stirlingS2 == null) {
  374.             // the cache has never been initialized, compute the first numbers
  375.             // by direct recurrence relation

  377.             // as S(26,9) = 11201516780955125625 is larger than Long.MAX_VALUE
  378.             // we must stop computation at row 26
  379.             final int maxIndex = 26;
  380.             stirlingS2 = new long[maxIndex][];
  381.             stirlingS2[0] = new long[] { 1l };
  382.             for (int i = 1; i < stirlingS2.length; ++i) {
  383.                 stirlingS2[i] = new long[i + 1];
  384.                 stirlingS2[i][0] = 0;
  385.                 stirlingS2[i][1] = 1;
  386.                 stirlingS2[i][i] = 1;
  387.                 for (int j = 2; j < i; ++j) {
  388.                     stirlingS2[i][j] = j * stirlingS2[i - 1][j] + stirlingS2[i - 1][j - 1];
  389.                 }
  390.             }

  392.             // atomically save the cache
  393.             STIRLING_S2.compareAndSet(null, stirlingS2);

  395.         }

  397.         if (n < stirlingS2.length) {
  398.             // the number is in the small cache
  399.             return stirlingS2[n][k];
  400.         } else {
  401.             // use explicit formula to compute the number without caching it
  402.             if (k == 0) {
  403.                 return 0;
  404.             } else if (k == 1 || k == n) {
  405.                 return 1;
  406.             } else if (k == 2) {
  407.                 return (1l << (n - 1)) - 1l;
  408.             } else if (k == n - 1) {
  409.                 return binomialCoefficient(n, 2);
  410.             } else {
  411.                 // definition formula: note that this may trigger some overflow
  412.                 long sum = 0;
  413.                 long sign = ((k & 0x1) == 0) ? 1 : -1;
  414.                 for (int j = 1; j <= k; ++j) {
  415.                     sign = -sign;
  416.                     sum += sign * binomialCoefficient(k, j) * ArithmeticUtils.pow(j, n);
  417.                     if (sum < 0) {
  418.                         // there was an overflow somewhere
  419.                         throw new MathRuntimeException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  420.                                                           n, 0, stirlingS2.length - 1);
  421.                     }
  422.                 }
  423.                 return sum / factorial(k);
  424.             }
  425.         }

  427.     }

  429.     /**
  430.      * Returns an iterator whose range is the k-element subsets of {0, ..., n - 1}
  431.      * represented as {@code int[]} arrays.
  432.      * <p>
  433.      * The arrays returned by the iterator are sorted in descending order and
  434.      * they are visited in lexicographic order with significance from right to
  435.      * left. For example, combinationsIterator(4, 2) returns an Iterator that
  436.      * will generate the following sequence of arrays on successive calls to
  437.      * {@code next()}:</p><p>
  438.      * {@code [0, 1], [0, 2], [1, 2], [0, 3], [1, 3], [2, 3]}
  439.      * </p><p>
  440.      * If {@code k == 0} an Iterator containing an empty array is returned and
  441.      * if {@code k == n} an Iterator containing [0, ..., n -1] is returned.</p>
  442.      *
  443.      * @param n Size of the set from which subsets are selected.
  444.      * @param k Size of the subsets to be enumerated.
  445.      * @return an {@link Iterator iterator} over the k-sets in n.
  446.      * @throws MathIllegalArgumentException if {@code n < 0}.
  447.      * @throws MathIllegalArgumentException if {@code k > n}.
  448.      */
  449.     public static Iterator<int[]> combinationsIterator(int n, int k) {
  450.         return new Combinations(n, k).iterator();
  451.     }

  453.     /**
  454.      * Check binomial preconditions.
  455.      *
  456.      * @param n Size of the set.
  457.      * @param k Size of the subsets to be counted.
  458.      * @throws MathIllegalArgumentException if {@code n < 0}.
  459.      * @throws MathIllegalArgumentException if {@code k > n}.
  460.      */
  461.     public static void checkBinomial(final int n,
  462.                                      final int k)
  463.         throws MathIllegalArgumentException {
  464.         if (n < k) {
  465.             throw new MathIllegalArgumentException(LocalizedCoreFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
  466.                                                 k, n, true);
  467.         }
  468.         if (n < 0) {
  469.             throw new MathIllegalArgumentException(LocalizedCoreFormats.BINOMIAL_NEGATIVE_PARAMETER, n);
  470.         }
  471.     }

  473.     /** Compute the Bell number (number of partitions of a set).
  474.      * @param n number of elements of the set
  475.      * @return Bell number Bₙ
  476.      * @since 2.2
  477.      */
  478.     public static long bellNumber(final int n) {

  480.         if (n < 0) {
  481.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL, n, 0);
  482.         }
  483.         if (n > MAX_BELL) {
  484.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE, n, MAX_BELL);
  485.         }

  487.         long[] bell = BELL.get();

  489.         if (bell == null) {

  491.             // the cache has never been initialized, compute the numbers using the Bell triangle
  492.             // storage for one line of the Bell triangle
  493.             bell = new long[MAX_BELL];
  494.             bell[0] = 1l;

  496.             final long[] row = new long[bell.length];
  497.             for (int k = 1; k < row.length; ++k) {

  499.                 // first row, with one element
  500.                 row[0] = 1l;

  502.                 // iterative computation of rows
  503.                 for (int i = 1; i < k; ++i) {
  504.                     long previous = row[0];
  505.                     row[0] = row[i - 1];
  506.                     for (int j = 1; j <= i; ++j) {
  507.                         long rj = row[j - 1] + previous;
  508.                         previous = row[j];
  509.                         row[j] = rj;
  510.                     }
  511.                 }

  513.                 bell[k] = row[k - 1];

  515.             }

  517.             BELL.compareAndSet(null, bell);

  519.         }

  521.         return bell[n];

  523.     }

  525.     /** Generate a stream of partitions of a list.
  526.      * <p>
  527.      * This method implements the iterative algorithm described in
  528.      * <a href="https://academic.oup.com/comjnl/article/32/3/281/331557">Short Note:
  529.      * A Fast Iterative Algorithm for Generating Set Partitions</a>
  530.      * by B. Djokić, M. Miyakawa, S. Sekiguchi, I. Semba, and I. Stojmenović
  531.      * (The Computer Journal, Volume 32, Issue 3, 1989, Pages 281–282,
  532.      * <a href="https://doi.org/10.1093/comjnl/32.3.281">https://doi.org/10.1093/comjnl/32.3.281</a>
  533.      * </p>
  534.      * @param <T> type of the list elements
  535.      * @param list list to partition
  536.      * @return stream of partitions of the list, each partition is an array or parts
  537.      * and each part is a list of elements
  538.      * @since 2.2
  539.      */
  540.     public static <T> Stream<List<T>[]> partitions(final List<T> list) {

  542.         // handle special cases of empty and singleton lists
  543.         if (list.size() < 2) {
  544.             @SuppressWarnings("unchecked")
  545.             final List<T>[] partition = (List<T>[]) Array.newInstance(List.class, 1);
  546.             partition[0] = list;
  547.             final Stream.Builder<List<T>[]> builder = Stream.builder();
  548.             return builder.add(partition).build();
  549.         }

  551.         return StreamSupport.stream(Spliterators.spliteratorUnknownSize(new PartitionsIterator<>(list),
  552.                                                                         Spliterator.DISTINCT | Spliterator.NONNULL |
  553.                                                                         Spliterator.IMMUTABLE | Spliterator.ORDERED),
  554.                                     false);

  556.     }

  558.     /** Generate a stream of permutations of a list.
  559.      * <p>
  560.      * This method implements the Steinhaus–Johnson–Trotter algorithm
  561.      * with Even's speedup
  562.      * <a href="https://en.wikipedia.org/wiki/Steinhaus%E2%80%93Johnson%E2%80%93Trotter_algorithm">Steinhaus–Johnson–Trotter algorithm</a>
  563.      * @param <T> type of the list elements
  564.      * @param list list to permute
  565.      * @return stream of permutations of the list
  566.      * @since 2.2
  567.      */
  568.     public static <T> Stream<List<T>> permutations(final List<T> list) {

  570.         // handle special cases of empty and singleton lists
  571.         if (list.size() < 2) {
  572.             return Stream.of(list);
  573.         }

  575.         return StreamSupport.stream(Spliterators.spliteratorUnknownSize(new PermutationsIterator<>(list),
  576.                                                                         Spliterator.DISTINCT | Spliterator.NONNULL |
  577.                                                                         Spliterator.IMMUTABLE | Spliterator.ORDERED),
  578.                                     false);

  580.     }

  582.     /**
  583.      * Class for computing the natural logarithm of the factorial of {@code n}.
  584.      * It allows to allocate a cache of precomputed values.
  585.      * In case of cache miss, computation is preformed by a call to
  586.      * {@link Gamma#logGamma(double)}.
  587.      */
  588.     public static final class FactorialLog {
  589.         /**
  590.          * Precomputed values of the function:
  591.          * {@code LOG_FACTORIALS[i] = log(i!)}.
  592.          */
  593.         private final double[] LOG_FACTORIALS;

  595.         /**
  596.          * Creates an instance, reusing the already computed values if available.
  597.          *
  598.          * @param numValues Number of values of the function to compute.
  599.          * @param cache Existing cache.
  600.          * @throw MathIllegalArgumentException if {@code n < 0}.
  601.          */
  602.         private FactorialLog(int numValues,
  603.                              double[] cache) {
  604.             if (numValues < 0) {
  605.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL, numValues, 0);
  606.             }

  608.             LOG_FACTORIALS = new double[numValues];

  610.             final int beginCopy = 2;
  611.             final int endCopy = cache == null || cache.length <= beginCopy ?
  612.                 beginCopy : cache.length <= numValues ?
  613.                 cache.length : numValues;

  615.             // Copy available values.
  616.             for (int i = beginCopy; i < endCopy; i++) {
  617.                 LOG_FACTORIALS[i] = cache[i];
  618.             }

  620.             // Precompute.
  621.             for (int i = endCopy; i < numValues; i++) {
  622.                 LOG_FACTORIALS[i] = LOG_FACTORIALS[i - 1] + FastMath.log(i);
  623.             }
  624.         }

  626.         /**
  627.          * Creates an instance with no precomputed values.
  628.          * @return instance with no precomputed values
  629.          */
  630.         public static FactorialLog create() {
  631.             return new FactorialLog(0, null);
  632.         }

  634.         /**
  635.          * Creates an instance with the specified cache size.
  636.          *
  637.          * @param cacheSize Number of precomputed values of the function.
  638.          * @return a new instance where {@code cacheSize} values have been
  639.          * precomputed.
  640.          * @throws MathIllegalArgumentException if {@code n < 0}.
  641.          */
  642.         public FactorialLog withCache(final int cacheSize) {
  643.             return new FactorialLog(cacheSize, LOG_FACTORIALS);
  644.         }

  646.         /**
  647.          * Computes {@code log(n!)}.
  648.          *
  649.          * @param n Argument.
  650.          * @return {@code log(n!)}.
  651.          * @throws MathIllegalArgumentException if {@code n < 0}.
  652.          */
  653.         public double value(final int n) {
  654.             if (n < 0) {
  655.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.FACTORIAL_NEGATIVE_PARAMETER,
  656.                                                n);
  657.             }

  659.             // Use cache of precomputed values.
  660.             if (n < LOG_FACTORIALS.length) {
  661.                 return LOG_FACTORIALS[n];
  662.             }

  664.             // Use cache of precomputed factorial values.
  665.             if (n < FACTORIALS.length) {
  666.                 return FastMath.log(FACTORIALS[n]);
  667.             }

  669.             // Delegate.
  670.             return Gamma.logGamma(n + 1);
  671.         }
  672.     }
  673. }
Created with JaCoCo 0.8.12.202403310830