Test generic hamming coding routines (SECDEC). These can fix 1 bit errors and detect two bit errors.

gistfile1.cpp

/* hamming_test.cpp
Test generic hamming coding routines (SECDEC). These can fix 1 bit errors and detect two bit errors.
(C) 2015 Niall Douglas http://www.nedprod.com

On my 3.9Ghz Ivy Bridge Intel 3770K:

For no intrinsics:
Calculating 4Kb: 12.06Mb/sec
     Fixing 4Kb: 14.98Mb/sec

With intrinsics (bitscan for set bits):
Calculating 4Kb: 73.78Mb/sec
     Fixing 4Kb: 64.15Mb/sec

With tables (bitscan for set bits):
Calculating 4Kb: 82.41Mb/sec
     Fixing 4Kb: 72.58Mb/sec
     
With OpenMP (fixed bit masking):
Calculating 4Kb: 592Mb/sec
     Fixing 4Kb: 189Mb/sec
*/

#include <iostream>
#include <random>
#include <vector>
#include <chrono>
#include <string.h>

#define USE_INTRINSICS

static size_t ecc_calculate(const char *buffer, size_t length);
__attribute__((noinline)) size_t test_ecc_calculate_4096(const char *buffer)
{
  return ecc_calculate(buffer, 4096);
}

// Called before any other usage
static size_t ecc_twospowers_[8*sizeof(size_t)];
static std::vector<size_t> ecc_table;
static size_t ecc_twospowers(size_t i)
{
#if 1
  return ecc_twospowers_[i];
#else
  return ((size_t)1<<i);
#endif
}
static void ecc_initialise(size_t length)
{
  for(size_t n=0; n<8*sizeof(size_t); n++)
    ecc_twospowers_[n]=((size_t)1<<n);
  ecc_table.clear();
  length*=8;
  ecc_table.resize(length);
  for(size_t i=0; i<length; i++)
  {
    // Make a code bit
    size_t b=i+1;
#ifdef USE_INTRINSICS
    size_t topbit=8*sizeof(size_t)-__builtin_clzl(b);
    b+=topbit;
    if(b>=ecc_twospowers(topbit)) b++;
    //while(b>ecc_twospowers(_topbit+1)) _topbit++;
    //b+=_topbit;
    //if(b>=ecc_twospowers(_topbit)) b++;
#else
    for(size_t p=0; ecc_twospowers(p)<(b+1); p++)
      b++;
#endif
    ecc_table[i]=b;
  }
}

// Return how many bits of ECC are needed for length bytes
static size_t ecc_bit_length(size_t length)
{
  length*=8;
  // This is (data bits + parity bits + 1) <= 2^(parity bits)
  for(size_t p=1; p<sizeof(size_t)*8; p++)
    if((length+p+1)<=ecc_twospowers(p))
      return p;
  return 0;
}

static size_t ecc_calculate(const char *buffer, size_t length)
{
  size_t ecc=0;
#ifdef USE_INTRINSICS
  size_t _topbit=0, _temp=1, i=0, _i=0, _bitset=0, *_buffer=(size_t *) buffer, *_end=(size_t *)(buffer+length);
  length*=8;
  _i-=8*sizeof(size_t);
  for(;;)
  {
    // Clear the previous set bit and get next set bit
    _temp^=ecc_twospowers(_bitset);
    if(!_temp)
    {
      if(_buffer>=_end)
        break;
      do
      {
        _temp=*_buffer++;
        _i+=8*sizeof(size_t);
      } while(!_temp && _buffer<_end);
      if(!_temp) break;
    }
    _bitset=__builtin_ctzl(_temp); // FIXME Assumes little endian!
    //std::cout << "Offset " << _i << " is " << std::hex << _temp << std::dec << " next set bit is " << _bitset << std::endl;
    if((i=_i+_bitset)>=length)
      break;
#else
  length*=8;
  for(size_t i=0; i<length; i++)
  {
    // Find next set bit
    if(!(buffer[i/8] & ecc_twospowers(i % 8)))
      continue;
#endif
#if 1
    size_t b=ecc_table[i];
#else
    // Make a code bit
    size_t b=i+1;
#ifdef USE_INTRINSICS
    size_t topbit=8*sizeof(size_t)-__builtin_clzl(b);
    b+=topbit;
    if(b>=ecc_twospowers(topbit)) b++;
    //while(b>ecc_twospowers(_topbit+1)) _topbit++;
    //b+=_topbit;
    //if(b>=ecc_twospowers(_topbit)) b++;
#else
    for(size_t p=0; ecc_twospowers(p)<(b+1); p++)
      b++;
#endif
#endif
    ecc^=b;
  }
  return ecc;
}

// Increase parallelism by working in batches
static size_t ecc_calculate2(const char *buffer, size_t length)
{
  typedef unsigned char unit_t;
  size_t ecc=0;
  unit_t *_buffer=(unit_t *) buffer;
#pragma omp parallel for reduction(^:ecc)
  for(size_t i=0; i<length; i+=sizeof(unit_t))
  {
    unit_t c=_buffer[i/sizeof(unit_t)];
    if(!c)
      continue;
    char bitset[8*sizeof(unit_t)];
    for(size_t n=0; n<8*sizeof(unit_t); n++)
      bitset[n]=!!(c&((size_t)1<<n));
#if 1
    for(size_t n=0; n<8*sizeof(unit_t); n++)
    {
      if(bitset[n])
      {
        ecc^=ecc_table[i*8+n];
      }
    }
#else
    // Make code bits
    size_t b[8*sizeof(unit_t)];
    for(size_t n=0; n<8*sizeof(unit_t); n++)
      b[n]=i*8+n+1;
    size_t t[8*sizeof(unit_t)];
    for(size_t n=0; n<8*sizeof(unit_t); n++)
      t[n]=8*sizeof(size_t)-__builtin_clzl(b[n]);
    for(size_t n=0; n<8*sizeof(unit_t); n++)
      b[n]+=t[n];
    for(size_t n=0; n<8*sizeof(unit_t); n++)
      if(b[n]>=ecc_twospowers(t[n])) b[n]++;
    for(size_t n=0; n<8*sizeof(unit_t); n++)
    {
      if(bitset[n])
      {
        ecc^=b[n];
      }
    }
#endif
  }
  return ecc;
}

// Simple OpenMP
static size_t ecc_calculate3(const char *buffer, size_t length)
{
  typedef size_t unit_t;
  size_t ecc=0;
  unit_t *_buffer=(unit_t *) buffer;
#pragma omp parallel for reduction(^:ecc)
  for(size_t i=0; i<length; i+=sizeof(unit_t))
  {
    unit_t c=_buffer[i/sizeof(unit_t)];
    while(c)
    {
      size_t bitset=__builtin_ctzl(c);
      if(i+bitset/8>=length) break;
      c&=~ecc_twospowers(bitset);
      ecc^=ecc_table[i*8+bitset];
    }
  }
  return ecc;
}

static bool is_single_bit_set(size_t v)
{
#ifdef USE_INTRINSICS
  return __builtin_popcountl(v)==1;
#else
  size_t count=0;
  for(size_t i=0; i<8*sizeof(size_t); i++)
    if(v&ecc_twospowers(i))
      count++;
  return count==1;
#endif
}

static size_t ecc_find_badbit(const char *buffer, size_t length, size_t eccdiff)
{
  length*=8;
  // Is the ECC itself corrupt?
  if(is_single_bit_set(eccdiff))
    return (size_t)-1;
  for(size_t i=0, b=1; i<length; i++, b++)
  {
    // Skip parity bits
    while(is_single_bit_set(b))
      b++;
    if(b==eccdiff)
      return i;
  }
  return (size_t)-1;
}

#define ecc_calculate ecc_calculate2
int main(int argc, const char *argv[])
{
  int errcode=0;
  if(argc<2)
  {
    std::cerr << "Usage: " << argv[0] << " <bytes>" << std::endl;
    return 1;
  }
  size_t bytes=atoi(argv[1]);
  //bytes=4096;
  std::cout << "Filling a buffer with " << bytes << " of randomness" << std::endl;
  std::vector<char> buffer(bytes);
  {
    std::random_device r;
    std::random_device::result_type *b=(std::random_device::result_type *) buffer.data();
    for(size_t n=0; n<buffer.size()/sizeof(*b); n++)
      b[n]=r();
  }
  ecc_initialise(bytes);
  size_t eccbits=ecc_bit_length(bytes);
  std::cout << "ECC will be " << eccbits << " bits long" << std::endl;
  if(eccbits/8+1>8*sizeof(size_t))
  {
    std::cerr << "ECC would exceed the length of a size_t!" << std::endl;
    return 1;
  }
  size_t ecc=ecc_calculate(buffer.data(), bytes);
  std::cout << "ECC was calculated to be " << std::hex << ecc << std::dec << std::endl;
  if(ecc!=ecc_calculate2(buffer.data(), bytes))
  {
    std::cout << "ERROR: Batch ecc compute does not equal non-batch!" << std::endl;
    errcode=1;
  }

  std::cout << "Flipping every bit in the buffer to see if it is correctly detected ..." << std::endl;  
  auto begin=std::chrono::high_resolution_clock::now();  
  for(size_t toflip=0; toflip<bytes*8; toflip++)
  {
    buffer[toflip/8]^=((size_t)1<<(toflip%8));
    size_t newecc=ecc_calculate(buffer.data(), bytes);
    if(ecc==newecc)
    {
      std::cerr << "ERROR: Flipping bit " << toflip << " not detected!" << std::endl;
      errcode=2;
    }
    else if(is_single_bit_set(ecc ^ newecc))
    {
      std::cerr << "ERROR: He claims ECC is corrupt!" << std::endl;
      errcode=2;
    }
    else
    {
      size_t badbit=ecc_find_badbit(buffer.data(), bytes, ecc ^ newecc);
      if(badbit!=toflip)
      {
        std::cerr << "ERROR: Bad bit " << badbit << " is not the bit " << toflip << " we flipped!" << std::endl;
        errcode=2;
      }
//      else
//        std::cout << "SUCCESS: Bit flip " << toflip << " correctly detected" << std::endl;
    }
    buffer[toflip/8]^=((size_t)1<<(toflip%8));
  }
  auto end=std::chrono::high_resolution_clock::now();
  auto diff=std::chrono::duration_cast<std::chrono::duration<double, std::ratio<1, 1>>>(end-begin);
  std::cout << "Checking and fixing is approximately " << (bytes*10000/diff.count()/1024/1024) << " Mb/sec" << std::endl;
  
  std::cout << "\nFlipping two bits in the buffer to see if it is correctly detected ..." << std::endl;  
  buffer[0]^=1;
  begin=std::chrono::high_resolution_clock::now();  
  for(size_t toflip=1; toflip<bytes*8; toflip++)
  {
    buffer[toflip/8]^=((size_t)1<<(toflip%8));
    size_t newecc=ecc_calculate(buffer.data(), bytes);
    if(ecc==newecc)
    {
      std::cerr << "ERROR: Flipping bits 0 and " << toflip << " not detected!" << std::endl;
      errcode=3;
    }
//    else
//    {
//      std::cout << "SUCCESS: Flipping bits 0 and " << toflip << " correctly detected" << std::endl;
//    }
    buffer[toflip/8]^=((size_t)1<<(toflip%8));
  }
  end=std::chrono::high_resolution_clock::now();
  diff=std::chrono::duration_cast<std::chrono::duration<double, std::ratio<1, 1>>>(end-begin);
  std::cout << "Calculating is approximately " << (bytes*10000/diff.count()/1024/1024) << " Mb/sec" << std::endl;

  std::cout << "\nCalculating speeds ..." << std::endl;
  size_t foo=0;
  begin=std::chrono::high_resolution_clock::now();  
  for(size_t n=0; n<10000; n++)
  {
    buffer[0]=(char)n;
    foo+=ecc_calculate(buffer.data(), bytes);
  }
  end=std::chrono::high_resolution_clock::now();
  diff=std::chrono::duration_cast<std::chrono::duration<double, std::ratio<1, 1>>>(end-begin);
  if(foo)
    std::cout << "Calculating is approximately " << (bytes*10000/diff.count()/1024/1024) << " Mb/sec" << std::endl;
  foo=0;
  begin=std::chrono::high_resolution_clock::now();  
  for(size_t n=0; n<10000; n++)
  {
    buffer[0]=(char)n;
    foo+=ecc_find_badbit(buffer.data(), bytes, ecc ^ ecc_calculate(buffer.data(), bytes));
  }
  end=std::chrono::high_resolution_clock::now();
  diff=std::chrono::duration_cast<std::chrono::duration<double, std::ratio<1, 1>>>(end-begin);
  if(foo)
    std::cout << "Checking and fixing is approximately " << (bytes*10000/diff.count()/1024/1024) << " Mb/sec" << std::endl;
  return errcode;
}