cfd78e200cb30f0a716843f9ae800ec10b152caa
[strongswan.git] / src / libstrongswan / plugins / ntru / ntru_crypto / ntru_crypto_ntru_poly.c
1 /******************************************************************************
2 * NTRU Cryptography Reference Source Code
3 * Copyright (c) 2009-2013, by Security Innovation, Inc. All rights reserved.
4 *
5 * ntru_crypto_ntru_poly.c is a component of ntru-crypto.
6 *
7 * Copyright (C) 2009-2013 Security Innovation
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 *
23 *****************************************************************************/
24
25 #include <stdlib.h>
26 #include <string.h>
27 #include "ntru_crypto_ntru_poly.h"
28
29 /* ntru_poly_check_min_weight
30 *
31 * Checks that the number of 0, +1, and -1 trinary ring elements meet or exceed
32 * a minimum weight.
33 */
34
35 bool
36 ntru_poly_check_min_weight(
37 uint16_t num_els, /* in - degree of polynomial */
38 uint8_t *ringels, /* in - pointer to trinary ring elements */
39 uint16_t min_wt) /* in - minimum weight */
40 {
41 uint16_t wt[3];
42 uint16_t i;
43
44 wt[0] = wt[1] = wt[2] = 0;
45 for (i = 0; i < num_els; i++) {
46 ++wt[ringels[i]];
47 }
48 if ((wt[0] < min_wt) || (wt[1] < min_wt) || (wt[2] < min_wt)) {
49 return FALSE;
50 }
51 return TRUE;
52 }
53
54
55 /* ntru_ring_mult_indices
56 *
57 * Multiplies ring element (polynomial) "a" by ring element (polynomial) "b"
58 * to produce ring element (polynomial) "c" in (Z/qZ)[X]/(X^N - 1).
59 * This is a convolution operation.
60 *
61 * Ring element "b" is a sparse trinary polynomial with coefficients -1, 0,
62 * and 1. It is specified by a list, bi, of its nonzero indices containing
63 * indices for the bi_P1_len +1 coefficients followed by the indices for the
64 * bi_M1_len -1 coefficients.
65 * The indices are in the range [0,N).
66 *
67 * The result array "c" may share the same memory space as input array "a",
68 * input array "b", or temp array "t".
69 *
70 * This assumes q is 2^r where 8 < r < 16, so that overflow of the sum
71 * beyond 16 bits does not matter.
72 */
73
74 void
75 ntru_ring_mult_indices(
76 uint16_t const *a, /* in - pointer to ring element a */
77 uint16_t bi_P1_len, /* in - no. of +1 coefficients in b */
78 uint16_t bi_M1_len, /* in - no. of -1 coefficients in b */
79 uint16_t const *bi, /* in - pointer to the list of nonzero
80 indices of ring element b,
81 containing indices for the +1
82 coefficients followed by the
83 indices for -1 coefficients */
84 uint16_t N, /* in - no. of coefficients in a, b, c */
85 uint16_t q, /* in - large modulus */
86 uint16_t *t, /* in - temp buffer of N elements */
87 uint16_t *c) /* out - address for polynomial c */
88 {
89 uint16_t mod_q_mask = q - 1;
90 uint16_t i, j, k;
91
92 /* t[(i+k)%N] = sum i=0 through N-1 of a[i], for b[k] = -1 */
93
94 for (k = 0; k < N; k++)
95 t[k] = 0;
96 for (j = bi_P1_len; j < bi_P1_len + bi_M1_len; j++) {
97 k = bi[j];
98 for (i = 0; k < N; ++i, ++k)
99 t[k] = t[k] + a[i];
100 for (k = 0; i < N; ++i, ++k)
101 t[k] = t[k] + a[i];
102 }
103
104 /* t[(i+k)%N] = -(sum i=0 through N-1 of a[i] for b[k] = -1) */
105
106 for (k = 0; k < N; k++)
107 t[k] = -t[k];
108
109 /* t[(i+k)%N] += sum i=0 through N-1 of a[i] for b[k] = +1 */
110
111 for (j = 0; j < bi_P1_len; j++) {
112 k = bi[j];
113 for (i = 0; k < N; ++i, ++k)
114 t[k] = t[k] + a[i];
115 for (k = 0; i < N; ++i, ++k)
116 t[k] = t[k] + a[i];
117 }
118
119 /* c = (a * b) mod q */
120
121 for (k = 0; k < N; k++)
122 c[k] = t[k] & mod_q_mask;
123 }
124
125
126 /* ntru_ring_mult_product_indices
127 *
128 * Multiplies ring element (polynomial) "a" by ring element (polynomial) "b"
129 * to produce ring element (polynomial) "c" in (Z/qZ)[X]/(X^N - 1).
130 * This is a convolution operation.
131 *
132 * Ring element "b" is represented by the product form b1 * b2 + b3, where
133 * b1, b2, and b3 are each a sparse trinary polynomial with coefficients -1,
134 * 0, and 1. It is specified by a list, bi, of the nonzero indices of b1, b2,
135 * and b3, containing the indices for the +1 coefficients followed by the
136 * indices for the -1 coefficients for each polynomial in that order.
137 * The indices are in the range [0,N).
138 *
139 * The result array "c" may share the same memory space as input array "a",
140 * or input array "b".
141 *
142 * This assumes q is 2^r where 8 < r < 16, so that overflow of the sum
143 * beyond 16 bits does not matter.
144 */
145
146 void
147 ntru_ring_mult_product_indices(
148 uint16_t *a, /* in - pointer to ring element a */
149 uint16_t b1i_len, /* in - no. of +1 or -1 coefficients in b1 */
150 uint16_t b2i_len, /* in - no. of +1 or -1 coefficients in b2 */
151 uint16_t b3i_len, /* in - no. of +1 or -1 coefficients in b3 */
152 uint16_t const *bi, /* in - pointer to the list of nonzero
153 indices of polynomials b1, b2, b3,
154 containing indices for the +1
155 coefficients followed by the
156 indices for -1 coefficients for
157 each polynomial */
158 uint16_t N, /* in - no. of coefficients in a, b, c */
159 uint16_t q, /* in - large modulus */
160 uint16_t *t, /* in - temp buffer of 2N elements */
161 uint16_t *c) /* out - address for polynomial c */
162 {
163 uint16_t *t2 = t + N;
164 uint16_t mod_q_mask = q - 1;
165 uint16_t i;
166
167
168 /* t2 = a * b1 */
169 ntru_ring_mult_indices(a, b1i_len, b1i_len, bi, N, q, t, t2);
170
171 /* t2 = (a * b1) * b2 */
172 ntru_ring_mult_indices(t2, b2i_len, b2i_len, bi + (b1i_len << 1), N, q,
173 t, t2);
174
175 /* t = a * b3 */
176 ntru_ring_mult_indices(a, b3i_len, b3i_len,
177 bi + ((b1i_len + b2i_len) << 1), N, q, t, t);
178
179 /* c = (a * b1 * b2) + (a * b3) */
180 for (i = 0; i < N; i++)
181 c[i] = (t2[i] + t[i]) & mod_q_mask;
182 }
183
184
185 /* ntru_ring_mult_coefficients
186 *
187 * Multiplies ring element (polynomial) "a" by ring element (polynomial) "b"
188 * to produce ring element (polynomial) "c" in (Z/qZ)[X]/(X^N - 1).
189 * This is a convolution operation.
190 *
191 * Ring element "b" has coefficients in the range [0,N).
192 *
193 * This assumes q is 2^r where 8 < r < 16, so that overflow of the sum
194 * beyond 16 bits does not matter.
195 */
196
197 void
198 ntru_ring_mult_coefficients(
199 uint16_t const *a, /* in - pointer to polynomial a */
200 uint16_t const *b, /* in - pointer to polynomial b */
201 uint16_t N, /* in - no. of coefficients in a, b, c */
202 uint16_t q, /* in - large modulus */
203 uint16_t *c) /* out - address for polynomial c */
204 {
205 uint16_t const *bptr = b;
206 uint16_t mod_q_mask = q - 1;
207 uint16_t i, k;
208
209 /* c[k] = sum(a[i] * b[k-i]) mod q */
210 memset(c, 0, N * sizeof(uint16_t));
211 for (k = 0; k < N; k++) {
212 i = 0;
213 while (i <= k)
214 c[k] += a[i++] * *bptr--;
215 bptr += N;
216 while (i < N)
217 c[k] += a[i++] * *bptr--;
218 c[k] &= mod_q_mask;
219 ++bptr;
220 }
221 }
222
223
224 /* ntru_ring_inv
225 *
226 * Finds the inverse of a polynomial, a, in (Z/2^rZ)[X]/(X^N - 1).
227 *
228 * This assumes q is 2^r where 8 < r < 16, so that operations mod q can
229 * wait until the end, and only 16-bit arrays need to be used.
230 */
231
232 bool
233 ntru_ring_inv(
234 uint16_t *a, /* in - pointer to polynomial a */
235 uint16_t N, /* in - no. of coefficients in a */
236 uint16_t q, /* in - large modulus */
237 uint16_t *t, /* in - temp buffer of 2N elements */
238 uint16_t *a_inv) /* out - address for polynomial a^-1 */
239 {
240 uint8_t *b = (uint8_t *)t; /* b cannot be in a_inv since it must be
241 rotated and copied there as a^-1 mod 2 */
242 uint8_t *c = b + N; /* c cannot be in a_inv since it exchanges
243 with b, and b cannot be in a_inv */
244 uint8_t *f = c + N;
245 uint8_t *g = (uint8_t *)a_inv; /* g needs N + 1 bytes */
246 uint16_t *t2 = t + N;
247 uint16_t deg_b;
248 uint16_t deg_c;
249 uint16_t deg_f;
250 uint16_t deg_g;
251 uint16_t k = 0;
252 bool done = FALSE;
253 uint16_t i, j;
254
255 /* form a^-1 in (Z/2Z)[X]/X^N - 1) */
256 memset(b, 0, (N << 1)); /* clear to init b, c */
257
258 /* b(X) = 1 */
259 b[0] = 1;
260 deg_b = 0;
261
262 /* c(X) = 0 (cleared above) */
263 deg_c = 0;
264
265 /* f(X) = a(X) mod 2 */
266 for (i = 0; i < N; i++)
267 f[i] = (uint8_t)(a[i] & 1);
268 deg_f = N - 1;
269
270 /* g(X) = X^N - 1 */
271 g[0] = 1;
272 memset(g + 1, 0, N - 1);
273 g[N] = 1;
274 deg_g = N;
275
276 /* until f(X) = 1 */
277
278 while (!done)
279 {
280
281 /* while f[0] = 0, f(X) /= X, c(X) *= X, k++ */
282
283 for (i = 0; (i <= deg_f) && (f[i] == 0); ++i);
284 if (i > deg_f)
285 return FALSE;
286 if (i) {
287 f = f + i;
288 deg_f = deg_f - i;
289 deg_c = deg_c + i;
290 for (j = deg_c; j >= i; j--)
291 c[j] = c[j-i];
292 for (j = 0; j < i; j++)
293 c[j] = 0;
294 k = k + i;
295 }
296
297 /* adjust degree of f(X) if the highest coefficients are zero
298 * Note: f[0] = 1 from above so the loop will terminate.
299 */
300
301 while (f[deg_f] == 0)
302 --deg_f;
303
304 /* if f(X) = 1, done
305 * Note: f[0] = 1 from above, so only check the x term and up
306 */
307
308 for (i = 1; (i <= deg_f) && (f[i] == 0); ++i);
309 if (i > deg_f) {
310 done = TRUE;
311 break;
312 }
313
314 /* if deg_f < deg_g, f <-> g, b <-> c */
315
316 if (deg_f < deg_g) {
317 uint8_t *x;
318
319 x = f;
320 f = g;
321 g = x;
322 deg_f ^= deg_g;
323 deg_g ^= deg_f;
324 deg_f ^= deg_g;
325 x = b;
326 b = c;
327 c = x;
328 deg_b ^= deg_c;
329 deg_c ^= deg_b;
330 deg_b ^= deg_c;
331 }
332
333 /* f(X) += g(X), b(X) += c(X) */
334
335 for (i = 0; i <= deg_g; i++)
336 f[i] ^= g[i];
337
338 if (deg_c > deg_b)
339 deg_b = deg_c;
340 for (i = 0; i <= deg_c; i++)
341 b[i] ^= c[i];
342 }
343
344 /* a^-1 in (Z/2Z)[X]/(X^N - 1) = b(X) shifted left k coefficients */
345
346 j = 0;
347 if (k >= N)
348 k = k - N;
349 for (i = k; i < N; i++)
350 a_inv[j++] = (uint16_t)(b[i]);
351 for (i = 0; i < k; i++)
352 a_inv[j++] = (uint16_t)(b[i]);
353
354 /* lift a^-1 in (Z/2Z)[X]/(X^N - 1) to a^-1 in (Z/qZ)[X]/(X^N -1) */
355
356 for (j = 0; j < 4; ++j) { /* assumes 256 < q <= 65536 */
357
358 /* a^-1 = a^-1 * (2 - a * a^-1) mod q */
359
360 memcpy(t2, a_inv, N * sizeof(uint16_t));
361 ntru_ring_mult_coefficients(a, t2, N, q, t);
362 for (i = 0; i < N; ++i)
363 t[i] = q - t[i];
364 t[0] = t[0] + 2;
365 ntru_ring_mult_coefficients(t2, t, N, q, a_inv);
366 }
367
368 return TRUE;
369
370
371 }
372
373