diff --git a/src/main/java/org/apache/commons/math4/analysis/differentiation/DSCompiler.java b/src/main/java/org/apache/commons/math4/analysis/differentiation/DSCompiler.java
index 3c1838dc2..dd857de21 100644
--- a/src/main/java/org/apache/commons/math4/analysis/differentiation/DSCompiler.java
+++ b/src/main/java/org/apache/commons/math4/analysis/differentiation/DSCompiler.java
@@ -26,7 +26,7 @@ import org.apache.commons.math4.exception.MathArithmeticException;
 import org.apache.commons.math4.exception.MathInternalError;
 import org.apache.commons.math4.exception.NotPositiveException;
 import org.apache.commons.math4.exception.NumberIsTooLargeException;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.FactorialDouble;
 import org.apache.commons.math4.util.FastMath;
 import org.apache.commons.math4.util.MathArrays;
 
@@ -123,6 +123,8 @@ import org.apache.commons.math4.util.MathArrays;
  * @since 3.1
  */
 public class DSCompiler {
+    /** Cache for factorials. */
+    private static FactorialDouble FACTORIAL = FactorialDouble.create().withCache(30);
 
     /** Array of all compilers created so far. */
     private static AtomicReference<DSCompiler[][]> compilers =
@@ -1803,8 +1805,7 @@ public class DSCompiler {
             for (int k = 0; k < orders.length; ++k) {
                 if (orders[k] > 0) {
                     try {
-                        term *= FastMath.pow(delta[k], orders[k]) /
-                        CombinatoricsUtils.factorial(orders[k]);
+                        term *= FastMath.pow(delta[k], orders[k]) / FACTORIAL.value(orders[k]);
                     } catch (NotPositiveException e) {
                         // this cannot happen
                         throw new MathInternalError(e);
diff --git a/src/main/java/org/apache/commons/math4/analysis/interpolation/HermiteInterpolator.java b/src/main/java/org/apache/commons/math4/analysis/interpolation/HermiteInterpolator.java
index 5a3081506..439eb1cc6 100644
--- a/src/main/java/org/apache/commons/math4/analysis/interpolation/HermiteInterpolator.java
+++ b/src/main/java/org/apache/commons/math4/analysis/interpolation/HermiteInterpolator.java
@@ -27,7 +27,7 @@ import org.apache.commons.math4.exception.MathArithmeticException;
 import org.apache.commons.math4.exception.NoDataException;
 import org.apache.commons.math4.exception.ZeroException;
 import org.apache.commons.math4.exception.util.LocalizedFormats;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.Factorial;
 
 /** Polynomial interpolator using both sample values and sample derivatives.
  * <p>
@@ -86,11 +86,13 @@ public class HermiteInterpolator implements UnivariateDifferentiableVectorFuncti
     public void addSamplePoint(final double x, final double[] ... value)
         throws ZeroException, MathArithmeticException {
 
+        if (value.length > 20) {
+            throw new MathArithmeticException(LocalizedFormats.NUMBER_TOO_LARGE, value.length, 20);
+        }
         for (int i = 0; i < value.length; ++i) {
-
             final double[] y = value[i].clone();
             if (i > 1) {
-                double inv = 1.0 / CombinatoricsUtils.factorial(i);
+                double inv = 1.0 / Factorial.value(i);
                 for (int j = 0; j < y.length; ++j) {
                     y[j] *= inv;
                 }
diff --git a/src/main/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtils.java b/src/main/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtils.java
index 92e9e6d9a..018df8126 100644
--- a/src/main/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtils.java
+++ b/src/main/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtils.java
@@ -22,7 +22,7 @@ import java.util.List;
 import java.util.Map;
 
 import org.apache.commons.math4.fraction.BigFraction;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
 import org.apache.commons.math4.util.FastMath;
 
 /**
@@ -331,7 +331,7 @@ public class PolynomialsUtils {
         final int[][] coeff = new int[dp1][dp1];
         for (int i = 0; i < dp1; i++){
             for(int j = 0; j <= i; j++){
-                coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j);
+                coeff[i][j] = (int) BinomialCoefficient.value(i, j);
             }
         }
 
diff --git a/src/main/java/org/apache/commons/math4/distribution/PascalDistribution.java b/src/main/java/org/apache/commons/math4/distribution/PascalDistribution.java
index f1768657c..9a693d193 100644
--- a/src/main/java/org/apache/commons/math4/distribution/PascalDistribution.java
+++ b/src/main/java/org/apache/commons/math4/distribution/PascalDistribution.java
@@ -19,8 +19,9 @@ package org.apache.commons.math4.distribution;
 import org.apache.commons.math4.exception.NotStrictlyPositiveException;
 import org.apache.commons.math4.exception.OutOfRangeException;
 import org.apache.commons.math4.exception.util.LocalizedFormats;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficientDouble;
+import org.apache.commons.numbers.combinatorics.LogBinomialCoefficient;
 import org.apache.commons.numbers.gamma.RegularizedBeta;
-import org.apache.commons.math4.util.CombinatoricsUtils;
 import org.apache.commons.math4.util.FastMath;
 
 /**
@@ -125,7 +126,7 @@ public class PascalDistribution extends AbstractIntegerDistribution {
         if (x < 0) {
             ret = 0.0;
         } else {
-            ret = CombinatoricsUtils.binomialCoefficientDouble(x +
+            ret = BinomialCoefficientDouble.value(x +
                   numberOfSuccesses - 1, numberOfSuccesses - 1) *
                   FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
                   FastMath.pow(1.0 - probabilityOfSuccess, x);
@@ -140,7 +141,7 @@ public class PascalDistribution extends AbstractIntegerDistribution {
         if (x < 0) {
             ret = Double.NEGATIVE_INFINITY;
         } else {
-            ret = CombinatoricsUtils.binomialCoefficientLog(x +
+            ret = LogBinomialCoefficient.value(x +
                   numberOfSuccesses - 1, numberOfSuccesses - 1) +
                   logProbabilityOfSuccess * numberOfSuccesses +
                   log1mProbabilityOfSuccess * x;
diff --git a/src/main/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTest.java b/src/main/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTest.java
index ef459176f..40d9d91ff 100644
--- a/src/main/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTest.java
+++ b/src/main/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTest.java
@@ -22,6 +22,7 @@ import java.util.Arrays;
 
 import org.apache.commons.rng.simple.RandomSource;
 import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficientDouble;
 import org.apache.commons.math4.distribution.EnumeratedRealDistribution;
 import org.apache.commons.math4.distribution.RealDistribution;
 import org.apache.commons.math4.distribution.AbstractRealDistribution;
@@ -41,7 +42,6 @@ import org.apache.commons.math4.linear.Array2DRowFieldMatrix;
 import org.apache.commons.math4.linear.FieldMatrix;
 import org.apache.commons.math4.linear.MatrixUtils;
 import org.apache.commons.math4.linear.RealMatrix;
-import org.apache.commons.math4.util.CombinatoricsUtils;
 import org.apache.commons.math4.util.FastMath;
 import org.apache.commons.math4.util.MathArrays;
 import org.apache.commons.math4.util.MathUtils;
@@ -959,7 +959,7 @@ public class KolmogorovSmirnovTest {
      */
     public double exactP(double d, int n, int m, boolean strict) {
        return 1 - n(m, n, m, n, calculateIntegralD(d, m, n, strict), strict) /
-               CombinatoricsUtils.binomialCoefficientDouble(n + m, m);
+           BinomialCoefficientDouble.value(n + m, m);
     }
 
     /**
diff --git a/src/test/java/org/apache/commons/math4/analysis/differentiation/DSCompilerTest.java b/src/test/java/org/apache/commons/math4/analysis/differentiation/DSCompilerTest.java
index 428450541..c466a5c38 100644
--- a/src/test/java/org/apache/commons/math4/analysis/differentiation/DSCompilerTest.java
+++ b/src/test/java/org/apache/commons/math4/analysis/differentiation/DSCompilerTest.java
@@ -23,7 +23,7 @@ import java.util.Map;
 
 import org.apache.commons.math4.analysis.differentiation.DSCompiler;
 import org.apache.commons.math4.exception.DimensionMismatchException;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
 import org.junit.Assert;
 import org.junit.Test;
 
@@ -37,7 +37,7 @@ public class DSCompilerTest {
     public void testSize() {
         for (int i = 0; i < 6; ++i) {
             for (int j = 0; j < 6; ++j) {
-                long expected = CombinatoricsUtils.binomialCoefficient(i + j, i);
+                long expected = BinomialCoefficient.value(i + j, i);
                 Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
                 Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
             }
diff --git a/src/test/java/org/apache/commons/math4/analysis/differentiation/DerivativeStructureTest.java b/src/test/java/org/apache/commons/math4/analysis/differentiation/DerivativeStructureTest.java
index 7d2a1be11..824b3e433 100644
--- a/src/test/java/org/apache/commons/math4/analysis/differentiation/DerivativeStructureTest.java
+++ b/src/test/java/org/apache/commons/math4/analysis/differentiation/DerivativeStructureTest.java
@@ -29,7 +29,7 @@ import org.apache.commons.math4.exception.NumberIsTooLargeException;
 import org.apache.commons.rng.UniformRandomProvider;
 import org.apache.commons.rng.simple.RandomSource;
 import org.apache.commons.numbers.core.ArithmeticUtils;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.Factorial;
 import org.apache.commons.math4.util.FastMath;
 import org.apache.commons.numbers.core.Precision;
 import org.junit.Assert;
@@ -188,7 +188,7 @@ public class DerivativeStructureTest extends ExtendedFieldElementAbstractTest<De
             DerivativeStructure r = new DerivativeStructure(1, 6, 0, x).reciprocal();
             Assert.assertEquals(1 / x, r.getValue(), 1.0e-15);
             for (int i = 1; i < r.getOrder(); ++i) {
-                double expected = ArithmeticUtils.pow(-1, i) * CombinatoricsUtils.factorial(i) /
+                double expected = ArithmeticUtils.pow(-1, i) * Factorial.value(i) /
                                   FastMath.pow(x, i + 1);
                 Assert.assertEquals(expected, r.getPartialDerivative(i), 1.0e-15 * FastMath.abs(expected));
             }
@@ -797,7 +797,7 @@ public class DerivativeStructureTest extends ExtendedFieldElementAbstractTest<De
                 DerivativeStructure log = new DerivativeStructure(1, maxOrder, 0, x).log();
                 Assert.assertEquals(FastMath.log(x), log.getValue(), epsilon[0]);
                 for (int n = 1; n <= maxOrder; ++n) {
-                    double refDer = -CombinatoricsUtils.factorial(n - 1) / FastMath.pow(-x, n);
+                    double refDer = -Factorial.value(n - 1) / FastMath.pow(-x, n);
                     Assert.assertEquals(refDer, log.getPartialDerivative(n), epsilon[n]);
                 }
             }
diff --git a/src/test/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtilsTest.java b/src/test/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtilsTest.java
index e9d2d09fe..9387e3492 100644
--- a/src/test/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtilsTest.java
+++ b/src/test/java/org/apache/commons/math4/analysis/polynomials/PolynomialsUtilsTest.java
@@ -18,7 +18,7 @@ package org.apache.commons.math4.analysis.polynomials;
 
 import org.apache.commons.math4.analysis.UnivariateFunction;
 import org.apache.commons.math4.analysis.integration.IterativeLegendreGaussIntegrator;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
 import org.apache.commons.math4.util.FastMath;
 import org.apache.commons.numbers.core.Precision;
 import org.junit.Assert;
@@ -292,7 +292,7 @@ public class PolynomialsUtilsTest {
             for (int w = 0; w < 10; ++w) {
                 for (int i = 0; i < 10; ++i) {
                     PolynomialFunction jacobi = PolynomialsUtils.createJacobiPolynomial(i, v, w);
-                    double binomial = CombinatoricsUtils.binomialCoefficient(v + i, i);
+                    double binomial = BinomialCoefficient.value(v + i, i);
                     Assert.assertTrue(Precision.equals(binomial, jacobi.value(1.0), 1));
                 }
             }
diff --git a/src/test/java/org/apache/commons/math4/linear/InverseHilbertMatrix.java b/src/test/java/org/apache/commons/math4/linear/InverseHilbertMatrix.java
index 141f34b17..3a44e9d08 100644
--- a/src/test/java/org/apache/commons/math4/linear/InverseHilbertMatrix.java
+++ b/src/test/java/org/apache/commons/math4/linear/InverseHilbertMatrix.java
@@ -21,7 +21,7 @@ import org.apache.commons.math4.linear.ArrayRealVector;
 import org.apache.commons.math4.linear.RealLinearOperator;
 import org.apache.commons.math4.linear.RealVector;
 import org.apache.commons.numbers.core.ArithmeticUtils;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
 
 /**
  * This class implements inverses of Hilbert Matrices as
@@ -58,11 +58,11 @@ public class InverseHilbertMatrix
      */
     public long getEntry(final int i, final int j) {
         long val = i + j + 1;
-        long aux = CombinatoricsUtils.binomialCoefficient(n + i, n - j - 1);
+        long aux = BinomialCoefficient.value(n + i, n - j - 1);
         val = ArithmeticUtils.mulAndCheck(val, aux);
-        aux = CombinatoricsUtils.binomialCoefficient(n + j, n - i - 1);
+        aux = BinomialCoefficient.value(n + j, n - i - 1);
         val = ArithmeticUtils.mulAndCheck(val, aux);
-        aux = CombinatoricsUtils.binomialCoefficient(i + j, i);
+        aux = BinomialCoefficient.value(i + j, i);
         val = ArithmeticUtils.mulAndCheck(val, aux);
         val = ArithmeticUtils.mulAndCheck(val, aux);
         return ((i + j) & 1) == 0 ? val : -val;
diff --git a/src/test/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTestTest.java b/src/test/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTestTest.java
index d6635685a..66cc38756 100644
--- a/src/test/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTestTest.java
+++ b/src/test/java/org/apache/commons/math4/stat/inference/KolmogorovSmirnovTestTest.java
@@ -25,7 +25,7 @@ import org.apache.commons.math4.distribution.NormalDistribution;
 import org.apache.commons.math4.distribution.UniformRealDistribution;
 import org.apache.commons.rng.simple.RandomSource;
 import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.math4.util.CombinatoricsUtils;
+import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
 import org.apache.commons.math4.util.FastMath;
 import org.apache.commons.math4.util.MathArrays;
 import org.apache.commons.math4.exception.NotANumberException;
@@ -644,7 +644,7 @@ public class KolmogorovSmirnovTestTest {
                 counts[x] += 1;
             }
 
-            final int numCombinations = (int) CombinatoricsUtils.binomialCoefficient(arraySize, numberOfTrueValues);
+            final int numCombinations = (int) BinomialCoefficient.value(arraySize, numberOfTrueValues);
 
             final long[] observed = new long[numCombinations];
             final double[] expected = new double[numCombinations];