double a;
        int b;
        final int[] istack = new int[NSTACK];
        ir = n - 1;

        for (; ; ) {
            // Insertion sort when subarray is small enough
            if (ir - l < M) {
                for (j = l + 1; j <= ir; j++) {
                    a = array[j + fromIndex];
                    b = indices[j + fromIndex];
                    for (i = j - 1; i >= l; i--) {
                        if (array[i + fromIndex] <= a) {
                            break;
                        }
                        array[i + 1 + fromIndex] = array[i + fromIndex];
                        indices[i + 1 + fromIndex] = indices[i + fromIndex];
                    }
                    array[i + 1 + fromIndex] = a;
                    indices[i + 1 + fromIndex] = b;
                }
                if (jstack < 0) {
                    break;
                }
                // Pop stack and begin a new round of partitioning
                ir = istack[jstack--];
                l = istack[jstack--];
            } else {
                // Choose median of left, center, and right elements as
                // partitioning element "a". Also rearrange so that a(l) <= a(l+1)
                // <= a(ir)
                k = (l + ir) >> 1;
                swap(array, k + fromIndex, l + 1 + fromIndex);
                swapIndices(indices, k + fromIndex, l + 1 + fromIndex);
                if (array[l + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + fromIndex, ir + fromIndex);
                    swapIndices(indices, l + fromIndex, ir + fromIndex);
                }
                if (array[l + 1 + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + 1 + fromIndex, ir + fromIndex);
                    swapIndices(indices, l + 1 + fromIndex, ir + fromIndex);
                }
                if (array[l + fromIndex] > array[l + 1 + fromIndex]) {
                    swap(array, l + fromIndex, l + 1 + fromIndex);
                    swapIndices(indices, l + fromIndex, l + 1 + fromIndex);
                }
                // Initialize pointers for partitioning
                i = l + 1;
                j = ir;
                // Partitioning element
                a = array[l + 1 + fromIndex];
                b = indices[l + 1 + fromIndex];
                // Beginning of innermost loop
                for (; ; ) {
                    // Scan up to find element > a
                    do {
                        i++;
                    } while (array[i + fromIndex] < a);
                    // Scan down to find element < a
                    do {
                        j--;
                    } while (array[j + fromIndex] > a);
                    // Pointers crossed. Partitioning complete
                    if (j < i) {
                        break;
                    }
                    // Exchange elements
                    swap(array, i + fromIndex, j + fromIndex);
                    swapIndices(indices, i + fromIndex, j + fromIndex);
                    // End of innermost loop
                }
                // Insert partitioning element
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;
                indices[l + 1 + fromIndex] = indices[j + fromIndex];
                indices[j + fromIndex] = b;
                jstack += 2;
                // NSTACK too small in sort
                if (jstack >= NSTACK) {
                    throw new SortingException();
                }
                // Push pointers to larger subarray on stack; process smaller
                // subarray immediately
                if (ir - i + 1 >= j - l) {
                    istack[jstack] = ir;
                    istack[jstack - 1] = i;
                    ir = j - 1;
                } else {
                    istack[jstack] = j - 1;
                    istack[jstack - 1] = l;
                    l = i;
                }
            }
        }

        return indices;
    }

