Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- //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()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement