lfsr.cpp

Go to the documentation of this file.

/* lsfr.cpp :  Linear shift feedback register
 * Copyright (C) 2017  Daniel Haley
 * 
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "atomprobe/helper/maths/lfsr.h"
#include "atomprobe/helper/aptAssert.h"

namespace AtomProbe 
{

//For the table to work, we need the sizeof(size_T) at preprocess time
#ifndef SIZEOF_SIZE_T
#error sizeof(size_t) macro is undefined... At time of writing, this is usually 4 (32 bit) or 8. You can work it out from a simple C++ program which prints out sizeof(size_t). This cant be done automatically due to preprocessor behaviour.
#endif

//Maximum period linear shift register values (computed by
//other program for Galois polynomial)
//Unless otherwise noted, all these entries have been verified using the
//verifyTable routine. 
//
//If you don't trust me, (who doesn't trust some random person on the internet?) 
//you can re-run the verification routine. 
//
//Note that verification time *doubles* with every entry, so full 64-bit verification
//is computationally intensive. I achieved 40 bits in half a day. 48 bits took over a week.
const size_t maximumLinearTable[] = {
      0x03,
      0x06,
      0x0C,
      0x14,
      0x30,
      0x60,
      0xb8,
    0x0110,
    0x0240,
    0x0500,
    0x0e08,
    0x1c80,
    0x3802,
    0x6000,
    0xb400,
    0x12000,
    0x20400,
    0x72000,
    0x90000,
    0x140000,
    0x300000,
    0x420000,
    0xD80000,
    0x1200000,
    0x3880000,
    0x7200000,
    0x9000000,
    0x14000000,
    0x32800000,
    0x48000000,

#if (SIZEOF_SIZE_T > 4)
    0xA3000000,
    0x100080000,
    0x262000000,
    0x500000000,
    0x801000000,
    0x1940000000,
    0x3180000000,
    0x4400000000,
    0x9C00000000,
    0x12000000000,
    0x29400000000,
    0x63000000000,
    0xA6000000000,
    0x1B0000000000,
    0x20E000000000,
    0x420000000000,
    0x894000000000,
    //Maximal linear table entries below line are unverified 
    //Verifying the full table might take a few months of computing time
    //But who needs to count beyond 2^49-1 (~10^14) anyway??
    0x1008000000000,

    //Really, there are more entries beyond this, but I consider it pretty much not worth the effort.
#endif
};

void LinearFeedbackShiftReg::setMaskPeriod(unsigned int newMask)
{
    //Don't fall off the table
    ASSERT((newMask-3) < sizeof(maximumLinearTable)/sizeof(size_t));

    maskVal=maximumLinearTable[newMask-3];

    //Set the mask to be all ones
    totalMask=0;
    for(size_t ui=0;ui<newMask;ui++)
        totalMask|= (size_t)(1)<<ui;


}

bool LinearFeedbackShiftReg::verifyTable(size_t maxLen)
{
    size_t tableLen =  sizeof(maximumLinearTable)/sizeof(size_t);
    
    //check each one is actually the full 2^n-1 period
    if(maxLen)
    {
        ASSERT(maxLen < tableLen);
        tableLen=maxLen;
    }

    //For the 32 bit table, this works pretty quickly.
    //for the 64  bit table, this takes a month or so
    //The test here is a little weak, all it tests is that the return to
    // 1 is such that there are enough visitations to visit each number exactly once
    // no more no less. However, it doesn't check each number is actually visited, 
    // as this would be memory intensive
    for(size_t n=3;n<tableLen+3;n++)
    {
        size_t period;
        setState(1);
        setMaskPeriod(n);
        period=0;
        do
        {
            clock();
            period++;
        }
        while(lfsr!=1);


        //we should have counted every bit position in the polynomial (except 0)
        //otherwise, this is not the maximal linear sequence
        if(period != ((size_t)(1)<<(n-(size_t)1)) -(size_t)(1))
            return false;
    }
    return true;
}

size_t LinearFeedbackShiftReg::clock()
{
    typedef size_t ull;

    lfsr = (lfsr >> 1) ^  ( (-(lfsr & (ull)(1u))) & maskVal ); 
    lfsr&=totalMask;
    if( lfsr == 0u)
        lfsr=1u;

    return lfsr;
}

}