blob: 61fe2483f4fd7760d2dcfee24f11badaae678ef4 [file] [log] [blame]
markus@openbsd.org5f004622019-01-30 19:51:15 +00001/* $OpenBSD: sntrup4591761.c,v 1.3 2019/01/30 19:51:15 markus Exp $ */
2
3/*
4 * Public Domain, Authors:
5 * - Daniel J. Bernstein
6 * - Chitchanok Chuengsatiansup
7 * - Tanja Lange
8 * - Christine van Vredendaal
9 */
10
Darren Tuckera0ca4002019-04-01 20:07:23 +110011#include "includes.h"
12
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000013#include <string.h>
14#include "crypto_api.h"
15
djm@openbsd.org533cfb02019-01-21 22:18:24 +000016/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/int32_sort.h */
17#ifndef int32_sort_h
18#define int32_sort_h
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000019
20
djm@openbsd.org533cfb02019-01-21 22:18:24 +000021static void int32_sort(crypto_int32 *,int);
22
23#endif
24
25/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/int32_sort.c */
26/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
27
28
29static void minmax(crypto_int32 *x,crypto_int32 *y)
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000030{
djm@openbsd.org533cfb02019-01-21 22:18:24 +000031 crypto_uint32 xi = *x;
32 crypto_uint32 yi = *y;
33 crypto_uint32 xy = xi ^ yi;
34 crypto_uint32 c = yi - xi;
35 c ^= xy & (c ^ yi);
36 c >>= 31;
37 c = -c;
38 c &= xy;
39 *x = xi ^ c;
40 *y = yi ^ c;
41}
42
43static void int32_sort(crypto_int32 *x,int n)
44{
45 int top,p,q,i;
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000046
47 if (n < 2) return;
48 top = 1;
49 while (top < n - top) top += top;
50
51 for (p = top;p > 0;p >>= 1) {
52 for (i = 0;i < n - p;++i)
53 if (!(i & p))
djm@openbsd.org533cfb02019-01-21 22:18:24 +000054 minmax(x + i,x + i + p);
55 for (q = top;q > p;q >>= 1)
56 for (i = 0;i < n - q;++i)
57 if (!(i & p))
58 minmax(x + i + p,x + i + q);
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000059 }
60}
61
djm@openbsd.org533cfb02019-01-21 22:18:24 +000062/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/small.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000063#ifndef small_h
64#define small_h
65
66
67typedef crypto_int8 small;
68
69static void small_encode(unsigned char *,const small *);
70
71static void small_decode(small *,const unsigned char *);
72
73
74static void small_random(small *);
75
76static void small_random_weightw(small *);
77
78#endif
79
djm@openbsd.org533cfb02019-01-21 22:18:24 +000080/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/mod3.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +000081#ifndef mod3_h
82#define mod3_h
83
84
85/* -1 if x is nonzero, 0 otherwise */
86static inline int mod3_nonzero_mask(small x)
87{
88 return -x*x;
89}
90
91/* input between -100000 and 100000 */
92/* output between -1 and 1 */
93static inline small mod3_freeze(crypto_int32 a)
94{
95 a -= 3 * ((10923 * a) >> 15);
96 a -= 3 * ((89478485 * a + 134217728) >> 28);
97 return a;
98}
99
100static inline small mod3_minusproduct(small a,small b,small c)
101{
102 crypto_int32 A = a;
103 crypto_int32 B = b;
104 crypto_int32 C = c;
105 return mod3_freeze(A - B * C);
106}
107
108static inline small mod3_plusproduct(small a,small b,small c)
109{
110 crypto_int32 A = a;
111 crypto_int32 B = b;
112 crypto_int32 C = c;
113 return mod3_freeze(A + B * C);
114}
115
116static inline small mod3_product(small a,small b)
117{
118 return a * b;
119}
120
121static inline small mod3_sum(small a,small b)
122{
123 crypto_int32 A = a;
124 crypto_int32 B = b;
125 return mod3_freeze(A + B);
126}
127
128static inline small mod3_reciprocal(small a1)
129{
130 return a1;
131}
132
133static inline small mod3_quotient(small num,small den)
134{
135 return mod3_product(num,mod3_reciprocal(den));
136}
137
138#endif
139
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000140/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/modq.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000141#ifndef modq_h
142#define modq_h
143
144
145typedef crypto_int16 modq;
146
147/* -1 if x is nonzero, 0 otherwise */
148static inline int modq_nonzero_mask(modq x)
149{
150 crypto_int32 r = (crypto_uint16) x;
151 r = -r;
152 r >>= 30;
153 return r;
154}
155
156/* input between -9000000 and 9000000 */
157/* output between -2295 and 2295 */
158static inline modq modq_freeze(crypto_int32 a)
159{
160 a -= 4591 * ((228 * a) >> 20);
161 a -= 4591 * ((58470 * a + 134217728) >> 28);
162 return a;
163}
164
165static inline modq modq_minusproduct(modq a,modq b,modq c)
166{
167 crypto_int32 A = a;
168 crypto_int32 B = b;
169 crypto_int32 C = c;
170 return modq_freeze(A - B * C);
171}
172
173static inline modq modq_plusproduct(modq a,modq b,modq c)
174{
175 crypto_int32 A = a;
176 crypto_int32 B = b;
177 crypto_int32 C = c;
178 return modq_freeze(A + B * C);
179}
180
181static inline modq modq_product(modq a,modq b)
182{
183 crypto_int32 A = a;
184 crypto_int32 B = b;
185 return modq_freeze(A * B);
186}
187
188static inline modq modq_square(modq a)
189{
190 crypto_int32 A = a;
191 return modq_freeze(A * A);
192}
193
194static inline modq modq_sum(modq a,modq b)
195{
196 crypto_int32 A = a;
197 crypto_int32 B = b;
198 return modq_freeze(A + B);
199}
200
201static inline modq modq_reciprocal(modq a1)
202{
203 modq a2 = modq_square(a1);
204 modq a3 = modq_product(a2,a1);
205 modq a4 = modq_square(a2);
206 modq a8 = modq_square(a4);
207 modq a16 = modq_square(a8);
208 modq a32 = modq_square(a16);
209 modq a35 = modq_product(a32,a3);
210 modq a70 = modq_square(a35);
211 modq a140 = modq_square(a70);
212 modq a143 = modq_product(a140,a3);
213 modq a286 = modq_square(a143);
214 modq a572 = modq_square(a286);
215 modq a1144 = modq_square(a572);
216 modq a1147 = modq_product(a1144,a3);
217 modq a2294 = modq_square(a1147);
218 modq a4588 = modq_square(a2294);
219 modq a4589 = modq_product(a4588,a1);
220 return a4589;
221}
222
223static inline modq modq_quotient(modq num,modq den)
224{
225 return modq_product(num,modq_reciprocal(den));
226}
227
228#endif
229
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000230/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/params.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000231#ifndef params_h
232#define params_h
233
234#define q 4591
235/* XXX: also built into modq in various ways */
236
237#define qshift 2295
238#define p 761
239#define w 286
240
241#define rq_encode_len 1218
242#define small_encode_len 191
243
244#endif
245
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000246/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/r3.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000247#ifndef r3_h
248#define r3_h
249
250
251static void r3_mult(small *,const small *,const small *);
252
253extern int r3_recip(small *,const small *);
254
255#endif
256
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000257/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000258#ifndef rq_h
259#define rq_h
260
261
262static void rq_encode(unsigned char *,const modq *);
263
264static void rq_decode(modq *,const unsigned char *);
265
266static void rq_encoderounded(unsigned char *,const modq *);
267
268static void rq_decoderounded(modq *,const unsigned char *);
269
270static void rq_round3(modq *,const modq *);
271
272static void rq_mult(modq *,const modq *,const small *);
273
274int rq_recip3(modq *,const small *);
275
276#endif
277
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000278/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/swap.h */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000279#ifndef swap_h
280#define swap_h
281
282static void swap(void *,void *,int,int);
283
284#endif
285
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000286/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/dec.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000287/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
288
289#ifdef KAT
290#endif
291
292
293int crypto_kem_sntrup4591761_dec(
294 unsigned char *k,
295 const unsigned char *cstr,
296 const unsigned char *sk
297)
298{
299 small f[p];
300 modq h[p];
301 small grecip[p];
302 modq c[p];
303 modq t[p];
304 small t3[p];
305 small r[p];
306 modq hr[p];
307 unsigned char rstr[small_encode_len];
308 unsigned char hash[64];
309 int i;
310 int result = 0;
311 int weight;
312
313 small_decode(f,sk);
314 small_decode(grecip,sk + small_encode_len);
315 rq_decode(h,sk + 2 * small_encode_len);
316
317 rq_decoderounded(c,cstr + 32);
318
319 rq_mult(t,c,f);
320 for (i = 0;i < p;++i) t3[i] = mod3_freeze(modq_freeze(3*t[i]));
321
322 r3_mult(r,t3,grecip);
323
324#ifdef KAT
325 {
326 int j;
327 printf("decrypt r:");
328 for (j = 0;j < p;++j)
329 if (r[j] == 1) printf(" +%d",j);
330 else if (r[j] == -1) printf(" -%d",j);
331 printf("\n");
332 }
333#endif
334
335 weight = 0;
336 for (i = 0;i < p;++i) weight += (1 & r[i]);
337 weight -= w;
338 result |= modq_nonzero_mask(weight); /* XXX: puts limit on p */
339
340 rq_mult(hr,h,r);
341 rq_round3(hr,hr);
342 for (i = 0;i < p;++i) result |= modq_nonzero_mask(hr[i] - c[i]);
343
344 small_encode(rstr,r);
345 crypto_hash_sha512(hash,rstr,sizeof rstr);
346 result |= crypto_verify_32(hash,cstr);
347
348 for (i = 0;i < 32;++i) k[i] = (hash[32 + i] & ~result);
349 return result;
350}
351
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000352/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/enc.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000353/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
354
355#ifdef KAT
356#endif
357
358
359int crypto_kem_sntrup4591761_enc(
360 unsigned char *cstr,
361 unsigned char *k,
362 const unsigned char *pk
363)
364{
365 small r[p];
366 modq h[p];
367 modq c[p];
368 unsigned char rstr[small_encode_len];
369 unsigned char hash[64];
370
371 small_random_weightw(r);
372
373#ifdef KAT
374 {
375 int i;
376 printf("encrypt r:");
377 for (i = 0;i < p;++i)
378 if (r[i] == 1) printf(" +%d",i);
379 else if (r[i] == -1) printf(" -%d",i);
380 printf("\n");
381 }
382#endif
383
384 small_encode(rstr,r);
385 crypto_hash_sha512(hash,rstr,sizeof rstr);
386
387 rq_decode(h,pk);
388 rq_mult(c,h,r);
389 rq_round3(c,c);
390
391 memcpy(k,hash + 32,32);
392 memcpy(cstr,hash,32);
393 rq_encoderounded(cstr + 32,c);
394
395 return 0;
396}
397
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000398/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/keypair.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000399/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
400
401
402#if crypto_kem_sntrup4591761_PUBLICKEYBYTES != rq_encode_len
403#error "crypto_kem_sntrup4591761_PUBLICKEYBYTES must match rq_encode_len"
404#endif
405#if crypto_kem_sntrup4591761_SECRETKEYBYTES != rq_encode_len + 2 * small_encode_len
406#error "crypto_kem_sntrup4591761_SECRETKEYBYTES must match rq_encode_len + 2 * small_encode_len"
407#endif
408
409int crypto_kem_sntrup4591761_keypair(unsigned char *pk,unsigned char *sk)
410{
411 small g[p];
412 small grecip[p];
413 small f[p];
414 modq f3recip[p];
415 modq h[p];
416
417 do
418 small_random(g);
419 while (r3_recip(grecip,g) != 0);
420
421 small_random_weightw(f);
422 rq_recip3(f3recip,f);
423
424 rq_mult(h,f3recip,g);
425
426 rq_encode(pk,h);
427 small_encode(sk,f);
428 small_encode(sk + small_encode_len,grecip);
429 memcpy(sk + 2 * small_encode_len,pk,rq_encode_len);
430
431 return 0;
432}
433
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000434/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/r3_mult.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000435/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
436
437
438static void r3_mult(small *h,const small *f,const small *g)
439{
440 small fg[p + p - 1];
441 small result;
442 int i, j;
443
444 for (i = 0;i < p;++i) {
445 result = 0;
446 for (j = 0;j <= i;++j)
447 result = mod3_plusproduct(result,f[j],g[i - j]);
448 fg[i] = result;
449 }
450 for (i = p;i < p + p - 1;++i) {
451 result = 0;
452 for (j = i - p + 1;j < p;++j)
453 result = mod3_plusproduct(result,f[j],g[i - j]);
454 fg[i] = result;
455 }
456
457 for (i = p + p - 2;i >= p;--i) {
458 fg[i - p] = mod3_sum(fg[i - p],fg[i]);
459 fg[i - p + 1] = mod3_sum(fg[i - p + 1],fg[i]);
460 }
461
462 for (i = 0;i < p;++i)
463 h[i] = fg[i];
464}
465
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000466/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/r3_recip.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000467/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
468
469
470/* caller must ensure that x-y does not overflow */
471static int smaller_mask_r3_recip(int x,int y)
472{
473 return (x - y) >> 31;
474}
475
476static void vectormod3_product(small *z,int len,const small *x,const small c)
477{
478 int i;
479 for (i = 0;i < len;++i) z[i] = mod3_product(x[i],c);
480}
481
482static void vectormod3_minusproduct(small *z,int len,const small *x,const small *y,const small c)
483{
484 int i;
485 for (i = 0;i < len;++i) z[i] = mod3_minusproduct(x[i],y[i],c);
486}
487
488static void vectormod3_shift(small *z,int len)
489{
490 int i;
491 for (i = len - 1;i > 0;--i) z[i] = z[i - 1];
492 z[0] = 0;
493}
494
495/*
496r = s^(-1) mod m, returning 0, if s is invertible mod m
497or returning -1 if s is not invertible mod m
498r,s are polys of degree <p
499m is x^p-x-1
500*/
501int r3_recip(small *r,const small *s)
502{
503 const int loops = 2*p + 1;
504 int loop;
505 small f[p + 1];
506 small g[p + 1];
Tim Rice00991152019-03-31 22:14:22 -0700507 small u[2*p + 2];
508 small v[2*p + 2];
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000509 small c;
510 int i;
511 int d = p;
512 int e = p;
513 int swapmask;
514
515 for (i = 2;i < p;++i) f[i] = 0;
516 f[0] = -1;
517 f[1] = -1;
518 f[p] = 1;
519 /* generalization: can initialize f to any polynomial m */
520 /* requirements: m has degree exactly p, nonzero constant coefficient */
521
522 for (i = 0;i < p;++i) g[i] = s[i];
523 g[p] = 0;
524
525 for (i = 0;i <= loops;++i) u[i] = 0;
526
527 v[0] = 1;
528 for (i = 1;i <= loops;++i) v[i] = 0;
529
530 loop = 0;
531 for (;;) {
532 /* e == -1 or d + e + loop <= 2*p */
533
534 /* f has degree p: i.e., f[p]!=0 */
535 /* f[i]==0 for i < p-d */
536
537 /* g has degree <=p (so it fits in p+1 coefficients) */
538 /* g[i]==0 for i < p-e */
539
540 /* u has degree <=loop (so it fits in loop+1 coefficients) */
541 /* u[i]==0 for i < p-d */
542 /* if invertible: u[i]==0 for i < loop-p (so can look at just p+1 coefficients) */
543
544 /* v has degree <=loop (so it fits in loop+1 coefficients) */
545 /* v[i]==0 for i < p-e */
546 /* v[i]==0 for i < loop-p (so can look at just p+1 coefficients) */
547
548 if (loop >= loops) break;
549
550 c = mod3_quotient(g[p],f[p]);
551
552 vectormod3_minusproduct(g,p + 1,g,f,c);
553 vectormod3_shift(g,p + 1);
554
555#ifdef SIMPLER
556 vectormod3_minusproduct(v,loops + 1,v,u,c);
557 vectormod3_shift(v,loops + 1);
558#else
559 if (loop < p) {
560 vectormod3_minusproduct(v,loop + 1,v,u,c);
561 vectormod3_shift(v,loop + 2);
562 } else {
563 vectormod3_minusproduct(v + loop - p,p + 1,v + loop - p,u + loop - p,c);
564 vectormod3_shift(v + loop - p,p + 2);
565 }
566#endif
567
568 e -= 1;
569
570 ++loop;
571
572 swapmask = smaller_mask_r3_recip(e,d) & mod3_nonzero_mask(g[p]);
573 swap(&e,&d,sizeof e,swapmask);
574 swap(f,g,(p + 1) * sizeof(small),swapmask);
575
576#ifdef SIMPLER
577 swap(u,v,(loops + 1) * sizeof(small),swapmask);
578#else
579 if (loop < p) {
580 swap(u,v,(loop + 1) * sizeof(small),swapmask);
581 } else {
582 swap(u + loop - p,v + loop - p,(p + 1) * sizeof(small),swapmask);
583 }
584#endif
585 }
586
587 c = mod3_reciprocal(f[p]);
588 vectormod3_product(r,p,u + p,c);
589 return smaller_mask_r3_recip(0,d);
590}
591
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000592/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/randomsmall.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000593/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
594
595
596static void small_random(small *g)
597{
598 int i;
599
600 for (i = 0;i < p;++i) {
601 crypto_uint32 r = small_random32();
602 g[i] = (small) (((1073741823 & r) * 3) >> 30) - 1;
603 }
604}
605
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000606/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/randomweightw.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000607/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
608
609
610static void small_random_weightw(small *f)
611{
612 crypto_int32 r[p];
613 int i;
614
615 for (i = 0;i < p;++i) r[i] = small_random32();
616 for (i = 0;i < w;++i) r[i] &= -2;
617 for (i = w;i < p;++i) r[i] = (r[i] & -3) | 1;
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000618 int32_sort(r,p);
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000619 for (i = 0;i < p;++i) f[i] = ((small) (r[i] & 3)) - 1;
620}
621
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000622/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000623/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
624
625
626static void rq_encode(unsigned char *c,const modq *f)
627{
628 crypto_int32 f0, f1, f2, f3, f4;
629 int i;
630
631 for (i = 0;i < p/5;++i) {
632 f0 = *f++ + qshift;
633 f1 = *f++ + qshift;
634 f2 = *f++ + qshift;
635 f3 = *f++ + qshift;
636 f4 = *f++ + qshift;
637 /* now want f0 + 6144*f1 + ... as a 64-bit integer */
638 f1 *= 3;
639 f2 *= 9;
640 f3 *= 27;
641 f4 *= 81;
642 /* now want f0 + f1<<11 + f2<<22 + f3<<33 + f4<<44 */
643 f0 += f1 << 11;
644 *c++ = f0; f0 >>= 8;
645 *c++ = f0; f0 >>= 8;
646 f0 += f2 << 6;
647 *c++ = f0; f0 >>= 8;
648 *c++ = f0; f0 >>= 8;
649 f0 += f3 << 1;
650 *c++ = f0; f0 >>= 8;
651 f0 += f4 << 4;
652 *c++ = f0; f0 >>= 8;
653 *c++ = f0; f0 >>= 8;
654 *c++ = f0;
655 }
656 /* XXX: using p mod 5 = 1 */
657 f0 = *f++ + qshift;
658 *c++ = f0; f0 >>= 8;
659 *c++ = f0;
660}
661
662static void rq_decode(modq *f,const unsigned char *c)
663{
664 crypto_uint32 c0, c1, c2, c3, c4, c5, c6, c7;
665 crypto_uint32 f0, f1, f2, f3, f4;
666 int i;
667
668 for (i = 0;i < p/5;++i) {
669 c0 = *c++;
670 c1 = *c++;
671 c2 = *c++;
672 c3 = *c++;
673 c4 = *c++;
674 c5 = *c++;
675 c6 = *c++;
676 c7 = *c++;
677
678 /* f0 + f1*6144 + f2*6144^2 + f3*6144^3 + f4*6144^4 */
679 /* = c0 + c1*256 + ... + c6*256^6 + c7*256^7 */
680 /* with each f between 0 and 4590 */
681
682 c6 += c7 << 8;
683 /* c6 <= 23241 = floor(4591*6144^4/2^48) */
684 /* f4 = (16/81)c6 + (1/1296)(c5+[0,1]) - [0,0.75] */
685 /* claim: 2^19 f4 < x < 2^19(f4+1) */
686 /* where x = 103564 c6 + 405(c5+1) */
687 /* proof: x - 2^19 f4 = (76/81)c6 + (37/81)c5 + 405 - (32768/81)[0,1] + 2^19[0,0.75] */
688 /* at least 405 - 32768/81 > 0 */
689 /* at most (76/81)23241 + (37/81)255 + 405 + 2^19 0.75 < 2^19 */
690 f4 = (103564*c6 + 405*(c5+1)) >> 19;
691
692 c5 += c6 << 8;
693 c5 -= (f4 * 81) << 4;
694 c4 += c5 << 8;
695
696 /* f0 + f1*6144 + f2*6144^2 + f3*6144^3 */
697 /* = c0 + c1*256 + c2*256^2 + c3*256^3 + c4*256^4 */
698 /* c4 <= 247914 = floor(4591*6144^3/2^32) */
699 /* f3 = (1/54)(c4+[0,1]) - [0,0.75] */
700 /* claim: 2^19 f3 < x < 2^19(f3+1) */
701 /* where x = 9709(c4+2) */
702 /* proof: x - 2^19 f3 = 19418 - (1/27)c4 - (262144/27)[0,1] + 2^19[0,0.75] */
703 /* at least 19418 - 247914/27 - 262144/27 > 0 */
704 /* at most 19418 + 2^19 0.75 < 2^19 */
705 f3 = (9709*(c4+2)) >> 19;
706
707 c4 -= (f3 * 27) << 1;
708 c3 += c4 << 8;
709 /* f0 + f1*6144 + f2*6144^2 */
710 /* = c0 + c1*256 + c2*256^2 + c3*256^3 */
711 /* c3 <= 10329 = floor(4591*6144^2/2^24) */
712 /* f2 = (4/9)c3 + (1/576)c2 + (1/147456)c1 + (1/37748736)c0 - [0,0.75] */
713 /* claim: 2^19 f2 < x < 2^19(f2+1) */
714 /* where x = 233017 c3 + 910(c2+2) */
715 /* proof: x - 2^19 f2 = 1820 + (1/9)c3 - (2/9)c2 - (32/9)c1 - (1/72)c0 + 2^19[0,0.75] */
716 /* at least 1820 - (2/9)255 - (32/9)255 - (1/72)255 > 0 */
717 /* at most 1820 + (1/9)10329 + 2^19 0.75 < 2^19 */
718 f2 = (233017*c3 + 910*(c2+2)) >> 19;
719
720 c2 += c3 << 8;
721 c2 -= (f2 * 9) << 6;
722 c1 += c2 << 8;
723 /* f0 + f1*6144 */
724 /* = c0 + c1*256 */
725 /* c1 <= 110184 = floor(4591*6144/2^8) */
726 /* f1 = (1/24)c1 + (1/6144)c0 - (1/6144)f0 */
727 /* claim: 2^19 f1 < x < 2^19(f1+1) */
728 /* where x = 21845(c1+2) + 85 c0 */
729 /* proof: x - 2^19 f1 = 43690 - (1/3)c1 - (1/3)c0 + 2^19 [0,0.75] */
730 /* at least 43690 - (1/3)110184 - (1/3)255 > 0 */
731 /* at most 43690 + 2^19 0.75 < 2^19 */
732 f1 = (21845*(c1+2) + 85*c0) >> 19;
733
734 c1 -= (f1 * 3) << 3;
735 c0 += c1 << 8;
736 f0 = c0;
737
738 *f++ = modq_freeze(f0 + q - qshift);
739 *f++ = modq_freeze(f1 + q - qshift);
740 *f++ = modq_freeze(f2 + q - qshift);
741 *f++ = modq_freeze(f3 + q - qshift);
742 *f++ = modq_freeze(f4 + q - qshift);
743 }
744
745 c0 = *c++;
746 c1 = *c++;
747 c0 += c1 << 8;
748 *f++ = modq_freeze(c0 + q - qshift);
749}
750
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000751/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_mult.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000752/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
753
754
755static void rq_mult(modq *h,const modq *f,const small *g)
756{
757 modq fg[p + p - 1];
758 modq result;
759 int i, j;
760
761 for (i = 0;i < p;++i) {
762 result = 0;
763 for (j = 0;j <= i;++j)
764 result = modq_plusproduct(result,f[j],g[i - j]);
765 fg[i] = result;
766 }
767 for (i = p;i < p + p - 1;++i) {
768 result = 0;
769 for (j = i - p + 1;j < p;++j)
770 result = modq_plusproduct(result,f[j],g[i - j]);
771 fg[i] = result;
772 }
773
774 for (i = p + p - 2;i >= p;--i) {
775 fg[i - p] = modq_sum(fg[i - p],fg[i]);
776 fg[i - p + 1] = modq_sum(fg[i - p + 1],fg[i]);
777 }
778
779 for (i = 0;i < p;++i)
780 h[i] = fg[i];
781}
782
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000783/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_recip3.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000784/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
785
786
787/* caller must ensure that x-y does not overflow */
788static int smaller_mask_rq_recip3(int x,int y)
789{
790 return (x - y) >> 31;
791}
792
793static void vectormodq_product(modq *z,int len,const modq *x,const modq c)
794{
795 int i;
796 for (i = 0;i < len;++i) z[i] = modq_product(x[i],c);
797}
798
799static void vectormodq_minusproduct(modq *z,int len,const modq *x,const modq *y,const modq c)
800{
801 int i;
802 for (i = 0;i < len;++i) z[i] = modq_minusproduct(x[i],y[i],c);
803}
804
805static void vectormodq_shift(modq *z,int len)
806{
807 int i;
808 for (i = len - 1;i > 0;--i) z[i] = z[i - 1];
809 z[0] = 0;
810}
811
812/*
813r = (3s)^(-1) mod m, returning 0, if s is invertible mod m
814or returning -1 if s is not invertible mod m
815r,s are polys of degree <p
816m is x^p-x-1
817*/
818int rq_recip3(modq *r,const small *s)
819{
820 const int loops = 2*p + 1;
821 int loop;
822 modq f[p + 1];
823 modq g[p + 1];
Tim Rice00991152019-03-31 22:14:22 -0700824 modq u[2*p + 2];
825 modq v[2*p + 2];
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000826 modq c;
827 int i;
828 int d = p;
829 int e = p;
830 int swapmask;
831
832 for (i = 2;i < p;++i) f[i] = 0;
833 f[0] = -1;
834 f[1] = -1;
835 f[p] = 1;
836 /* generalization: can initialize f to any polynomial m */
837 /* requirements: m has degree exactly p, nonzero constant coefficient */
838
839 for (i = 0;i < p;++i) g[i] = 3 * s[i];
840 g[p] = 0;
841
842 for (i = 0;i <= loops;++i) u[i] = 0;
843
844 v[0] = 1;
845 for (i = 1;i <= loops;++i) v[i] = 0;
846
847 loop = 0;
848 for (;;) {
849 /* e == -1 or d + e + loop <= 2*p */
850
851 /* f has degree p: i.e., f[p]!=0 */
852 /* f[i]==0 for i < p-d */
853
854 /* g has degree <=p (so it fits in p+1 coefficients) */
855 /* g[i]==0 for i < p-e */
856
857 /* u has degree <=loop (so it fits in loop+1 coefficients) */
858 /* u[i]==0 for i < p-d */
859 /* if invertible: u[i]==0 for i < loop-p (so can look at just p+1 coefficients) */
860
861 /* v has degree <=loop (so it fits in loop+1 coefficients) */
862 /* v[i]==0 for i < p-e */
863 /* v[i]==0 for i < loop-p (so can look at just p+1 coefficients) */
864
865 if (loop >= loops) break;
866
867 c = modq_quotient(g[p],f[p]);
868
869 vectormodq_minusproduct(g,p + 1,g,f,c);
870 vectormodq_shift(g,p + 1);
871
872#ifdef SIMPLER
873 vectormodq_minusproduct(v,loops + 1,v,u,c);
874 vectormodq_shift(v,loops + 1);
875#else
876 if (loop < p) {
877 vectormodq_minusproduct(v,loop + 1,v,u,c);
878 vectormodq_shift(v,loop + 2);
879 } else {
880 vectormodq_minusproduct(v + loop - p,p + 1,v + loop - p,u + loop - p,c);
881 vectormodq_shift(v + loop - p,p + 2);
882 }
883#endif
884
885 e -= 1;
886
887 ++loop;
888
889 swapmask = smaller_mask_rq_recip3(e,d) & modq_nonzero_mask(g[p]);
890 swap(&e,&d,sizeof e,swapmask);
891 swap(f,g,(p + 1) * sizeof(modq),swapmask);
892
893#ifdef SIMPLER
894 swap(u,v,(loops + 1) * sizeof(modq),swapmask);
895#else
896 if (loop < p) {
897 swap(u,v,(loop + 1) * sizeof(modq),swapmask);
898 } else {
899 swap(u + loop - p,v + loop - p,(p + 1) * sizeof(modq),swapmask);
900 }
901#endif
902 }
903
904 c = modq_reciprocal(f[p]);
905 vectormodq_product(r,p,u + p,c);
906 return smaller_mask_rq_recip3(0,d);
907}
908
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000909/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_round3.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000910/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
911
912
913static void rq_round3(modq *h,const modq *f)
914{
915 int i;
916
917 for (i = 0;i < p;++i)
918 h[i] = ((21846 * (f[i] + 2295) + 32768) >> 16) * 3 - 2295;
919}
920
djm@openbsd.org533cfb02019-01-21 22:18:24 +0000921/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_rounded.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +0000922/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
923
924
925static void rq_encoderounded(unsigned char *c,const modq *f)
926{
927 crypto_int32 f0, f1, f2;
928 int i;
929
930 for (i = 0;i < p/3;++i) {
931 f0 = *f++ + qshift;
932 f1 = *f++ + qshift;
933 f2 = *f++ + qshift;
934 f0 = (21846 * f0) >> 16;
935 f1 = (21846 * f1) >> 16;
936 f2 = (21846 * f2) >> 16;
937 /* now want f0 + f1*1536 + f2*1536^2 as a 32-bit integer */
938 f2 *= 3;
939 f1 += f2 << 9;
940 f1 *= 3;
941 f0 += f1 << 9;
942 *c++ = f0; f0 >>= 8;
943 *c++ = f0; f0 >>= 8;
944 *c++ = f0; f0 >>= 8;
945 *c++ = f0;
946 }
947 /* XXX: using p mod 3 = 2 */
948 f0 = *f++ + qshift;
949 f1 = *f++ + qshift;
950 f0 = (21846 * f0) >> 16;
951 f1 = (21846 * f1) >> 16;
952 f1 *= 3;
953 f0 += f1 << 9;
954 *c++ = f0; f0 >>= 8;
955 *c++ = f0; f0 >>= 8;
956 *c++ = f0;
957}
958
959static void rq_decoderounded(modq *f,const unsigned char *c)
960{
961 crypto_uint32 c0, c1, c2, c3;
962 crypto_uint32 f0, f1, f2;
963 int i;
964
965 for (i = 0;i < p/3;++i) {
966 c0 = *c++;
967 c1 = *c++;
968 c2 = *c++;
969 c3 = *c++;
970
971 /* f0 + f1*1536 + f2*1536^2 */
972 /* = c0 + c1*256 + c2*256^2 + c3*256^3 */
973 /* with each f between 0 and 1530 */
974
975 /* f2 = (64/9)c3 + (1/36)c2 + (1/9216)c1 + (1/2359296)c0 - [0,0.99675] */
976 /* claim: 2^21 f2 < x < 2^21(f2+1) */
977 /* where x = 14913081*c3 + 58254*c2 + 228*(c1+2) */
978 /* proof: x - 2^21 f2 = 456 - (8/9)c0 + (4/9)c1 - (2/9)c2 + (1/9)c3 + 2^21 [0,0.99675] */
979 /* at least 456 - (8/9)255 - (2/9)255 > 0 */
980 /* at most 456 + (4/9)255 + (1/9)255 + 2^21 0.99675 < 2^21 */
981 f2 = (14913081*c3 + 58254*c2 + 228*(c1+2)) >> 21;
982
983 c2 += c3 << 8;
984 c2 -= (f2 * 9) << 2;
985 /* f0 + f1*1536 */
986 /* = c0 + c1*256 + c2*256^2 */
987 /* c2 <= 35 = floor((1530+1530*1536)/256^2) */
988 /* f1 = (128/3)c2 + (1/6)c1 + (1/1536)c0 - (1/1536)f0 */
989 /* claim: 2^21 f1 < x < 2^21(f1+1) */
990 /* where x = 89478485*c2 + 349525*c1 + 1365*(c0+1) */
991 /* proof: x - 2^21 f1 = 1365 - (1/3)c2 - (1/3)c1 - (1/3)c0 + (4096/3)f0 */
992 /* at least 1365 - (1/3)35 - (1/3)255 - (1/3)255 > 0 */
993 /* at most 1365 + (4096/3)1530 < 2^21 */
994 f1 = (89478485*c2 + 349525*c1 + 1365*(c0+1)) >> 21;
995
996 c1 += c2 << 8;
997 c1 -= (f1 * 3) << 1;
998
999 c0 += c1 << 8;
1000 f0 = c0;
1001
1002 *f++ = modq_freeze(f0 * 3 + q - qshift);
1003 *f++ = modq_freeze(f1 * 3 + q - qshift);
1004 *f++ = modq_freeze(f2 * 3 + q - qshift);
1005 }
1006
1007 c0 = *c++;
1008 c1 = *c++;
1009 c2 = *c++;
1010
1011 f1 = (89478485*c2 + 349525*c1 + 1365*(c0+1)) >> 21;
1012
1013 c1 += c2 << 8;
1014 c1 -= (f1 * 3) << 1;
1015
1016 c0 += c1 << 8;
1017 f0 = c0;
1018
1019 *f++ = modq_freeze(f0 * 3 + q - qshift);
1020 *f++ = modq_freeze(f1 * 3 + q - qshift);
1021}
1022
djm@openbsd.org533cfb02019-01-21 22:18:24 +00001023/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/small.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +00001024/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
1025
1026
1027/* XXX: these functions rely on p mod 4 = 1 */
1028
1029/* all coefficients in -1, 0, 1 */
1030static void small_encode(unsigned char *c,const small *f)
1031{
1032 small c0;
1033 int i;
1034
1035 for (i = 0;i < p/4;++i) {
1036 c0 = *f++ + 1;
1037 c0 += (*f++ + 1) << 2;
1038 c0 += (*f++ + 1) << 4;
1039 c0 += (*f++ + 1) << 6;
1040 *c++ = c0;
1041 }
1042 c0 = *f++ + 1;
1043 *c++ = c0;
1044}
1045
1046static void small_decode(small *f,const unsigned char *c)
1047{
1048 unsigned char c0;
1049 int i;
1050
1051 for (i = 0;i < p/4;++i) {
1052 c0 = *c++;
1053 *f++ = ((small) (c0 & 3)) - 1; c0 >>= 2;
1054 *f++ = ((small) (c0 & 3)) - 1; c0 >>= 2;
1055 *f++ = ((small) (c0 & 3)) - 1; c0 >>= 2;
1056 *f++ = ((small) (c0 & 3)) - 1;
1057 }
1058 c0 = *c++;
1059 *f++ = ((small) (c0 & 3)) - 1;
1060}
1061
djm@openbsd.org533cfb02019-01-21 22:18:24 +00001062/* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/swap.c */
djm@openbsd.orgdfd59162019-01-21 10:20:12 +00001063/* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
1064
1065
1066static void swap(void *x,void *y,int bytes,int mask)
1067{
1068 int i;
1069 char xi, yi, c, t;
1070
1071 c = mask;
1072
1073 for (i = 0;i < bytes;++i) {
1074 xi = i[(char *) x];
1075 yi = i[(char *) y];
1076 t = c & (xi ^ yi);
1077 xi ^= t;
1078 yi ^= t;
1079 i[(char *) x] = xi;
1080 i[(char *) y] = yi;
1081 }
1082}
1083