I'm submitting the median paper to http://2017.programmingconference.org/track/programming-2017-papers. (It was rejected by ALENEX... the short story is there was no open-source implementation and they simply didn't buy the results.)

Now I have a C++ implementation available. I wanted to ask for your help with lower/upper median of 4, which I coded like this (based on the median of 5 posted here by user "tn" aka Teppo Niinimäki):

template <bool leanRight, class T>
void medianOf(T* r, size_t a, size_t b, size_t c, size_t d)
{
    assert(a != b && a != c && a != d && b != c && b != d
        && c != d);
    if (leanRight)
    {
        // In the median of 5 algorithm, consider r[e] infinite
        if (r[c] < r[a]) {
            std::swap(r[a], r[c]);
        }
        if (r[d] < r[b]) {
            std::swap(r[b], r[d]);
        }
        if (r[d] < r[c]) {
            std::swap(r[c], r[d]);
            std::swap(r[a], r[b]);
        }
        if (r[c] < r[b]) {
            std::swap(r[b], r[c]);
        }
        assert(r[a] <= r[c] && r[b] <= r[c] && r[c] <= r[d]);
    }
    else
    {
        // In the median of 5 algorithm consider r[a] infinitely small, then
        // change b->a. c->b, d->c, e->d
        if (r[c] < r[a]) {
            std::swap(r[a], r[c]);
        }
        if (r[c] < r[b]) {
            std::swap(r[b], r[c]);
        }
        if (r[d] < r[a]) {
            std::swap(r[a], r[d]);
        }
        if (r[d] < r[b]) {
            std::swap(r[b], r[d]);
        } else {
            if (r[b] < r[a]) {
                std::swap(r[a], r[b]);
            }
        }
        assert(r[a] <= r[b] && r[b] <= r[c] && r[b] <= r[d]);
    }
}

I started with Teppo's median-of-5 algorithm and eliminated the first/last element. Do you think you can do better?


Thanks,

Andrei

On 26.11.2016 18:17, Andrei Alexandrescu wrote:
> I'm submitting the median paper to
> http://2017.programmingconference.org/track/programming-2017-papers. (It
> was rejected by ALENEX... the short story is there was no open-source
> implementation and they simply didn't buy the results.)
>
> Now I have a C++ implementation available. I wanted to ask for your help
> with lower/upper median of 4, which I coded like this (based on the
> median of 5 posted here by user "tn" aka Teppo Niinimäki):
>
> template <bool leanRight, class T>
> void medianOf(T* r, size_t a, size_t b, size_t c, size_t d)
> {
>     assert(a != b && a != c && a != d && b != c && b != d
>         && c != d);
>     if (leanRight)
>     {
>         // In the median of 5 algorithm, consider r[e] infinite
>         if (r[c] < r[a]) {
>             std::swap(r[a], r[c]);
>         }
>         if (r[d] < r[b]) {
>             std::swap(r[b], r[d]);
>         }
>         if (r[d] < r[c]) {
>             std::swap(r[c], r[d]);
>             std::swap(r[a], r[b]);
>         }
>         if (r[c] < r[b]) {
>             std::swap(r[b], r[c]);
>         }
>         assert(r[a] <= r[c] && r[b] <= r[c] && r[c] <= r[d]);
>     }
>     else
>     {
>         // In the median of 5 algorithm consider r[a] infinitely small,
> then
>         // change b->a. c->b, d->c, e->d
>         if (r[c] < r[a]) {
>             std::swap(r[a], r[c]);
>         }
>         if (r[c] < r[b]) {
>             std::swap(r[b], r[c]);
>         }
>         if (r[d] < r[a]) {
>             std::swap(r[a], r[d]);
>         }
>         if (r[d] < r[b]) {
>             std::swap(r[b], r[d]);
>         } else {
>             if (r[b] < r[a]) {
>                 std::swap(r[a], r[b]);
>             }
>         }
>         assert(r[a] <= r[b] && r[b] <= r[c] && r[b] <= r[d]);
>     }
> }
>
> I started with Teppo's median-of-5 algorithm and eliminated the
> first/last element. Do you think you can do better?
>
>
> Thanks,
>
> Andrei

The following code (which I found using a quick brute-force search) uses at most 4 comparisons.


int median(bool leanRight)(int[4] a){
    static if(leanRight){
        return a[a[0]<a[1]?a[2]<a[3]?a[1]<a[3]?a[1]<a[2]?2:1:a[0]<a[3]?3:0:a[1]<a[2]?a[1]<a[3]?3:1:a[0]<a[2]?2:0:a[2]<a[3]?a[0]<a[3]?a[0]<a[2]?2:0:a[1]<a[3]?3:1:a[0]<a[2]?a[0]<a[3]?3:0:a[1]<a[2]?2:1];
    }else{
        return a[a[0]<a[1]?a[2]<a[3]?a[0]<a[2]?a[1]<a[2]?1:2:a[0]<a[3]?0:3:a[0]<a[3]?a[1]<a[3]?1:3:a[0]<a[2]?0:2:a[2]<a[3]?a[1]<a[2]?a[0]<a[2]?0:2:a[1]<a[3]?1:3:a[1]<a[3]?a[0]<a[3]?0:3:a[1]<a[2]?1:2];
    }
}

