001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.util;
018
019import java.util.concurrent.atomic.AtomicReference;
020
021import org.apache.commons.numbers.core.ArithmeticUtils;
022import org.apache.commons.math4.exception.MathArithmeticException;
023import org.apache.commons.math4.exception.NotPositiveException;
024import org.apache.commons.math4.exception.NumberIsTooLargeException;
025import org.apache.commons.math4.exception.util.LocalizedFormats;
026import org.apache.commons.numbers.combinatorics.Factorial;
027import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
028
029/**
030 * Combinatorial utilities.
031 *
032 * @since 3.3
033 */
034public final class CombinatoricsUtils {
035    /** Stirling numbers of the second kind. */
036    static final AtomicReference<long[][]> STIRLING_S2 = new AtomicReference<> (null);
037
038    /** Private constructor (class contains only static methods). */
039    private CombinatoricsUtils() {}
040
041    /**
042     * Returns the <a
043     * href="http://mathworld.wolfram.com/StirlingNumberoftheSecondKind.html">
044     * Stirling number of the second kind</a>, "{@code S(n,k)}", the number of
045     * ways of partitioning an {@code n}-element set into {@code k} non-empty
046     * subsets.
047     * <p>
048     * The preconditions are {@code 0 <= k <= n } (otherwise
049     * {@code NotPositiveException} is thrown)
050     * </p>
051     * @param n the size of the set
052     * @param k the number of non-empty subsets
053     * @return {@code S(n,k)}
054     * @throws NotPositiveException if {@code k < 0}.
055     * @throws NumberIsTooLargeException if {@code k > n}.
056     * @throws MathArithmeticException if some overflow happens, typically for n exceeding 25 and
057     * k between 20 and n-2 (S(n,n-1) is handled specifically and does not overflow)
058     * @since 3.1
059     */
060    public static long stirlingS2(final int n, final int k)
061        throws NotPositiveException, NumberIsTooLargeException, MathArithmeticException {
062        if (k < 0) {
063            throw new NotPositiveException(k);
064        }
065        if (k > n) {
066            throw new NumberIsTooLargeException(k, n, true);
067        }
068
069        long[][] stirlingS2 = STIRLING_S2.get();
070
071        if (stirlingS2 == null) {
072            // the cache has never been initialized, compute the first numbers
073            // by direct recurrence relation
074
075            // as S(26,9) = 11201516780955125625 is larger than Long.MAX_VALUE
076            // we must stop computation at row 26
077            final int maxIndex = 26;
078            stirlingS2 = new long[maxIndex][];
079            stirlingS2[0] = new long[] { 1l };
080            for (int i = 1; i < stirlingS2.length; ++i) {
081                stirlingS2[i] = new long[i + 1];
082                stirlingS2[i][0] = 0;
083                stirlingS2[i][1] = 1;
084                stirlingS2[i][i] = 1;
085                for (int j = 2; j < i; ++j) {
086                    stirlingS2[i][j] = j * stirlingS2[i - 1][j] + stirlingS2[i - 1][j - 1];
087                }
088            }
089
090            // atomically save the cache
091            STIRLING_S2.compareAndSet(null, stirlingS2);
092
093        }
094
095        if (n < stirlingS2.length) {
096            // the number is in the small cache
097            return stirlingS2[n][k];
098        } else {
099            // use explicit formula to compute the number without caching it
100            if (k == 0) {
101                return 0;
102            } else if (k == 1 || k == n) {
103                return 1;
104            } else if (k == 2) {
105                return (1l << (n - 1)) - 1l;
106            } else if (k == n - 1) {
107                return BinomialCoefficient.value(n, 2);
108            } else {
109                // definition formula: note that this may trigger some overflow
110                long sum = 0;
111                long sign = ((k & 0x1) == 0) ? 1 : -1;
112                for (int j = 1; j <= k; ++j) {
113                    sign = -sign;
114                    sum += sign * BinomialCoefficient.value(k, j) * ArithmeticUtils.pow(j, n);
115                    if (sum < 0) {
116                        // there was an overflow somewhere
117                        throw new MathArithmeticException(LocalizedFormats.ARGUMENT_OUTSIDE_DOMAIN,
118                                                          n, 0, stirlingS2.length - 1);
119                    }
120                }
121                return sum / Factorial.value(k);
122            }
123        }
124    }
125}