Count nucleotides in 2-bit encoding (Slower than bitmath, but still fun.)

nuccount.cpp

INLINE u32 nuccount(u64 kmer, unsigned k) {
    kmer ^= XOR_MASK;
    u32 counts(0);
#if 1
    static const u32 lut4 [] {
        67108864, 50397184, 50331904, 50331649, 50397184, 33685504, 33620224, 33619969, 50331904, 33620224, 33554944, 33554689, 50331649, 33619969, 33554689, 33554434,
        50397184, 33685504, 33620224, 33619969, 33685504, 16973824, 16908544, 16908289, 33620224, 16908544, 16843264, 16843009, 33619969, 16908289, 16843009, 16842754,
        50331904, 33620224, 33554944, 33554689, 33620224, 16908544, 16843264, 16843009, 33554944, 16843264, 16777984, 16777729, 33554689, 16843009, 16777729, 16777474,
        50331649, 33619969, 33554689, 33554434, 33619969, 16908289, 16843009, 16842754, 33554689, 16843009, 16777729, 16777474, 33554434, 16842754, 16777474, 16777219,
        50397184, 33685504, 33620224, 33619969, 33685504, 16973824, 16908544, 16908289, 33620224, 16908544, 16843264, 16843009, 33619969, 16908289, 16843009, 16842754,
        33685504, 16973824, 16908544, 16908289, 16973824, 262144, 196864, 196609, 16908544, 196864, 131584, 131329, 16908289, 196609, 131329, 131074,
        33620224, 16908544, 16843264, 16843009, 16908544, 196864, 131584, 131329, 16843264, 131584, 66304, 66049, 16843009, 131329, 66049, 65794,
        33619969, 16908289, 16843009, 16842754, 16908289, 196609, 131329, 131074, 16843009, 131329, 66049, 65794, 16842754, 131074, 65794, 65539,
        50331904, 33620224, 33554944, 33554689, 33620224, 16908544, 16843264, 16843009, 33554944, 16843264, 16777984, 16777729, 33554689, 16843009, 16777729, 16777474,
        33620224, 16908544, 16843264, 16843009, 16908544, 196864, 131584, 131329, 16843264, 131584, 66304, 66049, 16843009, 131329, 66049, 65794,
        33554944, 16843264, 16777984, 16777729, 16843264, 131584, 66304, 66049, 16777984, 66304, 1024, 769, 16777729, 66049, 769, 514,
        33554689, 16843009, 16777729, 16777474, 16843009, 131329, 66049, 65794, 16777729, 66049, 769, 514, 16777474, 65794, 514, 259,
        50331649, 33619969, 33554689, 33554434, 33619969, 16908289, 16843009, 16842754, 33554689, 16843009, 16777729, 16777474, 33554434, 16842754, 16777474, 16777219,
        33619969, 16908289, 16843009, 16842754, 16908289, 196609, 131329, 131074, 16843009, 131329, 66049, 65794, 16842754, 131074, 65794, 65539,
        33554689, 16843009, 16777729, 16777474, 16843009, 131329, 66049, 65794, 16777729, 66049, 769, 514, 16777474, 65794, 514, 259,
        33554434, 16842754, 16777474, 16777219, 16842754, 131074, 65794, 65539, 16777474, 65794, 514, 259, 16777219, 65539, 259, 4,
    };
    static const u32 lut3 [] {
        50331648, 33619968, 33554688, 33554433, 33619968, 16908288, 16843008, 16842753, 33554688, 16843008, 16777728, 16777473, 33554433, 16842753, 16777473, 16777218,
        33619968, 16908288, 16843008, 16842753, 16908288, 196608, 131328, 131073, 16843008, 131328, 66048, 65793, 16842753, 131073, 65793, 65538,
        33554688, 16843008, 16777728, 16777473, 16843008, 131328, 66048, 65793, 16777728, 66048, 768, 513, 16777473, 65793, 513, 258,
        33554433, 16842753, 16777473, 16777218, 16842753, 131073, 65793, 65538, 16777473, 65793, 513, 258, 16777218, 65538, 258, 3,
    };
#endif
    static const u32 lut2 [] {
        33554432, 16842752, 16777472, 16777217, 16842752, 131072, 65792, 65537, 16777472, 65792, 512, 257, 16777217, 65537, 257, 2,
    };
    static const u32 lut1 [] {16777216, 65536, 256, 1};
    static const u32 lut0 [] {0};
#if 1
    static const u32 *luts[] {lut0, lut1, lut2, lut3};
    k >>= 2;
    switch(k) {
        case 8: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 7: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 6: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 5: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 4: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 3: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 2: counts += lut4[kmer & 0xFF], kmer >>= 8, --k; [[fallthrough]]
        case 1: counts += lut4[kmer & 0xFF], kmer >>= 8, --k;
    }
#else
    static const u32 *luts[] {lut0, lut1};
    k >>= 1;
    switch(k) {
        case 16: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 15: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 14: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 13: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 12: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 11: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 10: counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 9:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 8:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 7:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 6:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 5:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 4:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 3:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 2:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k; [[fallthrough]]
        case 1:  counts += lut2[kmer & 0xFF], kmer >>= 4, --k;
    }
#endif
    return counts += luts[k][kmer];
}

This is inferior to the constant-time bit-bashing technique used in Bowtie2. I'm mostly glad it worked at all.