I guess you can also adapt the leanRight case to the leanLeft case by inverting all comparisons, then you should also get an algorithm that uses at most 4 comparisons.

On 11/26/2016 03:19 PM, Timon Gehr wrote:
>
> The following code (which I found using a quick brute-force search) uses
> at most 4 comparisons.
>
>
> int median(bool leanRight)(int[4] a){
>     static if(leanRight){
>         return
> a[a[0]<a[1]?a[2]<a[3]?a[1]<a[3]?a[1]<a[2]?2:1:a[0]<a[3]?3:0:a[1]<a[2]?a[1]<a[3]?3:1:a[0]<a[2]?2:0:a[2]<a[3]?a[0]<a[3]?a[0]<a[2]?2:0:a[1]<a[3]?3:1:a[0]<a[2]?a[0]<a[3]?3:0:a[1]<a[2]?2:1];
>
>     }else{
>         return
> a[a[0]<a[1]?a[2]<a[3]?a[0]<a[2]?a[1]<a[2]?1:2:a[0]<a[3]?0:3:a[0]<a[3]?a[1]<a[3]?1:3:a[0]<a[2]?0:2:a[2]<a[3]?a[1]<a[2]?a[0]<a[2]?0:2:a[1]<a[3]?1:3:a[1]<a[3]?a[0]<a[3]?0:3:a[1]<a[2]?1:2];
>
>     }
> }

Thanks! Sorry, I forgot to mention partitioning around the lower/upper median is also needed. It seems that changes the tradeoff space. -- Andrei

On Saturday, 26 November 2016 at 20:19:38 UTC, Timon Gehr wrote:
> The following code (which I found using a quick brute-force search) uses at most 4 comparisons.

 Which is probably the way to go. I was considering doing something similar for small sorting and the like, letting it get an optimized sort of say 8 items or less, and saving the results for quick compiling later.

On 26.11.2016 22:11, Andrei Alexandrescu wrote:
> On 11/26/2016 03:19 PM, Timon Gehr wrote:
>>
>> The following code (which I found using a quick brute-force search) uses
>> at most 4 comparisons.
>>
>>
>> int median(bool leanRight)(int[4] a){
>>     static if(leanRight){
>>         return
>> a[a[0]<a[1]?a[2]<a[3]?a[1]<a[3]?a[1]<a[2]?2:1:a[0]<a[3]?3:0:a[1]<a[2]?a[1]<a[3]?3:1:a[0]<a[2]?2:0:a[2]<a[3]?a[0]<a[3]?a[0]<a[2]?2:0:a[1]<a[3]?3:1:a[0]<a[2]?a[0]<a[3]?3:0:a[1]<a[2]?2:1];
>>
>>
>>     }else{
>>         return
>> a[a[0]<a[1]?a[2]<a[3]?a[0]<a[2]?a[1]<a[2]?1:2:a[0]<a[3]?0:3:a[0]<a[3]?a[1]<a[3]?1:3:a[0]<a[2]?0:2:a[2]<a[3]?a[1]<a[2]?a[0]<a[2]?0:2:a[1]<a[3]?1:3:a[1]<a[3]?a[0]<a[3]?0:3:a[1]<a[2]?1:2];
>>
>>
>>     }
>> }
>
> Thanks! Sorry, I forgot to mention partitioning around the lower/upper
> median is also needed. It seems that changes the tradeoff space. -- Andrei

Did you see the suggestion to make the leanLeft/leanRight cases symmetric, such that both use at most 4 branches?

On 26.11.2016 22:36, Timon Gehr wrote:
>>>
>>
>> Thanks! Sorry, I forgot to mention partitioning around the lower/upper
>> median is also needed. It seems that changes the tradeoff space. --
>> Andrei
>
> Did you see the suggestion to make the leanLeft/leanRight cases
> symmetric, such that both use at most 4 branches?

Code:

void medianOf(bool leanRight,T)(ref T[4] r, size_t a, size_t b, size_t c, size_t d){
    assert(a!=b&&a!=c&&a!=d&&b!=c&&b!=d&&c!=d);
    if(r[c]<r[a]) swap(r[a],r[c]);
    if(r[d]<r[b]) swap(r[b],r[d]);
    if(leanRight?r[d]<r[c]:r[b]<r[a]){
        swap(r[c],r[d]);
        swap(r[a],r[b]);
    }
    if(r[c]<r[b]) swap(r[b],r[c]);
    if(leanRight) assert(r[a]<=r[c]&&r[b]<=r[c]&&r[c]<=r[d]);
    else assert(r[a]<=r[b]&&r[b]<=r[c]&&r[b]<=r[d]);
}

