ROOT::Math C++ PRBS generator

[ferdymercury]

Describe the solution you'd like

It would be nice if ROOT provided a C++ function to generate a sequence of pseudo-random numbers that are frequently used in electronics testing, see:
https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence

Describe alternatives you've considered

• Perl https://metacpan.org/pod/Math::PRBS
• Matlab https://es.mathworks.com/help/serdes/ref/prbs.html
• Python https://www.reddit.com/r/askscience/comments/63a10q/for_prbs3_with_clock_input_on_each_gate_how_can/
• online examples https://www3.advantest.com/documents/11348/33b24c8a-c8cb-40b8-a2a7-37515ba4abc8 https://gist.github.com/mattbierner/d6d989bf26a7e54e7135

Additional context

Pulse generators usually generate this random bits. To know whether things are working in your data acquisition, you have to compare with the expected values.
https://cryptography.fandom.com/wiki/Linear_feedback_shift_register

There are also some snippets in ROOT using it, like here: https://root.cern/doc/master/civetweb_8c_source.html

[ferdymercury]

A possible implementation, that could go in some ROOT::Math classes:

#include <bitset>
#include <cmath>
#include <array>
#include <set>
#include <iostream>

using std::bitset;
using std::array;
using std::vector;
using std::set;
using std::cout;
using std::endl;

/**
 * @brief Generation of pseudo-random bits using a linear feedback shift register (LFSR), until a register value is repeated (or maxPeriod is reached)
 * @tparam k the length of the LFSR, or the order of the monic polynomial PRBS-k (last exponent)
 * @tparam nTaps the number of taps
 * @param start the start value of the LFSR
 * @param taps the taps that will be XOR-ed to calculate the new bit. They are the exponents of the monic polynomial. Ordering is unimportant.
 * @param wrapping if true, allow wrapping until maxPeriod or start is reached
 * @param verbose if true, print the lfsr and new bit on each step
 * @note Shift direction of the register is to the left <<, the newBit is set at bit position 0 (right)
 * @return the array of pseudo random bits, or an empty array if input was incorrect
 * @see https://en.wikipedia.org/wiki/Monic_polynomial
 * @see https://en.wikipedia.org/wiki/Linear-feedback_shift_register
 * @see https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence
 */
template <size_t k, size_t nTaps>
vector<bool>
LFSR(const bitset<k> start, const array<uint16_t, nTaps> taps, const bool wrapping = false, const bool verbose = false)
{
    vector<bool> result;
    //Sanity-checks
    if(taps.size()<2)
    {
        cerr << "At least two taps are needed" << endl;
        return result;
    }
    for(auto tap : taps)
    {
        if(tap > k)
        {
            cerr << "Out of range tap " << tap << endl;
            return result;
        }
    }
    if(start.none())
    {
        cerr << "A non-zero start value is needed" << endl;
        return result;
    }

    const uint32_t maxPeriod = pow(2,k) - 1;
    result.reserve(maxPeriod);

    set<uint32_t> lfsrHistory;
    bitset<k> lfsr(start);
    uint32_t i = 0;
    do
    {
        // XOR of all the tapped bits. We use -1 because an exponent N in the monic polynomial corresponds to bitset index N-1
        bool newBit = lfsr[taps.at(0) - 1];
        for(uint16_t j = 1; j < taps.size() ; ++j)
        {
            newBit ^= lfsr[taps.at(j) - 1];
        }

        lfsr <<= 1;
        lfsr[0] = newBit;
        if(verbose)
        {
            cout << i << " " << newBit << " " << lfsr << endl;
        }
        result.emplace_back(newBit);
        ++i;

        if(!wrapping)
        {
            if(lfsrHistory.count(lfsr.to_ulong()))
            {
                break;
            }
            lfsrHistory.insert(lfsr.to_ulong());
        }
    }
    while(lfsr != start && i < maxPeriod);
    result.shrink_to_fit();//only some special taps, namely the PRBS-k, will lead to the maxPeriod, others will stop earlier
    return result;
}

void test_prbs()
{
    //PRBS3
    array<uint16_t, 2> taps3 = {2, 3}; // Exponents of the monic polynomial
    auto prbs3 = LFSR(bitset<3>().flip(), taps3);// Start value all high

    //PRBS4
    array<uint16_t, 2> taps4 = {3, 4}; // Exponents of the monic polynomial
    auto prbs4 = LFSR(bitset<4>().flip(), taps4);// Start value all high

    //PRBS7
    array<uint16_t, 2> taps7 = {7, 6}; // Exponents of the monic polynomial
    auto prbs7 = LFSR(bitset<7>().flip(), taps7);// Start value all high

    for(auto prbs : {prbs3, prbs4, prbs7})
    {
        cout << prbs.size() << ":\t";
        for(auto p : prbs)
        {
            cout << p << " ";
        }
        cout << endl;
    }
    return 0;
}

[AdvaitDhingra]

Have you heard of the Bifurcation equation? It's quite old and I believe it was the first method of generating pseudo random numbers. Maybe this can be added as a feature by itself (i.e. multiple pseudo random number generator methods )

[ferdymercury]

Interesting. No, I hadn't heard about it. But I think PRBS is more about generating [0,1] bits for testing electronics, using easy/fast to implement LFSR, rather than generating floating point numbers.

[lmoneta]

Hi,

I am not expert of PRBS generators, but can you not use one of the current ROOT pseudo-random number generators and extract the bits from the generated numbers ?
If not, If you think is a useful addition, you can contribute with a PR which should include also a test and an example program (tutorial) showing its use

Cheers

Lorenzo

[ferdymercury]

The goal is the following. In electronics, the chips almost always generate PRBS sequences. When you acquire data in your system, you want to compare those bits, to see if any data was lost corrupted during the data acquisition chain. The goal is not to generate random sequences, for that, ROOT has already better generator tools. It is just to be able to compare the chip-generated ones with the expectation. And the ones in the chips are almost always of PRBS type.

Do you think the code snippet I wrote above would be ready for a PR ?
In what class or library would you suggest that I insert it?

[lmoneta]

Hi,
I understand it now. I would suggest you to make a new class that you can add to the MathCore library. The code snipped is probably a good starting point.
Thanks for this