diff --git a/commons-statistics-descriptive/pom.xml b/commons-statistics-descriptive/pom.xml
index ecb8f81..e53e4ca 100644
--- a/commons-statistics-descriptive/pom.xml
+++ b/commons-statistics-descriptive/pom.xml
@@ -66,12 +66,6 @@
             <scope>test</scope>
         </dependency>
 
-        <dependency>
-            <groupId>org.apache.commons</groupId>
-            <artifactId>commons-math3</artifactId>
-            <scope>test</scope>
-        </dependency>
-
         <!-- Required for DoubleEquivalence -->
         <dependency>
           <groupId>org.apache.commons</groupId>
diff --git a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/DoubleMath.java b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Interpolation.java
similarity index 91%
rename from commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/DoubleMath.java
rename to commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Interpolation.java
index 579ebcb..d2bf6dd 100644
--- a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/DoubleMath.java
+++ b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Interpolation.java
@@ -17,13 +17,13 @@
 package org.apache.commons.statistics.descriptive;
 
 /**
- * Support class for double math.
+ * Support class for interpolation.
  *
  * @since 1.1
  */
-final class DoubleMath {
+final class Interpolation {
     /** No instances. */
-    private DoubleMath() {}
+    private Interpolation() {}
 
     /**
      * Compute the arithmetic mean of the two values taking care to avoid overflow.
@@ -42,6 +42,18 @@
     }
 
     /**
+     * Compute the arithmetic mean of the two values.
+     *
+     * @param x Value.
+     * @param y Value.
+     * @return the mean
+     */
+    static double mean(int x, int y) {
+        // long arithmetic handles a 32-bit signed integer
+        return ((long) x + y) * 0.5;
+    }
+
+    /**
      * Linear interpolation between sorted values {@code a <= b} using the
      * interpolant {@code t} taking care to avoid overflow.
      *
diff --git a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Median.java b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Median.java
index 03c1869..1946ec8 100644
--- a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Median.java
+++ b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Median.java
@@ -36,6 +36,8 @@
  *
  * <p>The {@link NaNPolicy} can be used to change the behaviour on {@code NaN} values.
  *
+ * <p>Instances of this class are immutable and thread-safe.
+ *
  * @see #with(NaNPolicy)
  * @see <a href="https://en.wikipedia.org/wiki/Median">Median (Wikipedia)</a>
  * @since 1.1
@@ -141,14 +143,14 @@
         if (n <= 2) {
             switch (n) {
             case 2:
-                // Sort data handling NaN and signed zeros.
-                // Matches the default behaviour of Quantile using p=0.5.
+                // Sorting the array matches the behaviour of Quantile for n==2
+                // Handle NaN and signed zeros
                 if (Double.compare(x[1], x[0]) < 0) {
                     final double t = x[0];
                     x[0] = x[1];
                     x[1] = t;
                 }
-                return DoubleMath.mean(x[0], x[1]);
+                return Interpolation.mean(x[0], x[1]);
             case 1:
                 return x[0];
             default:
@@ -164,6 +166,47 @@
         }
         // Even: require (m-1, m)
         Selection.select(x, 0, n, new int[] {m - 1, m});
-        return DoubleMath.mean(x[m - 1], x[m]);
+        return Interpolation.mean(x[m - 1], x[m]);
+    }
+
+    /**
+     * Evaluate the median.
+     *
+     * <p>Note: This method may partially sort the input values if not configured to
+     * {@link #withCopy(boolean) copy} the input data.
+     *
+     * @param values Values.
+     * @return the median
+     */
+    public double evaluate(int[] values) {
+        final int[] x = copy ? values.clone() : values;
+        final int n = values.length;
+        // Special cases
+        if (n <= 2) {
+            switch (n) {
+            case 2:
+                // Sorting the array matches the behaviour of Quantile for n==2
+                if (x[1] < x[0]) {
+                    final int t = x[0];
+                    x[0] = x[1];
+                    x[1] = t;
+                }
+                return Interpolation.mean(x[0], x[1]);
+            case 1:
+                return x[0];
+            default:
+                return Double.NaN;
+            }
+        }
+        // Median index
+        final int m = n >>> 1;
+        // Odd
+        if ((n & 0x1) == 1) {
+            Selection.select(x, 0, n, m);
+            return x[m];
+        }
+        // Even: require (m-1, m)
+        Selection.select(x, 0, n, new int[] {m - 1, m});
+        return Interpolation.mean(x[m - 1], x[m]);
     }
 }
diff --git a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Quantile.java b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Quantile.java
index 29117b9..f974d1b 100644
--- a/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Quantile.java
+++ b/commons-statistics-descriptive/src/main/java/org/apache/commons/statistics/descriptive/Quantile.java
@@ -42,6 +42,8 @@
  *
  * <p>The {@link NaNPolicy} can be used to change the behaviour on {@code NaN} values.
  *
+ * <p>Instances of this class are immutable and thread-safe.
+ *
  * @see #with(NaNPolicy)
  * @see <a href="http://en.wikipedia.org/wiki/Quantile">Quantile (Wikipedia)</a>
  * @since 1.1
@@ -237,7 +239,7 @@
         // Partition and compute
         if (pos > i) {
             Selection.select(x, 0, n, new int[] {i, i + 1});
-            return DoubleMath.interpolate(x[i], x[i + 1], pos - i);
+            return Interpolation.interpolate(x[i], x[i + 1], pos - i);
         }
         Selection.select(x, 0, n, i);
         return x[i];
@@ -269,27 +271,95 @@
         }
 
         // Collect interpolation positions. We use the output q as storage.
-        final int[] indices = new int[p.length << 1];
-        int count = 0;
-        for (int k = 0; k < p.length; k++) {
-            final double pos = estimationType.index(p[k], n);
-            q[k] = pos;
-            final int i = (int) pos;
-            indices[count++] = i;
-            if (pos > i) {
-                // Require the next index for interpolation
-                indices[count++] = i + 1;
-            }
-        }
+        final int[] indices = computeIndices(n, p, q);
 
         // Partition
-        Selection.select(x, 0, n, truncate(indices, count));
+        Selection.select(x, 0, n, indices);
 
         // Compute
         for (int k = 0; k < p.length; k++) {
             final int i = (int) q[k];
             if (q[k] > i) {
-                q[k] = DoubleMath.interpolate(x[i], x[i + 1], q[k] - i);
+                q[k] = Interpolation.interpolate(x[i], x[i + 1], q[k] - i);
+            } else {
+                q[k] = x[i];
+            }
+        }
+        return q;
+    }
+
+    /**
+     * Evaluate the {@code p}-th quantile of the values.
+     *
+     * <p>Note: This method may partially sort the input values if not configured to
+     * {@link #withCopy(boolean) copy} the input data.
+     *
+     * <p><strong>Performance</strong>
+     *
+     * <p>It is not recommended to use this method for repeat calls for different quantiles
+     * within the same values. The {@link #evaluate(int[], double...)} method should be used
+     * which provides better performance.
+     *
+     * @param values Values.
+     * @param p Probability for the quantile to compute.
+     * @return the quantile
+     * @throws IllegalArgumentException if the probability {@code p} is not in the range {@code [0, 1]}
+     * @see #evaluate(int[], double...)
+     */
+    public double evaluate(int[] values, double p) {
+        checkProbability(p);
+        final int n = values.length;
+        // Special cases
+        if (n <= 1) {
+            return n == 0 ? Double.NaN : values[0];
+        }
+        final double pos = estimationType.index(p, n);
+        final int i = (int) pos;
+
+        // Partition and compute
+        final int[] x = copy ? values.clone() : values;
+        if (pos > i) {
+            Selection.select(x, 0, n, new int[] {i, i + 1});
+            return Interpolation.interpolate(x[i], x[i + 1], pos - i);
+        }
+        Selection.select(x, 0, n, i);
+        return x[i];
+    }
+
+    /**
+     * Evaluate the {@code p}-th quantiles of the values.
+     *
+     * <p>Note: This method may partially sort the input values if not configured to
+     * {@link #withCopy(boolean) copy} the input data.
+     *
+     * @param values Values.
+     * @param p Probabilities for the quantiles to compute.
+     * @return the quantiles
+     * @throws IllegalArgumentException if any probability {@code p} is not in the range {@code [0, 1]};
+     * or no probabilities are specified.
+     */
+    public double[] evaluate(int[] values, double... p) {
+        checkProbabilities(p);
+        final int n = values.length;
+        // Special cases
+        final double[] q = new double[p.length];
+        if (n <= 1) {
+            Arrays.fill(q, n == 0 ? Double.NaN : values[0]);
+            return q;
+        }
+
+        // Collect interpolation positions. We use the output q as storage.
+        final int[] indices = computeIndices(n, p, q);
+
+        // Partition
+        final int[] x = copy ? values.clone() : values;
+        Selection.select(x, 0, n, indices);
+
+        // Compute
+        for (int k = 0; k < p.length; k++) {
+            final int i = (int) q[k];
+            if (q[k] > i) {
+                q[k] = Interpolation.interpolate(x[i], x[i + 1], q[k] - i);
             } else {
                 q[k] = x[i];
             }
@@ -327,7 +397,7 @@
         final double v1 = values.applyAsDouble(i);
         if (pos > i) {
             final double v2 = values.applyAsDouble(i + 1);
-            return DoubleMath.interpolate(v1, v2, pos - i);
+            return Interpolation.interpolate(v1, v2, pos - i);
         }
         return v1;
     }
@@ -365,7 +435,7 @@
             final double v1 = values.applyAsDouble(i);
             if (pos > i) {
                 final double v2 = values.applyAsDouble(i + 1);
-                q[k] = DoubleMath.interpolate(v1, v2, pos - i);
+                q[k] = Interpolation.interpolate(v1, v2, pos - i);
             } else {
                 q[k] = v1;
             }
@@ -427,18 +497,33 @@
     }
 
     /**
-     * Truncate the array {@code a} to the given {@code size}. If the array length matches
-     * the {@code size} the original array is returned.
+     * Compute the indices required for quantile interpolation.
      *
-     * @param a Array.
-     * @param size Size.
-     * @return the array of the specified size
+     * <p>The zero-based interpolation index in {@code [0, n)} is
+     * saved into the working array {@code q} for each {@code p}.
+     *
+     * @param n Size of the data.
+     * @param p Probabilities for the quantiles to compute.
+     * @param q Working array for quantiles.
+     * @return the indices
      */
