fast pairwise

pairwiseFast.d

import std.range, std.traits;
import core.bitop : bsf;
    
import std.stdio;

/++
$(LUCKY Pairwise summation) algorithm. Range must be a finite range.
+/
F sumPairwise(F, R)(R data)
if(isInputRange!R && !isInfinite!R)
{
    version(LDC) pragma(LDC_never_inline);
    // Works for r with at least length < 2^^(64 + log2(16)), in keeping with the use of size_t
    // elsewhere in std.algorithm and std.range on 64 bit platforms.
    F[64] store = void;
    size_t idx = 0;

    static if (hasLength!R)
    {
        foreach(k; 0 .. data.length / 16)
        {
            static if (isRandomAccessRange!R && hasSlicing!R)
            {
                store[idx] = sumPairwise16!F(data);
                data = data[16 .. $];
            }
            else store[idx] = sumPairwiseN!(16, false, F)(data);
            foreach(_; 0 .. cast(uint)bsf(k+1))
            {
                store[idx -1] += store[idx];
                --idx;
            }
            ++idx;
        }    
        
        size_t i = 0;
        foreach (el; data)
        {
            store[idx] = el;
            foreach (_; 0 .. cast(uint)bsf(i+1))
            {
                store[idx - 1] += store[idx];
                --idx;
            }
            ++idx;
            ++i;
        }
    }
    else
    {
        size_t k = 0;
        while (!data.empty)
        {
            store[idx] = sumPairwiseN!(16, true, F)(data);
            foreach(_; 0 .. cast(uint)bsf(k+1))
            {
                store[idx -1] += store[idx];
                --idx;
            }
            ++idx;
            ++k;
        }
    }
   
    F s = store[idx - 1];
    foreach_reverse (j; 0 .. idx - 1)
        s += store[j];
    
    return s;
}

private auto sumPairwise8(F, R)(R r)
if(isRandomAccessRange!R)
{
    return (((cast(F)r[0] + r[1]) + (cast(F)r[2] + r[3]))
          + ((cast(F)r[4] + r[5]) + (cast(F)r[6] + r[7])));
}

private auto sumPairwise16(F, R)(R r)
if(isRandomAccessRange!R)
{
    return (((cast(F)r[ 0] + r[ 1]) + (cast(F)r[ 2] + r[ 3]))
          + ((cast(F)r[ 4] + r[ 5]) + (cast(F)r[ 6] + r[ 7])))
         + (((cast(F)r[ 8] + r[ 9]) + (cast(F)r[10] + r[11]))
          + ((cast(F)r[12] + r[13]) + (cast(F)r[14] + r[15])));
}

private auto sumPair(bool needEmptyChecks, F, R)(ref R r)
if(isForwardRange!R && !isRandomAccessRange!R)
{
    static if(needEmptyChecks) if(r.empty) return F(0);
    F s0 = r.front;
    r.popFront();
    static if(needEmptyChecks) if(r.empty) return s0;
    s0 += r.front;
    r.popFront();
    return s0;
}

private auto sumPairwiseN(size_t N, bool needEmptyChecks, F, R)(ref R r)
if(isForwardRange!R && !isRandomAccessRange!R)
{
    static assert(!(N & (N-1))); //isPow2
    static if(N == 2) return sumPair!(needEmptyChecks, F)(r);
    else return sumPairwiseN!(N/2, needEmptyChecks, F)(r)
        + sumPairwiseN!(N/2, needEmptyChecks, F)(r);
}