GNU Radio's DVBS2RX Package
gf.h
Go to the documentation of this file.
1/* -*- c++ -*- */
2/*
3 * Copyright (c) 2022 Igor Freire.
4 *
5 * This file is part of gr-dvbs2rx.
6 *
7 * SPDX-License-Identifier: GPL-3.0-or-later
8 */
9
10#ifndef INCLUDED_DVBS2RX_GF_H
11#define INCLUDED_DVBS2RX_GF_H
12
14#include <unordered_map>
15#include <bitset>
16#include <cstdint>
17#include <limits>
18#include <set>
19#include <stdexcept>
20#include <vector>
21
22namespace gr {
23namespace dvbs2rx {
24
25// Non-int types for storing GF(2) coefficients or GF(2^m) elements:
26typedef std::bitset<256> bitset256_t;
27
28template <typename T>
30
31template <typename T>
33
34/**
35 * @brief Test if bit is set
36 *
37 * @param x Bit register.
38 * @param i_bit Target bit index.
39 * @return true if bit is 1 and false otherwise.
40 */
41template <typename T>
42inline bool is_bit_set(const T& x, int i_bit)
43{
44 return x & (static_cast<T>(1) << i_bit);
45}
46
47/**
48 * @overload
49 * @note Template specialization for T = bitset256_t.
50 */
51template <>
52inline bool is_bit_set(const bitset256_t& x, int i_bit)
53{
54 return x.test(i_bit);
55}
56
57/**
58 * @brief Galois Field GF(2^m).
59 *
60 * @tparam T Base type for the field elements.
61 * @note See the reference implementation at
62 * https://github.com/igorauad/bch/blob/master/gf.py.
63 */
64template <typename T>
66{
67private:
68 const uint8_t m_m; // dimension of the GF(2^m) field
69 const uint32_t m_two_to_m_minus_one; // shortcut for (2^m - 1)
70 std::vector<T> m_table; // field elements
71 std::vector<T> m_table_nonzero; // non-zero field elements
72 std::unordered_map<T, uint32_t> m_exp_table; // map element alpha^i to its exponent i
73
74public:
75 /**
76 * @brief Construct a new Galois field object.
77 *
78 * @param prim_poly Primitive polynomial that generates the field.
79 */
80 galois_field(const gf2_poly<T>& prim_poly);
81
82 /**
83 * @brief Get the GF(2^m) element at a given index on the elements table.
84 *
85 * @param index Target index.
86 * @return T GF(2^m) element.
87 */
88 T operator[](uint32_t index) const { return m_table[index]; }
89
90 /**
91 * @brief Get the i-th power of the primitive element (alpha^i).
92 *
93 * @param i Exponent i of the target element alpha^i for i from 0 to "2**m - 2".
94 * @return T Element alpha^i.
95 * @note This function cannot obtain the zero (additive identity) element, given the
96 * zero element cannot be expressed as a power of the primitive element. To access the
97 * zero element, use the operator[] instead.
98 */
99 T get_alpha_i(uint32_t i) const { return m_table_nonzero[i % m_two_to_m_minus_one]; }
100
101 /**
102 * @brief Get the exponent i of a given element beta = alpha^i.
103 *
104 * @param beta Element beta that is a power of the primitive element alpha.
105 * @return T Exponent i.
106 * @note This function cannot obtain the exponent of the zero (additive identity)
107 * element, given the zero element cannot be expressed as a power of the primitive
108 * element alpha. A runtime error exception is raised if beta is the zero element.
109 */
110 uint32_t get_exponent(const T& beta) const { return m_exp_table.at(beta); }
111
112 /**
113 * @brief Multiply two elements from GF(2^m).
114 *
115 * @param a First multiplicand.
116 * @param b Second multiplicand.
117 * @return T Product a*b.
118 */
119 T multiply(const T& a, const T& b) const;
120
121 /**
122 * @brief Get the inverse beta^-1 from a GF(2^m) element beta.
123 *
124 * @param beta Element to invert.
125 * @return T Inverse beta^-1.
126 */
127 T inverse(const T& beta) const;
128
129 /**
130 * @brief Get the inverse from a GF(2^m) element alpha^i given by its exponent i.
131 *
132 * @param i Exponent i of the element beta = alpha^i to be inverted.
133 * @return T Inverse beta^-1.
134 */
135 T inverse_by_exp(uint32_t i) const;
136
137 /**
138 * @brief Divide two elements from GF(2^m).
139 *
140 * @param a Dividend.
141 * @param b Divisor.
142 * @return T Quotient a/b.
143 */
144 T divide(const T& a, const T& b) const;
145
146 /**
147 * @brief Get the conjugates of element beta.
148 *
149 * @param beta The element whose conjugates are to be computed.
150 * @return std::set<T> The set of distinct conjugates alpha^(i^(2^l)) in GF(2^m)
151 * associated with element beta=alpha^i.
152 */
153 std::set<T> get_conjugates(const T& beta) const;
154
155 /**
156 * @brief Compute the minimal polynomial associated with element beta.
157 *
158 * Computes the polynomial phi(x) over GF(2) of smallest degree having beta as root.
159 * See the notes in the reference Python implementation.
160 *
161 * @param beta GF(2^m) element.
162 * @return gf2_poly Minimal polynomial of beta as a polynomial over GF(2).
163 */
164 gf2_poly<T> get_min_poly(const T& beta) const;
165
166 /**
167 * @brief Get the dimension m of the GF(2^m) field.
168 *
169 * @return uint8_t Dimension m.
170 */
171 uint8_t get_m() const { return m_m; };
172};
173
174/**
175 * @brief Get the maximum degree a GF2 polynomial can have with type T.
176 *
177 * @return int Maximum degree.
178 */
179template <typename T>
180inline constexpr size_t get_max_gf2_poly_degree()
181{
182 return (sizeof(T) * 8 - 1);
183}
184
185/**
186 * @overload for T = bitset256_t.
187 */
188template <>
190{
191 return bitset256_t().size() - 1;
192}
193
194/**
195 * @brief Polynomial over GF(2).
196 *
197 * A polynomial whose coefficients are elements from GF(2), i.e., binary.
198 *
199 * @tparam T Type whose bits represent the binary polynomial coefficients.
200 */
201template <typename T>
203{
204private:
205 static constexpr int m_max_degree = get_max_gf2_poly_degree<T>();
206
207 T m_poly; // Polynomial coefficients
208 // NOTE: the LSB (bit 0) has the zero-degree coefficient, bit 1 has the coefficient of
209 // x, bit 2 of x^2, and so on. For instance, a polynomial "x^4 + x + 1" should be
210 // stored as "10011"/
211 int m_degree; // Polynomial degree
212 // NOTE: by convention, the zero polynomial has degree -1.
213
214public:
215 /**
216 * @brief Construct a new GF(2) polynomial object.
217 *
218 * @param coefs Binary coefficients.
219 */
220 gf2_poly(const T& coefs);
221
222 /**
223 * @brief GF(2) polynomial addition.
224 *
225 * @param x Polynomial to add.
226 * @return gf2_poly Addition result.
227 */
228 gf2_poly<T> operator+(const gf2_poly<T>& x) const { return m_poly ^ x.get_poly(); }
229
230 /**
231 * @brief Multiplication by a GF(2) scalar.
232 *
233 * @param x Binary scalar to multiply.
234 * @return gf2_poly Multiplication result.
235 */
236 gf2_poly<T> operator*(bool x) const { return x ? m_poly : 0; }
237
238 /**
239 * @brief Multiplication by another GF(2) polynomial.
240 *
241 * @param x Polynomial to multiply.
242 * @return gf2_poly Multiplication result.
243 */
245
246 /**
247 * @brief Equal comparator.
248 *
249 * @param x The other GF(2^m) polynomial.
250 * @return bool Whether they are equal.
251 */
252 bool operator==(const gf2_poly<T>& x) const { return m_poly == x.get_poly(); }
253
254 /**
255 * @brief Get the polynomial coefficients.
256 *
257 * @return T storage of polynomial coefficients.
258 */
259 const T& get_poly() const { return m_poly; }
260
261 /**
262 * @brief Get the polynomial degree.
263 *
264 * @return int Degree.
265 * @note By convention, the zero polynomial has degree -1.
266 */
267 int degree() const { return m_degree; }
268
269 /**
270 * @brief Test if the polynomial is the zero polynomial.
271 *
272 * @return true if the polynomial is zero and false otherwise.
273 */
274 bool is_zero() const { return m_degree == -1; }
275};
276
277/**
278 * @brief Compute the remainder of the division between two GF(2) polynomials.
279 *
280 * Computes the remainder of polynomial over GF(2) a(x) divided by another polynomial over
281 * GF(2) b(x), i.e., computes a(x) % b(x). The result has degree up to the degree of b(x)
282 * minus one. Hence, the result necessarily fits within the type used to store b(x).
283 *
284 * @tparam Ta Type of the dividend polynomial.
285 * @tparam Tb Type of the divisor polynomial.
286 * @param a Dividend polynomial.
287 * @param b Divisor polynomial.
288 * @return gf2_poly<Tb> Remainder result.
289 * @throws std::runtime_error if b(x) is a zero polynomial.
290 * @note The implementation requires the divisor type Tb to be larger than or equal to the
291 * dividend type Ta. Otherwise, a static assertion is raised.
292 */
293template <typename Ta, typename Tb>
295{
296 check_rem_types(a, b); // ensure the types are
297 if (b.degree() == -1) // zero divisor
298 throw std::runtime_error("Remainder of division by a zero polynomial");
299 if (a.degree() == -1) // zero dividend
300 return gf2_poly<Tb>(0);
301 if (a.degree() < b.degree()) // remainder is the dividend polynomial a(x) itself
302 return gf2_poly<Tb>(a.get_poly());
303
304 int b_degree = b.degree();
305 const Tb b_coefs = b.get_poly();
306 Tb remainder = a.get_poly(); // here type Tb must be large enough to store a(x)
307 for (int i = a.degree(); i >= b_degree; i--) {
308 if (is_bit_set(remainder, i))
309 remainder ^= b_coefs << (i - b_degree);
310 }
311 return remainder;
312}
313
314template <typename Ta, typename Tb>
315inline void check_rem_types(const gf2_poly<Ta> a, const gf2_poly<Tb> b)
316{
317 static_assert(sizeof(Tb) >= sizeof(Ta),
318 "Divisor type must be larger or equal than dividend type");
319}
320
321template <typename Tb>
323{
324 static_assert(sizeof(Tb) * 8 >= bitset256_t().size(),
325 "Divisor type must be larger or equal than dividend type");
326}
327
328template <typename Ta>
330{
331 static_assert(bitset256_t().size() >= sizeof(Ta) * 8,
332 "Divisor type must be larger or equal than dividend type");
333}
334
336{
337 static_assert(true); // the two types are the same
338}
339
340/**
341 * @brief Polynomial over GF(2^m).
342 *
343 * A polynomial whose coefficients are elements from a GF(2^m) extension field.
344 *
345 * @tparam T Type used to represent each polynomial coefficient as a GF(2^m) element.
346 */
347template <typename T>
349{
350private:
351 const galois_field<T>* m_gf; // Galois field
352 std::vector<T> m_poly; // Polynomial coefficients
353 // NOTE: index 0 has the zero-degree coefficient, index 1 has the coefficient of x,
354 // index 2 of x^2, and so on.
355 std::vector<uint32_t> m_nonzero_coef_idx; // Indexes (degrees) of the non-zero
356 // polynomial coefficients
357 std::vector<uint32_t> m_nonzero_coef_exp; // Exponents of the non-zero polynomial
358 // coefficients
359 size_t m_n_nonzero_coef; // Number of non-zero coefficients
360 // Example: a polynomial "alpha^4 x^3 + alpha^2 x + 1" would have the following data
361 // in the vectors: m_poly = [1, alpha^2, 0, alpha^4], m_nonzero_coef_idx = [0, 1, 3],
362 // m_nonzero_coef_exp = [0, 2, 4].
363 int m_degree; // Polynomial degree
364
365 /**
366 * @brief Set the polynomial degree.
367 */
368 void set_degree();
369
370 /**
371 * @brief Fill the indexes and exponents of the non-zero polynomial coefficients.
372 *
373 * Each coefficient can be represented as a GF(2^m) element alpha^j, where j is the
374 * exponent and alpha is the primitive element. To speed up computations, this
375 * function caches the exponents associated with the non-zero coefficients on a
376 * dedicated vector. Also, it caches the degrees of the terms associated with such
377 * non-zero coefficients.
378 */
379 void set_coef_exponents();
380
381public:
382 /**
383 * @brief Construct a new polynomial over GF(2^m).
384 *
385 * @param gf Reference Galois field.
386 * @param coefs Polynomial coefficients.
387 */
388 gf2m_poly(const galois_field<T>* const gf, std::vector<T>&& coefs);
389
390 /**
391 * @brief Construct a polynomial over GF(2^m) from a polynomial over GF(2).
392 *
393 * @param gf Reference Galois field.
394 * @param gf2_poly Reference polynomial over GF(2).
395 */
396 template <typename Y>
397 gf2m_poly(const galois_field<T>* const gf, const gf2_poly<Y>& gf2_poly) : m_gf(gf)
398 {
399 const Y& coefs = gf2_poly.get_poly();
400 for (int i = 0; i <= gf2_poly.degree(); i++) {
401 if (is_bit_set(coefs, i)) {
402 m_poly.push_back((*gf)[1]);
403 } else {
404 m_poly.push_back((*gf)[0]);
405 }
406 }
407 set_degree();
408 set_coef_exponents();
409 };
410
411 /**
412 * @brief GF(2^m) polynomial addition.
413 *
414 * @param x Polynomial to add.
415 * @return gf2m_poly<T> Addition result.
416 */
418
419 /**
420 * @brief Multiplication by a scalar.
421 *
422 * @param x Scalar to multiply.
423 * @return gf2m_poly<T> Multiplication result.
424 */
426
427 /**
428 * @brief Multiplication by another GF(2^m) polynomial.
429 *
430 * @param x Polynomial to multiply.
431 * @return gf2m_poly<T> Multiplication result.
432 */
434
435 /**
436 * @brief Equal comparator.
437 *
438 * @param x The other GF(2^m) polynomial.
439 * @return bool Whether they are equal.
440 */
441 bool operator==(const gf2m_poly<T>& x) const { return m_poly == x.get_poly(); }
442
443 /**
444 * @brief Access a polynomial coefficient.
445 *
446 * @param index Index of the target coefficient.
447 * @return T Polynomial coefficient.
448 */
449 T operator[](uint32_t index) const { return m_poly[index]; }
450
451 /**
452 * @brief Evaluate the polynomial for a given GF(2^m) element.
453 *
454 * Assuming the underlying polynomial is p(x), this function evaluates p(x) for a
455 * given x from GF(2^m).
456 *
457 * @param x GF(2^m) value for which the polynomial should be evaluated.
458 * @return T Evaluation result within GF(2^m).
459 */
460 T eval(T x) const;
461
462 /**
463 * @brief Evaluate the polynomial for an element alpha^i given by its exponent i.
464 *
465 * Assuming the underlying polynomial is p(x), this function evaluates p(x) for a
466 * given x from GF(2^m). The difference to the operator() is that this function takes
467 * x by its exponent i. If x = alpha^i, this function takes i as input instead of x.
468 *
469 * @param i Exponent of the GF(2^m) value for which the polynomial should be
470 * evaluated.
471 * @return T Evaluation result within GF(2^m).
472 * @note Since the zero element (additive identity) cannot be represented as a power
473 * of the primitive element, there is no exponent i for the zero element. Hence, this
474 * function can only evaluate the polynomial for non-zero elements.
475 */
476 T eval_by_exp(uint32_t i) const;
477
478 /**
479 * @brief Search roots by evaluating the polynomial for elements in a given range.
480 *
481 * Evaluates the polynomial for all alpha^i with i varying in [i_start, i_end] and
482 * returns the elements for which the polynomial evaluates to zero (the roots) given
483 * by their exponents. For instance, if alpha^2 and alpha^3 are roots of the
484 * polynomial, this function returns the vector [2, 3].
485 *
486 * The implementation does not simply rely on the eval() or eval_by_exp() functions.
487 * Instead, it is optimized by leveraging the fact that the evaluation is for a
488 * contiguous range of exponents. Hence, when searching for polynomial roots in
489 * GF(2^m), it is preferable to use this function instead of manually calling the
490 * eval() or eval_by_exp() functions.
491 *
492 * @param i_start Exponent of element alpha^i_start at the start of the range.
493 * @param i_end Exponent of element alpha^i_end at the end of the range.
494 * @param max_roots Maximum number of roots to be returned. When defined, the search
495 * is stopped earlier as soon as this number of roots is found.
496 * @return std::vector<uint32_t> Vector with the exponents associated with the GF(2^m)
497 * roots found in the range.
498 */
499 std::vector<uint32_t> search_roots_in_exp_range(
500 uint32_t i_start,
501 uint32_t i_end,
502 uint32_t max_roots = std::numeric_limits<uint32_t>::max()) const;
503
504 /**
505 * @brief Get the polynomial coefficients.
506 *
507 * @return const std::vector<T>& Reference to vector of polynomial coefficients.
508 */
509 const std::vector<T>& get_poly() const { return m_poly; }
510
511 /**
512 * @brief Get the polynomial degree.
513 *
514 * @return int Degree.
515 * @note By convention, the zero polynomial has degree -1.
516 */
517 int degree() const { return m_degree; }
518
519 /**
520 * @brief Convert the polynomial to a GF(2) polynomial.
521 *
522 * Works when all coefficients of the local polynomial are either unit or zero such
523 * that it can be reduced to a polynomial over GF(2). Otherwise, throws runtime error.
524 *
525 * @return gf2_poly<T> Polynomial over GF(2).
526 */
528};
529
530// Type definitions
535
536} // namespace dvbs2rx
537} // namespace gr
538
539#endif // INCLUDED_DVBS2RX_GF_H
Galois Field GF(2^m).
Definition gf.h:66
gf2_poly< T > get_min_poly(const T &beta) const
Compute the minimal polynomial associated with element beta.
uint8_t get_m() const
Get the dimension m of the GF(2^m) field.
Definition gf.h:171
T inverse(const T &beta) const
Get the inverse beta^-1 from a GF(2^m) element beta.
T divide(const T &a, const T &b) const
Divide two elements from GF(2^m).
T get_alpha_i(uint32_t i) const
Get the i-th power of the primitive element (alpha^i).
Definition gf.h:99
T multiply(const T &a, const T &b) const
Multiply two elements from GF(2^m).
std::set< T > get_conjugates(const T &beta) const
Get the conjugates of element beta.
T inverse_by_exp(uint32_t i) const
Get the inverse from a GF(2^m) element alpha^i given by its exponent i.
T operator[](uint32_t index) const
Get the GF(2^m) element at a given index on the elements table.
Definition gf.h:88
galois_field(const gf2_poly< T > &prim_poly)
Construct a new Galois field object.
uint32_t get_exponent(const T &beta) const
Get the exponent i of a given element beta = alpha^i.
Definition gf.h:110
Polynomial over GF(2).
Definition gf.h:203
gf2_poly< T > operator*(bool x) const
Multiplication by a GF(2) scalar.
Definition gf.h:236
gf2_poly< T > operator*(const gf2_poly< T > &x) const
Multiplication by another GF(2) polynomial.
int degree() const
Get the polynomial degree.
Definition gf.h:267
bool operator==(const gf2_poly< T > &x) const
Equal comparator.
Definition gf.h:252
const T & get_poly() const
Get the polynomial coefficients.
Definition gf.h:259
gf2_poly(const T &coefs)
Construct a new GF(2) polynomial object.
gf2_poly< T > operator+(const gf2_poly< T > &x) const
GF(2) polynomial addition.
Definition gf.h:228
bool is_zero() const
Test if the polynomial is the zero polynomial.
Definition gf.h:274
Polynomial over GF(2^m).
Definition gf.h:349
T eval_by_exp(uint32_t i) const
Evaluate the polynomial for an element alpha^i given by its exponent i.
const std::vector< T > & get_poly() const
Get the polynomial coefficients.
Definition gf.h:509
int degree() const
Get the polynomial degree.
Definition gf.h:517
gf2m_poly< T > operator*(const gf2m_poly< T > &x) const
Multiplication by another GF(2^m) polynomial.
gf2m_poly< T > operator*(T x) const
Multiplication by a scalar.
gf2m_poly(const galois_field< T > *const gf, std::vector< T > &&coefs)
Construct a new polynomial over GF(2^m).
gf2_poly< T > to_gf2_poly() const
Convert the polynomial to a GF(2) polynomial.
T eval(T x) const
Evaluate the polynomial for a given GF(2^m) element.
gf2m_poly(const galois_field< T > *const gf, const gf2_poly< Y > &gf2_poly)
Construct a polynomial over GF(2^m) from a polynomial over GF(2).
Definition gf.h:397
T operator[](uint32_t index) const
Access a polynomial coefficient.
Definition gf.h:449
gf2m_poly< T > operator+(const gf2m_poly< T > &x) const
GF(2^m) polynomial addition.
bool operator==(const gf2m_poly< T > &x) const
Equal comparator.
Definition gf.h:441
std::vector< uint32_t > search_roots_in_exp_range(uint32_t i_start, uint32_t i_end, uint32_t max_roots=std::numeric_limits< uint32_t >::max()) const
Search roots by evaluating the polynomial for elements in a given range.
#define DVBS2RX_API
Definition include/gnuradio/dvbs2rx/api.h:19
class DVBS2RX_API gf2_poly
Definition gf.h:29
std::bitset< 256 > bitset256_t
Definition gf.h:26
gf2_poly< Tb > operator%(const gf2_poly< Ta > a, const gf2_poly< Tb > b)
Compute the remainder of the division between two GF(2) polynomials.
Definition gf.h:294
bool is_bit_set(const T &x, int i_bit)
Test if bit is set.
Definition gf.h:42
void check_rem_types(const gf2_poly< Ta > a, const gf2_poly< Tb > b)
Definition gf.h:315
gf2_poly< uint64_t > gf2_poly_u64
Definition gf.h:533
constexpr size_t get_max_gf2_poly_degree()
Get the maximum degree a GF2 polynomial can have with type T.
Definition gf.h:180
constexpr size_t get_max_gf2_poly_degree< bitset256_t >()
Definition gf.h:189
gf2_poly< bitset256_t > gf2_poly_b256
Definition gf.h:534
gf2_poly< uint32_t > gf2_poly_u32
Definition gf.h:532
gf2_poly< uint16_t > gf2_poly_u16
Definition gf.h:531
Fixed-length double-ended queue with contiguous volk-aligned elements.
Definition gr_bch.h:22