-    private static int[] truncate(int[] a, int size) {
-        if (size < a.length) {
-            return Arrays.copyOf(a, size);
+    private int[] computeIndices(int n, double[] p, double[] q) {
+        final int[] indices = new int[p.length << 1];
+        int count = 0;
+        for (int k = 0; k < p.length; k++) {
+            final double pos = estimationType.index(p[k], n);
+            q[k] = pos;
+            final int i = (int) pos;
+            indices[count++] = i;
+            if (pos > i) {
+                // Require the next index for interpolation
+                indices[count++] = i + 1;
+            }
         }
-        return a;
+        if (count < indices.length) {
+            return Arrays.copyOf(indices, count);
+        }
+        return indices;
     }
 
     /**
diff --git a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/DoubleMathTest.java b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/InterpolationTest.java
similarity index 74%
rename from commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/DoubleMathTest.java
rename to commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/InterpolationTest.java
index d695bd4..8cf32d4 100644
--- a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/DoubleMathTest.java
+++ b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/InterpolationTest.java
@@ -27,13 +27,13 @@
 import org.junit.jupiter.params.provider.MethodSource;
 
 /**
- * Test for {@link DoubleMath}.
+ * Test for {@link Interpolation}.
  */
-class DoubleMathTest {
+class InterpolationTest {
     @ParameterizedTest
     @MethodSource
     void testMean(double x, double y, double expected) {
-        Assertions.assertEquals(expected, DoubleMath.mean(x, y));
+        Assertions.assertEquals(expected, Interpolation.mean(x, y));
     }
 
     static Stream<Arguments> testMean() {
@@ -41,6 +41,7 @@
         builder.add(Arguments.of(2, 3, 2.5));
         builder.add(Arguments.of(-4, 3, -0.5));
         builder.add(Arguments.of(-4, 4, 0));
+        builder.add(Arguments.of(-4, 5, 0.5));
         builder.add(Arguments.of(-0.0, -0.0, -0.0));
         builder.add(Arguments.of(-Double.MAX_VALUE, Double.MAX_VALUE, 0));
         builder.add(Arguments.of(Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_VALUE));
@@ -55,8 +56,30 @@
 
     @ParameterizedTest
     @MethodSource
+    void testIntMean(int x, int y, double expected) {
+        Assertions.assertEquals(expected, Interpolation.mean(x, y));
+    }
+
+    static Stream<Arguments> testIntMean() {
+        final Stream.Builder<Arguments> builder = Stream.builder();
+        builder.add(Arguments.of(2, 3, 2.5));
+        builder.add(Arguments.of(-4, 3, -0.5));
+        builder.add(Arguments.of(-4, 4, 0));
+        builder.add(Arguments.of(-4, 5, 0.5));
+        builder.add(Arguments.of(0, 0, 0));
+        builder.add(Arguments.of(Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE));
+        builder.add(Arguments.of(-Integer.MAX_VALUE, Integer.MAX_VALUE, 0));
+        builder.add(Arguments.of(Integer.MAX_VALUE, -Integer.MAX_VALUE, 0));
+        builder.add(Arguments.of(Integer.MIN_VALUE, Integer.MIN_VALUE, Integer.MIN_VALUE));
+        builder.add(Arguments.of(Integer.MIN_VALUE, Integer.MAX_VALUE, -0.5));
+        builder.add(Arguments.of(Integer.MAX_VALUE, Integer.MIN_VALUE, -0.5));
+        return builder.build();
+    }
+
+    @ParameterizedTest
+    @MethodSource
     void testInterpolate(double a, double b, double t, double expected) {
-        Assertions.assertEquals(expected, DoubleMath.interpolate(a, b, t));
+        Assertions.assertEquals(expected, Interpolation.interpolate(a, b, t));
     }
 
     static Stream<Arguments> testInterpolate() {
@@ -120,21 +143,21 @@
         // Interpolation between a==-0.0 and b>=0.0 is 0.0.
         // This is known behaviour and ignored as usage has t in (0, 1).
         if (Double.compare(a, -0.0) == 0 && Double.compare(b, 0.0) >= 0) {
-            Assertions.assertEquals(a, DoubleMath.interpolate(a, b, 0), 0.0, "exactness");
+            Assertions.assertEquals(a, Interpolation.interpolate(a, b, 0), 0.0, "exactness");
         } else {
             // result is -0.0 if a=b=-0.0
-            Assertions.assertEquals(a, DoubleMath.interpolate(a, b, 0), "exactness");
+            Assertions.assertEquals(a, Interpolation.interpolate(a, b, 0), "exactness");
         }
         // Not supported
         //Assertions.assertEquals(b, DoubleMath.interpolate(a, b, 1));
         Assertions.assertTrue(
-            Double.compare(DoubleMath.interpolate(a, b, t2), DoubleMath.interpolate(a, b, t1)) *
+            Double.compare(Interpolation.interpolate(a, b, t2), Interpolation.interpolate(a, b, t1)) *
                 Double.compare(t2, t1) * Double.compare(b, a) >= 0, "monotonicity");
         final double m = Double.MAX_VALUE;
-        Assertions.assertTrue(Double.isFinite(DoubleMath.interpolate(a * m, b * m, t1)), "boundedness");
-        Assertions.assertTrue(Double.isFinite(DoubleMath.interpolate(a * m, b * m, t2)), "boundedness");
-        Assertions.assertEquals(a, DoubleMath.interpolate(a, a, t1), "consistency");
-        Assertions.assertEquals(b, DoubleMath.interpolate(b, b, t1), "consistency");
+        Assertions.assertTrue(Double.isFinite(Interpolation.interpolate(a * m, b * m, t1)), "boundedness");
+        Assertions.assertTrue(Double.isFinite(Interpolation.interpolate(a * m, b * m, t2)), "boundedness");
+        Assertions.assertEquals(a, Interpolation.interpolate(a, a, t1), "consistency");
+        Assertions.assertEquals(b, Interpolation.interpolate(b, b, t1), "consistency");
     }
 
     static Stream<Arguments> testInterpolateProperties() {
@@ -154,26 +177,47 @@
     }
 
     @ParameterizedTest
-    @MethodSource
+    @MethodSource(value = {"anyFiniteDouble"})
     void testInterpolateEdge(double a, double b) {
         if (b < a) {
             testInterpolateEdge(b, a);
             return;
         }
         // Test the extremes of t in (0, 1)
-        Assertions.assertTrue(DoubleMath.interpolate(a, b, 0.0) >= a, () -> "t=0 a=" + a + " b=" + b);
-        Assertions.assertTrue(DoubleMath.interpolate(a, b, Double.MIN_VALUE) >= a, () -> "t=min a=" + a + " b=" + b);
-        Assertions.assertTrue(DoubleMath.interpolate(a, b, Math.nextDown(1.0)) <= b, () -> "t=0.999... a=" + a + " b=" + b);
+        Assertions.assertTrue(Interpolation.interpolate(a, b, 0.0) >= a, () -> "t=0 a=" + a + " b=" + b);
+        Assertions.assertTrue(Interpolation.interpolate(a, b, Double.MIN_VALUE) >= a, () -> "t=min a=" + a + " b=" + b);
+        Assertions.assertTrue(Interpolation.interpolate(a, b, Math.nextDown(1.0)) <= b, () -> "t=0.999... a=" + a + " b=" + b);
         // This fails with the current implementation as it assumes interpolation never uses t=1
         //Assertions.assertTrue(DoubleMath.interpolate(a, b, 1.0) <= b, () -> "t=1 a=" + a + " b=" + b);
     }
 
-    static Stream<Arguments> testInterpolateEdge() {
+    @ParameterizedTest
+    @MethodSource(value = {"anyFiniteDouble"})
+    void testMeanVsInterpolate(double a, double b) {
+        if (b < a) {
+            testMeanVsInterpolate(b, a);
+            return;
+        }
+        // The mean should be the exact double
+        // but the interpolation has more float additions.
+        // The multiplication by the interpolant 0.5 is exact for normal numbers.
+        double expected = a + b;
+        if (!Double.isFinite(expected) || Math.abs(expected) < Double.MIN_NORMAL) {
+            return;
+        }
+        expected *= 0.5;
+        Assertions.assertEquals(expected, Interpolation.mean(a, b));
+        Assertions.assertEquals(expected, Interpolation.interpolate(a, b, 0.5), Math.ulp(expected),
+            () -> a + ", " + b);
+    }
+
+    static Stream<Arguments> anyFiniteDouble() {
         final Stream.Builder<Arguments> builder = Stream.builder();
         final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
         final long infBits = Double.doubleToRawLongBits(Double.POSITIVE_INFINITY);
         final long signBit = Long.MIN_VALUE;
         for (int i = 0; i < 50; i++) {
+            // doubles in [-1, 1)
             builder.add(Arguments.of(signedDouble(rng), signedDouble(rng)));
             // Any finite double
             final long m = rng.nextLong(infBits);
diff --git a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/MedianTest.java b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/MedianTest.java
index f832f5b..c691be3 100644
--- a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/MedianTest.java
+++ b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/MedianTest.java
@@ -19,8 +19,6 @@
 
 import java.util.Arrays;
 import java.util.stream.Stream;
-import org.apache.commons.math3.stat.descriptive.rank.Percentile;
-import org.apache.commons.math3.stat.ranking.NaNStrategy;
 import org.apache.commons.rng.UniformRandomProvider;
 import org.apache.commons.rng.simple.RandomSource;
 import org.junit.jupiter.api.Assertions;
@@ -41,18 +39,24 @@
     }
 
     @ParameterizedTest
-    @MethodSource(value = {"testMedian"})
-    void testMedian(double[] values, double expected) {
+    @MethodSource(value = {"testDoubleMedian"})
+    void testDoubleMedian(double[] values, double expected) {
         final double[] copy = values.clone();
         Assertions.assertEquals(expected, Median.withDefaults().evaluate(values));
-        // Test the result and data (modified in-place) match the quantile implementation
-        Assertions.assertEquals(expected, Quantile.withDefaults().evaluate(copy, 0.5));
+        // Test the result and data (modified in-place) match the quantile implementation.
+        // Interpolation may be different by a single ulp.
+        final double q = Quantile.withDefaults().evaluate(copy, 0.5);
+        if (Double.isNaN(expected)) {
+            Assertions.assertEquals(expected, q);
+        } else {
+            Assertions.assertEquals(expected, q, Math.ulp(expected));
+        }
         Assertions.assertArrayEquals(values, copy);
     }
 
     @ParameterizedTest
-    @MethodSource(value = {"testMedian"})
-    void testMedianExcludeNaN(double[] values, double expected) {
+    @MethodSource(value = {"testDoubleMedian"})
+    void testDoubleMedianExcludeNaN(double[] values, double expected) {
         // If NaN is present then the result will change from expected so ignore this
         Assumptions.assumeTrue(Arrays.stream(values).filter(Double::isNaN).count() == 0);
         // Note: Use copy here. This checks that the copy of the data
@@ -71,11 +75,8 @@
         }
     }
 
-    static Stream<Arguments> testMedian() {
+    static Stream<Arguments> testDoubleMedian() {
         final Stream.Builder<Arguments> builder = Stream.builder();
-        final Percentile p = new Percentile(50).withNaNStrategy(NaNStrategy.FIXED);
-        // Note: Cannot use CM when NaN is adjacent to the middle of an odd length
-        // as it always interpolates pairs and uses: low + 0.0 * (NaN - low)
         for (final double[] x : new double[][] {
             {1},
             {1, 2},
@@ -87,20 +88,25 @@
             {Double.NaN, Double.NaN, 1, 2, 3, 4},
             {Double.MAX_VALUE, Double.MAX_VALUE},
             {-Double.MAX_VALUE, -Double.MAX_VALUE / 2},
+            // Cases where quantile interpolation using p=0.5 is different.
+            // Fail cases taken from adaption of InterpolationTest.testMeanVsInterpolate.
+            {6.125850131710258E31, 2.11712251424532992E17},
+            {3.550291387137841E117, 7.941355536862782E127},
+            {9.026950581570208E-93, 1.5840864549779843E-77},
         }) {
-            builder.add(Arguments.of(x, p.evaluate(x)));
+            builder.add(Arguments.of(x, evaluate(x)));
         }
-        // Cases where CM Percentile returns NaN
-        builder.add(Arguments.of(new double[]{1, 2, Double.NaN}, 2));
-        builder.add(Arguments.of(new double[]{Double.NaN, 1, 2, 3, Double.NaN}, 3));
+        // Cases where Commons Math Percentile returns NaN (because it always interpolates
+        // x[i] + g * (x[i+1] - x[i]) even when g==0)
+        builder.add(Arguments.of(new double[] {1, 2, Double.NaN}, 2));
+        builder.add(Arguments.of(new double[] {Double.NaN, 1, 2, 3, Double.NaN}, 3));
 
-        // Test against the percentile can fail at 1 ULP so used a fixed seed
-        final UniformRandomProvider rng = RandomSource.XO_SHI_RO_128_PP.create(26378461823L);
+        final UniformRandomProvider rng = RandomSource.XO_SHI_RO_128_PP.create();
         // Sizes above and below the threshold for partitioning
         double[] x;
         for (final int size : new int[] {5, 6, 50, 51}) {
             final double[] values = rng.doubles(size, -4.5, 1.5).toArray();
-            final double expected = p.evaluate(values);
+            final double expected = evaluate(values);
             for (int i = 0; i < 20; i++) {
                 x = TestHelper.shuffle(rng, values.clone());
                 builder.add(Arguments.of(x, expected));
@@ -126,13 +132,118 @@
         return builder.build();
     }
 
+    /**
+     * Evaluate the median using a full sort on a copy of the data.
+     *
+     * @param values Value.
+     * @return the median
+     */
+    private static double evaluate(double[] values) {
+        final double[] x = values.clone();
+        Arrays.sort(x);
+        final int m = x.length >> 1;
+        if ((x.length & 0x1) == 1) {
+            // odd
+            return x[m];
+        }
+        return Interpolation.mean(x[m - 1], x[m]);
+    }
+
     @Test
-    void testMedianWithCopy() {
-        final double[] values = {3, 4, 2, 1, 0};
+    void testDoubleMedianWithCopy() {
+        assertMedianWithCopy(new double[] {2, 1}, 1.5);
+        assertMedianWithCopy(new double[] {3, 2, 1}, 2);
+        assertMedianWithCopy(new double[] {4, 3, 2, 1}, 2.5);
+        assertMedianWithCopy(new double[] {5, 4, 3, 2, 1}, 3);
+        // Special case for 2 values with signed zeros (must be unordered)
+        assertMedianWithCopy(new double[] {0.0, -0.0}, 0.0);
+    }
+
+    private static void assertMedianWithCopy(double[] values, double expected) {
         final double[] original = values.clone();
-        Assertions.assertEquals(2, Median.withDefaults().withCopy(true).evaluate(values));
+        Assertions.assertEquals(expected, Median.withDefaults().withCopy(true).evaluate(values));
         Assertions.assertArrayEquals(original, values);
-        Assertions.assertEquals(2, Median.withDefaults().withCopy(false).evaluate(values));
+        Assertions.assertEquals(expected, Median.withDefaults().withCopy(false).evaluate(values));
+        Assertions.assertFalse(Arrays.equals(original, values));
+    }
+
+    @ParameterizedTest
+    @MethodSource(value = {"testIntMedian"})
+    void testIntMedian(int[] values, double expected) {
+        final int[] copy = values.clone();
+        Assertions.assertEquals(expected, Median.withDefaults().evaluate(values));
+        // Test the result and data (modified in-place) match the quantile implementation.
+        // Results are either integer or half-integer so this should be exact.
+        Assertions.assertEquals(expected, Quantile.withDefaults().evaluate(copy, 0.5));
+        Assertions.assertArrayEquals(values, copy);
+    }
+
+    static Stream<Arguments> testIntMedian() {
+        final Stream.Builder<Arguments> builder = Stream.builder();
+        for (final int[] x : new int[][] {
+            {1},
+            {1, 2},
+            {2, 1},
+            {1, 2, 3, 4},
+            {Integer.MAX_VALUE, Integer.MAX_VALUE / 2},
+            {Integer.MIN_VALUE, Integer.MIN_VALUE / 2},
+        }) {
+            builder.add(Arguments.of(x, evaluate(x)));
+        }
+
+        final UniformRandomProvider rng = RandomSource.XO_SHI_RO_128_PP.create();
+        // Sizes above and below the threshold for partitioning
+        int[] x;
+        for (final int size : new int[] {5, 6, 50, 51}) {
+            final int[] values = rng.ints(size, -4500, 1500).toArray();
+            final double expected = evaluate(values);
+            for (int i = 0; i < 20; i++) {
+                x = TestHelper.shuffle(rng, values.clone());
+                builder.add(Arguments.of(x, expected));
+            }
+            // Special values
+            for (final int y : new int[] {0, 1, Integer.MAX_VALUE, Integer.MIN_VALUE}) {
+                x = new int[size];
+                Arrays.fill(x, y);
+                builder.add(Arguments.of(x, y));
+            }
+        }
+        // Special cases
+        builder.add(Arguments.of(new int[] {}, Double.NaN));
+        builder.add(Arguments.of(new int[] {-Integer.MAX_VALUE, Integer.MAX_VALUE}, 0));
+        return builder.build();
+    }
+
+    /**
+     * Evaluate the median using a full sort on a copy of the data.
+     *
+     * @param values Value.
+     * @return the median
+     */
+    private static double evaluate(int[] values) {
+        final int[] x = values.clone();
+        Arrays.sort(x);
+        final int m = x.length >> 1;
+        if ((x.length & 0x1) == 1) {
+            // odd
+            return x[m];
+        }
+        return Interpolation.mean(x[m - 1], x[m]);
+    }
+
+    @Test
+    void testIntMedianWithCopy() {
+        assertMedianWithCopy(new int[] {2, 1}, 1.5);
+        assertMedianWithCopy(new int[] {3, 2, 1}, 2);
+        assertMedianWithCopy(new int[] {4, 3, 2, 1}, 2.5);
+        assertMedianWithCopy(new int[] {5, 4, 3, 2, 1}, 3);
+    }
+
+    private static void assertMedianWithCopy(int[] values, double expected) {
+        final int[] original = values.clone();
+        Assertions.assertEquals(expected, Median.withDefaults().withCopy(true).evaluate(values));
+        Assertions.assertArrayEquals(original, values);
+        Assertions.assertEquals(expected, Median.withDefaults().withCopy(false).evaluate(values));
         Assertions.assertFalse(Arrays.equals(original, values));
     }
 }
diff --git a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/QuantileTest.java b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/QuantileTest.java
index 417dc08..cf4fae7 100644
--- a/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/QuantileTest.java
+++ b/commons-statistics-descriptive/src/test/java/org/apache/commons/statistics/descriptive/QuantileTest.java
@@ -35,20 +35,6 @@
     /** Estimation types to test. */
     private static final EstimationMethod[] TYPES = EstimationMethod.values();
 
-    /**
-     * Interface to test the quantile for a single probability.
-     */
-    interface QuantileFunction {
-        double evaluate(Quantile m, double[] values, double p);
-    }
-
-    /**
-     * Interface to test the quantiles for a multiple probabilities.
-     */
-    interface QuantileFunction2 {
-        double[] evaluate(Quantile m, double[] values, double[] p);
-    }
-
     @Test
     void testNullPropertyThrows() {
         final Quantile m = Quantile.withDefaults();
@@ -92,11 +78,14 @@
 
     @Test
     void testBadQuantileThrows() {
-        final double[] values = {3, 4, 2, 1, 0};
+        final double[] values1 = {3, 4, 2, 1, 0};
+        final int[] values2 = {3, 4, 2, 1, 0};
         final Quantile m = Quantile.withDefaults();
         for (final double p : new double[] {-0.5, 1.2, Double.NaN}) {
-            Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values, p));
-            Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values, new double[] {p}));
+            Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values1, p));
+            Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values1, new double[] {p}));
+            Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values2, p));
+            Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values2, new double[] {p}));
             Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(10, i -> 1, p));
             Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(10, i -> 1, new double[] {p}));
         }
