|
| 1 | +// @(#)root/mathcore:$Id$ |
| 2 | +// Author: Fernando Hueso-González 04/08/2021 |
| 3 | + |
| 4 | +#ifndef ROOT_Math_LFSR |
| 5 | +#define ROOT_Math_LFSR |
| 6 | + |
| 7 | +#include <array> |
| 8 | +#include <bitset> |
| 9 | +#include <vector> |
| 10 | +#include <set> |
| 11 | +#include <cmath> |
| 12 | +#include "TError.h" |
| 13 | + |
| 14 | +/// Pseudo Random Binary Sequence (PRBS) generator namespace with functions based |
| 15 | +/// on linear feedback shift registers (LFSR) with a periodicity of 2^n-1 |
| 16 | +/// |
| 17 | +/// @note It should NOT be used for general-purpose random number generation or any |
| 18 | +/// statistical study, for those cases see e.g. std::mt19937 instead. |
| 19 | +/// |
| 20 | +/// The goal is to generate binary bit sequences with the same algorithm as the ones usually implemented |
| 21 | +/// in electronic chips, so that the theoretically expected ones can be compared with the acquired sequences. |
| 22 | +/// |
| 23 | +/// The main ingredients of a PRBS generator are a monic polynomial of maximum degree \f$n\f$, with coefficients |
| 24 | +/// either 0 or 1, and a <a href="https://www.nayuki.io/page/galois-linear-feedback-shift-register">Galois</a> |
| 25 | +/// linear-feedback shift register with a non-zero seed. When the monic polynomial exponents are chosen appropriately, |
| 26 | +/// the period of the resulting bit sequence (0s and 1s) yields \f$2^n - 1\f$. |
| 27 | +/// |
| 28 | +/// @sa https://gist.github.com/mattbierner/d6d989bf26a7e54e7135, https://root.cern/doc/master/civetweb_8c_source.html#l06030, |
| 29 | +/// https://cryptography.fandom.com/wiki/Linear_feedback_shift_register, https://www3.advantest.com/documents/11348/33b24c8a-c8cb-40b8-a2a7-37515ba4abc8, |
| 30 | +/// https://www.reddit.com/r/askscience/comments/63a10q/for_prbs3_with_clock_input_on_each_gate_how_can/, https://es.mathworks.com/help/serdes/ref/prbs.html, |
| 31 | +/// https://metacpan.org/pod/Math::PRBS, https://ez.analog.com/data_converters/high-speed_adcs/f/q-a/545335/ad9689-pn9-and-pn23 |
| 32 | + |
| 33 | + |
| 34 | + |
| 35 | +namespace ROOT::Math::LFSR { |
| 36 | + |
| 37 | + /** |
| 38 | + * @brief Generate the next pseudo-random bit using the current state of a linear feedback shift register (LFSR) and update it |
| 39 | + * @tparam k the length of the LFSR, usually also the order of the monic polynomial PRBS-k (last exponent) |
| 40 | + * @tparam nTaps the number of taps |
| 41 | + * @param lfsr the current value of the LFSR. Passed by reference, it will be updated with the next value |
| 42 | + * @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. |
| 43 | + * @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) |
| 44 | + * @return the new random bit |
| 45 | + * @throw an exception is thrown if taps are out of the range [1, k] |
| 46 | + * @see https://en.wikipedia.org/wiki/Monic_polynomial, https://en.wikipedia.org/wiki/Linear-feedback_shift_register, https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence |
| 47 | + */ |
| 48 | +template <size_t k, size_t nTaps> |
| 49 | +static bool NextLFSR(std::bitset<k> &lfsr, std::array<uint16_t, nTaps> taps, bool left = true) |
| 50 | +{ |
| 51 | + static_assert(k <= 32, "For the moment, only supported until k == 32."); |
| 52 | + static_assert(k > 0, "Non-zero degree is needed for the LFSR."); |
| 53 | + static_assert(nTaps > 0, "At least one tap is needed for the LFSR."); |
| 54 | + static_assert(nTaps <= k, "Cannot use more taps than polynomial order"); |
| 55 | + for (uint16_t j = 0; j < nTaps; ++j) { |
| 56 | + static_assert((taps[j] - 1) <= k && (taps[j] - 1) > 0, "Tap value is out of range [1,k]"); |
| 57 | + } |
| 58 | + |
| 59 | + // First, calculate the XOR (^) of all selected bits (marked by the taps) |
| 60 | + bool newBit = lfsr[taps[0] - 1]; // the exponent E of the polynomial correspond to index E - 1 in the bitset |
| 61 | + for (uint16_t j = 1; j < nTaps; ++j) { |
| 62 | + newBit ^= lfsr[taps[j] - 1]; |
| 63 | + } |
| 64 | + |
| 65 | + // Apply the shift to the register in the right direction, and overwrite the empty one with newBit |
| 66 | + if (left) { |
| 67 | + lfsr <<= 1; |
| 68 | + lfsr[0] = newBit; |
| 69 | + } else { |
| 70 | + lfsr >>= 1; |
| 71 | + lfsr[k - 1] = newBit; |
| 72 | + } |
| 73 | + |
| 74 | + return newBit; |
| 75 | +} |
| 76 | + |
| 77 | +/** |
| 78 | + * @brief Generation of a sequence of pseudo-random bits using a linear feedback shift register (LFSR), until a |
| 79 | + * register value is repeated (or maxPeriod is reached) |
| 80 | + * @tparam k the length of the LFSR, usually also the order of the monic polynomial PRBS-k (last exponent) |
| 81 | + * @tparam nTaps the number of taps |
| 82 | + * @tparam Output the type of the container where the bit result (0 or 1) is stored (e.g. char, bool). It's unsigned |
| 83 | + * char by default, use bool instead if you want to save memory |
| 84 | + * @param start the start value (seed) of the LFSR |
| 85 | + * @param taps the taps that will be XOR-ed to calculate the new bit. They are the exponents of the monic polynomial. |
| 86 | + * Ordering is unimportant. Note that an exponent E in the polynom maps to bit index E-1 in the LFSR. |
| 87 | + * @param[out] iterator container storing the array of pseudo random bits (iterator is left untouched if input |
| 88 | + * parameters were incorrect) Pass for example an iterator to std::vector<unsigned char> (or to bool if memory saving |
| 89 | + * is important) |
| 90 | + * @param left if true, the direction of the register shift is to the left <<, the newBit is set on lfsr at bit |
| 91 | + * position 0 (right). If false, shift is to the right and the newBit is stored at bit position (k-1) |
| 92 | + * @param wrapping if true, allow repetition of values in the LFSRhistory, until maxPeriod is reached or the repeated |
| 93 | + * value == start. Enabling this option saves memory as no history is kept |
| 94 | + * @param oppositeBit if true, use the high/low bit of the LFSR to store output (for left=true/false, respectively) |
| 95 | + * instead of the newBit returned by ::NextLFSR |
| 96 | + * @return the array of pseudo random bits, or an empty array if input was incorrect |
| 97 | + * @see https://en.wikipedia.org/wiki/Monic_polynomial, https://en.wikipedia.org/wiki/Linear-feedback_shift_register, |
| 98 | + * https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence |
| 99 | + */ |
| 100 | +template <size_t k, size_t nTaps, typename Output = unsigned char> |
| 101 | +static std::vector<Output> GenerateSequence(std::bitset<k> start, std::array<uint16_t, nTaps> taps, bool left = true, |
| 102 | + bool wrapping = false, bool oppositeBit = false) |
| 103 | +{ |
| 104 | + std::vector<bool> result; // Store result here |
| 105 | + |
| 106 | + // Sanity-checks |
| 107 | + static_assert(k <= 32, "For the moment, only supported until k == 32."); |
| 108 | + static_assert(k > 0, "Non-zero degree is needed for the LFSR."); |
| 109 | + static_assert(nTaps >= 2, "At least two taps are needed for a proper sequence"); |
| 110 | + static_assert(nTaps <= k, "Cannot use more taps than polynomial order"); |
| 111 | + for (auto tap : taps) { |
| 112 | + if (tap > k || tap == 0) { |
| 113 | + Error("ROOT::Math::LFSR", "Tap %u is out of range [1,%lu]", tap, k); |
| 114 | + return result; |
| 115 | + } |
| 116 | + } |
| 117 | + if (start.none()) { |
| 118 | + Error("ROOT::Math::LFSR", "A non-zero start value is needed"); |
| 119 | + return result; |
| 120 | + } |
| 121 | + |
| 122 | + // Calculate maximum period and pre-allocate space in result |
| 123 | + const uint32_t maxPeriod = pow(2, k) - 1; |
| 124 | + result.reserve(maxPeriod); |
| 125 | + |
| 126 | + std::set<uint32_t> lfsrHistory; // a placeholder to store the history of all different values of the LFSR |
| 127 | + std::bitset<k> lfsr(start); // a variable storing the current value of the LFSR |
| 128 | + uint32_t i = 0; // a loop counter |
| 129 | + if (oppositeBit) // if oppositeBit enabled, first value is already started with the seed |
| 130 | + result.emplace_back(left ? lfsr[k - 1] : lfsr[0]); |
| 131 | + |
| 132 | + // Loop now until maxPeriod or a lfsr value is repeated. If wrapping enabled, allow repeated values if not equal |
| 133 | + // to seed |
| 134 | + do { |
| 135 | + bool newBit = NextLFSR(lfsr, taps, left); |
| 136 | + |
| 137 | + if (!oppositeBit) |
| 138 | + result.emplace_back(newBit); |
| 139 | + else |
| 140 | + result.emplace_back(left ? lfsr[k - 1] : lfsr[0]); |
| 141 | + |
| 142 | + ++i; |
| 143 | + |
| 144 | + if (!wrapping) // If wrapping not allowed, break the loop once a repeated value is encountered |
| 145 | + { |
| 146 | + if (lfsrHistory.count(lfsr.to_ulong())) |
| 147 | + break; |
| 148 | + |
| 149 | + lfsrHistory.insert(lfsr.to_ulong()); // Add to the history |
| 150 | + } |
| 151 | + } while (lfsr != start && i < maxPeriod); |
| 152 | + |
| 153 | + if (oppositeBit) |
| 154 | + result.pop_back(); // remove last element, as we already pushed the one from the seed above the while loop |
| 155 | + |
| 156 | + result.shrink_to_fit(); // only some special taps will lead to the maxPeriod, others will stop earlier |
| 157 | + |
| 158 | + return; |
| 159 | +} |
| 160 | +} |
| 161 | + |
| 162 | +#endif |
0 commit comments