double a;
        final int[] istack = new int[NSTACK];
        ir = n - 1;

        for (; ; ) {
            // Insertion sort when subarray is small enough
            if (ir - l < M) {
                for (j = l + 1; j <= ir; j++) {
                    a = array[j + fromIndex];
                    for (i = j - 1; i >= l; i--) {
                        if (array[i + fromIndex] <= a) {
                            break;
                        }
                        array[i + 1 + fromIndex] = array[i + fromIndex];
                    }
                    array[i + 1 + fromIndex] = a;
                }
                if (jstack < 0) {
                    break;
                }
                // Pop stack and begin a new round of partitioning
                ir = istack[jstack--];
                l = istack[jstack--];
            } else {
                // Choose median of left, center, and right elements as
                // partitioning element "a". Also rearrange so that a(l) <= a(l+1)
                // <= a(ir)
                k = (l + ir) >> 1;
                swap(array, k + fromIndex, l + 1 + fromIndex);
                if (array[l + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + fromIndex, ir + fromIndex);
                }
                if (array[l + 1 + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + 1 + fromIndex, ir + fromIndex);
                }
                if (array[l + fromIndex] > array[l + 1 + fromIndex]) {
                    swap(array, l + fromIndex, l + 1 + fromIndex);
                }
                // Initialize pointers for partitioning
                i = l + 1;
                j = ir;
                // Partitioning element
                a = array[l + 1 + fromIndex];
                // Beginning of innermost loop
                for (; ; ) {
                    // Scan up to find element > a
                    do {
                        i++;
                    } while (array[i + fromIndex] < a);
                    // Scan down to find element < a
                    do {
                        j--;
                    } while (array[j + fromIndex] > a);
                    // Pointers crossed. Partitioning complete
                    if (j < i) {
                        break;
                    }
                    // Exchange elements
                    swap(array, i + fromIndex, j + fromIndex);
                    // End of innermost loop
                }
                // Insert partitioning element
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;
                jstack += 2;
                // NSTACK too small in sort
                if (jstack >= NSTACK) {
                    throw new SortingException();
                }
                // Push pointers to larger subarray on stack; process smaller
                // subarray immediately
                if (ir - i + 1 >= j - l) {
                    istack[jstack] = ir;
                    istack[jstack - 1] = i;
                    ir = j - 1;
                } else {
                    istack[jstack] = j - 1;
                    istack[jstack - 1] = l;
                    l = i;
                }
            }
        }
    }

    /**
     * Sorts provided array in ascending order so that {@code
     * array[i - 1] < array[i]} for any valid i.
     * This method modifies provided array so that
     * after execution of this method array elements are ordered.
     * An array containing the original indices where elements were
     * located is returned so that other arrays or collections can be kept
     * in the same order.
     *
     * @param array     Array to be sorted. After execution of this method
     *                  elements in array between fromIndex (inclusive) and toIndex
     *                  (exclusive) are modified so that they are on ascending order.
     * @param fromIndex Index were sorting starts (inclusive).
     * @param toIndex   Index were sorting stops (exclusive).
     * @return Array containing original location of elements that have been
     * sorted. Only elements between fromIndex (inclusive) and toIndex
     * (exclusive) are modified, the remaining ones are kept in natural
     * order.
     * @throws SortingException               If for some reason sorting fails.
     * @throws IllegalArgumentException       If {@code fromIndex > toIndex}.
     * @throws ArrayIndexOutOfBoundsException if {@code fromIndex < 0} or
     *                                        {@code toIndex > array.length}.
     */
    @Override
    public int[] sortWithIndices(final double[] array, final int fromIndex, final int toIndex)

