static int quickSelectAdaptive(double[] a, int left, int right, int ka, int kb,
            int[] bounds, int flags) {
        int l = left;
        int r = right;
        int m = flags;
        while (true) {
            // Select when ka and kb are close to the same end
            // |l|-----|ka|kkkkkkkk|kb|------|r|
            if (Math.min(kb - l, r - ka) < LINEAR_SORTSELECT_SIZE) {
                sortSelect(a, l, r, ka, kb);
                bounds[0] = kb;
                return ka;
            }

            // Only target ka; kb is assumed to be close
            int p0;
            final int n = r - l;
            // f in [0, 1]
            final double f = (double) (ka - l) / n;
            // Record the larger margin (start at 1/4) to create the estimated size.
            // step        L     R
            // far left    1/12  1/3   (use 1/4 + 1/32 + 1/64 ~ 0.328)
            // left        1/6   1/4
            // middle      2/9   2/9   (use 1/4 - 1/32 ~ 0.219)
            int margin = n >> 2;
            if (m < MODE_SAMPLING && r - l > FR_SAMPLING_SIZE) {
                // Floyd-Rivest sample step uses the same margins
                p0 = sampleStep(a, l, r, ka, bounds);
                if (f <= STEP_FAR_LEFT || f >= STEP_FAR_RIGHT) {
                    margin += (n >> 5) + (n >> 6);
                } else if (f > STEP_LEFT && f < STEP_RIGHT) {
                    margin -= n >> 5;
                }
            } else if (f <= STEP_LEFT) {
                if (f <= STEP_FAR_LEFT) {
                    margin += (n >> 5) + (n >> 6);
                    p0 = repeatedStepFarLeft(a, l, r, ka, bounds, m);
                } else {
                    p0 = repeatedStepLeft(a, l, r, ka, bounds, m);
                }
            } else if (f >= STEP_RIGHT) {
                if (f >= STEP_FAR_RIGHT) {
                    margin += (n >> 5) + (n >> 6);
                    p0 = repeatedStepFarRight(a, l, r, ka, bounds, m);
                } else {
                    p0 = repeatedStepRight(a, l, r, ka, bounds, m);
                }
            } else {
                margin -= n >> 5;
                p0 = repeatedStep(a, l, r, ka, bounds, m);
            }

            // Note: Here we expect [ka, kb] to be small and splitting is unlikely.
            //                   p0 p1
            // |l|--|ka|kkkk|kb|--|P|-------------------|r|
            // |l|----------------|P|--|ka|kkk|kb|------|r|
            // |l|-----------|ka|k|P|k|kb|--------------|r|
            final int p1 = bounds[0];
            if (kb < p0) {
                // Entirely on left side
                r = p0 - 1;
            } else if (ka > p1) {
                // Entirely on right side
                l = p1 + 1;
            } else {
                // Pivot splits [ka, kb]. Expect ends to be close to the pivot and finish.
                // Here we set the bounds for use after median-of-medians pivot selection.
                // In the event there are many equal values this allows collecting those
                // known to be equal together when moving around the medians sample.
                if (kb > p1) {
                    sortSelectLeft(a, p1 + 1, r, kb);
                    bounds[0] = kb;
                }
                if (ka < p0) {
                    sortSelectRight(a, l, p0 - 1, ka);
                    p0 = ka;
                }
                return p0;
            }
            // Update mode based on target partition size
            if (r - l > n - margin) {
                m++;
            }
        }
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Partitions a Floyd-Rivest sample around a pivot offset so that the input {@code k} will
     * fall in the smaller partition when the entire range is partitioned.
     *
     * <p>Assumes the range {@code r - l} is large.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int sampleStep(double[] a, int l, int r, int k, int[] upper) {

static void dualPivotQuickSelect(double[] a, int left, int right, UpdatingInterval k, int flags) {
        // If partitioning splits the interval then recursion is used for the left-most side(s)
        // and the right-most side remains within this function. If partitioning does
        // not split the interval then it remains within this function.
        int l = left;
        int r = right;
        int f = flags;
        int ka = k.left();
        int kb = k.right();
        final int[] upper = {0, 0, 0};
        while (true) {
            // Select when ka and kb are close to the same end,
            // or the entire range is small
            // |l|-----|ka|--------|kb|------|r|
            final int n = r - l;
            if (Math.min(kb - l, r - ka) < DP_SORTSELECT_SIZE ||
                n < (f & SORTSELECT_MASK)) {
                sortSelect(a, l, r, ka, kb);
                return;
            }
            if (kb - ka < DP_SORTSELECT_SIZE) {
                // Switch to single-pivot mode with Floyd-Rivest sub-sampling
                quickSelectAdaptive(a, l, r, ka, kb, upper, MODE_FR_SAMPLING);
                return;
            }
            if (f < 0) {
                // Excess recursion, switch to heap select
                heapSelect(a, l, r, ka, kb);
                return;
            }

            // Dual-pivot partitioning
            final int p0 = partition(a, l, r, upper);
            final int p1 = upper[0];

            // Recursion to max depth
            // Note: Here we possibly branch left, middle and right with multiple keys.
            // It is possible that the partition has split the keys
            // and the recursion proceeds with a reduced set in each region.
            //                   p0 p1               p2 p3
            // |l|--|ka|--k----k--|P|------k--|kb|----|P|----|r|
            //                 kb  |      ka
            f += RECURSION_INCREMENT;
            // Recurse left side if required
            if (ka < p0) {
                if (kb <= p1) {
                    // Entirely on left side
                    r = p0 - 1;
                    if (r < kb) {
                        kb = k.updateRight(r);
                    }
                    continue;
                }
                dualPivotQuickSelect(a, l, p0 - 1, k.splitLeft(p0, p1), f);
                // Here we must process middle and/or right
                ka = k.left();
            } else if (kb <= p1) {
                // No middle/right side
                return;
            } else if (ka <= p1) {
                // Advance lower bound
                ka = k.updateLeft(p1 + 1);
            }
            // Recurse middle if required
            final int p2 = upper[1];
            final int p3 = upper[2];
            if (ka < p2) {
                l = p1 + 1;
                if (kb <= p3) {
                    // Entirely in middle
                    r = p2 - 1;
                    if (r < kb) {
                        kb = k.updateRight(r);
                    }
                    continue;
                }
                dualPivotQuickSelect(a, l, p2 - 1, k.splitLeft(p2, p3), f);
                ka = k.left();
            } else if (kb <= p3) {
                // No right side
                return;
            } else if (ka <= p3) {
                ka = k.updateLeft(p3 + 1);
            }
            // Continue right
            l = p3 + 1;
        }
    }

    /**
     * Partition an array slice around 2 pivots. Partitioning exchanges array elements
     * such that all elements smaller than pivot are before it and all elements larger
     * than pivot are after it.
     *
     * <p>This method returns 4 points describing the pivot ranges of equal values.
     *
     * <pre>{@code
     *         |k0  k1|                |k2  k3|
     * |   <P  | ==P1 |  <P1 && <P2    | ==P2 |   >P   |
     * }</pre>
     *
     * <ul>
     * <li>k0: lower pivot1 point
     * <li>k1: upper pivot1 point (inclusive)
     * <li>k2: lower pivot2 point
     * <li>k3: upper pivot2 point (inclusive)
     * </ul>
     *
     * <p>Bounds are set so {@code i < k0}, {@code i > k3} and {@code k1 < i < k2} are
     * unsorted. When the range {@code [k0, k3]} contains fully sorted elements the result
     * is set to {@code k1 = k3; k2 == k0}. This can occur if
     * {@code P1 == P2} or there are zero or one value between the pivots
     * {@code P1 < v < P2}. Any sort/partition of ranges [left, k0-1], [k1+1, k2-1] and
     * [k3+1, right] must check the length is {@code > 1}.
     *
     * @param a Data array.
     * @param left Lower bound (inclusive).
     * @param right Upper bound (inclusive).
     * @param bounds Points [k1, k2, k3].
     * @return Lower bound (inclusive) of the pivot range [k0].
     */
    private static int partition(double[] a, int left, int right, int[] bounds) {

private static int repeatedStepLeft(double[] a, int l, int r, int k, int[] upper, int flags) {
        // Adapted from Alexandrescu (2016), algorithm 9.
        final int fp;
        final int s;
        int p;
        if (flags <= MODE_SAMPLING) {
            // Median into a 12th-tile
            fp = (r - l + 1) / 12;
            // Position the sample around the target k
            // Avoid bounds error due to rounding as (k-l)/(r-l) -> 1/12
            s = Math.max(k - mapDistance(k - l, l, r, fp), l + fp);
            p = k;
        } else {
            // i in 2nd quartile
            final int f = (r - l + 1) >> 2;
            final int f2 = f + f;
            final int end = l + f2;
            for (int i = l + f; i < end; i++) {
                Sorting.lowerMedian4(a, i - f, i, i + f, i + f2);
            }
            // i in 5th 12th-tile
            fp = f / 3;
            s = l + f + fp;
            // No adaption uses the middle to enforce strict margins
            p = s + (flags == MODE_ADAPTION ? mapDistance(k - l, l, r, fp) : (fp >>> 1));
        }
        final int e = s + fp - 1;
        for (int i = s; i <= e; i++) {
            Sorting.sort3(a, i - fp, i, i + fp);
        }
        p = quickSelectAdaptive(a, s, e, p, p, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, s, e, p, upper[0], upper);
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Assumes the range {@code r - l >= 11}; the caller is responsible for selection on a smaller
     * range.
     *
     * <p>Uses the Chen and Dumitrescu repeated step median-of-medians-of-medians algorithm
     * with the upper median of 4 then either median of 3 with the final sample placed in the
     * 8th 12th-tile, or max of 3 with the final sample in the 9th 12th-tile;
     * the pivot chosen from the sample is adaptive using the input {@code k}.
     *
     * <p>Given a pivot in the middle of the sample this has margins of 1/4 and 1/6.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @param flags Control flags.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int repeatedStepRight(double[] a, int l, int r, int k, int[] upper, int flags) {

private static int repeatedStepRight(double[] a, int l, int r, int k, int[] upper, int flags) {
        // Mirror image repeatedStepLeft using upper median into 3rd quartile
        final int fp;
        final int e;
        int p;
        if (flags <= MODE_SAMPLING) {
            // Median into a 12th-tile
            fp = (r - l + 1) / 12;
            // Position the sample around the target k
            // Avoid bounds error due to rounding as (r-k)/(r-l) -> 11/12
            e = Math.min(k + mapDistance(r - k, l, r, fp), r - fp);
            p = k;
        } else {
            // i in 3rd quartile
            final int f = (r - l + 1) >> 2;
            final int f2 = f + f;
            final int end = r - f2;
            for (int i = r - f; i > end; i--) {
                Sorting.upperMedian4(a, i - f2, i - f, i, i + f);
            }
            // i in 8th 12th-tile
            fp = f / 3;
            e = r - f - fp;
            // No adaption uses the middle to enforce strict margins
            p = e - (flags == MODE_ADAPTION ? mapDistance(r - k, l, r, fp) : (fp >>> 1));
        }
        final int s = e - fp + 1;
        for (int i = s; i <= e; i++) {
            Sorting.sort3(a, i - fp, i, i + fp);
        }
        p = quickSelectAdaptive(a, s, e, p, p, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, s, e, p, upper[0], upper);
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Assumes the range {@code r - l >= 11}; the caller is responsible for selection on a smaller
     * range.
     *
     * <p>Uses the Chen and Dumitrescu repeated step median-of-medians-of-medians algorithm
     * with the minimum of 4 then median of 3; the final sample is placed in the
     * 2nd 12th-tile; the pivot chosen from the sample is adaptive using the input {@code k}.
     *
     * <p>Given a pivot in the middle of the sample this has margins of 1/12 and 1/3.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @param flags Control flags.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int repeatedStepFarLeft(double[] a, int l, int r, int k, int[] upper, int flags) {

private static int repeatedStep(double[] a, int l, int r, int k, int[] upper, int flags) {
        // Adapted from Alexandrescu (2016), algorithm 8.
        final int fp;
        final int s;
        int p;
        if (flags <= MODE_SAMPLING) {
            // Median into a 12th-tile
            fp = (r - l + 1) / 12;
            // Position the sample around the target k
            s = k - mapDistance(k - l, l, r, fp);
            p = k;
        } else {
            // i in tertile [3f':6f')
            fp = (r - l + 1) / 9;
            final int f3 = 3 * fp;
            final int end = l + (f3 << 1);
            for (int i = l + f3; i < end; i++) {
                Sorting.sort3(a, i - f3, i, i + f3);
            }
            // 5th 9th-tile: [4f':5f')
            s = l + (fp << 2);
            // No adaption uses the middle to enforce strict margins
            p = s + (flags == MODE_ADAPTION ? mapDistance(k - l, l, r, fp) : (fp >>> 1));
        }
        final int e = s + fp - 1;
        for (int i = s; i <= e; i++) {
            Sorting.sort3(a, i - fp, i, i + fp);
        }
        p = quickSelectAdaptive(a, s, e, p, p, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, s, e, p, upper[0], upper);
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Assumes the range {@code r - l >= 11}; the caller is responsible for selection on a smaller
     * range.
     *
     * <p>Uses the Chen and Dumitrescu repeated step median-of-medians-of-medians algorithm
     * with the lower median of 4 then either median of 3 with the final sample placed in the
     * 5th 12th-tile, or min of 3 with the final sample in the 4th 12th-tile;
     * the pivot chosen from the sample is adaptive using the input {@code k}.
     *
     * <p>Given a pivot in the middle of the sample this has margins of 1/6 and 1/4.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @param flags Control flags.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int repeatedStepLeft(double[] a, int l, int r, int k, int[] upper, int flags) {

double vj = a[end];
        a[start] = a[left];
        a[end] = a[right];
        a[left] = vj;
        a[right] = vi;

        int i = start + 1;
        int j = end - 1;

        // Positioned for pre-in/decrement to write to pivot region
        int p0 = pivot0 == start ? i : pivot0;
        int p1 = pivot1 == end ? j : pivot1;

        while (true) {
            do {
                --i;
            } while (a[i] < v);
            do {
                ++j;
            } while (a[j] > v);
            vj = a[i];
            vi = a[j];
            a[i] = vi;
            a[j] = vj;
            // Move the equal values to pivot region
            if (vi == v) {
                a[i] = a[--p0];
                a[p0] = v;
            }
            if (vj == v) {
                a[j] = a[++p1];
                a[p1] = v;
            }
            // Termination check and finishing loops.
            // Note: This works even if pivot region is zero length (p1 == p0-1 due to
            // length 1 pivot region at either start/end) because we pre-inc/decrement
            // one side and post-inc/decrement the other side.
            if (i == left) {
                while (j < right) {
                    do {
                        ++j;
                    } while (a[j] > v);
                    final double w = a[j];

private static int sampleStep(double[] a, int l, int r, int k, int[] upper) {
        // Floyd-Rivest: use SELECT recursively on a sample of size S to get an estimate
        // for the (k-l+1)-th smallest element into a[k], biased slightly so that the
        // (k-l+1)-th element is expected to lie in the smaller set after partitioning.
        final int n = r - l + 1;
        final int ith = k - l + 1;
        final double z = Math.log(n);
        // sample size = 0.5 * n^(2/3)
        final double s = 0.5 * Math.exp(0.6666666666666666 * z);
        final double sd = 0.5 * Math.sqrt(z * s * (n - s) / n) * Integer.signum(ith - (n >> 1));
        final int ll = Math.max(l, (int) (k - ith * s / n + sd));
        final int rr = Math.min(r, (int) (k + (n - ith) * s / n + sd));
        // Sample recursion restarts from [ll, rr]
        final int p = quickSelectAdaptive(a, ll, rr, k, k, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, ll, rr, p, upper[0], upper);
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Assumes the range {@code r - l >= 8}; the caller is responsible for selection on a smaller
     * range. If using a 12th-tile for sampling then assumes {@code r - l >= 11}.
     *
     * <p>Uses the Chen and Dumitrescu repeated step median-of-medians-of-medians algorithm
     * with the median of 3 then median of 3; the final sample is placed in the
     * 5th 9th-tile; the pivot chosen from the sample is adaptive using the input {@code k}.
     *
     * <p>Given a pivot in the middle of the sample this has margins of 2/9 and 2/9.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @param flags Control flags.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int repeatedStep(double[] a, int l, int r, int k, int[] upper, int flags) {

final double w = a[great];
                a[great] = v;
                great--;
                // delay a[k] = w
                if (w < v1) {
                    a[k] = a[less];
                    a[less] = w;
                    less++;
                } else {
                    a[k] = w;
                }
            }
        }

        // Change to inclusive ends : a[less] < P1 : a[great] > P2
        less--;
        great++;
        // Move the pivots to correct locations
        a[left] = a[less];
        a[less] = v1;
        a[right] = a[great];
        a[great] = v2;

        // Record the pivot locations
        final int lower = less;
        bounds[2] = great;

        // equal elements
        // Original paper: If middle partition is bigger than a threshold
        // then check for equal elements.

        // Note: This is extra work. When performing partitioning the region of interest
        // may be entirely above or below the central region and this can be skipped.

        // Here we look for equal elements if the centre is more than 5/8 the length.
        // 5/8 = 1/2 + 1/8. Pivots must be different.
        if ((great - less) > (n >>> 1) + (n >>> 3) && v1 != v2) {

            // Fast-forward to reduce swaps. Changes inclusive ends to exclusive ends.
            // Since v1 != v2 these act as sentinels to prevent overrun.
            do {
                ++less;
            } while (a[less] == v1);
            do {
                --great;
            } while (a[great] == v2);

            // This copies the logic in the sorting loop using == comparisons
            EQUAL:
            for (int k = less - 1; ++k <= great;) {
                final double v = a[k];

final double u = a[i];
                    a[i] = a[i - f];
                    a[i - f] = u;
                }
            }
            // 2nd 12th-tile
            fp = f / 3;
            s = l + fp;
            // Lower margin has 2(d+1) elements; d == (position in sample) - s
            // Force k into the lower margin
            p = s + ((k - l) >>> 1);
        }
        final int e = s + fp - 1;
        for (int i = s; i <= e; i++) {
            Sorting.sort3(a, i - fp, i, i + fp);
        }
        p = quickSelectAdaptive(a, s, e, p, p, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, s, e, p, upper[0], upper);
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Assumes the range {@code r - l >= 11}; the caller is responsible for selection on a smaller
     * range.
     *
     * <p>Uses the Chen and Dumitrescu repeated step median-of-medians-of-medians algorithm
     * with the maximum of 4 then median of 3; the final sample is placed in the
     * 11th 12th-tile; the pivot chosen from the sample is adaptive using the input {@code k}.
     *
     * <p>Given a pivot in the middle of the sample this has margins of 1/3 and 1/12.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @param flags Control flags.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int repeatedStepFarRight(double[] a, int l, int r, int k, int[] upper, int flags) {

final double u = a[i];
                    a[i] = a[i + f];
                    a[i + f] = u;
                }
            }
            // 11th 12th-tile
            fp = f / 3;
            e = r - fp;
            // Upper margin has 2(d+1) elements; d == e - (position in sample)
            // Force k into the upper margin
            p = e - ((r - k) >>> 1);
        }
        final int s = e - fp + 1;
        for (int i = s; i <= e; i++) {
            Sorting.sort3(a, i - fp, i, i + fp);
        }
        p = quickSelectAdaptive(a, s, e, p, p, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, s, e, p, upper[0], upper);
    }

    /**
     * Expand a partition around a single pivot. Partitioning exchanges array
     * elements such that all elements smaller than pivot are before it and all
     * elements larger than pivot are after it. The central region is already
     * partitioned.
     *
     * <pre>{@code
     * |l             |s   |p0 p1|   e|                r|
     * |    ???       | <P | ==P | >P |        ???      |
     * }</pre>
     *
     * <p>This requires that {@code start != end}. However it handles
     * {@code left == start} and/or {@code end == right}.
     *
     * @param a Data array.
     * @param left Lower bound (inclusive).
     * @param right Upper bound (inclusive).
     * @param start Start of the partition range (inclusive).
     * @param end End of the partitioned range (inclusive).
     * @param pivot0 Lower pivot location (inclusive).
     * @param pivot1 Upper pivot location (inclusive).
     * @param upper Upper bound (inclusive) of the pivot range [k1].
     * @return Lower bound (inclusive) of the pivot range [k0].
     */
    // package-private for testing
    static int expandPartition(double[] a, int left, int right, int start, int end,

private static int repeatedStepFarLeft(double[] a, int l, int r, int k, int[] upper, int flags) {
        // Far step has been changed from the Alexandrescu (2016) step of lower-median-of-4, min-of-3
        // into the 4th 12th-tile to a min-of-4, median-of-3 into the 2nd 12th-tile.
        // The differences are:
        // - The upper margin when not sampling is 8/24 vs. 9/24; the lower margin remains at 1/12.
        // - The position of the sample is closer to the expected location of k < |A| / 12.
        // - Sampling mode uses a median-of-3 with adaptive k, matching the other step methods.
        //   A min-of-3 sample can create a pivot too small if used with adaption of k leaving
        //   k in the larger parition and a wasted iteration.
        // - Adaption is adjusted to force use of the lower margin when not sampling.
        final int fp;
        final int s;
        int p;
        if (flags <= MODE_SAMPLING) {
            // 2nd 12th-tile
            fp = (r - l + 1) / 12;
            s = l + fp;
            // Use adaption
            p = s + mapDistance(k - l, l, r, fp);
        } else {
            // i in 2nd quartile; min into i-f (1st quartile)
            final int f = (r - l + 1) >> 2;
            final int f2 = f + f;
            final int end = l + f2;
            for (int i = l + f; i < end; i++) {
                if (a[i + f] < a[i - f]) {
                    final double u = a[i + f];

private static int repeatedStepFarRight(double[] a, int l, int r, int k, int[] upper, int flags) {
        // Mirror image repeatedStepFarLeft
        final int fp;
        final int e;
        int p;
        if (flags <= MODE_SAMPLING) {
            // 11th 12th-tile
            fp = (r - l + 1) / 12;
            e = r - fp;
            // Use adaption
            p = e - mapDistance(r - k, l, r, fp);
        } else {
            // i in 3rd quartile; max into i+f (4th quartile)
            final int f = (r - l + 1) >> 2;
            final int f2 = f + f;
            final int end = r - f2;
            for (int i = r - f; i > end; i--) {
                if (a[i - f] > a[i + f]) {
                    final double u = a[i - f];

final double w = x[c];
        if (w < u) {
            if (v < w) {
                x[a] = v;
                x[b] = w;
                x[c] = u;
                return;
            }
            if (u < v) {
                x[a] = w;
                x[b] = u;
                x[c] = v;
                return;
            }
            // w < v < u
            x[a] = w;
            x[c] = u;
            return;
        }
        if (v < u) {
            // v < u < w
            x[a] = v;
            x[b] = u;
            return;
        }
        if (w < v) {
            // u < w < v
            x[b] = w;
            x[c] = v;
        }
        // u < v < w
    }

    /**
     * Sorts the elements at the given distinct indices in an array.
     *
     * @param x Data array.
     * @param a Index.
     * @param b Index.
     * @param c Index.
     * @param d Index.
     * @param e Index.
     */
    static void sort5(double[] x, int a, int b, int c, int d, int e) {

p = s + (flags == MODE_ADAPTION ? mapDistance(k - l, l, r, fp) : (fp >>> 1));
        }
        final int e = s + fp - 1;
        for (int i = s; i <= e; i++) {
            Sorting.sort3(a, i - fp, i, i + fp);
        }
        p = quickSelectAdaptive(a, s, e, p, p, upper, MODE_FR_SAMPLING);
        return expandPartition(a, l, r, s, e, p, upper[0], upper);
    }

    /**
     * Partition an array slice around a pivot. Partitioning exchanges array elements such
     * that all elements smaller than pivot are before it and all elements larger than
     * pivot are after it.
     *
     * <p>Assumes the range {@code r - l >= 11}; the caller is responsible for selection on a smaller
     * range.
     *
     * <p>Uses the Chen and Dumitrescu repeated step median-of-medians-of-medians algorithm
     * with the lower median of 4 then either median of 3 with the final sample placed in the
     * 5th 12th-tile, or min of 3 with the final sample in the 4th 12th-tile;
     * the pivot chosen from the sample is adaptive using the input {@code k}.
     *
     * <p>Given a pivot in the middle of the sample this has margins of 1/6 and 1/4.
     *
     * @param a Data array.
     * @param l Lower bound (inclusive).
     * @param r Upper bound (inclusive).
     * @param k Target index.
     * @param upper Upper bound (inclusive) of the pivot range.
     * @param flags Control flags.
     * @return Lower bound (inclusive) of the pivot range.
     */
    private static int repeatedStepLeft(double[] a, int l, int r, int k, int[] upper, int flags) {

private static int partition(double[] a, int left, int right, int[] bounds) {
        // Pick 2 pivots from 5 approximately uniform through the range.
        // Spacing is ~ 1/7 made using shifts. Other strategies are equal or much
        // worse. 1/7 = 5/35 ~ 1/8 + 1/64 : 0.1429 ~ 0.1406
        // Ensure the value is above zero to choose different points!
        final int n = right - left;
        final int step = 1 + (n >>> 3) + (n >>> 6);
        final int i3 = left + (n >>> 1);
        final int i2 = i3 - step;
        final int i1 = i2 - step;
        final int i4 = i3 + step;
        final int i5 = i4 + step;
        Sorting.sort5(a, i1, i2, i3, i4, i5);

        // Partition data using pivots P1 and P2 into less-than, greater-than or between.
        // Pivot values P1 & P2 are placed at the end. If P1 < P2, P2 acts as a sentinel.
        // k traverses the unknown region ??? and values moved if less-than or
        // greater-than:
        //
        // left        less              k       great         right
        // |P1|  <P1   |   P1 <= & <= P2 |    ???    |    >P2   |P2|
        //
        // <P1            (left, lt)
        // P1 <= & <= P2  [lt, k)
        // >P2            (gt, right)
        //
        // At the end pivots are swapped back to behind the less and great pointers.
        //
        // |  <P1        |P1|     P1<= & <= P2    |P2|      >P2    |

        // Swap ends to the pivot locations.
        final double v1 = a[i2];

final double w = a[great];
                    a[great] = v;
                    great--;
                    if (w == v1) {
                        a[k] = a[less];
                        a[less] = w;
                        less++;
                    } else {
                        a[k] = w;
                    }
                }
            }

            // Change to inclusive ends
            less--;
            great++;
        }

        // Between pivots in (less, great)
        if (v1 != v2 && less < great - 1) {
            // Record the pivot end points
            bounds[0] = less;
            bounds[1] = great;
        } else {
            // No unsorted internal region (set k1 = k3; k2 = k0)
            bounds[0] = bounds[2];
            bounds[1] = lower;
        }

        return lower;
    }