collatz glide kernel with reduction

//This kernel file should define LUTSTEP
//count collatz glide steps without actually comparing the value to its original. With reduction within the work group
//assuming SIEVESTEP = 32
//-cl-finite-math-only and -cl-no-signed-zeros may be used for faster floating point calculations

//constants
#define MAXSTEP 0xfffffffu
#define BITDIFF (32-LUTSTEP)
#define TABLEMASK ((1u << LUTSTEP) - 1)
#define ODDSCORE 1.5849625f

//Table of 3^N
__constant uint power3[21] =
{1u,3u,9u,27u,81u,243u,729u,2187u,6561u,19683u,
59049u,177147u,531441u,1594323u,4782969u,14348907u,43046721u,129140163u,387420489u,1162261467u,
3486784401u
};

//mul64() : takes two 32-bit uint and returns a uint2 as the multiply result
inline uint2 mul64(const uint a, const uint b){
  return (uint2)(a*b, mul_hi(a,b));
}

//collatzGlide() : count collatz glide and store the result of maximum glide.
//-----parameters-----
//result : s0 : max glide step. s1 : 'offset' of max glide. s23 : sum of steps (64-bit). Length = global work size
//sieve : Numbers needed to compute for every 2^32 numbers. Length = global work size
//table : look-up table for jumping. s0: x of 3^x mulVal. s1 : transformed b. Length = 2^ (LUTSTEP)
//negNadir : negative nadir score during a 'jump' by look-up table. Score is an approximation of log2 (val' / val). Length = 2^ (LUTSTEP)
//start : starting value (128-bit uint) as uint4
//l_maxStep : length = local work size
//l_totalOffset : length = local work size
//l_sumStep : length = local work size
//offset32 : actual offset = offset32 * (2^32)

__kernel void collatzGlide(
    __global uint *maxStep,
    __global ulong *totalOffset,
    __global ulong *sumStep,
    __global const uint *sieve,
    __global const uint2 *table,
    __global const float *negNadir,
    uint4 val,
    __local uint *l_maxStep,
    __local ulong *l_totalOffset,
    __local uint *l_sumStep,
    const uint offset32
){
    uint id = get_global_id(0);
    uint2 lut = val.s23; //lut : look-up table item. temporary store higher 64-bits of start
    val.s0 = sieve[id];
    val.s1 += offset32;
    val.s2 += (val.s1 < offset32); //carry operation of val
    val.s3 += (val.s2 < lut.s0);
    uint2 val_h = (uint2)(val.s3 < lut.s1 , 0u ); //val_h : higher 64-bit of val

    //save total offset to shared memory
    id = get_local_id(0);
    l_totalOffset[id] = upsample(offset32, val.s0);

    uint carryIdx = val.s0 & TABLEMASK; //carryIdx : 1. carry bits, 2. look-up table index
    uint stepCnt = 0; //stepCnt : count glide steps
    uint mulVal = 0; //mulVal : multiply value
    //uint2 lut; //lut : look-up table item
    float score = 0.0f; //score : approximation of log2(val'/ val)
    //The main loop, jump (LUTSTEP) iterations at a time, until the score hit nadir => the value will be less than its original value during next jump
    //The first iteration, val will always greater than its original because the sieve is applied; do{}while is used.
    do{
        lut = table[carryIdx]; //retrieve 'jump' table item
        stepCnt += LUTSTEP + lut.s0; //accumulate glide step
        mulVal = power3[lut.s0]; //multiply value = 3^(# odd encountered the in jump table)
        carryIdx = lut.s1;
        score += convert_float(lut.s0)* ODDSCORE - convert_float(LUTSTEP); //accumulate score, may use fma()

        //Do n = (n>>LUTSTEP)*a + b
        lut = mul64((val.s0 >> LUTSTEP) | (val.s1 << BITDIFF), mulVal); //bit 0-31
        val.s0 = lut.s0 + carryIdx;
        carryIdx = lut.s1 + (val.s0 < lut.s0); //carry bits
        lut = mul64((val.s1 >> LUTSTEP) | (val.s2 << BITDIFF), mulVal);//bit 32-63
        val.s1 = lut.s0 + carryIdx;
        carryIdx = lut.s1 + (val.s1 < lut.s0);
        lut = mul64((val.s2 >> LUTSTEP) | (val.s3 << BITDIFF), mulVal);//bit 64-95
        val.s2 = lut.s0 + carryIdx;
        carryIdx = lut.s1 + (val.s2 < lut.s0);
        lut = mul64((val.s3 >> LUTSTEP) | (val_h.s0 << BITDIFF), mulVal);//bit 96-127
        val.s3 = lut.s0 + carryIdx;
        carryIdx = lut.s1 + (val.s3 < lut.s0);
        lut = mul64((val_h.s0 >> LUTSTEP) | (val_h.s1 << BITDIFF), mulVal);//bit 128-159
        val_h.s0 = lut.s0 + carryIdx;
        carryIdx = lut.s1 + (val_h.s0 < lut.s0);
        lut = mul64((val_h.s1 >> LUTSTEP), mulVal);//bit 160-191 and overflow detection
        val_h.s1 = lut.s0 + carryIdx;

        lut.s1 += (val_h.s1 < lut.s0); //lut.s1 act as overflow bits
        carryIdx = val.s0 & TABLEMASK; //next table index
        //the loop should continue when all of the below are true
        //1. next nadir + score > 0 (val will be always greater than its original in the next jump so we can safely jump 'LUTSTEP' iterations)
        //2. no overflow
        //3. step < maxStep (2^28 -1)
        mulVal = (score > negNadir[carryIdx]) && (lut.s1 == 0) && (stepCnt < MAXSTEP);
    }while(mulVal);

    // Phase2 : find the position where the val fall under its original.
    do{
        mulVal = carryIdx & 1u; //odd or even?
        stepCnt += 1 + mulVal;
        carryIdx += select((uint)0, 1 + (carryIdx << 1), mulVal);
        carryIdx >>= 1;
        score += convert_float(mulVal) * ODDSCORE - 1.0f;
    }while(signbit(score) == 0); //while(score >= 0)

    //Save step to shared memory
    l_maxStep[id] = stepCnt;
    l_sumStep[id] = stepCnt;

    //Do reduction within the work group
    for(uint offset = get_local_size(0)/2; offset > 0; offset>>=1){
        barrier(CLK_LOCAL_MEM_FENCE);
        if(id < offset){
            uint idOther = id + offset;
            l_sumStep[id] += l_sumStep[idOther];

            stepCnt = l_maxStep[id];
            uint stepOther = l_maxStep[idOther];
            ulong offsetThis = l_totalOffset[id];
            ulong offsetOther = l_totalOffset[idOther];
            l_maxStep[id] = max(stepCnt, stepOther);
            l_totalOffset[id] = select(offsetOther, offsetThis, convert_ulong(stepCnt > stepOther));
        } //if(id < offset)
    } //for(offset)

    //save to output array (by the first thread in the work group)
    if(id == 0){
        uint wid = get_group_id(0);
        sumStep[wid] += l_sumStep[0]; //accumulate total steps
        uint stepOther = maxStep[wid];
        ulong offsetOther = totalOffset[wid];
        stepCnt = l_maxStep[0];
        ulong offsetThis = l_totalOffset[0];
        maxStep[wid] = max(stepCnt, stepOther);
        totalOffset[wid] = select( offsetOther, offsetThis, convert_ulong(stepCnt > stepOther));
    } //if(id == 0)
} //collatzGlide()