Implemented ntru_trits class
[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 /******************************************************************************
26 *
27 * File: ntru_crypto_ntru_poly.c
28 *
29 * Contents: Routines for generating and operating on polynomials in the
30 * NTRU algorithm.
31 *
32 *****************************************************************************/
33
34
35 #include <stdlib.h>
36 #include <string.h>
37 #include "ntru_crypto_ntru_poly.h"
38
39 #include "ntru_mgf1.h"
40
41 #include <utils/debug.h>
42
43 /* ntru_gen_poly
44 *
45 * Generates polynomials by creating for each polynomial, a list of the
46 * indices of the +1 coefficients followed by a list of the indices of
47 * the -1 coefficients.
48 *
49 * If a single polynomial is generated (non-product form), indices_counts
50 * contains a single value of the total number of indices (for +1 and -1
51 * comefficients combined).
52 *
53 * If multiple polynomials are generated (for product form), their lists of
54 * indices are sequentially stored in the indices buffer. Each byte of
55 * indices_counts contains the total number of indices (for +1 and -1
56 * coefficients combined) for a single polynomial, beginning with the
57 * low-order byte for the first polynomial. The high-order byte is unused.
58 *
59 * Returns NTRU_OK if successful.
60 * Returns HASH_BAD_ALG if the algorithm is not supported.
61 *
62 */
63
64 uint32_t
65 ntru_gen_poly(
66 hash_algorithm_t hash_algid, /* in - hash algorithm ID for
67 IGF-2 */
68 uint8_t min_calls, /* in - minimum no. of hash
69 calls */
70 uint16_t seed_len, /* in - no. of octets in seed */
71 uint8_t *seed, /* in - pointer to seed */
72 uint8_t *buf, /* in - pointer to working
73 buffer */
74 uint16_t N, /* in - max index + 1 */
75 uint8_t c_bits, /* in - no. bits for candidate */
76 uint16_t limit, /* in - conversion to index
77 limit */
78 bool is_product_form, /* in - if generating multiple
79 polys */
80 uint32_t indices_counts, /* in - nos. of indices needed */
81 uint16_t *indices) /* out - address for indices */
82 {
83 uint8_t md_len;
84 uint8_t *octets;
85 uint8_t *used;
86 uint8_t num_polys;
87 uint16_t num_indices;
88 uint16_t octets_available;
89 uint16_t index_cnt = 0;
90 uint8_t left = 0;
91 uint8_t num_left = 0;
92 ntru_mgf1_t *mgf1;
93
94 /* generate minimum MGF1 output */
95 DBG2(DBG_LIB, "MGF1 is seeded with %u bytes", seed_len);
96 mgf1 = ntru_mgf1_create(hash_algid, chunk_create(seed, seed_len), TRUE);
97 if (!mgf1)
98 {
99 return NTRU_MGF1_FAIL;
100 }
101 md_len = mgf1->get_hash_size(mgf1);
102 octets = buf;
103 octets_available = min_calls * md_len;
104
105 /* init indices counts for number of polynomials being generated */
106 if (is_product_form) {
107
108 /* number of indices for poly1 is in low byte of indices_counts,
109 * number of indices for poly2 and poly3 are in next higher bytes
110 */
111
112 num_polys = 3;
113 num_indices = (uint16_t)(indices_counts & 0xff);
114 indices_counts >>= 8;
115
116 } else {
117
118 /* number of bytes for poly is in low 16 bits of indices_counts */
119
120 num_polys = 1;
121 num_indices = (uint16_t)indices_counts;
122 }
123
124 /* init used-index array */
125
126 used = buf + octets_available;
127 memset(used, 0, N);
128
129 /* generate indices (IGF-2) for all polynomials */
130 DBG2(DBG_LIB, "MGF1 generates %u octets for %u indices",
131 octets_available, num_indices);
132 if (!mgf1->get_mask(mgf1, octets_available, octets))
133 {
134 mgf1->destroy(mgf1);
135 return NTRU_MGF1_FAIL;
136 }
137
138 while (num_polys > 0) {
139
140 /* generate indices for a single polynomial */
141
142 while (index_cnt < num_indices) {
143 uint16_t index;
144 uint8_t num_needed;
145
146 /* form next index to convert to an index */
147
148 do {
149 /* use any leftover bits first */
150
151 if (num_left != 0) {
152 index = left << (c_bits - num_left);
153 } else {
154 index = 0;
155 }
156
157 /* get the rest of the bits needed from new octets */
158
159 num_needed = c_bits - num_left;
160 while (num_needed != 0)
161 {
162
163 /* get another octet */
164 if (octets_available == 0)
165 {
166 octets = buf;
167 octets_available = md_len;
168
169 DBG2(DBG_LIB, "MGF1 generates another %u octets for the "
170 "remaining %u indices", octets_available,
171 num_indices - index_cnt);
172 if (!mgf1->get_mask(mgf1, octets_available, octets))
173 {
174 mgf1->destroy(mgf1);
175 return NTRU_MGF1_FAIL;
176 }
177 }
178 left = *octets++;
179 --octets_available;
180
181 if (num_needed <= 8)
182 {
183
184 /* all bits needed to fill the index are in this octet */
185
186 index |= ((uint16_t)(left)) >> (8 - num_needed);
187 num_left = 8 - num_needed;
188 num_needed = 0;
189 left &= 0xff >> (8 - num_left);
190
191 } else {
192
193 /* another octet will be needed after using this
194 * whole octet
195 */
196
197 index |= ((uint16_t)left) << (num_needed - 8);
198 num_needed -= 8;
199 }
200 }
201 } while (index >= limit);
202
203 /* form index and check if unique */
204
205 index %= N;
206 if (!used[index])
207 {
208 used[index] = 1;
209 indices[index_cnt] = index;
210 ++index_cnt;
211 }
212 }
213 --num_polys;
214
215 /* init for next polynomial if another polynomial to be generated */
216
217 if (num_polys > 0)
218 {
219 memset(used, 0, N);
220 num_indices = num_indices +
221 (uint16_t)(indices_counts & 0xff);
222 indices_counts >>= 8;
223 }
224 }
225 mgf1->destroy(mgf1);
226
227 return NTRU_OK;
228 }
229
230
231 /* ntru_poly_check_min_weight
232 *
233 * Checks that the number of 0, +1, and -1 trinary ring elements meet or exceed
234 * a minimum weight.
235 */
236
237 bool
238 ntru_poly_check_min_weight(
239 uint16_t num_els, /* in - degree of polynomial */
240 uint8_t *ringels, /* in - pointer to trinary ring elements */
241 uint16_t min_wt) /* in - minimum weight */
242 {
243 uint16_t wt[3];
244 uint16_t i;
245
246 wt[0] = wt[1] = wt[2] = 0;
247 for (i = 0; i < num_els; i++) {
248 ++wt[ringels[i]];
249 }
250 if ((wt[0] < min_wt) || (wt[1] < min_wt) || (wt[2] < min_wt)) {
251 return FALSE;
252 }
253 return TRUE;
254 }
255
256
257 /* ntru_ring_mult_indices
258 *
259 * Multiplies ring element (polynomial) "a" by ring element (polynomial) "b"
260 * to produce ring element (polynomial) "c" in (Z/qZ)[X]/(X^N - 1).
261 * This is a convolution operation.
262 *
263 * Ring element "b" is a sparse trinary polynomial with coefficients -1, 0,
264 * and 1. It is specified by a list, bi, of its nonzero indices containing
265 * indices for the bi_P1_len +1 coefficients followed by the indices for the
266 * bi_M1_len -1 coefficients.
267 * The indices are in the range [0,N).
268 *
269 * The result array "c" may share the same memory space as input array "a",
270 * input array "b", or temp array "t".
271 *
272 * This assumes q is 2^r where 8 < r < 16, so that overflow of the sum
273 * beyond 16 bits does not matter.
274 */
275
276 void
277 ntru_ring_mult_indices(
278 uint16_t const *a, /* in - pointer to ring element a */
279 uint16_t bi_P1_len, /* in - no. of +1 coefficients in b */
280 uint16_t bi_M1_len, /* in - no. of -1 coefficients in b */
281 uint16_t const *bi, /* in - pointer to the list of nonzero
282 indices of ring element b,
283 containing indices for the +1
284 coefficients followed by the
285 indices for -1 coefficients */
286 uint16_t N, /* in - no. of coefficients in a, b, c */
287 uint16_t q, /* in - large modulus */
288 uint16_t *t, /* in - temp buffer of N elements */
289 uint16_t *c) /* out - address for polynomial c */
290 {
291 uint16_t mod_q_mask = q - 1;
292 uint16_t i, j, k;
293
294 /* t[(i+k)%N] = sum i=0 through N-1 of a[i], for b[k] = -1 */
295
296 for (k = 0; k < N; k++)
297 t[k] = 0;
298 for (j = bi_P1_len; j < bi_P1_len + bi_M1_len; j++) {
299 k = bi[j];
300 for (i = 0; k < N; ++i, ++k)
301 t[k] = t[k] + a[i];
302 for (k = 0; i < N; ++i, ++k)
303 t[k] = t[k] + a[i];
304 }
305
306 /* t[(i+k)%N] = -(sum i=0 through N-1 of a[i] for b[k] = -1) */
307
308 for (k = 0; k < N; k++)
309 t[k] = -t[k];
310
311 /* t[(i+k)%N] += sum i=0 through N-1 of a[i] for b[k] = +1 */
312
313 for (j = 0; j < bi_P1_len; j++) {
314 k = bi[j];
315 for (i = 0; k < N; ++i, ++k)
316 t[k] = t[k] + a[i];
317 for (k = 0; i < N; ++i, ++k)
318 t[k] = t[k] + a[i];
319 }
320
321 /* c = (a * b) mod q */
322
323 for (k = 0; k < N; k++)
324 c[k] = t[k] & mod_q_mask;
325 }
326
327
328 /* ntru_ring_mult_product_indices
329 *
330 * Multiplies ring element (polynomial) "a" by ring element (polynomial) "b"
331 * to produce ring element (polynomial) "c" in (Z/qZ)[X]/(X^N - 1).
332 * This is a convolution operation.
333 *
334 * Ring element "b" is represented by the product form b1 * b2 + b3, where
335 * b1, b2, and b3 are each a sparse trinary polynomial with coefficients -1,
336 * 0, and 1. It is specified by a list, bi, of the nonzero indices of b1, b2,
337 * and b3, containing the indices for the +1 coefficients followed by the
338 * indices for the -1 coefficients for each polynomial in that order.
339 * The indices are in the range [0,N).
340 *
341 * The result array "c" may share the same memory space as input array "a",
342 * or input array "b".
343 *
344 * This assumes q is 2^r where 8 < r < 16, so that overflow of the sum
345 * beyond 16 bits does not matter.
346 */
347
348 void
349 ntru_ring_mult_product_indices(
350 uint16_t *a, /* in - pointer to ring element a */
351 uint16_t b1i_len, /* in - no. of +1 or -1 coefficients in b1 */
352 uint16_t b2i_len, /* in - no. of +1 or -1 coefficients in b2 */
353 uint16_t b3i_len, /* in - no. of +1 or -1 coefficients in b3 */
354 uint16_t const *bi, /* in - pointer to the list of nonzero
355 indices of polynomials b1, b2, b3,
356 containing indices for the +1
357 coefficients followed by the
358 indices for -1 coefficients for
359 each polynomial */
360 uint16_t N, /* in - no. of coefficients in a, b, c */
361 uint16_t q, /* in - large modulus */
362 uint16_t *t, /* in - temp buffer of 2N elements */
363 uint16_t *c) /* out - address for polynomial c */
364 {
365 uint16_t *t2 = t + N;
366 uint16_t mod_q_mask = q - 1;
367 uint16_t i;
368
369
370 /* t2 = a * b1 */
371 ntru_ring_mult_indices(a, b1i_len, b1i_len, bi, N, q, t, t2);
372
373 /* t2 = (a * b1) * b2 */
374 ntru_ring_mult_indices(t2, b2i_len, b2i_len, bi + (b1i_len << 1), N, q,
375 t, t2);
376
377 /* t = a * b3 */
378 ntru_ring_mult_indices(a, b3i_len, b3i_len,
379 bi + ((b1i_len + b2i_len) << 1), N, q, t, t);
380
381 /* c = (a * b1 * b2) + (a * b3) */
382 for (i = 0; i < N; i++)
383 c[i] = (t2[i] + t[i]) & mod_q_mask;
384 }
385
386
387 /* ntru_ring_mult_coefficients
388 *
389 * Multiplies ring element (polynomial) "a" by ring element (polynomial) "b"
390 * to produce ring element (polynomial) "c" in (Z/qZ)[X]/(X^N - 1).
391 * This is a convolution operation.
392 *
393 * Ring element "b" has coefficients in the range [0,N).
394 *
395 * This assumes q is 2^r where 8 < r < 16, so that overflow of the sum
396 * beyond 16 bits does not matter.
397 */
398
399 void
400 ntru_ring_mult_coefficients(
401 uint16_t const *a, /* in - pointer to polynomial a */
402 uint16_t const *b, /* in - pointer to polynomial b */
403 uint16_t N, /* in - no. of coefficients in a, b, c */
404 uint16_t q, /* in - large modulus */
405 uint16_t *c) /* out - address for polynomial c */
406 {
407 uint16_t const *bptr = b;
408 uint16_t mod_q_mask = q - 1;
409 uint16_t i, k;
410
411 /* c[k] = sum(a[i] * b[k-i]) mod q */
412 memset(c, 0, N * sizeof(uint16_t));
413 for (k = 0; k < N; k++) {
414 i = 0;
415 while (i <= k)
416 c[k] += a[i++] * *bptr--;
417 bptr += N;
418 while (i < N)
419 c[k] += a[i++] * *bptr--;
420 c[k] &= mod_q_mask;
421 ++bptr;
422 }
423 }
424
425
426 /* ntru_ring_inv
427 *
428 * Finds the inverse of a polynomial, a, in (Z/2^rZ)[X]/(X^N - 1).
429 *
430 * This assumes q is 2^r where 8 < r < 16, so that operations mod q can
431 * wait until the end, and only 16-bit arrays need to be used.
432 */
433
434 bool
435 ntru_ring_inv(
436 uint16_t *a, /* in - pointer to polynomial a */
437 uint16_t N, /* in - no. of coefficients in a */
438 uint16_t q, /* in - large modulus */
439 uint16_t *t, /* in - temp buffer of 2N elements */
440 uint16_t *a_inv) /* out - address for polynomial a^-1 */
441 {
442 uint8_t *b = (uint8_t *)t; /* b cannot be in a_inv since it must be
443 rotated and copied there as a^-1 mod 2 */
444 uint8_t *c = b + N; /* c cannot be in a_inv since it exchanges
445 with b, and b cannot be in a_inv */
446 uint8_t *f = c + N;
447 uint8_t *g = (uint8_t *)a_inv; /* g needs N + 1 bytes */
448 uint16_t *t2 = t + N;
449 uint16_t deg_b;
450 uint16_t deg_c;
451 uint16_t deg_f;
452 uint16_t deg_g;
453 uint16_t k = 0;
454 bool done = FALSE;
455 uint16_t i, j;
456
457 /* form a^-1 in (Z/2Z)[X]/X^N - 1) */
458 memset(b, 0, (N << 1)); /* clear to init b, c */
459
460 /* b(X) = 1 */
461 b[0] = 1;
462 deg_b = 0;
463
464 /* c(X) = 0 (cleared above) */
465 deg_c = 0;
466
467 /* f(X) = a(X) mod 2 */
468 for (i = 0; i < N; i++)
469 f[i] = (uint8_t)(a[i] & 1);
470 deg_f = N - 1;
471
472 /* g(X) = X^N - 1 */
473 g[0] = 1;
474 memset(g + 1, 0, N - 1);
475 g[N] = 1;
476 deg_g = N;
477
478 /* until f(X) = 1 */
479
480 while (!done)
481 {
482
483 /* while f[0] = 0, f(X) /= X, c(X) *= X, k++ */
484
485 for (i = 0; (i <= deg_f) && (f[i] == 0); ++i);
486 if (i > deg_f)
487 return FALSE;
488 if (i) {
489 f = f + i;
490 deg_f = deg_f - i;
491 deg_c = deg_c + i;
492 for (j = deg_c; j >= i; j--)
493 c[j] = c[j-i];
494 for (j = 0; j < i; j++)
495 c[j] = 0;
496 k = k + i;
497 }
498
499 /* adjust degree of f(X) if the highest coefficients are zero
500 * Note: f[0] = 1 from above so the loop will terminate.
501 */
502
503 while (f[deg_f] == 0)
504 --deg_f;
505
506 /* if f(X) = 1, done
507 * Note: f[0] = 1 from above, so only check the x term and up
508 */
509
510 for (i = 1; (i <= deg_f) && (f[i] == 0); ++i);
511 if (i > deg_f) {
512 done = TRUE;
513 break;
514 }
515
516 /* if deg_f < deg_g, f <-> g, b <-> c */
517
518 if (deg_f < deg_g) {
519 uint8_t *x;
520
521 x = f;
522 f = g;
523 g = x;
524 deg_f ^= deg_g;
525 deg_g ^= deg_f;
526 deg_f ^= deg_g;
527 x = b;
528 b = c;
529 c = x;
530 deg_b ^= deg_c;
531 deg_c ^= deg_b;
532 deg_b ^= deg_c;
533 }
534
535 /* f(X) += g(X), b(X) += c(X) */
536
537 for (i = 0; i <= deg_g; i++)
538 f[i] ^= g[i];
539
540 if (deg_c > deg_b)
541 deg_b = deg_c;
542 for (i = 0; i <= deg_c; i++)
543 b[i] ^= c[i];
544 }
545
546 /* a^-1 in (Z/2Z)[X]/(X^N - 1) = b(X) shifted left k coefficients */
547
548 j = 0;
549 if (k >= N)
550 k = k - N;
551 for (i = k; i < N; i++)
552 a_inv[j++] = (uint16_t)(b[i]);
553 for (i = 0; i < k; i++)
554 a_inv[j++] = (uint16_t)(b[i]);
555
556 /* lift a^-1 in (Z/2Z)[X]/(X^N - 1) to a^-1 in (Z/qZ)[X]/(X^N -1) */
557
558 for (j = 0; j < 4; ++j) { /* assumes 256 < q <= 65536 */
559
560 /* a^-1 = a^-1 * (2 - a * a^-1) mod q */
561
562 memcpy(t2, a_inv, N * sizeof(uint16_t));
563 ntru_ring_mult_coefficients(a, t2, N, q, t);
564 for (i = 0; i < N; ++i)
565 t[i] = q - t[i];
566 t[0] = t[0] + 2;
567 ntru_ring_mult_coefficients(t2, t, N, q, a_inv);
568 }
569
570 return TRUE;
571
572
573 }
574
575