[BACK]Return to sntrup4591761.c CVS log [TXT][DIR] Up to [local] / src / usr.bin / ssh

Annotation of src/usr.bin/ssh/sntrup4591761.c, Revision 1.3

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