Снятие рекурсии в Ярославском Dual Pivot Quicksort

Question

Я начал работать над модифицированной версией Yaroslavskiys Dual Pivot Quicksort и в качестве первого шага попытался удалить рекурсию (код основан на рекурсивной версии в этой статье). Возможно, это немного старомодно, но я все же хотел бы попробовать. Вот код C, который у меня есть. Он работает, как задумано, но меня немного беспокоят патологические случаи, когда локальный стек может переполняться:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define QS_STACK_SIZE     96
#define QS_ARR_SIZE       10000000

#define INT_SWAP(x_,y_)\
do {\
    int tmp_ = (x_);\
    (x_) = (y_);\
    (y_) = tmp_;\
} while(0)

void qsort_yaro(int *a, int left, int right, int div);

int main(void)
{
    int i;
    static int a[QS_ARR_SIZE];

    srand (time (NULL));

    for (i = 0; i < QS_ARR_SIZE; i ++)
        a[i] = rand();

    printf ("Sorting %i elements ...\n", QS_ARR_SIZE);
    qsort_yaro (a, 0, QS_ARR_SIZE - 1, 3);

    for (i = 0; i < QS_ARR_SIZE - 1; i ++)
        if (a[i] > a[i + 1])
            {
            printf ("ERROR : Sequence not ascending!\n");
            break;
            }

    if (i == QS_ARR_SIZE - 1)
        printf ("OK : Sequence ascending!\n");

    return 0;
}

void qsort_yaro(int *a, int left, int right, int div)
{
    int i, j, k, third;
    int len;
    int m1, m2;
    int pivot1, pivot2;
    int less, great;
    int dist;
    int hi[QS_STACK_SIZE];
    int lo[QS_STACK_SIZE];
    int lvl = -1;
    int max_stack = 0;

    lo[++lvl] = left;
    hi[lvl] = right;

    while (lvl >= 0)
        {
        left = lo[lvl];
        right = hi[lvl--];

        while ((len = right - left) >= 27)
            { // insertion sort for tiny array
            third = len / div;
            //
            // "medians"
            m1 = left + third;
            m2 = right - third;

            if (m1 <= left)
                m1 = left + 1;
            if (m2 >= right)
                m2 = right - 1;
            if (a[m1] < a[m2])
                {
                INT_SWAP(a[m1], a[left]);
                INT_SWAP(a[m2], a[right]);
                }
            else
                {
                INT_SWAP(a[m1], a[right]);
                INT_SWAP(a[m2], a[left]);
                }

            // pivots
            pivot1 = a[left];
            pivot2 = a[right];
            // pointers
            less = left + 1;
            great = right - 1;
            // sorting
            for (k = less; k <= great; k++)
                {
                if (a[k] < pivot1)
                    {
                    INT_SWAP(a[k], a[less]);
                    less++;
                    }
                else if (a[k] > pivot2)
                    {
                    while (k < great && a[great] > pivot2)
                        great--;
                    INT_SWAP(a[k], a[great]);
                    great--;
                    if (a[k] < pivot1)
                        {
                        INT_SWAP(a[k], a[less]);
                        less++;
                        }
                    }
                }
            // swaps
            dist = great - less;
            if (dist < 13)
                div++;

            INT_SWAP(a[less - 1], a[left]);
            INT_SWAP(a[great + 1], a[right]);

//            lo[++lvl] = left;
//            hi[lvl]   = less - 2;
//            lo[++lvl] = great + 2;
//            hi[lvl]   = right;

            // Push the smaller interval on top for less stack usage
            // subarrays
            if (less - 2 - left > right - great - 2)
                {
                lo[++lvl] = left;
                hi[lvl]   = less - 2;
                lo[++lvl] = great + 2;
                hi[lvl]   = right;
                }
            else
                {
                lo[++lvl] = great + 2;
                hi[lvl]   = right;
                lo[++lvl] = left;
                hi[lvl]   = less - 2;
                }
            if (max_stack < lvl)
                max_stack = lvl;
            // equal elements
            if (dist > len - 13 && pivot1 != pivot2)
                {
                for (k = less; k <= great; k++)
                    {
                    if (a[k] == pivot1)
                        {
                        INT_SWAP(a[k], a[less]);
                        less++;
                        }
                    else if (a[k] == pivot2)
                        {
                        INT_SWAP(a[k], a[great]);
                        great--;
                        if (a[k] == pivot1)
                            {
                            INT_SWAP(a[k], a[less]);
                            less++;
                            }
                        }
                    }
                }
            // subarray
            if (pivot1 < pivot2)
                {
                left  = less;
                right = great;
                }
            else
                {
                left  = lo[lvl];
                right = hi[lvl--];
                }
            }

        for (i = left + 1; i <= right; i++)
            for (j = i; j > left && a[j] < a[j - 1]; j--)
                INT_SWAP (a[j], a[j - 1]);
        }
    printf ("Max Stack Size : %i\n", max_stack);
}

Я подправил оригинальную реализацию этим небольшим изменением

        if (less - 2 - left > right - great - 2)
            {
            lo[++lvl] = left;
            hi[lvl]   = less - 2;
            lo[++lvl] = great + 2;
            hi[lvl]   = right;
            }
        else
            {
            lo[++lvl] = great + 2;
            hi[lvl]   = right;
            lo[++lvl] = left;
            hi[lvl]   = less - 2;
            }

которая основана на той же идее, что и в обычной итеративной реализации 1-pivot, где большие части помещаются в стек, чтобы ограничить размер стека до log(n). Здесь меньшая часть помещается сверху стека, чтобы вернуть его как можно раньше.

Я провел несколько экспериментов со случайно сгенерированными массивами из 10.000.000 элементов и мог реально сохранить до 5 уровней стека этим методом (максимум ~45).

Теперь я не нахожу строгого аргумента, почему размер стека должен быть ограничен константой в этой реализации, поэтому я либо пропустил здесь математику, либо мне нужна лучшая идея и реализация.