double a;
        l = 0;
        ir = n - 1;
        for (; ; ) {
            if (ir <= l + 1) {
                if (ir == l + 1 && array[ir + fromIndex] < array[l + fromIndex]) {
                    swap(array, l + fromIndex, ir + fromIndex);
                }
                return array[k + fromIndex];
            } else {
                mid = (l + ir) >> 1;
                swap(array, mid + fromIndex, l + 1 + fromIndex);
                if (array[l + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + fromIndex, ir + fromIndex);
                }
                if (array[l + 1 + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + 1 + fromIndex, ir + fromIndex);
                }
                if (array[l + fromIndex] > array[l + 1 + fromIndex]) {
                    swap(array, l + fromIndex, l + 1 + fromIndex);
                }
                i = l + 1;
                j = ir;
                a = array[l + 1 + fromIndex];
                for (; ; ) {
                    do {
                        i++;
                    } while (array[i + fromIndex] < a);

                    do {
                        j--;
                    } while (array[j + fromIndex] > a);

                    if (j < i) {
                        break;
                    }

                    swap(array, i + fromIndex, j + fromIndex);
                }
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;
                if (j >= k) {
                    ir = j - 1;
                }
                if (j <= k) {
                    l = i;
                }
            }
        }
    }

    /**
     * Returns the k-th sorted element in provided array of {@link Comparable}
     * starting at fromIndex and finishing at toIndex, elements outside this
     * range are ignored.
     * Selecting an element is usually faster than sorting the whole
     * array, and for that reason, when only a few sorted elements are
     * required, this method should be used instead of sort.
     * Because array is passed by reference, after executing this method
     * array is modified so that in k location it contains the k-th sorted
     * elements and on locations array[fromIndex] ... array[k-1 + fromIndex]
     * contains unsorted elements smaller than sorted element
     * array[k + fromIndex],  and on locations
     * array[k+1 + fromIndex] ... array[toIndex-1] contains unsorted
     * elements greater than sorted element array[k + fromIndex].
     *
     * @param k         Position of sorted element to be retrieved.
     * @param array     Array to be used for retrieving k-th sorted element.
     *                  Provided array is passed by reference and modified upon execution of
     *                  this method so that k-th location contains k-th sorted element, and
     *                  array[fromIndex] ... array[k-1 + fromIndex] contains unsorted
     *                  elements smaller than sorted element array[k + fromIndex] and
     *                  array[k+1 + fromIndex] ... array[toIndex-1] contains elements
     *                  greater than sorted element array[k].
     * @param fromIndex Index were sorting starts (inclusive).
     * @param toIndex   Index were sorting stops (exclusive).
     * @return The k-th sorted element in provided array.
     * @throws IllegalArgumentException       if k &lt; (toIndex - fromIndex) or
     *                                        fromIndex &lt; toIndex.
     * @throws ArrayIndexOutOfBoundsException if fromIndex or toIndex are
     *                                        outside array boundaries.
     */
    public float select(final int k, final float[] array, final int fromIndex, final int toIndex) {

swap(array, k + fromIndex, l + 1 + fromIndex);
                if (comparator.compare(array[l + fromIndex], array[ir + fromIndex]) > 0) {
                    swap(array, l + fromIndex, ir + fromIndex);
                }
                if (comparator.compare(array[l + 1 + fromIndex], array[ir + fromIndex]) > 0) {
                    swap(array, l + 1 + fromIndex, ir + fromIndex);
                }
                if (comparator.compare(array[l + fromIndex], array[l + 1 + fromIndex]) > 0) {
                    swap(array, l + fromIndex, l + 1 + fromIndex);
                }
                // Initialize pointers for partitioning
                i = l + 1;
                j = ir;
                // Partitioning element
                a = array[l + 1 + fromIndex];
                // Beginning of innermost loop
                for (; ; ) {
                    // Scan up to find element > a
                    do {
                        i++;
                    } while (comparator.compare(array[i + fromIndex], a) < 0);
                    // Scan down to find element < a
                    do {
                        j--;
                    } while (comparator.compare(array[j + fromIndex], a) > 0);
                    // Pointers crossed. Partitioning complete
                    if (j < i) {
                        break;
                    }
                    // Exchange elements
                    swap(array, i + fromIndex, j + fromIndex);
                    // End of innermost loop
                }
                // Insert partitioning element
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;

swap(array, k + fromIndex, l + 1 + fromIndex);
                if (array[l + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + fromIndex, ir + fromIndex);
                }
                if (array[l + 1 + fromIndex] > array[ir + fromIndex]) {
                    swap(array, l + 1 + fromIndex, ir + fromIndex);
                }
                if (array[l + fromIndex] > array[l + 1 + fromIndex]) {
                    swap(array, l + fromIndex, l + 1 + fromIndex);
                }
                // Initialize pointers for partitioning
                i = l + 1;
                j = ir;
                // Partitioning element
                a = array[l + 1 + fromIndex];
                // Beginning of innermost loop
                for (; ; ) {
                    // Scan up to find element > a
                    do {
                        i++;
                    } while (array[i + fromIndex] < a);
                    // Scan down to find element < a
                    do {
                        j--;
                    } while (array[j + fromIndex] > a);
                    // Pointers crossed. Partitioning complete
                    if (j < i) {
                        break;
                    }
                    // Exchange elements
                    swap(array, i + fromIndex, j + fromIndex);
                    // End of innermost loop
                }
                // Insert partitioning element
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;

} while (comparator.compare(array[j + fromIndex], a) > 0);
                    // Pointers crossed. Partitioning complete
                    if (j < i) {
                        break;
                    }
                    // Exchange elements
                    swap(array, i + fromIndex, j + fromIndex);
                    swapIndices(indices, i + fromIndex, j + fromIndex);
                    // End of innermost loop
                }
                // Insert partitioning element
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;
                indices[l + 1 + fromIndex] = indices[j + fromIndex];
                indices[j + fromIndex] = b;
                jstack += 2;
                // NSTACK too small in sort
                if (jstack >= NSTACK) {
                    throw new SortingException();
                }
                // Push pointers to larger subarray on stack; process smaller
                // subarray immediately
                if (ir - i + 1 >= j - l) {
                    istack[jstack] = ir;
                    istack[jstack - 1] = i;
                    ir = j - 1;
                } else {
                    istack[jstack] = j - 1;
                    istack[jstack - 1] = l;
                    l = i;
                }
            }
        }

        return indices;
    }

    /**
     * Returns sorting method of this class.
     *
     * @return Sorting method.
     */
    @Override
    public SortingMethod getMethod() {

public int[] sortWithIndices(final double[] array, final int fromIndex, final int toIndex) {

        if (fromIndex > toIndex) {
            throw new IllegalArgumentException();
        }
        if (fromIndex < 0 || toIndex > array.length) {
            throw new ArrayIndexOutOfBoundsException();
        }

        final int[] indices = getInitialIndicesVector(array.length);
        if (fromIndex == toIndex) {
            return indices;
        }

        int i;
        final int n = toIndex - fromIndex;

        for (i = n / 2 - 1; i >= 0; i--) {
            siftDownWithIndices(array, indices, i, n - 1, fromIndex);
        }
        for (i = n - 1; i > 0; i--) {
            swap(array, fromIndex, i + fromIndex);
            swapIndices(indices, fromIndex, i + fromIndex);
            siftDownWithIndices(array, indices, 0, i - 1, fromIndex);
        }

        return indices;
    }

} while (comparator.compare(array[j + fromIndex], a) > 0);
                    // Pointers crossed. Partitioning complete
                    if (j < i) {
                        break;
                    }
                    // Exchange elements
                    swap(array, i + fromIndex, j + fromIndex);
                    swapIndices(indices, i + fromIndex, j + fromIndex);
                    // End of innermost loop
                }
                // Insert partitioning element
                array[l + 1 + fromIndex] = array[j + fromIndex];
                array[j + fromIndex] = a;
                indices[l + 1 + fromIndex] = indices[j + fromIndex];
                indices[j + fromIndex] = b;
                jstack += 2;
                // NSTACK too small in sort
                if (jstack >= NSTACK) {
                    throw new SortingException();
                }
                // Push pointers to larger subarray on stack; process smaller
                // subarray immediately
                if (ir - i + 1 >= j - l) {
                    istack[jstack] = ir;
                    istack[jstack - 1] = i;
                    ir = j - 1;
                } else {
                    istack[jstack] = j - 1;
                    istack[jstack - 1] = l;
                    l = i;
                }
            }
        }

        return indices;
    }

