|
| 1 | +// @(#)root/mathcore:$Id$ |
| 2 | +// Author: Fernando Hueso-González 04/08/2021 |
| 3 | + |
| 4 | +#ifndef ROOT_TRandomBinary |
| 5 | +#define ROOT_TRandomBinary |
| 6 | + |
| 7 | +////////////////////////////////////////////////////////////////////////// |
| 8 | +// // |
| 9 | +// TRandomBinary // |
| 10 | +// // |
| 11 | +// Pseudo Random Binary Sequence generator class (periodicity = 2**n-1) // |
| 12 | +// // |
| 13 | +////////////////////////////////////////////////////////////////////////// |
| 14 | + |
| 15 | +#include <array> |
| 16 | +#include <bitset> |
| 17 | +#include <vector> |
| 18 | +#include <set> |
| 19 | +#include <cmath> |
| 20 | +#include "Rtypes.h" |
| 21 | +#include "TError.h" |
| 22 | + |
| 23 | +class TRandomBinary { |
| 24 | +public: |
| 25 | + /** |
| 26 | + * @brief Generate the next pseudo-random bit using the current state of a linear feedback shift register (LFSR) and update it |
| 27 | + * @tparam k the length of the LFSR, usually also the order of the monic polynomial PRBS-k (last exponent) |
| 28 | + * @tparam nTaps the number of taps |
| 29 | + * @param lfsr the current value of the LFSR. Passed by reference, it will be updated with the next value |
| 30 | + * @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. Note that an exponent E in the polynom maps to bit index E-1 in the LFSR. |
| 31 | + * @param left if true, the direction of the register shift is to the left <<, the newBit is set on lfsr at bit position 0 (right). If false, shift is to the right and the newBit is stored at bit position (k-1) |
| 32 | + * @return the new random bit |
| 33 | + * @throw an exception is thrown if taps are out of the range [1, k] |
| 34 | + * @see https://en.wikipedia.org/wiki/Monic_polynomial |
| 35 | + * @see https://en.wikipedia.org/wiki/Linear-feedback_shift_register |
| 36 | + * @see https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence |
| 37 | + */ |
| 38 | + template <size_t k, size_t nTaps> |
| 39 | + static bool |
| 40 | + NextLFSR(std::bitset<k>& lfsr, const std::array<uint16_t, nTaps> taps, const bool left = true) |
| 41 | + { |
| 42 | + static_assert(k <= 32, "For the moment, only supported until k == 32."); |
| 43 | + static_assert(k > 0, "Non-zero degree is needed for the LFSR."); |
| 44 | + static_assert(nTaps > 0, "At least one tap is needed for the LFSR."); |
| 45 | + static_assert(nTaps <= k, "Cannot use more taps than polynomial order"); |
| 46 | + |
| 47 | + // First, calculate the XOR (^) of all selected bits (marked by the taps) |
| 48 | + bool newBit = lfsr[taps.at(0) - 1]; // the exponent E of the polynomial correspond to index E - 1 in the bitset |
| 49 | + for(uint16_t j = 1; j < nTaps ; ++j) |
| 50 | + { |
| 51 | + newBit ^= lfsr[taps.at(j) - 1]; |
| 52 | + } |
| 53 | + |
| 54 | + //Apply the shift to the register in the right direction, and overwrite the empty one with newBit |
| 55 | + if(left) |
| 56 | + { |
| 57 | + lfsr <<= 1; |
| 58 | + lfsr[0] = newBit; |
| 59 | + } |
| 60 | + else |
| 61 | + { |
| 62 | + lfsr >>= 1; |
| 63 | + lfsr[k-1] = newBit; |
| 64 | + } |
| 65 | + |
| 66 | + return newBit; |
| 67 | + } |
| 68 | + |
| 69 | + /** |
| 70 | + * @brief Generation of a sequence of pseudo-random bits using a linear feedback shift register (LFSR), until a register value is repeated (or maxPeriod is reached) |
| 71 | + * @tparam k the length of the LFSR, usually also the order of the monic polynomial PRBS-k (last exponent) |
| 72 | + * @tparam nTaps the number of taps |
| 73 | + * @param start the start value (seed) of the LFSR |
| 74 | + * @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. Note that an exponent E in the polynom maps to bit index E-1 in the LFSR. |
| 75 | + * @param left if true, the direction of the register shift is to the left <<, the newBit is set on lfsr at bit position 0 (right). If false, shift is to the right and the newBit is stored at bit position (k-1) |
| 76 | + * @param wrapping if true, allow repetition of values in the LFSRhistory, until maxPeriod is reached or the repeated value == start. Enabling this option saves memory as no history is kept |
| 77 | + * @param oppositeBit if true, use the high/low bit of the LFSR to store output (for left=true/false, respectively) instead of the newBit returned by ::NextLFSR |
| 78 | + * @return the array of pseudo random bits, or an empty array if input was incorrect |
| 79 | + * @see https://en.wikipedia.org/wiki/Monic_polynomial |
| 80 | + * @see https://en.wikipedia.org/wiki/Linear-feedback_shift_register |
| 81 | + * @see https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence |
| 82 | + */ |
| 83 | + template <size_t k, size_t nTaps> |
| 84 | + static std::vector<bool> |
| 85 | + GenerateSequence(const std::bitset<k> start, const std::array<uint16_t, nTaps> taps, const bool left = true, const bool wrapping = false, const bool oppositeBit = false) |
| 86 | + { |
| 87 | + std::vector<bool> result; // Store result here |
| 88 | + |
| 89 | + //Sanity-checks |
| 90 | + static_assert(k <= 32, "For the moment, only supported until k == 32."); |
| 91 | + static_assert(k > 0, "Non-zero degree is needed for the LFSR."); |
| 92 | + static_assert(nTaps >= 2, "At least two taps are needed for a proper sequence"); |
| 93 | + static_assert(nTaps <= k, "Cannot use more taps than polynomial order"); |
| 94 | + for(auto tap : taps) { |
| 95 | + if(tap > k || tap == 0) { |
| 96 | + Error("TRandomBinary", "Tap %u is out of range [1,%lu]", tap, k); |
| 97 | + return result; |
| 98 | + } |
| 99 | + } |
| 100 | + if(start.none()) { |
| 101 | + Error("TRandomBinary", "A non-zero start value is needed"); |
| 102 | + return result; |
| 103 | + } |
| 104 | + |
| 105 | + // Calculate maximum period and pre-allocate space in result |
| 106 | + const uint32_t maxPeriod = pow(2,k) - 1; |
| 107 | + result.reserve(maxPeriod); |
| 108 | + |
| 109 | + std::set<uint32_t> lfsrHistory; // a placeholder to store the history of all different values of the LFSR |
| 110 | + std::bitset<k> lfsr(start); // a variable storing the current value of the LFSR |
| 111 | + uint32_t i = 0; // a loop counter |
| 112 | + if(oppositeBit) // if oppositeBit enabled, first value is already started with the seed |
| 113 | + result.emplace_back(left ? lfsr[k-1] : lfsr[0]); |
| 114 | + |
| 115 | + //Loop now until maxPeriod or a lfsr value is repeated. If wrapping enabled, allow repeated values if not equal to seed |
| 116 | + do { |
| 117 | + bool newBit = NextLFSR(lfsr, taps, left); |
| 118 | + |
| 119 | + if(!oppositeBit) |
| 120 | + result.emplace_back(newBit); |
| 121 | + else |
| 122 | + result.emplace_back(left ? lfsr[k-1] : lfsr[0]); |
| 123 | + |
| 124 | + ++i; |
| 125 | + |
| 126 | + if(!wrapping) // If wrapping not allowed, break the loop once a repeated value is encountered |
| 127 | + { |
| 128 | + if(lfsrHistory.count(lfsr.to_ulong())) |
| 129 | + break; |
| 130 | + |
| 131 | + lfsrHistory.insert(lfsr.to_ulong()); // Add to the history |
| 132 | + } |
| 133 | + } |
| 134 | + while(lfsr != start && i < maxPeriod); |
| 135 | + |
| 136 | + if(oppositeBit) |
| 137 | + result.pop_back();// remove last element, as we already pushed the one from the seed above the while loop |
| 138 | + |
| 139 | + result.shrink_to_fit();//only some special taps will lead to the maxPeriod, others will stop earlier |
| 140 | + |
| 141 | + return result; |
| 142 | + } |
| 143 | + |
| 144 | + TRandomBinary() = default; |
| 145 | + virtual ~TRandomBinary() = default; |
| 146 | + |
| 147 | + ClassDef(TRandomBinary,0) //Pseudo Random Binary Sequence (periodicity = 2**n - 1) |
| 148 | +}; |
| 149 | + |
| 150 | +#endif |
0 commit comments