SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
discrmath.h
Go to the documentation of this file.
1 // This file is part of LeanSDR Copyright (C) 2018 <pabr@pabr.org>.
2 // See the toplevel README for more information.
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program. If not, see <http://www.gnu.org/licenses/>.
16 
17 #ifndef LEANSDR_DISCRMATH_H
18 #define LEANSDR_DISCRMATH_H
19 
20 #include <cstddef>
21 
22 namespace leansdr
23 {
24 
25 // Polynomials with N binary coefficients.
26 // This is designed to compile into trivial inline code.
27 
28 // Bits are packed into words of type T.
29 // Unused most-significant bits of the last word are 0.
30 // T must be an unsigned type.
31 
32 template <typename T, int N>
33 struct bitvect
34 {
35  static const size_t WSIZE = sizeof(T) * 8;
36  static const size_t NW = (N + WSIZE - 1) / WSIZE;
37  T v[NW];
38 
39  bitvect() {}
40  bitvect(T val)
41  {
42  v[0] = val;
43  for (int i = 1; i < NW; ++i)
44  v[i] = 0;
45  }
46 
47  // Assign from another bitvect, with truncation or padding
48  template <int M>
50  {
51  int nw = (a.NW < NW) ? a.NW : NW;
52  for (int i = 0; i < nw; ++i)
53  v[i] = a.v[i];
54  if (M < N)
55  for (int i = a.NW; i < NW; ++i)
56  v[i] = 0;
57  if (M > N)
58  truncate_to_N();
59  return *this;
60  }
61 
62  bool operator[](unsigned int i) const
63  {
64  return v[i / WSIZE] & ((T)1 << (i & (WSIZE - 1)));
65  }
66 
67  // Add (in GF(2))
68 
69  template <int M>
71  {
72  int nw = (a.NW < NW) ? a.NW : NW;
73  for (int i = 0; i < nw; ++i)
74  v[i] ^= a.v[i];
75  if (M > N)
76  truncate_to_N();
77  return *this;
78  }
79 
80  // Multiply by X^n
81 
82  bitvect<T, N> &operator<<=(unsigned int n)
83  {
84  if (n >= WSIZE)
85  fail("shift exceeds word width");
86  T mask = ~(((T)-1) << n);
87  for (int i = NW - 1; i > 0; --i)
88  v[i] = (v[i] << n) | ((v[i - 1] >> (WSIZE - n)) & mask);
89  v[0] <<= n;
90  if (N != NW * WSIZE)
91  truncate_to_N();
92  return *this;
93  }
94 
95  private:
96  // Call this whenever bits N .. NW*WSIZE-1 may have been set.
97  inline void truncate_to_N()
98  {
99  v[NW - 1] &= ((T)-1) >> (NW * WSIZE - N);
100  }
101 
102 }; // bitvect
103 
104 // Return m modulo (X^N+p).
105 // p is typically a generator polynomial of degree N
106 // with the highest-coefficient monomial omitted.
107 // m is packed into nm words of type Tm.
108 // The MSB of m[0] is the highest-order coefficient of m.
109 
110 template <typename T, int N, typename Tm>
111 bitvect<T, N> divmod(const Tm *m, size_t nm, const bitvect<T, N> &p)
112 {
113  bitvect<T, N> res = 0;
114  const Tm bitmask = (Tm)1 << (sizeof(Tm) * 8 - 1);
115  for (; nm--; ++m)
116  {
117  Tm mi = *m;
118  for (int bit = sizeof(Tm) * 8; bit--; mi <<= 1)
119  {
120  // Multiply by X, save outgoing coeff of degree N
121  bool resN = res[N - 1];
122  res <<= 1;
123  // Add m[i]
124  if (mi & bitmask)
125  res.v[0] ^= 1;
126  // Modulo X^N+p
127  if (resN)
128  res += p;
129  }
130  }
131  return res;
132 }
133 
134 // Return (m*X^N) modulo (X^N+p).
135 // Same as divmod(), slightly faster than appending N zeroes.
136 
137 template <typename T, int N, typename Tm>
138 bitvect<T, N> shiftdivmod(const Tm *m, size_t nm, const bitvect<T, N> &p,
139  T init = 0)
140 {
141  bitvect<T, N> res;
142  for (int i = 0; i < res.NW; ++i)
143  res.v[i] = init;
144  const Tm bitmask = (Tm)1 << (sizeof(Tm) * 8 - 1);
145  for (; nm--; ++m)
146  {
147  Tm mi = *m;
148  for (int bit = sizeof(Tm) * 8; bit--; mi <<= 1)
149  {
150  // Multiply by X, save outgoing coeff of degree N
151  bool resN = res[N - 1];
152  res <<= 1;
153  // Add m[i]*X^N
154  resN ^= (bool)(mi & bitmask);
155  // Modulo X^N+p
156  if (resN)
157  res += p;
158  }
159  }
160  return res;
161 }
162 
163 template <typename T, int N>
164 bool operator==(const bitvect<T, N> &a, const bitvect<T, N> &b)
165 {
166  for (int i = 0; i < a.NW; ++i)
167  if (a.v[i] != b.v[i])
168  return false;
169  return true;
170 }
171 
172 // Add (in GF(2))
173 
174 template <typename T, int N>
176 {
177  bitvect<T, N> res;
178  for (int i = 0; i < a.NW; ++i)
179  res.v[i] = a.v[i] ^ b.v[i];
180  return res;
181 }
182 
183 // Polynomial multiplication.
184 
185 template <typename T, int N, int NB>
187 {
188  bitvect<T, N> res = 0;
189  for (int i = 0; i < NB; ++i, a <<= 1)
190  if (b[i])
191  res += a;
192  // TBD If truncation is needed, do it only once at the end.
193  return res;
194 }
195 
196 // Finite group GF(2^N), for small N.
197 
198 // GF(2) is the ring ({0,1},+,*).
199 // GF(2)[X] is the ring of polynomials with coefficients in GF(2).
200 // P(X) is an irreducible polynomial of GF(2)[X].
201 // N is the degree of P(x).
202 // GF(2)[X]/(P) is GF(2)[X] modulo P(X).
203 // (GF(2)[X]/(P), +) is a group with 2^N elements.
204 // (GF(2)[X]/(P)*, *) is a group with 2^N-1 elements.
205 // (GF(2)[X]/(P), +, *) is a field with 2^N elements, noted GF(2^N).
206 
207 // Te is a C++ integer type for encoding elements of GF(2^N).
208 // Binary coefficients are packed, with degree 0 at LSB.
209 // (Te)0 is 0
210 // (Te)1 is 1
211 // (Te)2 is X
212 // (Te)3 is X+1
213 // (Te)4 is X^2
214 // TRUNCP is the encoding of the generator, with highest monomial omitted.
215 // ALPHA is a primitive element of GF(2^N). Usually "2"=[X] is chosen.
216 
217 template <typename Te, int N, Te ALPHA, Te TRUNCP>
218 struct gf2n
219 {
220  typedef Te element;
221  static const Te alpha = ALPHA;
223  {
224  if (ALPHA != 2)
225  fail("alpha!=2 not implemented");
226  // Precompute log and exp tables.
227  Te alpha_i = 1; // ALPHA^0
228  for (int i = 0; i < (1 << N); ++i)
229  {
230  lut_exp[i] = alpha_i; // ALPHA^i
231  lut_exp[((1 << N) - 1) + i] = alpha_i; // Wrap to avoid modulo 2^N-1
232  lut_log[alpha_i] = i;
233  bool overflow = alpha_i & (1 << (N - 1));
234  alpha_i <<= 1; // Multiply by alpha=[X] i.e. increase degrees
235  alpha_i &= ~((~(Te)0) << N); // In case Te is wider than N bits
236  if (overflow)
237  alpha_i ^= TRUNCP; // Modulo P iteratively
238  }
239  }
240  inline Te add(Te x, Te y) { return x ^ y; } // Addition modulo 2
241  inline Te sub(Te x, Te y) { return x ^ y; } // Subtraction modulo 2
242  inline Te mul(Te x, Te y)
243  {
244  if (!x || !y)
245  return 0;
246  return lut_exp[lut_log[x] + lut_log[y]];
247  }
248  inline Te div(Te x, Te y)
249  {
250  if (!x)
251  return 0;
252  return lut_exp[lut_log[x] + ((1 << N) - 1) - lut_log[y]];
253  }
254  inline Te inv(Te x)
255  {
256  return lut_exp[((1 << N) - 1) - lut_log[x]];
257  }
258  inline Te exp(Te x) { return lut_exp[x]; }
259  inline Te log(Te x) { return lut_log[x]; }
260 
261  private:
262  Te lut_exp[(1 << N) * 2]; // Extra room for wrapping modulo 2^N-1
263  Te lut_log[1 << N];
264 }; // gf2n
265 
266 } // namespace leansdr
267 
268 #endif // LEANSDR_DISCRMATH_H
bitvect< T, N > operator*(bitvect< T, N > a, const bitvect< T, NB > &b)
Definition: discrmath.h:186
Te add(Te x, Te y)
Definition: discrmath.h:240
Te exp(Te x)
Definition: discrmath.h:258
bitvect< T, N > & operator<<=(unsigned int n)
Definition: discrmath.h:82
void truncate_to_N()
Definition: discrmath.h:97
Te inv(Te x)
Definition: discrmath.h:254
static const size_t WSIZE
Definition: discrmath.h:35
bool operator[](unsigned int i) const
Definition: discrmath.h:62
Te div(Te x, Te y)
Definition: discrmath.h:248
void fail(const char *s)
Definition: framework.cpp:11
bitvect< T, N > & copy(const bitvect< T, M > &a)
Definition: discrmath.h:49
static const size_t NW
Definition: discrmath.h:36
int32_t i
Definition: decimators.h:244
Te log(Te x)
Definition: discrmath.h:259
bitvect< T, N > divmod(const Tm *m, size_t nm, const bitvect< T, N > &p)
Definition: discrmath.h:111
bool operator==(const bitvect< T, N > &a, const bitvect< T, N > &b)
Definition: discrmath.h:164
bitvect< T, N > & operator+=(const bitvect< T, M > &a)
Definition: discrmath.h:70
Te mul(Te x, Te y)
Definition: discrmath.h:242
Te sub(Te x, Te y)
Definition: discrmath.h:241
bitvect< T, N > operator+(const bitvect< T, N > &a, const bitvect< T, N > &b)
Definition: discrmath.h:175
bitvect< T, N > shiftdivmod(const Tm *m, size_t nm, const bitvect< T, N > &p, T init=0)
Definition: discrmath.h:138
bitvect(T val)
Definition: discrmath.h:40