@@ -104,10 +93,14 @@
 
     @Test
     void testNoQuantilesThrows() {
-        final double[] values = {3, 4, 2, 1, 0};
+        final double[] values1 = {3, 4, 2, 1, 0};
+        final int[] values2 = {3, 4, 2, 1, 0};
         final Quantile m = Quantile.withDefaults();
-        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values));
-        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values, new double[0]));
+        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values1));
+        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values1, new double[0]));
+        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values2));
+        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(values2, new double[0]));
+        Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(10, i -> 1));
         Assertions.assertThrows(IllegalArgumentException.class, () -> m.evaluate(10, i -> 1, new double[0]));
     }
 
@@ -120,16 +113,32 @@
         }
     }
 
+    // double[]
+
+    /**
+     * Interface to test the quantile for a single probability.
+     */
+    interface DoubleQuantileFunction {
+        double evaluate(Quantile m, double[] values, double p);
+    }
+
+    /**
+     * Interface to test the quantiles for a multiple probabilities.
+     */
+    interface DoubleQuantileFunctionN {
+        double[] evaluate(Quantile m, double[] values, double[] p);
+    }
+
     @ParameterizedTest
-    @MethodSource(value = {"testQuantile"})
-    void testQuantile(double[] values, double[] p, double[][] expected, double delta) {
+    @MethodSource(value = {"testDoubleQuantile"})
+    void testDoubleQuantile(double[] values, double[] p, double[][] expected, double delta) {
         assertQuantile(Quantile.withDefaults(), values, p, expected, delta,
             Quantile::evaluate, Quantile::evaluate);
     }
 
     @ParameterizedTest
-    @MethodSource(value = {"testQuantile"})
-    void testQuantileExcludeNaN(double[] values, double[] p, double[][] expected, double delta) {
+    @MethodSource(value = {"testDoubleQuantile"})
+    void testDoubleQuantileExcludeNaN(double[] values, double[] p, double[][] expected, double delta) {
         // If NaN is present then the result will change from expected so ignore this
         Assumptions.assumeTrue(Arrays.stream(values).filter(Double::isNaN).count() == 0);
         // Note: Use copy here. This checks that the copy of the data
@@ -150,7 +159,7 @@
     }
 
     @ParameterizedTest
-    @MethodSource(value = {"testQuantile"})
+    @MethodSource(value = {"testDoubleQuantile"})
     void testQuantileSorted(double[] values, double[] p, double[][] expected, double delta) {
         assertQuantile(Quantile.withDefaults(), values, p, expected, delta,
             (m, x, q) -> {
@@ -167,7 +176,7 @@
 
     private static void assertQuantile(Quantile m, double[] values, double[] p,
         double[][] expected, double delta,
-        QuantileFunction f1, QuantileFunction2 f2) {
+        DoubleQuantileFunction f1, DoubleQuantileFunctionN fn) {
         Assertions.assertEquals(expected.length, TYPES.length);
         for (int i = 0; i < TYPES.length; i++) {
             final EstimationMethod type = TYPES[i];
@@ -178,22 +187,22 @@
                     assertEqualsOrExactlyEqual(expected[i][j], f1.evaluate(m, values.clone(), p[j]), delta,
                         () -> type.toString());
                 }
-                assertEqualsOrExactlyEqual(expected[i][j], f2.evaluate(m, values.clone(), new double[] {p[j]})[0], delta,
+                assertEqualsOrExactlyEqual(expected[i][j], fn.evaluate(m, values.clone(), new double[] {p[j]})[0], delta,
                     () -> type.toString());
             }
             // Bulk quantiles
             if (delta < 0) {
-                Assertions.assertArrayEquals(expected[i], f2.evaluate(m, values.clone(), p),
+                Assertions.assertArrayEquals(expected[i], fn.evaluate(m, values.clone(), p),
                     () -> type.toString());
             } else {
-                Assertions.assertArrayEquals(expected[i], f2.evaluate(m, values.clone(), p), delta,
+                Assertions.assertArrayEquals(expected[i], fn.evaluate(m, values.clone(), p), delta,
                     () -> type.toString());
             }
         }
     }
 
     /**
-     * Assert that {@code expected} and {@code actual} are equal within the given{@code delta}.
+     * Assert that {@code expected} and {@code actual} are equal within the given {@code delta}.
      * If the {@code delta} is negative it is ignored and values must be exactly equal.
      *
      * @param expected Expected
@@ -210,79 +219,79 @@
         }
     }
 
-    static Stream<Arguments> testQuantile() {
+    static Stream<Arguments> testDoubleQuantile() {
         final Stream.Builder<Arguments> builder = Stream.builder();
         // Special cases
         final double nan = Double.NaN;
-        addQuantiles(builder, new double[] {}, new double[] {0.75}, 1e-5,
+        addDoubleQuantiles(builder, new double[] {}, new double[] {0.75}, 1e-5,
             new double[] {nan, nan, nan, nan, nan, nan, nan, nan, nan});
-        addQuantiles(builder, new double[] {42}, new double[] {0.75}, 1e-5,
+        addDoubleQuantiles(builder, new double[] {42}, new double[] {0.75}, 1e-5,
             new double[] {42, 42, 42, 42, 42, 42, 42, 42, 42});
         // Cases from Commons Math PercentileTest
-        addQuantiles(builder, new double[] {1, 2, 3}, new double[] {0.75}, 1e-5,
+        addDoubleQuantiles(builder, new double[] {1, 2, 3}, new double[] {0.75}, 1e-5,
             new double[] {3, 3, 2, 2.25, 2.75, 3, 2.5, 2.83333, 2.81250});
-        addQuantiles(builder, new double[] {0, 1}, new double[] {0.25}, 1e-5,
+        addDoubleQuantiles(builder, new double[] {0, 1}, new double[] {0.25}, 1e-5,
             new double[] {0, 0, 0, 0, 0, 0, 0.25, 0, 0});
         final double[] d = new double[] {1, 3, 2, 4};
-        addQuantiles(builder, d, new double[] {0.3, 0.25, 0.75, 0.5}, 1e-5,
+        addDoubleQuantiles(builder, d, new double[] {0.3, 0.25, 0.75, 0.5}, 1e-5,
             new double[] {2, 2, 1, 1.2, 1.7, 1.5, 1.9, 1.63333, 1.65},
             new double[] {1, 1.5, 1, 1, 1.5, 1.25, 1.75, 1.41667, 1.43750},
             new double[] {3, 3.5, 3, 3, 3.5, 3.75, 3.25, 3.58333, 3.56250},
             new double[] {2, 2.5, 2, 2, 2.5, 2.5, 2.5, 2.5, 2.5});
         // NIST example
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {95.1772, 95.1567, 95.1937, 95.1959, 95.1442, 95.0610, 95.1591, 95.1195, 95.1772, 95.0925,
                 95.1990, 95.1682},
-            new double[] {0.9}, 1e-4,
+            new double[] {0.9}, 1e-5,
             new double[] {95.19590, 95.19590, 95.19590, 95.19546, 95.19683, 95.19807, 95.19568, 95.19724, 95.19714});
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {12.5, 12.0, 11.8, 14.2, 14.9, 14.5, 21.0, 8.2, 10.3, 11.3, 14.1, 9.9, 12.2, 12.0, 12.1, 11.0,
                 19.8, 11.0, 10.0, 8.8, 9.0, 12.3},
             new double[] {0.05}, 1e-4,
             new double[] {8.8000, 8.8000, 8.2000, 8.2600, 8.5600, 8.2900, 8.8100, 8.4700, 8.4925});
         // Special values tests
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {nan},
             new double[] {0.5}, 1e-4,
             new double[] {nan, nan, nan, nan, nan, nan, nan, nan, nan});
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {nan, nan},
             new double[] {0.5}, 1e-4,
             new double[] {nan, nan, nan, nan, nan, nan, nan, nan, nan});
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {1, nan},
             new double[] {0.5}, 1e-4,
             new double[] {1, nan, 1, 1, nan, nan, nan, nan, nan});
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {1, 2, nan},
             new double[] {0.5}, 1e-4,
             new double[] {2, 2, 2, 1.5, 2, 2, 2, 2, 2});
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {1, 2, nan, nan},
             new double[] {0.5}, 1e-4,
             new double[] {2, nan, 2, 2, nan, nan, nan, nan, nan});
         // min/max test. This hits edge cases for bounds clipping when computing
         // the index in [0, n).
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {5, 4, 3, 2, 1},
-            new double[] {0.0, 1.0}, 0.0,
+            new double[] {0.0, 1.0}, -1,
             new double[] {1, 1, 1, 1, 1, 1, 1, 1, 1},
             new double[] {5, 5, 5, 5, 5, 5, 5, 5, 5});
         // Note: This tests interpolation between -0.0 and -0.0, and -0.0 and 0.0.
         // When the quantile requires interpolation, the sign should be maintained
         // if the upper bound is -0.0.
         // No interpolation
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {-0.0, 0.0, 0.0},
             new double[] {0.0}, -1,
             new double[] {-0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0});
         // Interpolation between negative zeros
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {-0.0, -0.0, -0.0, -0.0, -0.0},
             new double[] {0.45}, -1,
             new double[] {-0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0});
         // Interpolation between -0.0 and 0.0
-        addQuantiles(builder,
+        addDoubleQuantiles(builder,
             new double[] {-0.0, -0.0, 0.0, 0.0, 0.0},
             new double[] {0.45}, -1,
             // Here only HF3 rounds to k=2; other discrete methods to k=3;
@@ -354,7 +363,7 @@
      * @param p Quantiles to compute.
      * @param expected Expected result for each p for every estimation type.
      */
-    private static void addQuantiles(Stream.Builder<Arguments> builder,
+    private static void addDoubleQuantiles(Stream.Builder<Arguments> builder,
         double[] x, double[] p, double delta, double[]... expected) {
         Assertions.assertEquals(p.length, expected.length);
         for (final double[] e : expected) {
@@ -371,7 +380,7 @@
     }
 
     @Test
-    void testQuantileWithCopy() {
+    void testDoubleQuantileWithCopy() {
         final double[] values = {3, 4, 2, 1, 0};
         final double[] original = values.clone();
         Assertions.assertEquals(2, Quantile.withDefaults().withCopy(true).evaluate(values, 0.5));
@@ -381,7 +390,7 @@
     }
 
     @Test
-    void testQuantileWithCopy2() {
+    void testDoubleQuantileWithCopy2() {
         final double[] values = {3, 4, 2, 1, 0};
         final double[] original = values.clone();
         Assertions.assertEquals(2, Quantile.withDefaults().withCopy(true).evaluate(values, new double[] {0.5})[0]);
@@ -389,4 +398,211 @@
         Assertions.assertEquals(2, Quantile.withDefaults().withCopy(false).evaluate(values, new double[] {0.5})[0]);
         Assertions.assertFalse(Arrays.equals(original, values));
     }
+
+    // int[]
+
+    /**
+     * Interface to test the quantile for a single probability.
+     */
+    interface IntQuantileFunction {
+        double evaluate(Quantile m, int[] values, double p);
+    }
+
+    /**
+     * Interface to test the quantiles for a multiple probabilities.
+     */
+    interface IntQuantileFunctionN {
+        double[] evaluate(Quantile m, int[] values, double[] p);
+    }
+
+    @ParameterizedTest
+    @MethodSource(value = {"testIntQuantile"})
+    void testIntQuantile(int[] values, double[] p, double[][] expected, double delta) {
+        assertQuantile(Quantile.withDefaults(), values, p, expected, delta,
+            Quantile::evaluate, Quantile::evaluate);
+    }
+
+    @ParameterizedTest
+    @MethodSource(value = {"testIntQuantile"})
+    void testQuantileSorted(int[] values, double[] p, double[][] expected, double delta) {
+        assertQuantile(Quantile.withDefaults(), values, p, expected, delta,
+            (m, x, q) -> {
+                // No clone here as later calls with the same array will also sort it
+                Arrays.sort(x);
+                return m.evaluate(x.length, i -> x[i], q);
+            },
+            (m, x, q) -> {
+                // No clone here as later calls with the same array will also sort it
+                Arrays.sort(x);
+                return m.evaluate(x.length, i -> x[i], q);
+            });
+    }
+
+    private static void assertQuantile(Quantile m, int[] values, double[] p,
+        double[][] expected, double delta,
+        IntQuantileFunction f1, IntQuantileFunctionN fn) {
+        Assertions.assertEquals(expected.length, TYPES.length);
+        for (int i = 0; i < TYPES.length; i++) {
+            final EstimationMethod type = TYPES[i];
+            m = m.with(type);
+            // Single quantiles
+            for (int j = 0; j < p.length; j++) {
+                if (f1 != null) {
+                    assertEqualsOrExactlyEqual(expected[i][j], f1.evaluate(m, values.clone(), p[j]), delta,
+                        () -> type.toString());
+                }
+                assertEqualsOrExactlyEqual(expected[i][j], fn.evaluate(m, values.clone(), new double[] {p[j]})[0], delta,
+                    () -> type.toString());
+            }
+            // Bulk quantiles
+            if (delta < 0) {
+                Assertions.assertArrayEquals(expected[i], fn.evaluate(m, values.clone(), p),
+                    () -> type.toString());
+            } else {
+                Assertions.assertArrayEquals(expected[i], fn.evaluate(m, values.clone(), p), delta,
+                    () -> type.toString());
+            }
+        }
+    }
+
+    static Stream<Arguments> testIntQuantile() {
+        final Stream.Builder<Arguments> builder = Stream.builder();
+        // Special cases
+        final double nan = Double.NaN;
+        addIntQuantiles(builder, new int[] {}, new double[] {0.75}, 1e-5,
+            new double[] {nan, nan, nan, nan, nan, nan, nan, nan, nan});
+        addIntQuantiles(builder, new int[] {42}, new double[] {0.75}, 1e-5,
+            new double[] {42, 42, 42, 42, 42, 42, 42, 42, 42});
+        // Cases from Commons Math PercentileTest
+        addIntQuantiles(builder, new int[] {1, 2, 3}, new double[] {0.75}, 1e-5,
+            new double[] {3, 3, 2, 2.25, 2.75, 3, 2.5, 2.83333, 2.81250});
+        addIntQuantiles(builder, new int[] {0, 1}, new double[] {0.25}, 1e-5,
+            new double[] {0, 0, 0, 0, 0, 0, 0.25, 0, 0});
+        final int[] d = new int[] {1, 3, 2, 4};
+        addIntQuantiles(builder, d, new double[] {0.3, 0.25, 0.75, 0.5}, 1e-5,
+            new double[] {2, 2, 1, 1.2, 1.7, 1.5, 1.9, 1.63333, 1.65},
+            new double[] {1, 1.5, 1, 1, 1.5, 1.25, 1.75, 1.41667, 1.43750},
+            new double[] {3, 3.5, 3, 3, 3.5, 3.75, 3.25, 3.58333, 3.56250},
+            new double[] {2, 2.5, 2, 2, 2.5, 2.5, 2.5, 2.5, 2.5});
+        // NIST example
+        // Scale example 1 by 1e4
+        addIntQuantiles(builder,
+            new int[] {951772, 951567, 951937, 951959, 951442, 950610, 951591, 951195, 951772, 950925,
+                951990, 951682},
+            new double[] {0.9}, 1e-1,
+            new double[] {951959.0, 951959.0, 951959.0, 951954.6, 951968.3, 951980.7, 951956.8, 951972.4, 951971.4});
+        // Scale example 2 by 10
+        addIntQuantiles(builder,
+            new int[] {125, 120, 118, 142, 149, 145, 210, 82, 103, 113, 141, 99, 122, 120, 121, 110,
+                198, 110, 100, 88, 90, 123},
+            new double[] {0.05}, 1e-3,
+            new double[] {88.000, 88.000, 82.000, 82.600, 85.600, 82.900, 88.100, 84.700, 84.925});
+        // min/max test. This hits edge cases for bounds clipping when computing
+        // the index in [0, n).
+        addIntQuantiles(builder,
+            new int[] {5, 4, 3, 2, 1},
+            new double[] {0.0, 1.0}, -1,
+            new double[] {1, 1, 1, 1, 1, 1, 1, 1, 1},
+            new double[] {5, 5, 5, 5, 5, 5, 5, 5, 5});
+        // Interpolation between zeros
+        addIntQuantiles(builder,
+            new int[] {0, 0, 0, 0, 0},
+            new double[] {0.45}, -1,
+            new double[] {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
+
+        // Test data samples using R version 4.3.3
+        // require(stats)
+        // x = as.integer(c(1, 2, 3))   %% data
+        // p = c(0.1, 0.25, 0.5, 0.975) %% probabilities
+        // for (t in c(1:9)) { cat('{', paste(quantile(x, p, type=t), collapse=', '), "}, \n", sep="") }
+
+        final double[] p = {0.1, 0.25, 0.5, 0.975};
+
+        /** Poisson samples: scipy.stats.poisson.rvs(45, size=100). */
+        final int[] v4 = {42, 51, 38, 38, 49, 48, 42, 47, 51, 46, 45, 35, 39, 42, 49, 55, 53, 46, 49, 56,
+            42, 46, 42, 53, 43, 55, 49, 52, 51, 45, 40, 49, 39, 40, 46, 43, 46, 48, 36, 44, 40, 49, 49, 43, 45, 44, 41, 55,
+            52, 45, 57, 41, 43, 44, 38, 52, 44, 45, 43, 42, 38, 37, 47, 42, 47, 45, 70, 45, 50, 47, 46, 50, 47, 35, 43, 52,
+            51, 41, 45, 42, 45, 53, 46, 48, 51, 43, 63, 48, 49, 41, 58, 51, 59, 43, 39, 32, 35, 46, 50, 50};
+
+        builder.add(Arguments.of(v4, p, new double[][] {
+            {38, 42, 46, 59},
+            {38.5, 42, 46, 59},
+            {38, 42, 46, 59},
+            {38, 42, 46, 58.5},
+            {38.5, 42, 46, 59},
+            {38.1, 42, 46, 60.9},
+            {38.9, 42, 46, 58.525},
+            {38.3666666666667, 42, 46, 59.6333333333333},
+            {38.4, 42, 46, 59.475},
+        }, 1e-12));
+
+        // Discrete 'normal' samples: scipy.stats.norm.rvs(loc=3.4, scale=2250000, size=100).astype(int)
+        final int[] v = {659592, -849723, -1765944, 2610353, 2192883, -265784, 1126824, 1412069, 918175, -1066899,
+            -1289922, -1359925, 549577, 2891935, 5498383, 649055, -3774137, 3026349, 4317084, 16068, 1747179, -94833,
+            -891275, 146951, 659679, -1298483, -9717, -2372749, -213892, -956213, -2380241, -3265588, 515620, 334156,
+            2595489, -463102, -1490342, -700231, -245959, -1596233, 1702451, 1265231, 2338985, -1298796, -2493882,
+            -1679698, 251933, -3446511, 3437103, -2940127, -1996991, -695605, -3127437, 2895523, 1659299, 935033,
+            609115, -1245544, -1131839, 1645603, -1673754, -4318740, -163129, -1733950, 2609546, -536282, -2873472,
+            -204545, 872775, 448272, 1048969, -1781997, -3602571, -2346653, 2084923, 1289364, 2450980, -1809052,
+            -422631, 1895287, 72169, -4933595, 2602830, -1106753, 1295126, 1671634, 929809, -3094175, -787509, -769431,
+            -209387, 1517866, -1861080, 1863380, -699593, -174200, 2132930, 1957489, -340803, -2263330};
+
+        builder.add(Arguments.of(v, p, new double[][] {
+            {-2873472, -1490342, -204545, 3437103},
+            {-2683677, -1425133.5, -189372.5, 3437103},
+            {-2873472, -1490342, -204545, 3437103},
+            {-2873472, -1490342, -204545, 3231726},
+            {-2683677, -1425133.5, -189372.5, 3437103},
+            {-2835513, -1457737.75, -189372.5, 3855093.975},
+            {-2531841, -1392529.25, -189372.5, 3241994.85},
+            {-2734289, -1436001.58333333, -189372.5, 3576433.325},
+            {-2721636, -1433284.5625, -189372.5, 3541600.74374999},
+        }, 1e-8));
+
+        return builder.build();
+    }
+
+    /**
+     * Adds the quantiles.
+     *
+     * @param builder Builder.
+     * @param x Data.
+     * @param p Quantiles to compute.
+     * @param expected Expected result for each p for every estimation type.
+     */
+    private static void addIntQuantiles(Stream.Builder<Arguments> builder,
+        int[] x, double[] p, double delta, double[]... expected) {
+        Assertions.assertEquals(p.length, expected.length);
+        for (final double[] e : expected) {
+            Assertions.assertEquals(e.length, TYPES.length);
+        }
+        // Transpose
+        final double[][] t = new double[TYPES.length][p.length];
+        for (int i = 0; i < t.length; i++) {
+            for (int j = 0; j < p.length; j++) {
+                t[i][j] = expected[j][i];
+            }
+        }
+        builder.add(Arguments.of(x, p, t, delta));
+    }
+
+    @Test
+    void testIntQuantileWithCopy() {
+        final int[] values = {3, 4, 2, 1, 0};
+        final int[] original = values.clone();
+        Assertions.assertEquals(2, Quantile.withDefaults().withCopy(true).evaluate(values, 0.5));
+        Assertions.assertArrayEquals(original, values);
+        Assertions.assertEquals(2, Quantile.withDefaults().withCopy(false).evaluate(values, 0.5));
+        Assertions.assertFalse(Arrays.equals(original, values));
+    }
+
+    @Test
+    void testIntQuantileWithCopy2() {
+        final int[] values = {3, 4, 2, 1, 0};
+        final int[] original = values.clone();
+        Assertions.assertEquals(2, Quantile.withDefaults().withCopy(true).evaluate(values, new double[] {0.5})[0]);
+        Assertions.assertArrayEquals(original, values);
+        Assertions.assertEquals(2, Quantile.withDefaults().withCopy(false).evaluate(values, new double[] {0.5})[0]);
+        Assertions.assertFalse(Arrays.equals(original, values));
+    }
 }
diff --git a/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MedianPerformance.java b/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MedianPerformance.java
index 760b4ec..4622c80 100644
--- a/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MedianPerformance.java
+++ b/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/MedianPerformance.java
@@ -103,7 +103,7 @@
             // we set it anyway for completeness.
             Objects.requireNonNull(name);
             if (JDK.equals(name)) {
-                function = MedianPerformance::sortMedian;
+                function = DoubleFunctionSource::sortMedian;
             } else if (CM3.equals(name)) {
                 final org.apache.commons.math3.stat.descriptive.rank.Median m =
                     new org.apache.commons.math3.stat.descriptive.rank.Median();
@@ -118,37 +118,106 @@
                 throw new IllegalStateException("Unknown double[] function: " + name);
             }
         }
+
+        /**
+         * Sort the values and compute the median.
+         *
+         * @param values Values.
+         * @return the median
+         */
+        private static double sortMedian(double[] values) {
+            // Implicit NPE
+            final int n = values.length;
+            // Special cases
+            if (n <= 2) {
+                switch (n) {
+                case 2:
+                    return (values[0] + values[1]) * 0.5;
+                case 1:
+                    return values[0];
+                default:
+                    return Double.NaN;
+                }
+            }
+            // A sort is required
+            Arrays.sort(values);
+            final int k = n >>> 1;
+            // Odd
+            if ((n & 0x1) == 0x1) {
+                return values[k];
+            }
+            // Even
+            return (values[k - 1] + values[k]) * 0.5;
+        }
     }
 
     /**
-     * Sort the values and compute the median.
-     *
-     * @param values Values.
-     * @return the median
+     * Source of a {@link ToDoubleFunction} for a {@code int[]}.
      */
-    static double sortMedian(double[] values) {
-        // Implicit NPE
-        final int n = values.length;
-        // Special cases
-        if (n <= 2) {
-            switch (n) {
-            case 2:
-                return (values[0] + values[1]) * 0.5;
-            case 1:
-                return values[0];
-            default:
-                return Double.NaN;
+    @State(Scope.Benchmark)
+    public static class IntFunctionSource {
+        /** Name of the source. */
+        @Param({JDK, STATISTICS})
+        private String name;
+
+        /** The action. */
+        private ToDoubleFunction<int[]> function;
+
+        /**
+         * @return the function
+         */
+        public ToDoubleFunction<int[]> getFunction() {
+            return function;
+        }
+
+        /**
+         * Create the function.
+         */
+        @Setup
+        public void setup() {
+            // Note: Functions defensively copy the data by default
+            // Note: KeyStratgey does not matter for single / paired keys but
+            // we set it anyway for completeness.
+            Objects.requireNonNull(name);
+            if (JDK.equals(name)) {
+                function = IntFunctionSource::sortMedian;
+            } else if (STATISTICS.equals(name)) {
+                function = Median.withDefaults()::evaluate;
+            } else {
+                throw new IllegalStateException("Unknown int[] function: " + name);
             }
         }
-        // A sort is required
-        Arrays.sort(values);
-        final int k = n >>> 1;
-        // Odd
-        if ((n & 0x1) == 0x1) {
-            return values[k];
+
+        /**
+         * Sort the values and compute the median.
+         *
+         * @param values Values.
+         * @return the median
+         */
+        private static double sortMedian(int[] values) {
+            // Implicit NPE
+            final int n = values.length;
+            // Special cases
+            if (n <= 2) {
+                switch (n) {
+                case 2:
+                    return (values[0] + values[1]) * 0.5;
+                case 1:
+                    return values[0];
+                default:
+                    return Double.NaN;
+                }
+            }
+            // A sort is required
+            Arrays.sort(values);
+            final int k = n >>> 1;
+            // Odd
+            if ((n & 0x1) == 0x1) {
+                return values[k];
+            }
+            // Even
+            return (values[k - 1] + values[k]) * 0.5;
         }
-        // Even
-        return (values[k - 1] + values[k]) * 0.5;
     }
 
     /**
@@ -159,11 +228,27 @@
      * @param bh Data sink.
      */
     @Benchmark
-    public void median(DoubleFunctionSource function, DataSource source, Blackhole bh) {
+    public void doubleMedian(DoubleFunctionSource function, DataSource source, Blackhole bh) {
         final int size = source.size();
         final ToDoubleFunction<double[]> fun = function.getFunction();
         for (int j = -1; ++j < size;) {
             bh.consume(fun.applyAsDouble(source.getData(j)));
         }
     }
+
+    /**
+     * Create the statistic using an array.
+     *
+     * @param function Source of the function.
+     * @param source Source of the data.
+     * @param bh Data sink.
+     */
+    @Benchmark
+    public void intMedian(IntFunctionSource function, DataSource source, Blackhole bh) {
+        final int size = source.size();
+        final ToDoubleFunction<int[]> fun = function.getFunction();
+        for (int j = -1; ++j < size;) {
+            bh.consume(fun.applyAsDouble(source.getIntData(j)));
+        }
+    }
 }
diff --git a/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/QuantilePerformance.java b/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/QuantilePerformance.java
index a5a16e1..ede1c07 100644
--- a/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/QuantilePerformance.java
+++ b/commons-statistics-examples/examples-jmh/src/main/java/org/apache/commons/statistics/examples/jmh/descriptive/QuantilePerformance.java
@@ -25,6 +25,7 @@
 import java.util.Objects;
 import java.util.concurrent.TimeUnit;
 import java.util.concurrent.atomic.AtomicInteger;
+import java.util.function.BiFunction;
 import java.util.function.BinaryOperator;
 import java.util.logging.Logger;
 import org.apache.commons.rng.UniformRandomProvider;
@@ -248,6 +249,18 @@
         /**
          * Gets the sample for the given {@code index}.
          *
+         * <p>This is returned in a randomized order per iteration.
+         *
+         * @param index Index.
+         * @return the data sample
+         */
+        public int[] getIntData(int index) {
+            return getIntDataSample(order[index]);
+        }
+
+        /**
+         * Gets the sample for the given {@code index}.
+         *
          * @param index Index.
          * @return the data sample
          */
@@ -261,6 +274,22 @@
         }
 
         /**
+         * Gets the sample for the given {@code index}.
+         *
+         * @param index Index.
+         * @return the data sample
+         */
+        protected int[] getIntDataSample(int index) {
+            // For parity with other methods do not use data.clone()
+            final int[] a = data[index];
+            final int[] x = new int[a.length];
+            for (int i = -1; ++i < a.length;) {
+                x[i] = a[i];
+            }
+            return x;
+        }
+
+        /**
          * Gets the sample size for the given {@code index}.
          *
          * @param index Index.
@@ -908,7 +937,7 @@
      * Source of a {@link BinaryOperator} for a {@code double[]} and quantiles.
      */
     @State(Scope.Benchmark)
-    public static class QuantileFunctionSource {
+    public static class DoubleQuantileFunctionSource {
         /** Name of the source. */
         @Param({JDK, CM3, CM4, STATISTICS})
         private String name;
@@ -931,7 +960,7 @@
             // Note: Functions should not defensively copy the data
             // as a clone is passed in from the data source.
             if (JDK.equals(name)) {
-                function = QuantilePerformance::sortQuantile;
+                function = DoubleQuantileFunctionSource::sortQuantile;
             } else if (CM3.equals(name)) {
                 // No way to avoid a data copy here. CM does
                 // defensive copying for most array input.
@@ -968,48 +997,127 @@
                 throw new IllegalStateException("Unknown double[] function: " + name);
             }
         }
-    }
 
-    /**
-     * Sort the values and compute the median.
-     *
-     * @param values Values.
-     * @param p p-th quantiles to compute.
-     * @return the quantiles
-     */
-    static double[] sortQuantile(double[] values, double[] p) {
-        // Implicit NPE
-        final int n = values.length;
-        if (p.length == 0) {
-            throw new IllegalArgumentException("No quantiles specified");
-        }
-        for (final double pp : p) {
-            checkQuantile(pp);
-        }
-        // Special cases
-        final double[] q = new double[p.length];
-        if (n <= 1) {
-            Arrays.fill(q, n == 0 ? Double.NaN : values[0]);
+        /**
+         * Sort the values and compute the median.
+         *
+         * @param values Values.
+         * @param p p-th quantiles to compute.
+         * @return the quantiles
+         */
+        private static double[] sortQuantile(double[] values, double[] p) {
+            // Implicit NPE
+            final int n = values.length;
+            if (p.length == 0) {
+                throw new IllegalArgumentException("No quantiles specified for double[] data");
+            }
+            for (final double pp : p) {
+                checkQuantile(pp);
+            }
+            // Special cases
+            final double[] q = new double[p.length];
+            if (n <= 1) {
+                Arrays.fill(q, n == 0 ? Double.NaN : values[0]);
+                return q;
+            }
+            // A sort is required
+            Arrays.sort(values);
+            for (int i = 0; i < p.length; i++) {
+                // EstimationMethod.HF6 (as per the Apache Commons Math Percentile
+                // legacy implementation)
+                final double pos = p[i] * (n + 1);
+                final double fpos = Math.floor(pos);
+                final int j = (int) fpos;
+                final double g = pos - fpos;
+                if (j < 1) {
+                    q[i] = values[0];
+                } else if (j >= n) {
+                    q[i] = values[n - 1];
+                } else {
+                    q[i] = (1 - g) * values[j - 1] + g * values[j];
+                }
+            }
             return q;
         }
-        // A sort is required
-        Arrays.sort(values);
-        for (int i = 0; i < p.length; i++) {
-            // EstimationMethod.HF6 (as per the Apache Commons Math Percentile
-            // legacy implementation)
-            final double pos = p[i] * (n + 1);
-            final double fpos = Math.floor(pos);
-            final int j = (int) fpos;
-            final double g = pos - fpos;
-            if (j < 1) {
-                q[i] = values[0];
-            } else if (j >= n) {
-                q[i] = values[n - 1];
+    }
+
+
+    /**
+     * Source of a {@link BiFunction} for an {@code int[]} and quantiles.
+     */
+    @State(Scope.Benchmark)
+    public static class IntQuantileFunctionSource {
+        /** Name of the source. */
+        @Param({JDK, CM3, CM4, STATISTICS})
+        private String name;
+
+        /** The action. */
+        private BiFunction<int[], double[], double[]> function;
+
+        /**
+         * @return the function
+         */
+        public BiFunction<int[], double[], double[]> getFunction() {
+            return function;
+        }
+
+        /**
+         * Create the function.
+         */
+        @Setup
+        public void setup() {
+            // Note: Functions should not defensively copy the data
+            // as a clone is passed in from the data source.
+            if (JDK.equals(name)) {
+                function = IntQuantileFunctionSource::sortQuantile;
+            } else if (STATISTICS.equals(name)) {
+                function = Quantile.withDefaults()::evaluate;
             } else {
-                q[i] = (1 - g) * values[j - 1] + g * values[j];
+                throw new IllegalStateException("Unknown int[] function: " + name);
             }
         }
-        return q;
+
+        /**
+         * Sort the values and compute the median.
+         *
+         * @param values Values.
+         * @param p p-th quantiles to compute.
+         * @return the quantiles
+         */
+        private static double[] sortQuantile(int[] values, double[] p) {
+            // Implicit NPE
+            final int n = values.length;
+            if (p.length == 0) {
+                throw new IllegalArgumentException("No quantiles specified for int[] data");
+            }
+            for (final double pp : p) {
+                checkQuantile(pp);
+            }
+            // Special cases
+            final double[] q = new double[p.length];
+            if (n <= 1) {
+                Arrays.fill(q, n == 0 ? Double.NaN : values[0]);
+                return q;
+            }
+            // A sort is required
+            Arrays.sort(values);
+            for (int i = 0; i < p.length; i++) {
+                // EstimationMethod.HF6 (as per the Apache Commons Math Percentile
+                // legacy implementation)
+                final double pos = p[i] * (n + 1);
+                final double fpos = Math.floor(pos);
+                final int j = (int) fpos;
+                final double g = pos - fpos;
+                if (j < 1) {
+                    q[i] = values[0];
+                } else if (j >= n) {
+                    q[i] = values[n - 1];
+                } else {
+                    q[i] = (1 - g) * values[j - 1] + g * values[j];
+                }
+            }
+            return q;
+        }
     }
 
     /**
@@ -1033,7 +1141,7 @@
      * @param bh Data sink.
      */
     @Benchmark
-    public void quantiles(QuantileFunctionSource function, DataSource source,
+    public void doubleQuantiles(DoubleQuantileFunctionSource function, DataSource source,
             QuantileSource quantiles, Blackhole bh) {
         final int size = source.size();
         final double[] p = quantiles.getData();
@@ -1052,7 +1160,7 @@
      * @param bh Data sink.
      */
     @Benchmark
-    public void quantileRange(QuantileFunctionSource function, DataSource source,
+    public void doubleQuantileRange(DoubleQuantileFunctionSource function, DataSource source,
             QuantileRangeSource quantiles, Blackhole bh) {
         final int size = source.size();
         final double[] p = quantiles.getData();
@@ -1061,4 +1169,42 @@
             bh.consume(fun.apply(source.getData(j), p));
         }
     }
+
+    /**
+     * Create the statistic using an array and given quantiles.
+     *
+     * @param function Source of the function.
+     * @param source Source of the data.
+     * @param quantiles Source of the quantiles.
+     * @param bh Data sink.
+     */
+    @Benchmark
+    public void intQuantiles(IntQuantileFunctionSource function, DataSource source,
+            QuantileSource quantiles, Blackhole bh) {
+        final int size = source.size();
+        final double[] p = quantiles.getData();
+        final BiFunction<int[], double[], double[]> fun = function.getFunction();
+        for (int j = -1; ++j < size;) {
+            bh.consume(fun.apply(source.getIntData(j), p));
+        }
+    }
+
+    /**
+     * Create the statistic using an array and given quantiles.
+     *
+     * @param function Source of the function.
+     * @param source Source of the data.
+     * @param quantiles Source of the quantiles.
+     * @param bh Data sink.
+     */
+    @Benchmark
+    public void intQuantileRange(IntQuantileFunctionSource function, DataSource source,
+            QuantileRangeSource quantiles, Blackhole bh) {
+        final int size = source.size();
+        final double[] p = quantiles.getData();
+        final BiFunction<int[], double[], double[]> fun = function.getFunction();
+        for (int j = -1; ++j < size;) {
+            bh.consume(fun.apply(source.getIntData(j), p));
+        }
+    }
 }