public void sort(final double[] array, final int fromIndex, final int toIndex) {

        if (fromIndex > toIndex) {
            throw new IllegalArgumentException();
        }
        if (fromIndex < 0 || toIndex > array.length) {
            throw new ArrayIndexOutOfBoundsException();
        }
        if (fromIndex == toIndex) {
            return;
        }

        int i;
        final int n = toIndex - fromIndex;

        for (i = n / 2 - 1; i >= 0; i--) {
            siftDown(array, i, n - 1, fromIndex);
        }
        for (i = n - 1; i > 0; i--) {
            swap(array, fromIndex, i + fromIndex);
            siftDown(array, 0, i - 1, fromIndex);
        }
    }

    /**
     * Sorts provided array in ascending order so that {@code
     * array[i - 1] < array[i]} for any valid i.
     * This method modifies provided array so that
     * after execution of this method array elements are ordered.
     * An array containing the original indices where elements were
     * located is returned so that other arrays or collections can be kept
     * in the same order.
     *
     * @param array     Array to be sorted. After execution of this method
     *                  elements in array between fromIndex (inclusive) and toIndex
     *                  (exclusive) are modified so that they are on ascending order.
     * @param fromIndex Index were sorting starts (inclusive).
     * @param toIndex   Index were sorting stops (exclusive).
     * @return Array containing original location of elements that have been
     * sorted. Only elements between fromIndex (inclusive) and toIndex
     * (exclusive) are modified, the remaining ones are kept in natural
     * order.
     * @throws IllegalArgumentException       If {@code fromIndex > toIndex}.
     * @throws ArrayIndexOutOfBoundsException if {@code fromIndex < 0} or
     *                                        {@code toIndex > array.length}.
     */
    @Override
    public int[] sortWithIndices(final double[] array, final int fromIndex, final int toIndex) {

double v;

        do {
            inc *= INCREMENT_FACTOR;
            inc++;
        } while (inc <= n);

        // Loop over the partial sorts
        do {
            inc /= INCREMENT_FACTOR;
            // Outer loop of straight insertion
            for (int i = inc; i < n; i++) {
                v = array[i + fromIndex];
                b = indices[i + fromIndex];
                j = i;

                // Inner loop of straight insertion
                while (array[j - inc + fromIndex] > v) {
                    array[j + fromIndex] = array[j - inc + fromIndex];
                    indices[j + fromIndex] = indices[j - inc + fromIndex];
                    j -= inc;
                    if (j < inc) {
                        break;
                    }
                }
                array[j + fromIndex] = v;
                indices[j + fromIndex] = b;
            }
        } while (inc > MIN_INCREMENT);

        return indices;
    }

