/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  
  package org.apache.commons.numbers.examples.jmh.arrays;
  
  import java.util.Arrays;
  import java.util.BitSet;
  
  /**
   * A K<sup>th</sup> selector implementation to pick up the K<sup>th</sup> ordered element
   * from a data array containing the input numbers. Uses a partial sort of the data.
   *
   * <p>Note: The search may use a cache of pivots. A pivot is a K<sup>th</sup> index that
   * corresponds to a value correctly positioned within the equivalent fully sorted data array.
   * Each step of the algorithm partitions the array between a lower and upper bound into
   * values above or below a value chosen from the interval. The index of this value after the
   * partition is stored as a pivot. A subsequent search into the array can use the pivots to
   * quickly bracket the interval for the next target index.
   *
   * <p>The maximum supported cache size is 2^30 - 1. This ensures that any valid index i
   * into the cache can be increased for the next level using 2 * i + 2 without overflow.
   *
   * <p>The method for choosing the value within the interval to pivot around is specified by
   * the {@link PivotingStrategy}. Ideally this should provide a guess of the middle (median)
   * of the interval. The value must be within the interval, thus using for example the
   * mean of the end values is not an option. The default uses a guess for the median from
   * 3 values that are likely to be representative of the range of values.
   *
   * <p>This implementation provides the option to return the (K+1) ordered element along with
   * the K<sup>th</sup>. This uses the position of K and the most recent upper bound on the
   * bracket known to contain values greater than the K<sup>th</sup>. This prevents using a
   * second search into the array for the (K+1) element.
   *
   * @since 1.2
   */
  class KthSelector {
      /** Empty pivots array. */
      static final int[] NO_PIVOTS = {};
      /** Minimum selection size for insertion sort rather than selection.
       * Dual-pivot quicksort used 27 in the original paper. */
      private static final int MIN_SELECT_SIZE = 17;
  
      /** A {@link PivotingStrategy} used for pivoting. */
      private final PivotingStrategy pivotingStrategy;
  
      /** Minimum selection size for insertion sort rather than selection. */
      private final int minSelectSize;
  
      /**
       * Constructor with default {@link PivotingStrategy#MEDIAN_OF_3 median of 3} pivoting
       * strategy.
       */
      KthSelector() {
          this(PivotingStrategy.MEDIAN_OF_3);
      }
  
      /**
       * Constructor with specified pivoting strategy.
       *
       * @param pivotingStrategy Pivoting strategy to use.
       */
      KthSelector(PivotingStrategy pivotingStrategy) {
          this(pivotingStrategy, MIN_SELECT_SIZE);
      }
  
      /**
       * Constructor with specified pivoting strategy and select size.
       *
       * @param pivotingStrategy Pivoting strategy to use.
       * @param minSelectSize Minimum selection size for insertion sort rather than selection.
       */
      KthSelector(PivotingStrategy pivotingStrategy, int minSelectSize) {
          this.pivotingStrategy = pivotingStrategy;
          this.minSelectSize = minSelectSize;
      }
  
      /**
       * Select K<sup>th</sup> value in the array. Optionally select the next value after K.
       *
       * <p>Note: If K+1 is requested this method assumes it is a valid index into the array.
       *
       * <p>Uses a single-pivot partition method.
       *
       * @param data Data array to use to find out the K<sup>th</sup> value.
       * @param k Index whose value in the array is of interest.
      * @param kp1 K+1<sup>th</sup> value (if not null)
      * @return K<sup>th</sup> value
      */
     double selectSP(double[] data, int k, double[] kp1) {
         int begin = 0;
         int end = data.length;
         while (end - begin > minSelectSize) {
             // Select a pivot and partition data array around it
             final int pivot = partitionSP(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k));
             if (k == pivot) {
                 // The pivot was exactly the element we wanted
                 return finalSelection(data, k, kp1, end);
             } else if (k < pivot) {
                 // The element is in the left partition
                 end = pivot;
             } else {
                 // The element is in the right partition
                 begin = pivot + 1;
             }
         }
         sortRange(data, begin, end);
         if (kp1 != null) {
             // Either end == data.length and k+1 is sorted; or
             // end == pivot where data[k] <= data[pivot] <= data[pivot+j] for all j
             kp1[0] = data[k + 1];
         }
         return data[k];
     }
 
     /**
      * Select K<sup>th</sup> value in the array. Optionally select the next value after K.
      *
      * <p>Note: If K+1 is requested this method assumes it is a valid index into the array.
      *
      * <p>Uses a single-pivot partition method with special handling of NaN and signed zeros.
      * Correction of signed zeros requires a sweep across the entire range.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param k Index whose value in the array is of interest.
      * @param kp1 K+1<sup>th</sup> value (if not null)
      * @return K<sup>th</sup> value
      */
     double selectSPN(double[] data, int k, double[] kp1) {
         // Handle NaN
         final int length = sortNaN(data);
         if (k >= length) {
             if (kp1 != null) {
                 kp1[0] = Double.NaN;
             }
             return Double.NaN;
         }
 
         int begin = 0;
         int end = length;
         while (end - begin > minSelectSize) {
             // Select a pivot and partition data array around it
             final int pivot = partitionSPN(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k));
             if (k == pivot) {
                 // The pivot was exactly the element we wanted
                 if (data[k] == 0) {
                     orderSignedZeros(data, 0, length);
                 }
                 return finalSelection(data, k, kp1, end);
             } else if (k < pivot) {
                 // The element is in the left partition
                 end = pivot;
             } else {
                 // The element is in the right partition
                 begin = pivot + 1;
             }
         }
         insertionSort(data, begin, end, begin != 0);
         if (data[k] == 0) {
             orderSignedZeros(data, 0, length);
         }
         if (kp1 != null) {
             // Either end == data.length and k+1 is sorted; or
             // end == pivot where data[k] <= data[pivot] <= data[pivot+j] for all j
             kp1[0] = data[k + 1];
         }
         return data[k];
     }
 
     /**
      * Select K<sup>th</sup> value in the array. Optionally select the next value after K.
      *
      * <p>Note: If K+1 is requested this method assumes it is a valid index into the array.
      *
      * <p>Uses a single-pivot partition method with a heap cache.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots. Maximum supported heap size is 2^30 - 1.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param pivotsHeap Cached pivots heap that can be used for efficient estimation.
      * @param k Index whose value in the array is of interest.
      * @param kp1 K+1<sup>th</sup> value (if not null)
      * @return K<sup>th</sup> value
      */
     double selectSPH(double[] data, int[] pivotsHeap, int k, double[] kp1) {
         final int heapLength = pivotsHeap.length;
         if (heapLength == 0) {
             // No pivots
             return selectSP(data, k, kp1);
         }
         int begin = 0;
         int end = data.length;
         int node = 0;
         while (end - begin > minSelectSize) {
             int pivot;
 
             if (node < heapLength && pivotsHeap[node] >= 0) {
                 // The pivot has already been found in a previous call
                 // and the array has already been partitioned around it
                 pivot = pivotsHeap[node];
             } else {
                 // Select a pivot and partition data array around it
                 pivot = partitionSP(data, begin, end,
                     pivotingStrategy.pivotIndex(data, begin, end - 1, k));
                 if (node < heapLength) {
                     pivotsHeap[node] = pivot;
                 }
             }
 
             if (k == pivot) {
                 // The pivot was exactly the element we wanted
                 return finalSelection(data, k, kp1, end);
             } else if (k < pivot) {
                 // The element is in the left partition
                 end = pivot;
                 if (node < heapLength) {
                     node = Math.min((node << 1) + 1, heapLength);
                 }
             } else {
                 // The element is in the right partition
                 begin = pivot + 1;
                 if (node < heapLength) {
                     node = Math.min((node << 1) + 2, heapLength);
                 }
             }
         }
         sortRange(data, begin, end);
         if (kp1 != null) {
             // Either end == data.length and k+1 is sorted; or
             // end == pivot where data[k] <= data[pivot] <= data[pivot+j] for all j
             kp1[0] = data[k + 1];
         }
         return data[k];
     }
 
     /**
      * Select K<sup>th</sup> value in the array. Optionally select the next value after K.
      *
      * <p>Note: If K+1 is requested this method assumes it is a valid index into the array
      * (i.e. K is not the last index in the array).
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param k Index whose value in the array is of interest.
      * @param kp1 K+1<sup>th</sup> value (if not null)
      * @param end Upper bound (exclusive) of the interval containing K.
      * This should be either a pivot point {@code data[k] <= data[end]} or the length
      * of the data array.
      * @return K<sup>th</sup> value
      */
     private static double finalSelection(double[] data, int k, double[] kp1, int end) {
         if (kp1 != null) {
             // After partitioning all elements above k are greater than or equal to k.
             // Find the minimum of the elements above.
             // Set the k+1 limit as either a pivot or the end of the data.
             final int limit = Math.min(end, data.length - 1);
             double min = data[k + 1];
             for (int i = k + 2; i <= limit; i++) {
                 if (DoubleMath.lessThan(data[i], min)) {
                     min = data[i];
                 }
             }
             kp1[0] = min;
         }
         return data[k];
     }
 
     /**
      * Partition the array such that indices {@code k} correspond to their correctly
      * sorted value in the equivalent fully sorted array. For all indices {@code k}
      * and any index {@code i}:
      *
      * <pre>{@code
      * data[i < k] <= data[k] <= data[k < i]
      * }</pre>
      *
      * <p>Uses a single-pivot partition method.
      *
      * @param data Values.
      * @param k Indices.
      */
     void partitionSP(double[] data, int... k) {
         final int n = k.length;
         if (n <= 1) {
             if (n == 1) {
                 selectSP(data, k[0], null);
             }
             return;
         }
         // Multiple pivots
         final int length = data.length;
         final BitSet pivots = new BitSet(length);
 
         for (int i = 0; i < n; i++) {
             final int kk = k[i];
             if (pivots.get(kk)) {
                 // Already sorted
                 continue;
             }
             int begin;
             int end;
             if (i == 0) {
                 begin = 0;
                 end = length;
             } else {
                 // Start inclusive
                 begin = pivots.previousSetBit(kk) + 1;
                 end = pivots.nextSetBit(kk + 1);
                 if (end < 0) {
                     end = length;
                 }
             }
             partitionSP(data, begin, end, pivots, kk);
         }
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Uses a single-pivot partition method.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivots Cache of pivot points.
      * @param k Index whose value in the array is of interest.
      */
     private void partitionSP(double[] data, int begin, int end, BitSet pivots, int k) {
         // Find the unsorted range containing k
         while (end - begin > minSelectSize) {
             // Select a value and partition data array around it
             final int pivot = partitionSP(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k));
             pivots.set(pivot);
             if (k == pivot) {
                 // The pivot was exactly the element we wanted
                 return;
             } else if (k < pivot) {
                 // The element is in the left partition
                 end = pivot;
             } else {
                 // The element is in the right partition
                 begin = pivot + 1;
             }
         }
         setPivots(begin, end, pivots);
         sortRange(data, begin, end);
     }
 
     /**
      * 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>Uses a single-pivot partition method.
      *
      * @param data Data array.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivot Initial index of the pivot.
      * @return index of the pivot after partition
      */
     private static int partitionSP(double[] data, int begin, int end, int pivot) {
         final double value = data[pivot];
         data[pivot] = data[begin];
 
         int i = begin + 1;
         int j = end - 1;
         while (i < j) {
             while (i < j && DoubleMath.greaterThan(data[j], value)) {
                 --j;
             }
             while (i < j && DoubleMath.lessThan(data[i], value)) {
                 ++i;
             }
             if (i < j) {
                 final double tmp = data[i];
                 data[i++] = data[j];
                 data[j--] = tmp;
             }
         }
 
         if (i >= end || DoubleMath.greaterThan(data[i], value)) {
             --i;
         }
         data[begin] = data[i];
         data[i] = value;
         return i;
     }
 
     /**
      * Partition the array such that indices {@code k} correspond to their correctly
      * sorted value in the equivalent fully sorted array. For all indices {@code k}
      * and any index {@code i}:
      *
      * <pre>{@code
      * data[i < k] <= data[k] <= data[k < i]
      * }</pre>
      *
      * <p>Uses a single-pivot partition method.
      *
      * @param data Values.
      * @param k Indices.
      */
     void partitionSPN(double[] data, int... k) {
         final int n = k.length;
         if (n <= 1) {
             if (n == 1) {
                 selectSPN(data, k[0], null);
             }
             return;
         }
         // Multiple pivots
 
         // Handle NaN
         final int length = sortNaN(data);
         if (length < 1) {
             return;
         }
 
         final BitSet pivots = new BitSet(length);
 
         // Flag any pivots that are zero
         boolean zeros = false;
         for (int i = 0; i < n; i++) {
             final int kk = k[i];
             if (kk >= length || pivots.get(kk)) {
                 // Already sorted
                 continue;
             }
             int begin;
             int end;
             if (i == 0) {
                 begin = 0;
                 end = length;
             } else {
                 // Start inclusive
                 begin = pivots.previousSetBit(kk) + 1;
                 end = pivots.nextSetBit(kk + 1);
                 if (end < 0) {
                     end = length;
                 }
             }
             partitionSPN(data, begin, end, pivots, kk);
             zeros = zeros || data[kk] == 0;
         }
 
         // Handle signed zeros
         if (zeros) {
             orderSignedZeros(data, 0, length);
         }
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values. Does not partition
      * around signed zeros.
      *
      * <p>Uses a single-pivot partition method.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivots Cache of pivot points.
      * @param k Index whose value in the array is of interest.
      */
     private void partitionSPN(double[] data, int begin, int end, BitSet pivots, int k) {
         // Find the unsorted range containing k
         while (end - begin > minSelectSize) {
             // Select a value and partition data array around it
             final int pivot = partitionSPN(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k));
             pivots.set(pivot);
             if (k == pivot) {
                 // The pivot was exactly the element we wanted
                 return;
             } else if (k < pivot) {
                 // The element is in the left partition
                 end = pivot;
             } else {
                 // The element is in the right partition
                 begin = pivot + 1;
             }
         }
         setPivots(begin, end, pivots);
         sortRange(data, begin, end);
     }
 
     /**
      * 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>Note: Requires that the range contains no NaN values. Does not partition
      * around signed zeros.
      *
      * <p>Uses a single-pivot partition method.
      *
      * @param data Data array.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivot Initial index of the pivot.
      * @return index of the pivot after partition
      */
     private static int partitionSPN(double[] data, int begin, int end, int pivot) {
         final double value = data[pivot];
         data[pivot] = data[begin];
 
         int i = begin + 1;
         int j = end - 1;
         while (i < j) {
             while (i < j && data[j] > value) {
                 --j;
             }
             while (i < j && data[i] < value) {
                 ++i;
             }
             if (i < j) {
                 final double tmp = data[i];
                 data[i++] = data[j];
                 data[j--] = tmp;
             }
         }
 
         if (i >= end || data[i] > value) {
             --i;
         }
         data[begin] = data[i];
         data[i] = value;
         return i;
     }
 
     /**
      * Sort an array range.
      *
      * @param data Data array.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private static void sortRange(double[] data, int begin, int end) {
         Arrays.sort(data, begin, end);
     }
 
     /**
      * Sorts an array using an insertion sort.
      *
      * <p>Note: Requires that the range contains no NaN values. It does not respect the
      * order of signed zeros.
      *
      * <p>This method is fast up to approximately 40 - 80 values.
      *
      * <p>The {@code internal} flag indicates that the value at {@code data[begin - 1]}
      * is sorted.
      *
      * @param data Data array.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param internal Internal flag.
      */
     private static void insertionSort(double[] data, int begin, int end, boolean internal) {
         Sorting.sort(data, begin, end - 1, internal);
     }
 
     /**
      * Move NaN values to the end of the array.
      * This allows all other values to be compared using {@code <, ==, >} operators (with
      * the exception of signed zeros).
      *
      * @param data Values.
      * @return end of non-NaN values
      */
     private static int sortNaN(double[] data) {
         int end = data.length;
         // Avoid unnecessary moves
         while (--end > 0) {
             if (!Double.isNaN(data[end])) {
                 break;
             }
         }
         end++;
         for (int i = end; i > 0;) {
             final double v = data[--i];
             if (Double.isNaN(v)) {
                 data[i] = data[--end];
                 data[end] = v;
             }
         }
         return end;
     }
 
     /**
      * Detect and fix the sort order of signed zeros. Assumes the data may have been
      * partially ordered around zero.
      *
      * <p>Searches for zeros if {@code data[begin] <= 0} and {@code data[end - 1] >= 0}.
      * If zeros are discovered in the range then they are assumed to be continuous.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private static void fixSignedZeros(double[] data, int begin, int end) {
         int j;
         if (data[begin] <= 0 && data[end - 1] >= 0) {
             int i = begin;
             while (data[i] < 0) {
                 i++;
             }
             j = end - 1;
             while (data[j] > 0) {
                 j--;
             }
             sortZero(data, i, j + 1);
         }
     }
 
     /**
      * Count the number of signed zeros in the range and order them to be correctly
      * sorted. This checks all values in the range. It does not assume zeros are continuous.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private static void orderSignedZeros(double[] data, int begin, int end) {
         int c = countSignedZeros(data, begin, end);
         if (c != 0) {
             int i = begin - 1;
             while (++i < end) {
                 if (data[i] == 0) {
                     data[i] = -0.0;
                     if (--c == 0) {
                         break;
                     }
                 }
             }
             while (++i < end) {
                 if (data[i] == 0) {
                     data[i] = 0.0;
                 }
             }
         }
     }
 
     /**
      * Count the number of signed zeros (-0.0).
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @return the count
      */
     static int countSignedZeros(double[] data, int begin, int end) {
         // Count negative zeros
         int c = 0;
         for (int i = begin; i < end; i++) {
             if (data[i] == 0 && Double.doubleToRawLongBits(data[i]) < 0) {
                 c++;
             }
         }
         return c;
     }
 
     /**
      * Sort a range of all zero values.
      * This orders -0.0 before 0.0.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     static void sortZero(double[] data, int begin, int end) {
         // Count negative zeros
         int c = 0;
         for (int i = begin; i < end; i++) {
             if (Double.doubleToRawLongBits(data[i]) < 0) {
                 c++;
             }
         }
         // Replace
         if (c != 0) {
             int i = begin;
             while (c-- > 0) {
                 data[i++] = -0.0;
             }
             while (i < end) {
                 data[i++] = 0.0;
             }
         }
     }
 
     /**
      * Sets the pivots.
      *
      * @param from Start (inclusive)
      * @param to End (exclusive)
      * @param pivots the pivots
      */
     private static void setPivots(int from, int to, BitSet pivots) {
         if (from + 1 == to) {
             pivots.set(from);
         } else {
             pivots.set(from, to);
         }
     }
 
     /**
      * Partition the array such that indices {@code k} correspond to their correctly
      * sorted value in the equivalent fully sorted array. For all indices {@code k}
      * and any index {@code i}:
      *
      * <pre>{@code
      * data[i < k] <= data[k] <= data[k < i]
      * }</pre>
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method by Sedgewick.
      *
      * @param data Values.
      * @param k Indices.
      */
     void partitionSBM(double[] data, int... k) {
         final int n = k.length;
         if (n < 1) {
             return;
         }
 
         // Handle NaN
         final int length = sortNaN(data);
         if (length < 1) {
             return;
         }
 
         if (n == 1) {
             if (k[0] < length) {
                 partitionSBM(data, 0, length, k[0]);
             }
             return;
         }
         // Special case for partition around adjacent indices (for interpolation)
         if (n == 2 && k[0] + 1 == k[1]) {
             if (k[0] < length) {
                 final int p = partitionSBM(data, 0, length, k[0]);
                 if (p > k[1]) {
                     partitionMin(data, k[1], p);
                 }
             }
             return;
         }
 
         // To partition all k requires not moving any pivot k after it has been
         // processed. This is supported using two strategies:
         //
         // 1. Processing k in sorted order:
         // (k1, end), (k2, end), (k3, end), ... , k1 <= k2 <= k3
         // This can reorder each region during processing without destroying sorted k.
         //
         // 2. Processing unique k and visiting array regions only once:
         // Pre-process the pivots to make them unique and store the entire sorted
         // region between the end pivots (k1, kn) in a BitSet type structure:
         // |k1|......|k2|....|p|k3|k4|pppp|......|kn|
         // k can be processed in any order, e.g. k3. We use already sorted regions
         // |p| to bracket the search for each k, and skip k that are already sorted (k4).
         // Worst case storage cost is Order(N / 64).
         // The advantage is never visiting any part of the array twice. If the pivots
         // saturate the entire range then performance degrades to the speed of
         // the sort of the entire array.
 
         // Multiple pivots
         final BitSet pivots = new BitSet(length);
 
         for (int i = 0; i < n; i++) {
             final int kk = k[i];
             if (kk >= length || pivots.get(kk)) {
                 // Already sorted
                 continue;
             }
             int begin;
             int end;
             if (i == 0) {
                 begin = 0;
                 end = length;
             } else {
                 // Start inclusive
                 begin = pivots.previousSetBit(kk) + 1;
                 end = pivots.nextSetBit(kk + 1);
                 if (end < 0) {
                     end = length;
                 }
             }
             partitionSBM(data, begin, end, pivots, kk);
         }
     }
 
     /**
      * Move the minimum value to the start of the range.
      *
      * <p>Note: Requires that the range contains no NaN values.
      * Does not respect the ordering of signed zeros.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     static void partitionMin(double[] data, int begin, int end) {
         int i = begin;
         double min = data[i];
         int j = i;
         while (++i < end) {
             if (data[i] < min) {
                 min = data[i];
                 j = i;
             }
         }
         //swap(data, begin, j)
         data[j] = data[begin];
         data[begin] = min;
     }
 
     /**
      * Move the maximum value to the end of the range.
      *
      * <p>Note: Requires that the range contains no NaN values.
      * Does not respect the ordering of signed zeros.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     static void partitionMax(double[] data, int begin, int end) {
         int i = end - 1;
         double max = data[i];
         int j = i;
         while (--i >= begin) {
             if (data[i] > max) {
                 max = data[i];
                 j = i;
             }
         }
         //swap(data, end - 1, j)
         data[j] = data[end - 1];
         data[end - 1] = max;
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method by Sedgewick.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivots Cache of pivot points.
      * @param k Index whose value in the array is of interest.
      */
     private void partitionSBM(double[] data, int begin, int end, BitSet pivots, int k) {
         // Find the unsorted range containing k
         final int[] upper = {0};
         while (end - begin > minSelectSize) {
             // Select a value and partition data array around it
             final int from = partitionSBM(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k), upper);
             final int to = upper[0];
             setPivots(from, to, pivots);
             if (k >= to) {
                 // The element is in the right partition
                 begin = to;
             } else if (k < from) {
                 // The element is in the left partition
                 end = from;
             } else {
                 // The range contains the element we wanted
                 return;
             }
         }
         setPivots(begin, end, pivots);
         insertionSort(data, begin, end, begin != 0);
         fixSignedZeros(data, begin, end);
     }
 
     /**
      * 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>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method by Sedgewick.
      *
      * @param data Data array.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivot Initial index of the pivot.
      * @param upper Upper bound (exclusive) of the pivot range.
      * @return Lower bound (inclusive) of the pivot range.
      */
     private static int partitionSBM(double[] data, int begin, int end, int pivot, int[] upper) {
         // Single-pivot Bentley-McIlroy quicksort handling equal keys (Sedgewick's algorithm).
         //
         // Partition data using pivot P into less-than, greater-than or equal.
         // P is placed at the end to act as a sentinal.
         // k traverses the unknown region ??? and values moved if equal (l) or greater (g):
         //
         // left    p       i            j         q    right
         // |  ==P  |  <P   |     ???    |   >P    | ==P  |P|
         //
         // At the end P and additional equal values are swapped back to the centre.
         //
         // |         <P        | ==P |            >P        |
         //
         // Adapted from Sedgewick "Quicksort is optimal"
         // https://sedgewick.io/wp-content/themes/sedgewick/talks/2002QuicksortIsOptimal.pdf
         //
         // The algorithm has been changed so that:
         // - A pivot point must be provided.
         // - An edge case where the search meets in the middle is handled.
         // - Equal value data is not swapped to the end. Since the value is fixed then
         //   only the less than / greater than value must be moved from the end inwards.
         //   The end is then assumed to be the equal value. This would not work with
         //   object references. Equivalent swap calls are commented.
         // - Added a fast-forward over initial range containing the pivot.
 
         final int l = begin;
         final int r = end - 1;
 
         int p = l;
         int q = r;
 
         // Use the pivot index to set the upper sentinal value
         final double v = data[pivot];
         data[pivot] = data[r];
         data[r] = v;
 
         // Special case: count signed zeros
         int c = 0;
         if (v == 0) {
             c = countSignedZeros(data, begin, end);
         }
 
         // Fast-forward over equal regions to reduce swaps
         while (data[p] == v) {
             if (++p == q) {
                 // Edge-case: constant value
                 if (c != 0) {
                     sortZero(data, begin, end);
                 }
                 upper[0] = end;
                 return begin;
             }
         }
         // Cannot overrun as the prior scan using p stopped before the end
         while (data[q - 1] == v) {
             q--;
         }
 
         int i = p - 1;
         int j = q;
 
         for (;;) {
             do {
                 ++i;
             } while (data[i] < v);
             while (v < data[--j]) {
                 // Stop at l (not i) allows scan loops to be independent
                 if (j == l) {
                     break;
                 }
             }
             if (i >= j) {
                 // Edge-case if search met on an internal pivot value
                 // (not at the greater equal region, i.e. i < q).
                 // Move this to the lower-equal region.
                 if (i == j && v == data[i]) {
                     //swap(data, i++, p++)
                     //data[i++] = data[p++];
                     data[i++] = data[p];
                     data[p++] = v;
                 }
                 break;
             }
             //swap(data, i, j)
             final double vj = data[i];
             final double vi = data[j];
             data[i] = vi;
             data[j] = vj;
             if (vi == v) {
                 //swap(data, i, p++)
                 //data[i] = data[p++];
                 data[i] = data[p];
                 data[p++] = v;
             }
             if (vj == v) {
                 //swap(data, j, --q)
                 data[j] = data[--q];
                 data[q] = v;
             }
         }
         // i is at the end (exclusive) of the less-than region
 
         // Place pivot value in centre
         //swap(data, r, i)
         data[r] = data[i];
         data[i] = v;
 
         // Move equal regions to the centre.
         // Set the pivot range [j, i) and move this outward for equal values.
         j = i++;
 
         // less-equal:
         //   for (int k = l; k < p; k++):
         //     swap(data, k, --j)
         // greater-equal:
         //   for (int k = r; k-- > q; i++) {
         //     swap(data, k, i)
 
         // Move the minimum of less-equal or less-than
         int move = Math.min(p - l, j - p);
         final int lower = j - (p - l);
         for (int k = l; move-- > 0; k++) {
             data[k] = data[--j];
             data[j] = v;
         }
         // Move the minimum of greater-equal or greater-than
         move = Math.min(r - q, q - i);
         upper[0] = i + (r - q);
         for (int k = r; move-- > 0; i++) {
             data[--k] = data[i];
             data[i] = v;
         }
 
         // Special case: fixed signed zeros
         if (c != 0) {
             p = lower;
             while (c-- > 0) {
                 data[p++] = -0.0;
             }
             while (p < upper[0]) {
                 data[p++] = 0.0;
             }
         }
 
         return lower;
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method by Sedgewick.
      *
      * <p>Returns the last known pivot location adjacent to K.
      * If {@code p <= k} the range [p, min{k+2, data.length}) is sorted.
      * If {@code p > k} then p is a pivot.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param k Index whose value in the array is of interest.
      * @return the bound index
      */
     private int partitionSBM(double[] data, int begin, int end, int k) {
         // Find the unsorted range containing k
         final int[] upper = {0};
         while (end - begin > minSelectSize) {
             // Select a value and partition data array around it
             final int from = partitionSBM(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k), upper);
             final int to = upper[0];
             if (k >= to) {
                 // The element is in the right partition
                 begin = to;
             } else if (k < from) {
                 // The element is in the left partition
                 end = from;
             } else {
                 // The range contains the element we wanted
                 return end;
             }
         }
         insertionSort(data, begin, end, begin != 0);
         fixSignedZeros(data, begin, end);
         // Either end == data.length and k+1 is sorted; or
         // end == pivot and k+1 is sorted
         return begin;
     }
 
     /**
      * Partition the array such that indices {@code k} correspond to their correctly
      * sorted value in the equivalent fully sorted array. For all indices {@code k}
      * and any index {@code i}:
      *
      * <pre>{@code
      * data[i < k] <= data[k] <= data[k < i]
      * }</pre>
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method.
      *
      * @param data Values.
      * @param k Indices.
      */
     void partitionBM(double[] data, int... k) {
         final int n = k.length;
         if (n < 1) {
             return;
         }
 
         // Handle NaN
         final int length = sortNaN(data);
         if (length < 1) {
             return;
         }
 
         if (n == 1) {
             partitionBM(data, 0, length, k[0]);
             return;
         }
         // Special case for partition around adjacent indices (for interpolation)
         if (n == 2 && k[0] + 1 == k[1]) {
             final int p = partitionBM(data, 0, length, k[0]);
             if (p > k[1]) {
                 partitionMin(data, k[1], p);
             }
             return;
         }
 
         // Multiple pivots
         final BitSet pivots = new BitSet(length);
 
         for (int i = 0; i < n; i++) {
             final int kk = k[i];
             if (kk >= length || pivots.get(kk)) {
                 // Already sorted
                 continue;
             }
             int begin;
             int end;
             if (i == 0) {
                 begin = 0;
                 end = length;
             } else {
                 // Start inclusive
                 begin = pivots.previousSetBit(kk) + 1;
                 end = pivots.nextSetBit(kk + 1);
                 if (end < 0) {
                     end = length;
                 }
             }
             partitionBM(data, begin, end, pivots, kk);
         }
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivots Cache of pivot points.
      * @param k Index whose value in the array is of interest.
      */
     private void partitionBM(double[] data, int begin, int end, BitSet pivots, int k) {
         // Find the unsorted range containing k
         final int[] upper = {0};
         while (end - begin > minSelectSize) {
             // Select a value and partition data array around it
             final int from = partitionBM(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k), upper);
             final int to = upper[0];
             setPivots(from, to, pivots);
             if (k >= to) {
                 // The element is in the right partition
                 begin = to;
             } else if (k < from) {
                 // The element is in the left partition
                 end = from;
             } else {
                 // The range contains the element we wanted
                 return;
             }
         }
         setPivots(begin, end, pivots);
         insertionSort(data, begin, end, begin != 0);
         fixSignedZeros(data, begin, end);
     }
 
     /**
      * 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>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method.
      *
      * @param data Data array.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivot Initial index of the pivot.
      * @param upper Upper bound (exclusive) of the pivot range.
      * @return Lower bound (inclusive) of the pivot range.
      */
     private static int partitionBM(double[] data, int begin, int end, int pivot, int[] upper) {
         // Partition method handling equal keys, Bentley-McIlroy quicksort.
         //
         // Adapted from program 7 in Bentley-McIlroy (1993)
         // Engineering a sort function
         // SOFTWARE—PRACTICE AND EXPERIENCE, VOL.23(11), 1249–1265
         //
         // 3-way partition of the data using a pivot value into
         // less-than, equal or greater-than.
         //
         // First partition data into 4 reqions by scanning the unknown region from
         // left (i) and right (j) and moving equal values to the ends:
         //                  i->       <-j                  end
         // l        p       |           |         q       r|
         // | equal  | less  |  unknown  | greater | equal ||
         //
         //                    <-j                          end
         // l        p             i               q       r|
         // | equal  | less        |       greater | equal ||
         //
         // Then the equal values are copied from the ends to the centre:
         // | less        |        equal      |     greater |
 
         final int l = begin;
         final int r = end - 1;
 
         int i = l;
         int j = r;
         int p = l;
         int q = r;
 
         final double v = data[pivot];
 
         // Special case: count signed zeros
         int c = 0;
         if (v == 0) {
             c = countSignedZeros(data, begin, end);
         }
 
         for (;;) {
             while (i <= j && data[i] <= v) {
                 if (data[i] == v) {
                     //swap(data, i, p++)
                     data[i] = data[p];
                     data[p++] = v;
                 }
                 i++;
             }
             while (j >= i && data[j] >= v) {
                 if (v == data[j]) {
                     //swap(data, j, q--)
                     data[j] = data[q];
                     data[q--] = v;
                 }
                 j--;
             }
             if (i > j) {
                 break;
             }
             swap(data, i++, j--);
         }
 
         // Move equal regions to the centre.
         int s = Math.min(p - l, i - p);
         for (int k = l; s > 0; k++, s--) {
             //swap(data, k, i - s)
             data[k] = data[i - s];
             data[i - s] = v;
         }
         s = Math.min(q - j, r - q);
         for (int k = i; s > 0; k++, s--) {
             //swap(data, end - s, k)
             data[end - s] = data[k];
             data[k] = v;
         }
 
         // Set output range
         i = i - p + l;
         j = j - q + end;
         upper[0] = j;
 
         // Special case: fixed signed zeros
         if (c != 0) {
             p = i;
             while (c-- > 0) {
                 data[p++] = -0.0;
             }
             while (p < j) {
                 data[p++] = 0.0;
             }
         }
 
         return i;
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method.
      *
      * <p>Returns the last known pivot location adjacent to K.
      * If {@code p <= k} the range [p, min{k+2, data.length}) is sorted.
      * If {@code p > k} then p is a pivot.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param k Index whose value in the array is of interest.
      * @return the bound index
      */
     private int partitionBM(double[] data, int begin, int end, int k) {
         // Find the unsorted range containing k
         final int[] upper = {0};
         while (end - begin > minSelectSize) {
             // Select a value and partition data array around it
             final int from = partitionBM(data, begin, end,
                 pivotingStrategy.pivotIndex(data, begin, end - 1, k), upper);
             final int to = upper[0];
             if (k >= to) {
                 // The element is in the right partition
                 begin = to;
             } else if (k < from) {
                 // The element is in the left partition
                 end = from;
             } else {
                 // The range contains the element we wanted
                 return end;
             }
         }
         insertionSort(data, begin, end, begin != 0);
         fixSignedZeros(data, begin, end);
         // Either end == data.length and k+1 is sorted; or
         // end == pivot and k+1 is sorted
         return begin;
     }
 
     /**
      * Partition the array such that indices {@code k} correspond to their correctly
      * sorted value in the equivalent fully sorted array. For all indices {@code k}
      * and any index {@code i}:
      *
      * <pre>{@code
      * data[i < k] <= data[k] <= data[k < i]
      * }</pre>
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * @param data Values.
      * @param k Indices.
      */
     void partitionDP(double[] data, int... k) {
         final int n = k.length;
         if (n < 1) {
             return;
         }
 
         // Handle NaN
         final int length = sortNaN(data);
         if (length < 1) {
             return;
         }
 
         if (n == 1) {
             partitionDP(data, 0, length, (BitSet) null, k[0]);
             return;
         }
         // Special case for partition around adjacent indices (for interpolation)
         if (n == 2 && k[0] + 1 == k[1]) {
             final int p = partitionDP(data, 0, length, (BitSet) null, k[0]);
             if (p > k[1]) {
                 partitionMin(data, k[1], p);
             }
             return;
         }
 
         // Multiple pivots
         final BitSet pivots = new BitSet(length);
 
         for (int i = 0; i < n; i++) {
             final int kk = k[i];
             if (kk >= length || pivots.get(kk)) {
                 // Already sorted
                 continue;
             }
             int begin;
             int end;
             if (i == 0) {
                 begin = 0;
                 end = length;
             } else {
                 // Start inclusive
                 begin = pivots.previousSetBit(kk) + 1;
                 end = pivots.nextSetBit(kk + 1);
                 if (end < 0) {
                     end = length;
                 }
             }
             partitionDP(data, begin, end, pivots, kk);
         }
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * <p>Returns the pivot location adjacent to K to signal if K+1 is sorted.
      * If {@code p <= k} the range [p, min{k+2, data.length}) is sorted.
      * If {@code p > k} then p is a pivot.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivots Cache of pivot points.
      * @param k Index whose value in the array is of interest.
      * @return the bound index
      */
     private int partitionDP(double[] data, int begin, int end, BitSet pivots, int k) {
         // Find the unsorted range containing k
         final int[] bounds = new int[4];
         int div = 3;
         while (end - begin > minSelectSize) {
             div = partitionDP(data, begin, end, bounds, div);
             final int k0 = bounds[0];
             final int k1 = bounds[1];
             final int k2 = bounds[2];
             final int k3 = bounds[3];
             // sorted in [k0, k1) and (k2, k3]
             if (pivots != null) {
                 setPivots(k0, k1, pivots);
                 setPivots(k2 + 1, k3 + 1, pivots);
             }
             if (k > k3) {
                 // The element is in the right partition
                 begin = k3 + 1;
             } else if (k < k0) {
                 // The element is in the left partition
                 end = k0;
             } else if (k >= k1 && k <= k2) {
                 // Internal unsorted region
                 begin = k1;
                 end = k2 + 1;
             } else {
                 // The sorted ranges contain the element we wanted.
                 // Return a pivot (k0; k2+1) to signal if k+1 is sorted.
                 if (k + 1 < k1) {
                     return k0;
                 }
                 if (k + 1 < k3) {
                     return k2 + 1;
                 }
                 return end;
             }
         }
         if (pivots != null) {
             setPivots(begin, end, pivots);
         }
         insertionSort(data, begin, end, begin != 0);
         fixSignedZeros(data, begin, end);
         // Either end == data.length and k+1 is sorted; or
         // end == pivot and k+1 is sorted
         return begin;
     }
 
     /**
      * 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>Note: Requires that the range contains no NaN values. If the range contains
      * signed zeros and one is chosen as a pivot point the sort order of zeros is correct.
      *
      * <p>This method returns 4 points: the lower and upper pivots and bounds for
      * the internal range of unsorted values.
      * <ul>
      * <li>k0: lower pivot point: {@code a[k] < a[k0]} for {@code k < k0}.
      * <li>k1: the start (inclusive) of the unsorted range: {@code k0 < k1}.
      * <li>k2: the end (inclusive) of the unsorted range: {@code k2 <= k3}.
      * <li>k3: upper pivot point: {@code a[k3] < a[k]} for {@code k3 < k}.
      * </ul>
      *
      * <p>Bounds are set so {@code [k0, k1)} and {@code (k2, k3]} are fully sorted.
      * When the range {@code [k0, k3]} contains fully sorted elements
      * the result is set to {@code k1 = k3+1} and {@code k2 = k3}.
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * @param a Data array.
      * @param left Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param bounds Points [k0, k1, k2, k3].
      * @param div Divisor for the range to pick medians.
      * @return div
      */
     private int partitionDP(double[] a, int left, int end, int[] bounds, int div) {
         // Dual-pivot quicksort method by Vladimir Yaroslavskiy.
         //
         // 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 sentinal.
         // k traverses the unknown region ??? and values moved if less (l) or greater (g):
         //
         // left        l                k           g         right
         // |P1|  <P1   |   P1<= & <= P2 |    ???    |    >P2   |P2|
         //
         // At the end pivots are swapped back to behind the l and g pointers.
         //
         // |  <P1        |P1|     P1<= & <= P2    |P2|      >P2   |
         //
         // Adapted from Yaroslavskiy
         // http://codeblab.com/wp-content/uploads/2009/09/DualPivotQuicksort.pdf
         //
         // Modified to allow partial sorting (partitioning):
         // - Ignore insertion sort for tiny array (handled by calling code)
         // - Ignore recursive calls for a full sort (handled by calling code)
         // - Check for equal elements if a pivot is a signed zero
         // - Fix signed zeros within the region between pivots
         // - Change to fast-forward over initial ascending / descending runs
         // - Change to a single-pivot partition method if the pivots are equal
 
         final int right = end - 1;
         final int len = right - left;
 
         // Find pivots:
 
         // Original method: Guess medians using 1/3 and 2/3 of range
         final int third = len / div;
         int m1 = left + third;
         int m2 = right - third;
         if (m1 <= left) {
             m1 = left + 1;
         }
         if (m2 >= right) {
             m2 = right - 1;
         }
         if (a[m1] < a[m2]) {
             swap(a, m1, left);
             swap(a, m2, right);
         } else {
             swap(a, m1, right);
             swap(a, m2, left);
         }
         // pivots
         final double pivot1 = a[left];
         final double pivot2 = a[right];
 
         // Single pivot sort
         if (pivot1 == pivot2) {
             final int lower = partitionSBM(a, left, end, m1, bounds);
             final int upper = bounds[0];
             // Set dual pivot range
             bounds[0] = lower;
             bounds[3] = upper - 1;
             // Fully sorted internally
             bounds[1] = upper;
             bounds[2] = upper - 1;
             return div;
         }
 
         // Special case: Handle signed zeros
         int c = 0;
         if (pivot1 == 0 || pivot2 == 0) {
             c = countSignedZeros(a, left, end);
         }
 
         // pointers
         int less = left + 1;
         int great = right - 1;
 
         // Fast-forward ascending / descending runs to reduce swaps
         while (a[less] < pivot1) {
             less++;
         }
         while (a[great] > pivot2) {
             great--;
         }
 
         // sorting
         SORTING:
         for (int k = less; k <= great; k++) {
             final double v = a[k];
             if (v < pivot1) {
                 //swap(a, k, less++)
                 a[k] = a[less];
                 a[less] = v;
                 less++;
             } else if (v > pivot2) {
                 // Original
                 //while (k < great && a[great] > pivot2) {
                 //    great--;
                 //}
                 while (a[great] > pivot2) {
                     if (great-- == k) {
                         // Done
                         break SORTING;
                     }
                 }
                 // swap(a, k, great--)
                 // if (a[k] < pivot1)
                 //    swap(a, k, less++)
                 final double w = a[great];
                 a[great] = v;
                 great--;
                 // a[k] = w
                 if (w < pivot1) {
                     a[k] = a[less];
                     a[less] = w;
                     less++;
                 } else {
                     a[k] = w;
                 }
             }
         }
         // swaps
         final int dist = great - less;
         // Original paper: If middle partition (dist) is less than 13
         // then increase 'div' by 1. This means that the two outer partitions
         // contained most of the data and choosing medians should take
         // values closer to the edge. The middle will be sorted by quicksort.
         // 13 = 27 / 2 where 27 is the threshold for quicksort.
         if (dist < (minSelectSize >>> 1)) {
             div++;
         }
         //swap(a, less - 1, left)
         //swap(a, great + 1, right)
         a[left] = a[less - 1];
         a[less - 1] = pivot1;
         a[right] = a[great + 1];
         a[great + 1] = pivot2;
 
         // unsorted in [less, great]
 
         // Set the pivots
         bounds[0] = less - 1;
         bounds[3] = great + 1;
         //partitionDP(a, left, less - 2, div)
         //partitionDP(a, great + 2, right, div)
 
         // equal elements
         // Original paper: If middle partition (dist) is bigger
         // than (length - 13) then check for equal elements, i.e.
         // if the middle was very large there may be many repeated elements.
         // 13 = 27 / 2 where 27 is the threshold for quicksort.
         // We always do this if the pivots are signed zeros.
         if ((dist > len - (minSelectSize >>> 1) || c != 0) && pivot1 != pivot2) {
             // Fast-forward to reduce swaps
             while (a[less] == pivot1) {
                 less++;
             }
             while (a[great] == pivot2) {
                 great--;
             }
             // This copies the logic in the sorting loop using == comparisons
             EQUAL:
             for (int k = less; k <= great; k++) {
                 final double v = a[k];
                 if (v == pivot1) {
                     //swap(a, k, less++)
                     a[k] = a[less];
                     a[less] = v;
                     less++;
                 } else if (v == pivot2) {
                     while (a[great] == pivot2) {
                         if (great-- == k) {
                             // Done
                             break EQUAL;
                         }
                     }
                     final double w = a[great];
                     a[great] = v;
                     great--;
                     if (w == pivot1) {
                         a[k] = a[less];
                         a[less] = w;
                         less++;
                     } else {
                         a[k] = w;
                     }
                 }
             }
         }
         // unsorted in [less, great]
         if (pivot1 < pivot2 && less < great) {
             //partitionDP(a, less, great, div)
             bounds[1] = less;
             bounds[2] = great;
         } else {
             // Fully sorted
             bounds[1] = bounds[3] + 1;
             bounds[2] = bounds[3];
         }
 
         // Fix signed zeros
         if (c != 0) {
             int i;
             if (pivot1 == 0) {
                 i = bounds[0];
                 while (c-- > 0) {
                     a[i++] = -0.0;
                 }
                 while (i < end && a[i] == 0) {
                     a[i++] = 0.0;
                 }
             } else {
                 i = bounds[3];
                 while (a[i] == 0) {
                     a[i--] = 0.0;
                     if (i == left) {
                         break;
                     }
                 }
                 while (c-- > 0) {
                     a[++i] = -0.0;
                 }
             }
         }
 
         return div;
     }
 
     /**
      * Partition the array such that indices {@code k} correspond to their correctly
      * sorted value in the equivalent fully sorted array. For all indices {@code k}
      * and any index {@code i}:
      *
      * <pre>{@code
      * data[i < k] <= data[k] <= data[k < i]
      * }</pre>
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * @param data Values.
      * @param k Indices.
      */
     void partitionDP5(double[] data, int... k) {
         final int n = k.length;
         if (n < 1) {
             return;
         }
 
         // Handle NaN
         final int length = sortNaN(data);
         if (length < 1) {
             return;
         }
 
         if (n == 1) {
             partitionDP5(data, 0, length, (BitSet) null, k[0]);
             return;
         }
         // Special case for partition around adjacent indices (for interpolation)
         if (n == 2 && k[0] + 1 == k[1]) {
             final int p = partitionDP5(data, 0, length, (BitSet) null, k[0]);
             if (p > k[1]) {
                 partitionMin(data, k[1], p);
             }
             return;
         }
 
         // Multiple pivots
         final BitSet pivots = new BitSet(length);
 
         for (int i = 0; i < n; i++) {
             final int kk = k[i];
             if (kk >= length || pivots.get(kk)) {
                 // Already sorted
                 continue;
             }
             int begin;
             int end;
             if (i == 0) {
                 begin = 0;
                 end = length;
             } else {
                 // Start inclusive
                 begin = pivots.previousSetBit(kk) + 1;
                 end = pivots.nextSetBit(kk + 1);
                 if (end < 0) {
                     end = length;
                 }
             }
             partitionDP5(data, begin, end, pivots, kk);
         }
     }
 
     /**
      * Partition around the K<sup>th</sup> value in the array.
      *
      * <p>This method can be used for repeat calls to identify K<sup>th</sup> values in
      * the same array by caching locations of pivots (correctly sorted indices).
      *
      * <p>Note: Requires that the range contains no NaN values.
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * <p>Returns the pivot location adjacent to K to signal if K+1 is sorted.
      * If {@code p <= k} the range [p, min{k+2, data.length}) is sorted.
      * If {@code p > k} then p is a pivot.
      *
      * @param data Data array to use to find out the K<sup>th</sup> value.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param pivots Cache of pivot points.
      * @param k Index whose value in the array is of interest.
      * @return the bound index
      */
     private int partitionDP5(double[] data, int begin, int end, BitSet pivots, int k) {
         // Find the unsorted range containing k
         final int[] bounds = new int[4];
         while (end - begin > minSelectSize) {
             partitionDP5(data, begin, end, bounds);
             final int k0 = bounds[0];
             final int k1 = bounds[1];
             final int k2 = bounds[2];
             final int k3 = bounds[3];
             // sorted in [k0, k1) and (k2, k3]
             if (pivots != null) {
                 setPivots(k0, k1, pivots);
                 setPivots(k2 + 1, k3 + 1, pivots);
             }
             if (k > k3) {
                 // The element is in the right partition
                 begin = k3 + 1;
             } else if (k < k0) {
                 // The element is in the left partition
                 end = k0;
             } else if (k >= k1 && k <= k2) {
                 // Internal unsorted region
                 begin = k1;
                 end = k2 + 1;
             } else {
                 // The sorted ranges contain the element we wanted.
                 // Return a pivot (k0; k2+1) to signal if k+1 is sorted.
                 if (k + 1 < k1) {
                     return k0;
                 }
                 if (k + 1 < k3) {
                     return k2 + 1;
                 }
                 return end;
             }
         }
         if (pivots != null) {
             setPivots(begin, end, pivots);
         }
         insertionSort(data, begin, end, begin != 0);
         fixSignedZeros(data, begin, end);
         // Either end == data.length and k+1 is sorted; or
         // end == pivot and k+1 is sorted
         return begin;
     }
 
     /**
      * 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>Note: Requires that the range contains no NaN values. If the range contains
      * signed zeros and one is chosen as a pivot point the sort order of zeros is correct.
      *
      * <p>This method returns 4 points: the lower and upper pivots and bounds for
      * the internal range of unsorted values.
      * <ul>
      * <li>k0: lower pivot point: {@code a[k] < a[k0]} for {@code k < k0}.
      * <li>k1: the start (inclusive) of the unsorted range: {@code k0 < k1}.
      * <li>k2: the end (inclusive) of the unsorted range: {@code k2 <= k3}.
      * <li>k3: upper pivot point: {@code a[k3] < a[k]} for {@code k3 < k}.
      * </ul>
      *
      * <p>Bounds are set so {@code [k0, k1)} and {@code (k2, k3]} are fully sorted.
      * When the range {@code [k0, k3]} contains fully sorted elements
      * the result is set to {@code k1 = k3+1} and {@code k2 = k3}.
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * @param a Data array.
      * @param left Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param bounds Points [k0, k1, k2, k3].
      */
     private static void partitionDP5(double[] a, int left, int end, int[] bounds) {
         // Dual-pivot quicksort method by Vladimir Yaroslavskiy.
         //
         // Adapted from:
         //
         // http://codeblab.com/wp-content/uploads/2009/09/DualPivotQuicksort.pdf
         //
         // Modified to allow partial sorting (partitioning):
         // - Choose a pivot using 5 sorted points from the range.
         // - Ignore insertion sort for tiny array (handled by calling code)
         // - Ignore recursive calls for a full sort (handled by calling code)
         // - Check for equal elements if a pivot is a signed zero
         // - Fix signed zeros within the region between pivots
         // - Change to fast-forward over initial ascending / descending runs
         // - Change to a single-pivot partition method if the pivots are equal
 
         final int right = end - 1;
         final int len = right - left;
 
         // Find pivots:
 
         // Original method: Guess medians using 1/3 and 2/3 of range.
         // Here we sort 5 points and choose 2 and 4 as the pivots: 1/6, 1/3, 1/2, 2/3, 5/6
         // 1/6 ~ 1/8 + 1/32. Ensure the value is above zero to choose different points!
         // This is safe is len >= 4.
         final int sixth = 1 + (len >>> 3) + (len >>> 5);
         final int p3 = left + (len >>> 1);
         final int p2 = p3 - sixth;
         final int p1 = p2 - sixth;
         final int p4 = p3 + sixth;
         final int p5 = p4 + sixth;
         Sorting.sort5(a, p1, p2, p3, p4, p5);
 
         // For testing
         //p2 = DualPivotingStrategy.SORT_5.pivotIndex(a, left, end - 1, bounds);
         //p4 = bounds[0];
 
         final double pivot1 = a[p2];
         final double pivot2 = a[p4];
 
         // Add property to control this switch so we can benchmark not using it.
 
         if (pivot1 == pivot2) {
             // pivots == median !
             // Switch to a single pivot sort around the estimated median
             final int lower = partitionSBM(a, left, end, p3, bounds);
             final int upper = bounds[0];
             // Set dual pivot range
             bounds[0] = lower;
             bounds[3] = upper - 1;
             // No unsorted internal region
             bounds[1] = upper;
             bounds[2] = upper - 1;
             return;
         }
 
         // Special case: Handle signed zeros
         int c = 0;
         if (pivot1 == 0 || pivot2 == 0) {
             c = countSignedZeros(a, left, end);
         }
 
         // Move ends to the pivot locations.
         // After sorting the final pivot locations are overwritten.
         a[p2] = a[left];
         a[p4] = a[right];
         // It is assumed
         //a[left] = pivot1
         //a[right] = pivot2
 
         // pointers
         int less = left + 1;
         int great = right - 1;
 
         // Fast-forward ascending / descending runs to reduce swaps
         while (a[less] < pivot1) {
             less++;
         }
         while (a[great] > pivot2) {
             great--;
         }
 
         // sorting
         SORTING:
         for (int k = less; k <= great; k++) {
             final double v = a[k];
             if (v < pivot1) {
                 //swap(a, k, less++)
                 a[k] = a[less];
                 a[less] = v;
                 less++;
             } else if (v > pivot2) {
                 // Original
                 //while (k < great && a[great] > pivot2) {
                 //    great--;
                 //}
                 while (a[great] > pivot2) {
                     if (great-- == k) {
                         // Done
                         break SORTING;
                     }
                 }
                 // swap(a, k, great--)
                 // if (a[k] < pivot1)
                 //    swap(a, k, less++)
                 final double w = a[great];
                 a[great] = v;
                 great--;
                 // a[k] = w
                 if (w < pivot1) {
                     a[k] = a[less];
                     a[less] = w;
                     less++;
                 } else {
                     a[k] = w;
                 }
             }
         }
         // swaps
         //swap(a, less - 1, left)
         //swap(a, great + 1, right)
         a[left] = a[less - 1];
         a[less - 1] = pivot1;
         a[right] = a[great + 1];
         a[great + 1] = pivot2;
 
         // unsorted in [less, great]
 
         // Set the pivots
         bounds[0] = less - 1;
         bounds[3] = great + 1;
         //partitionDP5(a, left, less - 2)
         //partitionDP5(a, great + 2, right)
 
         // equal elements
         // Original paper: If middle partition (dist) is bigger
         // than (length - 13) then check for equal elements, i.e.
         // if the middle was very large there may be many repeated elements.
         // 13 = 27 / 2 where 27 is the threshold for quicksort.
 
         // Look for equal elements if the centre is more than 2/3 the length
         // We always do this if the pivots are signed zeros.
         if ((less < p1 && great > p5 || c != 0) && pivot1 != pivot2) {
 
             // Fast-forward to reduce swaps
             while (a[less] == pivot1) {
                 less++;
             }
             while (a[great] == pivot2) {
                 great--;
             }
 
             // This copies the logic in the sorting loop using == comparisons
             EQUAL:
             for (int k = less; k <= great; k++) {
                 final double v = a[k];
                 if (v == pivot1) {
                     //swap(a, k, less++)
                     a[k] = a[less];
                     a[less] = v;
                     less++;
                 } else if (v == pivot2) {
                     while (a[great] == pivot2) {
                         if (great-- == k) {
                             // Done
                             break EQUAL;
                         }
                     }
                     final double w = a[great];
                     a[great] = v;
                     great--;
                     if (w == pivot1) {
                         a[k] = a[less];
                         a[less] = w;
                         less++;
                     } else {
                         a[k] = w;
                     }
                 }
             }
         }
         // unsorted in [less, great]
         if (pivot1 < pivot2 && less < great) {
             //partitionDP5(a, less, great)
             bounds[1] = less;
             bounds[2] = great;
         } else {
             // Fully sorted
             bounds[1] = bounds[3] + 1;
             bounds[2] = bounds[3];
         }
 
         // Fix signed zeros
         if (c != 0) {
             int i;
             if (pivot1 == 0) {
                 i = bounds[0];
                 while (c-- > 0) {
                     a[i++] = -0.0;
                 }
                 while (i < end && a[i] == 0) {
                     a[i++] = 0.0;
                 }
             } else {
                 i = bounds[3];
                 while (a[i] == 0) {
                     a[i--] = 0.0;
                     if (i == left) {
                         break;
                     }
                 }
                 while (c-- > 0) {
                     a[++i] = -0.0;
                 }
             }
         }
     }
 
     /**
      * Swaps the two specified elements in the array.
      *
      * @param array Array.
      * @param i First index.
      * @param j Second index.
      */
     private static void swap(double[] array, int i, int j) {
         final double tmp = array[i];
         array[i] = array[j];
         array[j] = tmp;
     }
 
     /**
      * Sort the data.
      *
      * <p>Uses a single-pivot quicksort partition method.
      *
      * @param data Values.
      */
     void sortSP(double[] data) {
         sortSP(data, 0, data.length);
     }
 
     /**
      * Sort the data.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private void sortSP(double[] data, int begin, int end) {
         if (end - begin <= 1) {
             return;
         }
         final int i = partitionSP(data, begin, end,
             pivotingStrategy.pivotIndex(data, begin, end - 1, begin));
         sortSP(data, begin, i);
         sortSP(data, i + 1, end);
     }
 
     /**
      * Sort the data.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method.
      *
      * @param data Values.
      */
     void sortSBM(double[] data) {
         sortSBM(data, 0, sortNaN(data));
     }
 
     /**
      * Sort the data. Requires no NaN values in the range.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private void sortSBM(double[] data, int begin, int end) {
         if (end - begin <= minSelectSize) {
             insertionSort(data, begin, end, begin != 0);
             if (begin < end) {
                 fixSignedZeros(data, begin, end);
             }
             return;
         }
         final int[] to = {0};
         final int from = partitionSBM(data, begin, end,
             pivotingStrategy.pivotIndex(data, begin, end - 1, begin), to);
         sortSBM(data, begin, from);
         sortSBM(data, to[0], end);
     }
 
     /**
      * Sort the data.
      *
      * <p>Uses a Bentley-McIlroy quicksort partition method.
      *
      * @param data Values.
      */
     void sortBM(double[] data) {
         sortBM(data, 0, sortNaN(data));
     }
 
     /**
      * Sort the data. Requires no NaN values in the range.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private void sortBM(double[] data, int begin, int end) {
         if (end - begin <= minSelectSize) {
             insertionSort(data, begin, end, begin != 0);
             if (begin < end) {
                 fixSignedZeros(data, begin, end);
             }
             return;
         }
         final int[] to = {0};
         final int from = partitionBM(data, begin, end,
             pivotingStrategy.pivotIndex(data, begin, end - 1, begin), to);
         sortBM(data, begin, from);
         sortBM(data, to[0], end);
     }
 
     /**
      * Sort the data.
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy.
      *
      * @param data Values.
      */
     void sortDP(double[] data) {
         sortDP(data, 0, sortNaN(data), 3);
     }
 
     /**
      * Sort the data. Requires no NaN values in the range.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      * @param div Divisor for the range to pick medians.
      */
     private void sortDP(double[] data, int begin, int end, int div) {
         if (end - begin <= minSelectSize) {
             insertionSort(data, begin, end, begin != 0);
             if (begin < end) {
                 fixSignedZeros(data, begin, end);
             }
             return;
         }
         final int[] bounds = new int[4];
         div = partitionDP(data, begin, end, bounds, div);
         final int k0 = bounds[0];
         final int k1 = bounds[1];
         final int k2 = bounds[2];
         final int k3 = bounds[3];
         sortDP(data, begin, k0, div);
         sortDP(data, k3 + 1, end, div);
         sortDP(data, k1, k2 + 1, div);
     }
 
     /**
      * Sort the data.
      *
      * <p>Uses a dual-pivot quicksort method by Vladimir Yaroslavskiy optimised
      * to choose the pivots using 5 sorted points.
      *
      * @param data Values.
      */
     void sortDP5(double[] data) {
         sortDP5(data, 0, sortNaN(data));
     }
 
     /**
      * Sort the data. Requires no NaN values in the range.
      *
      * @param data Values.
      * @param begin Lower bound (inclusive).
      * @param end Upper bound (exclusive).
      */
     private void sortDP5(double[] data, int begin, int end) {
         if (end - begin <= minSelectSize) {
             insertionSort(data, begin, end, begin != 0);
             if (begin < end) {
                 fixSignedZeros(data, begin, end);
             }
             return;
         }
         final int[] bounds = new int[4];
         partitionDP5(data, begin, end, bounds);
         final int k0 = bounds[0];
         final int k1 = bounds[1];
         final int k2 = bounds[2];
         final int k3 = bounds[3];
         sortDP5(data, begin, k0);
         sortDP5(data, k3 + 1, end);
         sortDP5(data, k1, k2 + 1);
     }
 
     /**
      * Creates the pivots heap for a data array of the specified {@code length}.
      * If the array is too small to use the pivots heap then an empty array is returned.
      *
      * @param length Length.
      * @return the pivots heap
      */
     static int[] createPivotsHeap(int length) {
         if (length <= MIN_SELECT_SIZE) {
             return NO_PIVOTS;
         }
         // Size should be x^2 - 1, where x is the layers in the heap.
         // Do not create more pivots than the array length. When partitions are small
         // the pivots are no longer used so this does not have to contain all indices.
         // Default size in Commons Math Percentile class was 1023 (10 layers).
         final int n = nextPow2(length >>> 1);
         final int[] pivotsHeap = new int[Math.min(n, 1 << 10) - 1];
         Arrays.fill(pivotsHeap, -1);
         return pivotsHeap;
     }
 
     /**
      * Returns the closest power-of-two number greater than or equal to value.
      *
      * <p>Warning: This will return {@link Integer#MIN_VALUE} for any value above
      * {@code 1 << 30}. This is the next power of 2 as an unsigned integer.
      *
      * @param value the value (must be positive)
      * @return the closest power-of-two number greater than or equal to value
      */
     private static int nextPow2(int value) {
         // shift by -x is equal to shift by (32 - x) as only the low 5-bits are used.
         return 1 << -Integer.numberOfLeadingZeros(value - 1);
     }
 }