crypto/openssh/sntrup4591761.c

   1 /*  $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
  11 #include "includes.h"
  12
  13 #include <string.h>
  14 #include "crypto_api.h"
  15
  16 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/int32_sort.h */
  17 #ifndef int32_sort_h
  18 #define int32_sort_h
  19
  20
  21 static 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
  29 static void minmax(crypto_int32 *x,crypto_int32 *y)
  30 {
  31   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
  43 static void int32_sort(crypto_int32 *x,int n)
  44 {
  45   int top,p,q,i;
  46
  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))
  54         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);
  59   }
  60 }
  61
  62 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/small.h */
  63 #ifndef small_h
  64 #define small_h
  65
  66
  67 typedef crypto_int8 small;
  68
  69 static void small_encode(unsigned char *,const small *);
  70
  71 static void small_decode(small *,const unsigned char *);
  72
  73
  74 static void small_random(small *);
  75
  76 static void small_random_weightw(small *);
  77
  78 #endif
  79
  80 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/mod3.h */
  81 #ifndef mod3_h
  82 #define mod3_h
  83
  84
  85 /* -1 if x is nonzero, 0 otherwise */
  86 static 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 */
  93 static 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
 100 static 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
 108 static 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
 116 static inline small mod3_product(small a,small b)
 117 {
 118   return a * b;
 119 }
 120
 121 static 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
 128 static inline small mod3_reciprocal(small a1)
 129 {
 130   return a1;
 131 }
 132
 133 static inline small mod3_quotient(small num,small den)
 134 {
 135   return mod3_product(num,mod3_reciprocal(den));
 136 }
 137
 138 #endif
 139
 140 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/modq.h */
 141 #ifndef modq_h
 142 #define modq_h
 143
 144
 145 typedef crypto_int16 modq;
 146
 147 /* -1 if x is nonzero, 0 otherwise */
 148 static 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 */
 158 static 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
 165 static 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
 173 static 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
 181 static 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
 188 static inline modq modq_square(modq a)
 189 {
 190   crypto_int32 A = a;
 191   return modq_freeze(A * A);
 192 }
 193
 194 static 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
 201 static 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
 223 static inline modq modq_quotient(modq num,modq den)
 224 {
 225   return modq_product(num,modq_reciprocal(den));
 226 }
 227
 228 #endif
 229
 230 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/params.h */
 231 #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
 246 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/r3.h */
 247 #ifndef r3_h
 248 #define r3_h
 249
 250
 251 static void r3_mult(small *,const small *,const small *);
 252
 253 extern int r3_recip(small *,const small *);
 254
 255 #endif
 256
 257 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq.h */
 258 #ifndef rq_h
 259 #define rq_h
 260
 261
 262 static void rq_encode(unsigned char *,const modq *);
 263
 264 static void rq_decode(modq *,const unsigned char *);
 265
 266 static void rq_encoderounded(unsigned char *,const modq *);
 267
 268 static void rq_decoderounded(modq *,const unsigned char *);
 269
 270 static void rq_round3(modq *,const modq *);
 271
 272 static void rq_mult(modq *,const modq *,const small *);
 273
 274 int rq_recip3(modq *,const small *);
 275
 276 #endif
 277
 278 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/swap.h */
 279 #ifndef swap_h
 280 #define swap_h
 281
 282 static void swap(void *,void *,int,int);
 283
 284 #endif
 285
 286 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/dec.c */
 287 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 288
 289 #ifdef KAT
 290 #endif
 291
 292
 293 int 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
 352 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/enc.c */
 353 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 354
 355 #ifdef KAT
 356 #endif
 357
 358
 359 int 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
 398 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/keypair.c */
 399 /* 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
 409 int 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
 434 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/r3_mult.c */
 435 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 436
 437
 438 static 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
 466 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/r3_recip.c */
 467 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 468
 469
 470 /* caller must ensure that x-y does not overflow */
 471 static int smaller_mask_r3_recip(int x,int y)
 472 {
 473   return (x - y) >> 31;
 474 }
 475
 476 static 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
 482 static 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
 488 static 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 /*
 496 r = s^(-1) mod m, returning 0, if s is invertible mod m
 497 or returning -1 if s is not invertible mod m
 498 r,s are polys of degree <p
 499 m is x^p-x-1
 500 */
 501 int 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];
 507   small u[2*p + 2];
 508   small v[2*p + 2];
 509   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
 592 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/randomsmall.c */
 593 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 594
 595
 596 static 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
 606 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/randomweightw.c */
 607 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 608
 609
 610 static 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;
 618   int32_sort(r,p);
 619   for (i = 0;i < p;++i) f[i] = ((small) (r[i] & 3)) - 1;
 620 }
 621
 622 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq.c */
 623 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 624
 625
 626 static 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
 662 static 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
 751 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_mult.c */
 752 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 753
 754
 755 static 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
 783 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_recip3.c */
 784 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 785
 786
 787 /* caller must ensure that x-y does not overflow */
 788 static int smaller_mask_rq_recip3(int x,int y)
 789 {
 790   return (x - y) >> 31;
 791 }
 792
 793 static 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
 799 static 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
 805 static 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 /*
 813 r = (3s)^(-1) mod m, returning 0, if s is invertible mod m
 814 or returning -1 if s is not invertible mod m
 815 r,s are polys of degree <p
 816 m is x^p-x-1
 817 */
 818 int 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];
 824   modq u[2*p + 2];
 825   modq v[2*p + 2];
 826   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
 909 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_round3.c */
 910 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 911
 912
 913 static 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
 921 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/rq_rounded.c */
 922 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
 923
 924
 925 static 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
 959 static 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
1023 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/small.c */
1024 /* 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 */
1030 static 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
1046 static 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
1062 /* from libpqcrypto-20180314/crypto_kem/sntrup4591761/ref/swap.c */
1063 /* See https://ntruprime.cr.yp.to/software.html for detailed documentation. */
1064
1065
1066 static 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