On 26.11.2016 23:16, Timon Gehr wrote:
> On 26.11.2016 22:36, Timon Gehr wrote:
>>>>
>>>
>>> Thanks! Sorry, I forgot to mention partitioning around the lower/upper
>>> median is also needed. It seems that changes the tradeoff space. --
>>> Andrei
>>
>> Did you see the suggestion to make the leanLeft/leanRight cases
>> symmetric, such that both use at most 4 branches?
>
> Code:
>
> void medianOf(bool leanRight,T)(ref T[4] r, size_t a, size_t b, size_t
> c, size_t d){
>     assert(a!=b&&a!=c&&a!=d&&b!=c&&b!=d&&c!=d);
>     if(r[c]<r[a]) swap(r[a],r[c]);
>     if(r[d]<r[b]) swap(r[b],r[d]);
>     if(leanRight?r[d]<r[c]:r[b]<r[a]){
>         swap(r[c],r[d]);
>         swap(r[a],r[b]);
>     }
>     if(r[c]<r[b]) swap(r[b],r[c]);
>     if(leanRight) assert(r[a]<=r[c]&&r[b]<=r[c]&&r[c]<=r[d]);
>     else assert(r[a]<=r[b]&&r[b]<=r[c]&&r[b]<=r[d]);
> }
>

Slightly larger code, but uses fewer swaps:

void medianOf(bool leanRight,T)(ref T[4] r, size_t a, size_t b, size_t c, size_t d){
    assert(a!=b&&a!=c&&a!=d&&b!=c&&b!=d&&c!=d);
    if(r[c]<r[a]) swap(r[a],r[c]);
    if(r[d]<r[b]) swap(r[b],r[d]);
    static if(leanRight){
        if(r[d]<r[c]){
            swap(r[c],r[d]);
            if(r[c]<r[a]) swap(r[a],r[c]);
        }else if(r[c]<r[b]) swap(r[b],r[c]);
        assert(r[a]<=r[c]&&r[b]<=r[c]&&r[c]<=r[d]);
    }else{
        if(r[b]<r[a]){
            swap(r[a],r[b]);
            if(r[d]<r[b]) swap(r[b],r[d]);
        }else if(r[c]<r[b]) swap(r[b],r[c]);
        assert(r[a]<=r[b]&&r[b]<=r[c]&&r[b]<=r[d]);
    }
}

The best strategy is probably to exhaustively enumerate all Pareto-optimal algorithms and to benchmark them automatically within the larger algorithm.

Thanks! I'll experiment with these and let you know what was best. -- Andrei

On 11/26/2016 05:30 PM, Timon Gehr wrote:
>> Code:
>>
>> void medianOf(bool leanRight,T)(ref T[4] r, size_t a, size_t b, size_t
>> c, size_t d){
>>     assert(a!=b&&a!=c&&a!=d&&b!=c&&b!=d&&c!=d);
>>     if(r[c]<r[a]) swap(r[a],r[c]);
>>     if(r[d]<r[b]) swap(r[b],r[d]);
>>     if(leanRight?r[d]<r[c]:r[b]<r[a]){
>>         swap(r[c],r[d]);
>>         swap(r[a],r[b]);
>>     }
>>     if(r[c]<r[b]) swap(r[b],r[c]);
>>     if(leanRight) assert(r[a]<=r[c]&&r[b]<=r[c]&&r[c]<=r[d]);
>>     else assert(r[a]<=r[b]&&r[b]<=r[c]&&r[b]<=r[d]);
>> }
>>
>
> Slightly larger code, but uses fewer swaps:
>
> void medianOf(bool leanRight,T)(ref T[4] r, size_t a, size_t b, size_t
> c, size_t d){
>     assert(a!=b&&a!=c&&a!=d&&b!=c&&b!=d&&c!=d);
>     if(r[c]<r[a]) swap(r[a],r[c]);
>     if(r[d]<r[b]) swap(r[b],r[d]);
>     static if(leanRight){
>         if(r[d]<r[c]){
>             swap(r[c],r[d]);
>             if(r[c]<r[a]) swap(r[a],r[c]);
>         }else if(r[c]<r[b]) swap(r[b],r[c]);
>         assert(r[a]<=r[c]&&r[b]<=r[c]&&r[c]<=r[d]);
>     }else{
>         if(r[b]<r[a]){
>             swap(r[a],r[b]);
>             if(r[d]<r[b]) swap(r[b],r[d]);
>         }else if(r[c]<r[b]) swap(r[b],r[c]);
>         assert(r[a]<=r[b]&&r[b]<=r[c]&&r[b]<=r[d]);
>     }
> }
>
> The best strategy is probably to exhaustively enumerate all
> Pareto-optimal algorithms and to benchmark them automatically within the
> larger algorithm.