final double a = ra[l + fromIndex];
        final int b = rb[l + fromIndex];
        jold = l;
        j = 2 * l + 1;
        while (j <= r) {
            if (j < r && ra[j + fromIndex] < ra[j + 1 + fromIndex]) {
                j++;
            }
            if (a >= ra[j + fromIndex]) {
                break;
            }
            ra[jold + fromIndex] = ra[j + fromIndex];
            rb[jold + fromIndex] = rb[j + fromIndex];
            jold = j;
            j = 2 * j + 1;
        }
        ra[jold + fromIndex] = a;
        rb[jold + fromIndex] = b;
    }

    /**
     * Internal method to reorder sub-array ra.
     *
     * @param ra        sub-array ra.
     * @param l         l value.
     * @param r         r value.
     * @param fromIndex initial position.
     */
    private void siftDown(final float[] ra, final int l, final int r, final int fromIndex) {

final double a = ra[l + fromIndex];
        final int b = rb[l + fromIndex];
        jold = l;
        j = 2 * l + 1;
        while (j <= r) {
            if (j < r && ra[j + fromIndex] < ra[j + 1 + fromIndex]) {
                j++;
            }
            if (a >= ra[j + fromIndex]) {
                break;
            }
            ra[jold + fromIndex] = ra[j + fromIndex];
            rb[jold + fromIndex] = rb[j + fromIndex];
            jold = j;
            j = 2 * j + 1;
        }
        ra[jold + fromIndex] = a;
        rb[jold + fromIndex] = b;
    }

double v;

        do {
            inc *= INCREMENT_FACTOR;
            inc++;
        } while (inc <= n);

        // Loop over the partial sorts
        do {
            inc /= INCREMENT_FACTOR;
            // Outer loop of straight insertion
            for (int i = inc; i < n; i++) {
                v = array[i + fromIndex];
                j = i;

                // Inner loop of straight insertion
                while (array[j - inc + fromIndex] > v) {
                    array[j + fromIndex] = array[j - inc + fromIndex];
                    j -= inc;
                    if (j < inc) {
                        break;
                    }
                }
                array[j + fromIndex] = v;
            }
        } while (inc > MIN_INCREMENT);
    }

    /**
     * Sorts provided array in ascending order so that {@code
     * array[i - 1] < array[i]} for any valid i.
     * This method modifies provided array so that
     * after execution of this method array elements are ordered.
     * An array containing the original indices where elements were
     * located is returned so that other arrays or collections can be kept
     * in the same order.
     *
     * @param array     Array to be sorted. After execution of this method
     *                  elements in array between fromIndex (inclusive) and toIndex
     *                  (exclusive) are modified so that they are on ascending order.
     * @param fromIndex Index were sorting starts (inclusive).
     * @param toIndex   Index were sorting stops (exclusive).
     * @return Array containing original location of elements that have been
     * sorted. Only elements between fromIndex (inclusive) and toIndex
     * (exclusive) are modified, the remaining ones are kept in natural
     * order.
     * @throws IllegalArgumentException       If {@code fromIndex > toIndex}.
     * @throws ArrayIndexOutOfBoundsException if {@code fromIndex < 0} or
     *                                        {@code toIndex > array.length}.
     */
    @Override
    public int[] sortWithIndices(final double[] array, final int fromIndex, final int toIndex) {

final double a = ra[l + fromIndex];
        jold = l;
        j = 2 * l + 1;
        while (j <= r) {
            if (j < r && ra[j + fromIndex] < ra[j + 1 + fromIndex]) {
                j++;
            }
            if (a >= ra[j + fromIndex]) {
                break;
            }
            ra[jold + fromIndex] = ra[j + fromIndex];
            jold = j;
            j = 2 * j + 1;
        }
        ra[jold + fromIndex] = a;
    }

    /**
     * Internal method to reorder sub-array ra along with its corresponding
     * indices.
     *
     * @param ra        sub-array ra.
     * @param rb        sub-array rb.
     * @param l         l value.
     * @param r         r value.
     * @param fromIndex initial position.
     */
    private void siftDownWithIndices(final double[] ra, final int[] rb, final int l, final int r,