1 #include "fft.h"
2 #include "gf.h"
3 #include "parameters.h"
4 #include "parsing.h"
5 #include "reed_solomon.h"
6 #include <stdint.h>
7 #include <stdio.h>
8 #include <string.h>
9 /**
10 * @file reed_solomon.c
11 * Constant time implementation of Reed-Solomon codes
12 */
13
14
15 static void compute_syndromes(uint16_t *syndromes, uint8_t *cdw);
16 static uint16_t compute_elp(uint16_t *sigma, const uint16_t *syndromes);
17 static void compute_roots(uint8_t *error, uint16_t *sigma);
18 static void compute_z_poly(uint16_t *z, const uint16_t *sigma, uint16_t degree, const uint16_t *syndromes);
19 static void compute_error_values(uint16_t *error_values, const uint16_t *z, const uint8_t *error);
20 static void correct_errors(uint8_t *cdw, const uint16_t *error_values);
21
22 /**
23 * @brief Encodes a message message of PARAM_K bits to a Reed-Solomon codeword codeword of PARAM_N1 bytes
24 *
25 * Following @cite lin1983error (Chapter 4 - Cyclic Codes),
26 * We perform a systematic encoding using a linear (PARAM_N1 - PARAM_K)-stage shift register
27 * with feedback connections based on the generator polynomial PARAM_RS_POLY of the Reed-Solomon code.
28 *
29 * @param[out] cdw Array of size VEC_N1_SIZE_64 receiving the encoded message
30 * @param[in] msg Array of size VEC_K_SIZE_64 storing the message
31 */
PQCLEAN_HQCRMRS192_CLEAN_reed_solomon_encode(uint8_t * cdw,const uint8_t * msg)32 void PQCLEAN_HQCRMRS192_CLEAN_reed_solomon_encode(uint8_t *cdw, const uint8_t *msg) {
33 size_t i, j, k;
34 uint8_t gate_value = 0;
35
36 uint16_t tmp[PARAM_G] = {0};
37 uint16_t PARAM_RS_POLY [] = {RS_POLY_COEFS};
38 uint8_t prev, x;
39
40 for (i = 0; i < PARAM_N1; ++i) {
41 cdw[i] = 0;
42 }
43
44 for (i = 0; i < PARAM_K; ++i) {
45 gate_value = (uint8_t) (msg[PARAM_K - 1 - i] ^ cdw[PARAM_N1 - PARAM_K - 1]);
46
47 for (j = 0; j < PARAM_G; ++j) {
48 tmp[j] = PQCLEAN_HQCRMRS192_CLEAN_gf_mul(gate_value, PARAM_RS_POLY[j]);
49 }
50
51 prev = 0;
52 for (k = 0; k < PARAM_N1 - PARAM_K; k++) {
53 x = cdw[k];
54 cdw[k] = (uint8_t) (prev ^ tmp[k]);
55 prev = x;
56 }
57 }
58
59 memcpy(cdw + PARAM_N1 - PARAM_K, msg, PARAM_K);
60 }
61
62
63
64 /**
65 * @brief Computes 2 * PARAM_DELTA syndromes
66 *
67 * @param[out] syndromes Array of size 2 * PARAM_DELTA receiving the computed syndromes
68 * @param[in] cdw Array of size PARAM_N1 storing the received vector
69 */
compute_syndromes(uint16_t * syndromes,uint8_t * cdw)70 void compute_syndromes(uint16_t *syndromes, uint8_t *cdw) {
71 for (size_t i = 0; i < 2 * PARAM_DELTA; ++i) {
72 for (size_t j = 1; j < PARAM_N1; ++j) {
73 syndromes[i] ^= PQCLEAN_HQCRMRS192_CLEAN_gf_mul(cdw[j], alpha_ij_pow[i][j - 1]);
74 }
75 syndromes[i] ^= cdw[0];
76 }
77 }
78
79
80
81 /**
82 * @brief Computes the error locator polynomial (ELP) sigma
83 *
84 * This is a constant time implementation of Berlekamp's simplified algorithm (see @cite lin1983error (Chapter 6 - BCH Codes). <br>
85 * We use the letter p for rho which is initialized at -1. <br>
86 * The array X_sigma_p represents the polynomial X^(mu-rho)*sigma_p(X). <br>
87 * Instead of maintaining a list of sigmas, we update in place both sigma and X_sigma_p. <br>
88 * sigma_copy serves as a temporary save of sigma in case X_sigma_p needs to be updated. <br>
89 * We can properly correct only if the degree of sigma does not exceed PARAM_DELTA.
90 * This means only the first PARAM_DELTA + 1 coefficients of sigma are of value
91 * and we only need to save its first PARAM_DELTA - 1 coefficients.
92 *
93 * @returns the degree of the ELP sigma
94 * @param[out] sigma Array of size (at least) PARAM_DELTA receiving the ELP
95 * @param[in] syndromes Array of size (at least) 2*PARAM_DELTA storing the syndromes
96 */
compute_elp(uint16_t * sigma,const uint16_t * syndromes)97 static uint16_t compute_elp(uint16_t *sigma, const uint16_t *syndromes) {
98 uint16_t deg_sigma = 0;
99 uint16_t deg_sigma_p = 0;
100 uint16_t deg_sigma_copy = 0;
101 uint16_t sigma_copy[PARAM_DELTA + 1] = {0};
102 uint16_t X_sigma_p[PARAM_DELTA + 1] = {0, 1};
103 uint16_t pp = (uint16_t) -1; // 2*rho
104 uint16_t d_p = 1;
105 uint16_t d = syndromes[0];
106
107 uint16_t mask1, mask2, mask12;
108 uint16_t deg_X, deg_X_sigma_p;
109 uint16_t dd;
110 uint16_t mu;
111
112 uint16_t i;
113
114 sigma[0] = 1;
115 for (mu = 0; (mu < (2 * PARAM_DELTA)); ++mu) {
116 // Save sigma in case we need it to update X_sigma_p
117 memcpy(sigma_copy, sigma, 2 * (PARAM_DELTA));
118 deg_sigma_copy = deg_sigma;
119
120 dd = PQCLEAN_HQCRMRS192_CLEAN_gf_mul(d, PQCLEAN_HQCRMRS192_CLEAN_gf_inverse(d_p));
121
122 for (i = 1; (i <= mu + 1) && (i <= PARAM_DELTA); ++i) {
123 sigma[i] ^= PQCLEAN_HQCRMRS192_CLEAN_gf_mul(dd, X_sigma_p[i]);
124 }
125
126 deg_X = mu - pp;
127 deg_X_sigma_p = deg_X + deg_sigma_p;
128
129 // mask1 = 0xffff if(d != 0) and 0 otherwise
130 mask1 = -((uint16_t) - d >> 15);
131
132 // mask2 = 0xffff if(deg_X_sigma_p > deg_sigma) and 0 otherwise
133 mask2 = -((uint16_t) (deg_sigma - deg_X_sigma_p) >> 15);
134
135 // mask12 = 0xffff if the deg_sigma increased and 0 otherwise
136 mask12 = mask1 & mask2;
137 deg_sigma ^= mask12 & (deg_X_sigma_p ^ deg_sigma);
138
139 if (mu == (2 * PARAM_DELTA - 1)) {
140 break;
141 }
142
143 pp ^= mask12 & (mu ^ pp);
144 d_p ^= mask12 & (d ^ d_p);
145 for (i = PARAM_DELTA; i; --i) {
146 X_sigma_p[i] = (mask12 & sigma_copy[i - 1]) ^ (~mask12 & X_sigma_p[i - 1]);
147 }
148
149 deg_sigma_p ^= mask12 & (deg_sigma_copy ^ deg_sigma_p);
150 d = syndromes[mu + 1];
151
152 for (i = 1; (i <= mu + 1) && (i <= PARAM_DELTA); ++i) {
153 d ^= PQCLEAN_HQCRMRS192_CLEAN_gf_mul(sigma[i], syndromes[mu + 1 - i]);
154 }
155 }
156
157 return deg_sigma;
158 }
159
160
161
162 /**
163 * @brief Computes the error polynomial error from the error locator polynomial sigma
164 *
165 * See function PQCLEAN_HQCRMRS192_CLEAN_fft for more details.
166 *
167 * @param[out] error Array of 2^PARAM_M elements receiving the error polynomial
168 * @param[out] error_compact Array of PARAM_DELTA + PARAM_N1 elements receiving a compact representation of the vector error
169 * @param[in] sigma Array of 2^PARAM_FFT elements storing the error locator polynomial
170 */
compute_roots(uint8_t * error,uint16_t * sigma)171 static void compute_roots(uint8_t *error, uint16_t *sigma) {
172 uint16_t w[1 << PARAM_M] = {0};
173
174 PQCLEAN_HQCRMRS192_CLEAN_fft(w, sigma, PARAM_DELTA + 1);
175 PQCLEAN_HQCRMRS192_CLEAN_fft_retrieve_error_poly(error, w);
176 }
177
178
179
180 /**
181 * @brief Computes the polynomial z(x)
182 *
183 * See @cite lin1983error (Chapter 6 - BCH Codes) for more details.
184 *
185 * @param[out] z Array of PARAM_DELTA + 1 elements receiving the polynomial z(x)
186 * @param[in] sigma Array of 2^PARAM_FFT elements storing the error locator polynomial
187 * @param[in] degree Integer that is the degree of polynomial sigma
188 * @param[in] syndromes Array of 2 * PARAM_DELTA storing the syndromes
189 */
compute_z_poly(uint16_t * z,const uint16_t * sigma,uint16_t degree,const uint16_t * syndromes)190 static void compute_z_poly(uint16_t *z, const uint16_t *sigma, uint16_t degree, const uint16_t *syndromes) {
191 size_t i, j;
192 uint16_t mask;
193
194 z[0] = 1;
195
196 for (i = 1; i < PARAM_DELTA + 1; ++i) {
197 mask = -((uint16_t) (i - degree - 1) >> 15);
198 z[i] = mask & sigma[i];
199 }
200
201 z[1] ^= syndromes[0];
202
203 for (i = 2; i <= PARAM_DELTA; ++i) {
204 mask = -((uint16_t) (i - degree - 1) >> 15);
205 z[i] ^= mask & syndromes[i - 1];
206
207 for (j = 1; j < i; ++j) {
208 z[i] ^= mask & PQCLEAN_HQCRMRS192_CLEAN_gf_mul(sigma[j], syndromes[i - j - 1]);
209 }
210 }
211 }
212
213
214
215 /**
216 * @brief Computes the error values
217 *
218 * See @cite lin1983error (Chapter 6 - BCH Codes) for more details.
219 *
220 * @param[out] error_values Array of PARAM_DELTA elements receiving the error values
221 * @param[in] z Array of PARAM_DELTA + 1 elements storing the polynomial z(x)
222 * @param[in] z_degree Integer that is the degree of polynomial z(x)
223 * @param[in] error_compact Array of PARAM_DELTA + PARAM_N1 storing compact representation of the error
224 */
compute_error_values(uint16_t * error_values,const uint16_t * z,const uint8_t * error)225 static void compute_error_values(uint16_t *error_values, const uint16_t *z, const uint8_t *error) {
226 uint16_t beta_j[PARAM_DELTA] = {0};
227 uint16_t e_j[PARAM_DELTA] = {0};
228
229 uint16_t delta_counter;
230 uint16_t delta_real_value;
231 uint16_t found;
232 uint16_t mask1;
233 uint16_t mask2;
234 uint16_t tmp1;
235 uint16_t tmp2;
236 uint16_t inverse;
237 uint16_t inverse_power_j;
238
239 // Compute the beta_{j_i} page 31 of the documentation
240 delta_counter = 0;
241 for (size_t i = 0; i < PARAM_N1; i++) {
242 found = 0;
243 mask1 = (uint16_t) (-((int32_t)error[i]) >> 31); // error[i] != 0
244 for (size_t j = 0; j < PARAM_DELTA; j++) {
245 mask2 = ~((uint16_t) (-((int32_t) j ^ delta_counter) >> 31)); // j == delta_counter
246 beta_j[j] += mask1 & mask2 & gf_exp[i];
247 found += mask1 & mask2 & 1;
248 }
249 delta_counter += found;
250 }
251 delta_real_value = delta_counter;
252
253 // Compute the e_{j_i} page 31 of the documentation
254 for (size_t i = 0; i < PARAM_DELTA; ++i) {
255 tmp1 = 1;
256 tmp2 = 1;
257 inverse = PQCLEAN_HQCRMRS192_CLEAN_gf_inverse(beta_j[i]);
258 inverse_power_j = 1;
259
260 for (size_t j = 1; j <= PARAM_DELTA; ++j) {
261 inverse_power_j = PQCLEAN_HQCRMRS192_CLEAN_gf_mul(inverse_power_j, inverse);
262 tmp1 ^= PQCLEAN_HQCRMRS192_CLEAN_gf_mul(inverse_power_j, z[j]);
263 }
264 for (size_t k = 1; k < PARAM_DELTA; ++k) {
265 tmp2 = PQCLEAN_HQCRMRS192_CLEAN_gf_mul(tmp2, (1 ^ PQCLEAN_HQCRMRS192_CLEAN_gf_mul(inverse, beta_j[(i + k) % PARAM_DELTA])));
266 }
267 mask1 = (uint16_t) (((int16_t) i - delta_real_value) >> 15); // i < delta_real_value
268 e_j[i] = mask1 & PQCLEAN_HQCRMRS192_CLEAN_gf_mul(tmp1, PQCLEAN_HQCRMRS192_CLEAN_gf_inverse(tmp2));
269 }
270
271 // Place the delta e_{j_i} values at the right coordinates of the output vector
272 delta_counter = 0;
273 for (size_t i = 0; i < PARAM_N1; ++i) {
274 found = 0;
275 mask1 = (uint16_t) (-((int32_t)error[i]) >> 31); // error[i] != 0
276 for (size_t j = 0; j < PARAM_DELTA; j++) {
277 mask2 = ~((uint16_t) (-((int32_t) j ^ delta_counter) >> 31)); // j == delta_counter
278 error_values[i] += mask1 & mask2 & e_j[j];
279 found += mask1 & mask2 & 1;
280 }
281 delta_counter += found;
282 }
283 }
284
285
286
287 /**
288 * @brief Correct the errors
289 *
290 * @param[out] cdw Array of PARAM_N1 elements receiving the corrected vector
291 * @param[in] error Array of the error vector
292 * @param[in] error_values Array of PARAM_DELTA elements storing the error values
293 */
correct_errors(uint8_t * cdw,const uint16_t * error_values)294 static void correct_errors(uint8_t *cdw, const uint16_t *error_values) {
295 for (size_t i = 0; i < PARAM_N1; ++i) {
296 cdw[i] ^= error_values[i];
297 }
298 }
299
300
301
302 /**
303 * @brief Decodes the received word
304 *
305 * This function relies on six steps:
306 * <ol>
307 * <li> The first step, is the computation of the 2*PARAM_DELTA syndromes.
308 * <li> The second step is the computation of the error-locator polynomial sigma.
309 * <li> The third step, done by additive FFT, is finding the error-locator numbers by calculating the roots of the polynomial sigma and takings their inverses.
310 * <li> The fourth step, is the polynomial z(x).
311 * <li> The fifth step, is the computation of the error values.
312 * <li> The sixth step is the correction of the errors in the received polynomial.
313 * </ol>
314 * For a more complete picture on Reed-Solomon decoding, see Shu. Lin and Daniel J. Costello in Error Control Coding: Fundamentals and Applications @cite lin1983error
315 *
316 * @param[out] msg Array of size VEC_K_SIZE_64 receiving the decoded message
317 * @param[in] cdw Array of size VEC_N1_SIZE_64 storing the received word
318 */
PQCLEAN_HQCRMRS192_CLEAN_reed_solomon_decode(uint8_t * msg,uint8_t * cdw)319 void PQCLEAN_HQCRMRS192_CLEAN_reed_solomon_decode(uint8_t *msg, uint8_t *cdw) {
320 uint16_t syndromes[2 * PARAM_DELTA] = {0};
321 uint16_t sigma[1 << PARAM_FFT] = {0};
322 uint8_t error[1 << PARAM_M] = {0};
323 uint16_t z[PARAM_N1] = {0};
324 uint16_t error_values[PARAM_N1] = {0};
325 uint16_t deg;
326
327 // Calculate the 2*PARAM_DELTA syndromes
328 compute_syndromes(syndromes, cdw);
329
330 // Compute the error locator polynomial sigma
331 // Sigma's degree is at most PARAM_DELTA but the FFT requires the extra room
332 deg = compute_elp(sigma, syndromes);
333
334 // Compute the error polynomial error
335 compute_roots(error, sigma);
336
337 // Compute the polynomial z(x)
338 compute_z_poly(z, sigma, deg, syndromes);
339
340 // Compute the error values
341 compute_error_values(error_values, z, error);
342
343 // Correct the errors
344 correct_errors(cdw, error_values);
345
346 // Retrieve the message from the decoded codeword
347 memcpy(msg, cdw + (PARAM_G - 1), PARAM_K);
348
349 }
350