Packages
ntru_elixir
0.0.0
libntru wrapper for elixir. NTRU is a post quantom cryptography algorithm.
Current section
Files
Jump to
Current section
Files
libntru/src/poly.c
#include <stdlib.h>
#include <string.h>
#include <alloca.h>
#include "poly.h"
#include "poly_ssse3.h"
#include "poly_avx2.h"
#include "rand.h"
#include "err.h"
#include "arith.h"
#include "encparams.h"
#include "ntru_endian.h"
#define NTRU_KARATSUBA_THRESH_16 40
#define NTRU_KARATSUBA_THRESH_64 120
/***************************************
* NTRU Prime *
***************************************/
uint8_t ntruprime_mult_poly(NtruIntPoly *a, NtruIntPoly *b, NtruIntPoly *c, uint16_t modulus) {
uint16_t N = a->N;
if (N != b->N)
return 0;
c->N = N;
memset(&c->coeffs, 0, N * sizeof c->coeffs[0]);
uint16_t k;
for (k=0; k<N; k++) {
uint64_t ck1 = 0;
uint64_t i;
for (i=0; i<=k; i++)
ck1 += ((uint64_t)a->coeffs[i]) * b->coeffs[k-i];
uint64_t ck2 = 0;
for (i=k+1; i<N; i++)
ck2 += ((uint64_t)a->coeffs[i]) * b->coeffs[k+N-i];
uint64_t ck = c->coeffs[k] + ck1 + ck2;
c->coeffs[k] = ck % modulus;
if (k < N) {
uint64_t ck = c->coeffs[k+1] + ck2;
c->coeffs[k+1] = ck % modulus;
}
}
return 1;
}
uint8_t ntruprime_rand_tern(uint16_t N, NtruIntPoly *poly, NtruRandContext *rand_ctx) {
poly->N = N;
uint16_t i;
uint32_t arr[N];
if (ntru_rand_generate((uint8_t*)arr, N*sizeof(arr[0]), rand_ctx) != NTRU_SUCCESS)
return 0;
for (i=0; i<N; i++)
/* produces almost uniform output, there is a 1 in 2^32 chance of it being biased */
poly->coeffs[i] = arr[i] % 3;
return 1;
}
uint8_t ntruprime_rand_tern_t(uint16_t N, uint16_t t, NtruIntPoly *poly, NtruRandContext *rand_ctx) {
poly->N = N;
uint8_t rand_bits[(2*t+7)/8];
if (ntru_rand_generate((uint8_t*)rand_bits, sizeof rand_bits, rand_ctx) != NTRU_SUCCESS)
return 0;
uint16_t i;
for (i=0; i<2*t; i++) {
uint8_t r = (rand_bits[i/8] >> (i%8)) & 1;
poly->coeffs[i] = r ? 1 : 2;
}
for (; i<N; i++)
poly->coeffs[i] = 0;
uint32_t rand_indices[N];
if (ntru_rand_generate((uint8_t*)rand_indices, N*sizeof(rand_indices[0]), rand_ctx) != NTRU_SUCCESS)
return 0;
/* Fisher-Yates shuffle */
uint32_t rand_idx = 0;
for (i=N-1; i>0; i--) {
/* slightly biased; for N=739, P(poly contains a biased coeff) = 1/10798 */
uint32_t j = rand_indices[rand_idx] % (i+1);
uint16_t temp = poly->coeffs[i];
poly->coeffs[i] = poly->coeffs[j];
poly->coeffs[j] = temp;
rand_idx++;
}
return 1;
}
/* Zeros a polynomial and sets the number of coefficients */
void ntruprime_zero(NtruIntPoly *a, uint16_t N) {
a->N = N;
memset(&a->coeffs, 0, N * sizeof a->coeffs[0]);
}
/* Reduces a NtruIntPoly modulo x^N-x-1, where N = a->N. */
void ntruprime_reduce(NtruIntPoly *a, NtruIntPoly *b, uint16_t modulus) {
uint16_t N = a->N - 1;
memcpy(&b->coeffs, &a->coeffs, N * sizeof a->coeffs[0]);
b->coeffs[0] = (((uint64_t)b->coeffs[0]) + a->coeffs[N]) % modulus;
b->coeffs[1] = (((uint64_t)b->coeffs[1]) + a->coeffs[N]) % modulus;
b->N = N;
}
/* Copies a NtruIntPoly to another NtruIntPoly */
void ntruprime_copy(NtruIntPoly *a, NtruIntPoly *b) {
memcpy(&b->coeffs, &a->coeffs, a->N * sizeof a->coeffs[0]);
b->N = a->N;
}
/* Tests whether a polynomial has all zero coefficients */
uint8_t ntruprime_equals0(NtruIntPoly *a) {
uint16_t i;
for (i=0; i<a->N; i++)
if (a->coeffs[i] != 0)
return 0;
return 1;
}
/* Returns the degree of a polynomial. Returns 0 if the polynomial equals 0. */
uint16_t ntruprime_deg(NtruIntPoly *a) {
uint16_t i;
for (i=a->N-1; i>0; i--)
if (a->coeffs[i] != 0)
return i;
return 0;
}
/* Multiplies a polynomial by an integer */
void ntruprime_mult_mod(NtruIntPoly *a, uint16_t factor, uint16_t modulus) {
uint16_t i;
for (i=0; i<a->N; i++)
a->coeffs[i] = (((uint64_t)a->coeffs[i]) * factor) % modulus;
}
/** Subtracts b*u from a; a and b are assumed to be reduced mod modulus. */
void ntruprime_subtract_multiple(NtruIntPoly *a, NtruIntPoly *b, uint16_t u, uint16_t modulus) {
uint16_t N = a->N;
if (b->N > N)
N = b->N;
uint16_t i;
for (i=0; i<N; i++) {
uint64_t ai = a->coeffs[i];
ai += ((uint64_t)u) * (modulus-b->coeffs[i]);
a->coeffs[i] = ai % modulus;
}
}
/* Multiplies a polynomial by x^(-1) in (Z/qZ)[x][x^p-x-1] where p=a->N, q=modulus */
void ntruprime_div_x(NtruIntPoly *a, uint16_t modulus) {
uint16_t N = a->N;
uint64_t a0 = a->coeffs[0];
memmove(&a->coeffs[0], &a->coeffs[1], (N-1) * sizeof a->coeffs[0]);
a->coeffs[N-1] = a0;
a->coeffs[0] = (((uint64_t)a->coeffs[0])-a0+modulus) % modulus;
}
/* Calculates the multiplicative inverse of a mod modulus
using the extended Euclidean algorithm. */
uint16_t ntruprime_inv_int(uint16_t a, uint16_t modulus) {
int16_t x = 0;
int16_t lastx = 1;
int16_t y = 1;
int16_t lasty = 0;
int16_t b = modulus;
while (b != 0) {
int16_t quotient = a / b;
int16_t temp = a;
a = b;
b = temp % b;
temp = x;
x = lastx - quotient*x;
lastx = temp;
temp = y;
y = lasty - quotient*y;
lasty = temp;
}
if (lastx < 0)
lastx += modulus;
return lastx;
}
/* not constant time! */
uint8_t ntruprime_inv_poly(NtruIntPoly *a, NtruIntPoly *inv, uint16_t modulus) {
uint16_t N = a->N;
uint16_t k = 0;
NtruIntPoly *b = alloca(sizeof *b);
ntruprime_zero(b, N+1);
b->coeffs[0] = 1;
NtruIntPoly *c = alloca(sizeof *c);
ntruprime_zero(c, N+1);
/* f = a */
NtruIntPoly *f = alloca(sizeof *f);
f->N = N + 1;
ntruprime_copy(a, f);
f->coeffs[N] = 0;
/* g = x^p - x - 1 */
NtruIntPoly *g = alloca(sizeof *g);
ntruprime_zero(g, N+1);
g->coeffs[0] = modulus - 1;
g->coeffs[1] = modulus - 1;
g->coeffs[N] = 1;
for (;;) {
while (f->coeffs[0] == 0) {
uint16_t i;
/* f(x) = f(x) / x */
for (i=1; i<=N; i++)
f->coeffs[i-1] = f->coeffs[i];
f->coeffs[N] = 0;
/* c(x) = c(x) * x */
for (i=N-1; i>0; i--)
c->coeffs[i] = c->coeffs[i-1];
c->coeffs[0] = 0;
k++;
if (ntruprime_equals0(f)) /* not invertible */
return 0;
}
if (ntruprime_deg(f) == 0) {
uint16_t f0_inv = ntruprime_inv_int(f->coeffs[0], modulus);
ntruprime_mult_mod(b, f0_inv, modulus); /* b = b * f[0]^(-1) */
ntruprime_reduce(b, inv, modulus);
uint16_t i;
/* b = b * x^(-k) */
for (i=0; i<k; i++)
ntruprime_div_x(inv, modulus);
return 1;
}
if (ntruprime_deg(f) < ntruprime_deg(g)) {
/* exchange f and g */
NtruIntPoly *temp = f;
f = g;
g = temp;
/* exchange b and c */
temp = b;
b = c;
c = temp;
}
/* u = f[0] * g[0]^(-1) */
uint16_t g0_inv = ntruprime_inv_int(g->coeffs[0], modulus);
uint16_t u = (f->coeffs[0] * (uint64_t)g0_inv) % modulus;
ntruprime_subtract_multiple(f, g, u, modulus); /* f = f - u * g */
ntruprime_subtract_multiple(b, c, u, modulus); /* b = b - u * c */
}
}
/***************************************
* NTRUEncrypt *
***************************************/
uint8_t ntru_num_bits(uint16_t n) {
uint8_t b = 1;
while (n >>= 1)
b++;
return b;
}
uint8_t ntru_rand_tern(uint16_t N, uint16_t num_ones, uint16_t num_neg_ones, NtruTernPoly *poly, NtruRandContext *rand_ctx) {
int16_t coeffs[N];
memset(&coeffs, 0, N * sizeof coeffs[0]);
uint16_t rand_len = num_ones + num_neg_ones + 10; /* 10 more to avoid calling the RNG again, for up to 10 collisions */
uint16_t rand_data[rand_len];
if (ntru_rand_generate((uint8_t*)rand_data, rand_len*2, rand_ctx) != NTRU_SUCCESS)
return 0;
uint16_t r_idx = 0; /* index into rand_data */
uint16_t bits = ntru_num_bits(N);
uint16_t i = 0;
while (i < num_ones) {
uint16_t r = htole16(rand_data[r_idx]) >> (8*sizeof r - bits); /* 0 <= r < 2^bits */
r_idx++;
/* refill rand_data if we run out */
if (r_idx >= rand_len) {
if (ntru_rand_generate((uint8_t*)rand_data, rand_len*2, rand_ctx) != NTRU_SUCCESS)
return 0;
r_idx = 0;
}
if (r<N && !coeffs[r]) {
poly->ones[i] = r;
coeffs[r] = 1;
i++;
}
}
i = 0;
while (i < num_neg_ones) {
uint16_t r = htole16(rand_data[r_idx]) >> (8*sizeof r - bits); /* 0 <= r < 2^bits */
r_idx++;
/* refill rand_data if we run out */
if (r_idx >= rand_len) {
if (ntru_rand_generate((uint8_t*)rand_data, rand_len*2, rand_ctx) != NTRU_SUCCESS)
return 0;
r_idx = 0;
}
if (r<N && !coeffs[r]) {
poly->neg_ones[i] = r;
coeffs[r] = -1;
i++;
}
}
poly->N = N;
poly->num_ones = num_ones;
poly->num_neg_ones = num_neg_ones;
return 1;
}
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
uint8_t ntru_rand_prod(uint16_t N, uint16_t df1, uint16_t df2, uint16_t df3_ones, uint16_t df3_neg_ones, NtruProdPoly *poly, NtruRandContext *rand_ctx) {
poly->N = N;
uint8_t result = ntru_rand_tern(N, df1, df1, &poly->f1, rand_ctx);
result &= ntru_rand_tern(N, df2, df2, &poly->f2, rand_ctx);
result &= ntru_rand_tern(N, df3_ones, df3_neg_ones, &poly->f3, rand_ctx);
return result;
}
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
void ntru_add(NtruIntPoly *a, NtruIntPoly *b) {
uint16_t i;
for (i=0; i<b->N; i++)
a->coeffs[i] += b->coeffs[i];
}
void ntru_add_mod2_32(uint32_t *a, uint32_t *b, uint16_t len) {
uint16_t i;
for (i=0; i<len; i++)
a[i] ^= b[i];
}
void ntru_add_mod2_64(uint64_t *a, uint64_t *b, uint16_t len) {
uint16_t i;
for (i=0; i<len; i++)
a[i] ^= b[i];
}
void ntru_sub(NtruIntPoly *a, NtruIntPoly *b) {
uint16_t i;
for (i=0; i<b->N; i++)
a->coeffs[i] -= b->coeffs[i];
}
void ntru_neg_mod(NtruIntPoly *a, uint16_t modulus) {
uint16_t i;
for (i=0; i<a->N; i++)
a->coeffs[i] = modulus - a->coeffs[i];
}
void ntru_mult_int_16_base(int16_t *a, int16_t *b, int16_t *c, uint16_t len, uint16_t N, uint16_t mod_mask) {
memset(c, 0, 2*(2*len-1)); /* only needed if N < NTRU_KARATSUBA_THRESH_16 */
uint16_t c_idx = 0;
uint16_t k;
for (k=0; k<2*len-1; k++) {
int16_t ck = 0;
uint16_t i;
int16_t istart = k - len + 1;
if (istart < 0)
istart = 0;
int16_t iend = k + 1;
if (iend > len)
iend = len;
int16_t a_idx = k - istart;
for (i=istart; i<iend; i++) {
ck += b[i] * a[a_idx];
a_idx--;
if (a_idx < 0)
a_idx = len - 1;
}
c[c_idx] += ck;
c_idx++;
if (c_idx >= N)
c_idx = 0;
}
}
void ntru_mult_karatsuba_16(int16_t *a, int16_t *b, int16_t *c, uint16_t len, uint16_t N) {
if (len < NTRU_KARATSUBA_THRESH_16)
ntru_mult_int_16_base(a, b, c, len, N, -1);
else {
uint16_t len2 = len / 2;
int16_t z0[NTRU_INT_POLY_SIZE];
int16_t z1[NTRU_INT_POLY_SIZE];
int16_t z2[NTRU_INT_POLY_SIZE];
/* z0, z2 */
ntru_mult_karatsuba_16(a, b, z0, len2, N);
ntru_mult_karatsuba_16(a+len2, b+len2, z2, len-len2, N);
/* z1 */
int16_t lh1[NTRU_INT_POLY_SIZE];
int16_t lh2[NTRU_INT_POLY_SIZE];
uint16_t i;
for (i=0; i<len2; i++) {
lh1[i] = a[i] + a[len2+i];
lh2[i] = b[i] + b[len2+i];
}
if (len%2 != 0) {
lh1[len-len2-1] = a[len-1];
lh2[len-len2-1] = b[len-1];
}
ntru_mult_karatsuba_16(lh1, lh2, z1, len-len2, N);
for (i=0; i<2*len2-1; i++)
z1[i] -= z0[i];
z1[len] = 0;
for (i=0; i<2*(len-len2)-1; i++)
z1[i] -= z2[i];
/* c */
memset(c, 0, NTRU_INT_POLY_SIZE*2);
memcpy(c, z0, 2*(2*len2-1)); /* 2*len2-1 coefficients */
uint16_t c_idx = len2;
for (i=0; i<2*(len-len2)-1; i++) {
c[c_idx] += z1[i];
c_idx++;
if (c_idx >= N)
c_idx = 0;
}
c_idx = 2 * len2;
if (c_idx >= N)
c_idx = 0;
for (i=0; i<2*(len-len2)-1; i++) {
c[c_idx] += z2[i];
c_idx++;
if (c_idx >= N)
c_idx = 0;
}
}
}
uint8_t ntru_mult_int_16(NtruIntPoly *a, NtruIntPoly *b, NtruIntPoly *c, uint16_t mod_mask) {
uint16_t N = a->N;
if (N != b->N)
return 0;
c->N = N;
ntru_mult_karatsuba_16((int16_t*)&a->coeffs, (int16_t*)&b->coeffs, (int16_t*)&c->coeffs, N, N);
ntru_mod_mask(c, mod_mask);
return 1;
}
void ntru_mult_int_64_base(int16_t *a, int16_t *b, int16_t *c, uint16_t len, uint16_t N, uint16_t mod_mask) {
uint16_t N2 = (len+1) / 2;
uint64_t mod_mask_64 = mod_mask + (mod_mask<<25);
/* make 64-bit versions of a and b */
uint64_t a64[N2];
uint64_t b64[N2];
uint16_t i;
for (i=0; i<len/2; i++) {
a64[i] = (a[2*i]&mod_mask) + (((uint64_t)(a[2*i+1]&mod_mask))<<25);
b64[i] = (b[2*i]&mod_mask) + (((uint64_t)(b[2*i+1]&mod_mask))<<25);
}
if (len%2 == 1) {
a64[N2-1] = a[len-1] & mod_mask;
b64[N2-1] = b[len-1] & mod_mask;
}
/* multiply a64 by b64 */
uint16_t clen = 2 * N2; /* double capacity for intermediate result */
uint64_t c64[clen];
memset(&c64, 0, clen*8);
uint16_t overflow_ctr_start = (1<<(25-ntru_log2(mod_mask+1))) - 1;
uint16_t overflow_ctr_rem = overflow_ctr_start;
for (i=0; i<N2; i++) {
uint64_t j;
for (j=0; j<N2; j++) {
uint64_t ck = a64[i] * b64[j];
c64[i+j] += ck & mod_mask_64;
c64[i+j+1] += ck >> 50;
overflow_ctr_rem--;
if (!overflow_ctr_rem) {
uint16_t k;
for (k=0; k<clen; k++)
c64[k] &= mod_mask_64;
overflow_ctr_rem = overflow_ctr_start;
}
}
}
/* transform c64 into NtruIntPoly representation */
memset(c, 0, 2*(2*len-1));
uint16_t k = 0;
for (i=0; i<clen; i++) {
c[k] += c64[i];
if (++k >= N)
k = 0;
c[k] += c64[i] >> 25;
if (++k >= N)
k = 0;
}
}
void ntru_mult_karatsuba_64(int16_t *a, int16_t *b, int16_t *c, uint16_t len, uint16_t N, uint16_t mod_mask) {
if (len < NTRU_KARATSUBA_THRESH_64)
ntru_mult_int_64_base(a, b, c, len, N, mod_mask);
else {
uint16_t len2 = len / 2;
int16_t z0[NTRU_INT_POLY_SIZE];
int16_t z1[NTRU_INT_POLY_SIZE];
int16_t z2[NTRU_INT_POLY_SIZE];
/* z0, z2 */
ntru_mult_karatsuba_64(a, b, z0, len2, N, mod_mask);
ntru_mult_karatsuba_64(a+len2, b+len2, z2, len-len2, N, mod_mask);
/* z1 */
int16_t lh1[NTRU_INT_POLY_SIZE];
int16_t lh2[NTRU_INT_POLY_SIZE];
uint16_t i;
for (i=0; i<len2; i++) {
lh1[i] = a[i] + a[len2+i];
lh2[i] = b[i] + b[len2+i];
}
if (len%2 != 0) {
lh1[len-len2-1] = a[len-1];
lh2[len-len2-1] = b[len-1];
}
ntru_mult_karatsuba_64(lh1, lh2, z1, len-len2, N, mod_mask);
for (i=0; i<2*len2-1; i++)
z1[i] -= z0[i];
z1[len] = 0;
for (i=0; i<2*(len-len2)-1; i++)
z1[i] -= z2[i];
/* c */
memset(c, 0, NTRU_INT_POLY_SIZE*2);
memcpy(c, z0, 2*(2*len2-1)); /* 2*len2-1 coefficients */
uint16_t c_idx = len2;
for (i=0; i<2*(len-len2)-1; i++) {
c[c_idx] += z1[i];
c_idx++;
if (c_idx >= N)
c_idx = 0;
}
c_idx = 2 * len2;
if (c_idx >= N)
c_idx = 0;
for (i=0; i<2*(len-len2)-1; i++) {
c[c_idx] += z2[i];
c_idx++;
if (c_idx >= N)
c_idx = 0;
}
}
}
uint8_t ntru_mult_int_64(NtruIntPoly *a, NtruIntPoly *b, NtruIntPoly *c, uint16_t mod_mask) {
uint16_t N = a->N;
if (N != b->N)
return 0;
c->N = N;
ntru_mult_karatsuba_64((int16_t*)&a->coeffs, (int16_t*)&b->coeffs, (int16_t*)&c->coeffs, N, N, mod_mask);
ntru_mod_mask(c, mod_mask);
return 1;
}
uint8_t ntru_mult_tern_32(NtruIntPoly *a, NtruTernPoly *b, NtruIntPoly *c, uint16_t mod_mask) {
uint16_t N = a->N;
if (N != b->N)
return 0;
c->N = N;
memset(&c->coeffs, 0, N * sizeof c->coeffs[0]);
uint32_t mod_mask_32 = mod_mask;
mod_mask_32 += mod_mask_32 << 16;
typedef uint32_t __attribute__((__may_alias__)) uint32_t_alias;
/* make sure a.coeffs[i] <= mod_mask */
ntru_mod_mask(a, mod_mask_32);
uint16_t overflow_ctr_start = (1<<16)/(mod_mask+1) - 1;
uint16_t overflow_ctr_rem = overflow_ctr_start;
/* add coefficients that are multiplied by 1 */
uint16_t i;
for (i=0; i<b->num_ones; i++) {
int16_t j;
int16_t k = b->ones[i];
uint16_t j_end = N-2<b->ones[i] ? 0 : N-2-b->ones[i];
for (j=0; j<j_end; j+=2,k+=2)
*((uint32_t_alias*)&c->coeffs[k]) += *((uint32_t_alias*)&a->coeffs[j]);
for (; k<N; k++,j++)
c->coeffs[k] += a->coeffs[j];
for (k=0; j<N-2; j+=2,k+=2)
*((uint32_t_alias*)&c->coeffs[k]) += *((uint32_t_alias*)&a->coeffs[j]);
for (; j<N; j++,k++)
c->coeffs[k] += a->coeffs[j];
overflow_ctr_rem--;
if (!overflow_ctr_rem) {
ntru_mod_mask(c, mod_mask);
overflow_ctr_rem = overflow_ctr_start;
}
}
/* use inverse mask for subtraction */
mod_mask_32 = ~mod_mask_32;
for (i=0; i<N-2; i+=2)
*((uint32_t_alias*)&c->coeffs[i]) |= mod_mask_32;
for (; i<N; i++)
c->coeffs[i] |= mod_mask_32;
/* subtract coefficients that are multiplied by -1 */
overflow_ctr_rem = overflow_ctr_start;
for (i=0; i<b->num_neg_ones; i++) {
int16_t j;
int16_t k = b->neg_ones[i];
uint16_t j_end = N-2<b->neg_ones[i] ? 0 : N-2-b->neg_ones[i];
for (j=0; j<j_end; j+=2,k+=2)
*((uint32_t_alias*)&c->coeffs[k]) -= *((uint32_t_alias*)&a->coeffs[j]);
for (; k<N; k++,j++)
c->coeffs[k] -= a->coeffs[j];
for (k=0; j<N-2; j+=2,k+=2)
*((uint32_t_alias*)&c->coeffs[k]) -= *((uint32_t_alias*)&a->coeffs[j]);
for (; j<N; j++,k++)
c->coeffs[k] -= a->coeffs[j];
overflow_ctr_rem--;
if (!overflow_ctr_rem) {
for (j=0; j<N-2; j+=2)
*((uint32_t_alias*)&c->coeffs[j]) |= mod_mask_32;
for (; j<N; j++)
c->coeffs[j] |= mod_mask_32;
overflow_ctr_rem = overflow_ctr_start;
}
}
ntru_mod_mask(c, mod_mask);
return 1;
}
uint8_t ntru_mult_tern_64(NtruIntPoly *a, NtruTernPoly *b, NtruIntPoly *c, uint16_t mod_mask) {
uint16_t N = a->N;
if (N != b->N)
return 0;
memset(&c->coeffs, 0, N * sizeof c->coeffs[0]);
c->N = N;
uint64_t mod_mask_64 = mod_mask;
mod_mask_64 += mod_mask_64 << 16;
mod_mask_64 += mod_mask_64 << 32;
typedef uint64_t __attribute__((__may_alias__)) uint64_t_alias;
/* make sure a.coeffs[i] <= mod_mask */
ntru_mod_mask(a, mod_mask);
uint16_t overflow_ctr_start = (1<<16)/(mod_mask+1) - 1;
uint16_t overflow_ctr_rem = overflow_ctr_start;
/* add coefficients that are multiplied by 1 */
uint16_t i;
for (i=0; i<b->num_ones; i++) {
int16_t j;
int16_t k = b->ones[i];
uint16_t j_end = N-4<b->ones[i] ? 0 : N-4-b->ones[i];
for (j=0; j<j_end; j+=4,k+=4)
*((uint64_t_alias*)&c->coeffs[k]) += *((uint64_t_alias*)&a->coeffs[j]);
for (; k<N; k++,j++)
c->coeffs[k] += a->coeffs[j];
for (k=0; j<N-4; j+=4,k+=4)
*((uint64_t_alias*)&c->coeffs[k]) += *((uint64_t_alias*)&a->coeffs[j]);
for (; j<N; j++,k++)
c->coeffs[k] += a->coeffs[j];
overflow_ctr_rem--;
if (!overflow_ctr_rem) {
ntru_mod_mask(c, mod_mask);
overflow_ctr_rem = overflow_ctr_start;
}
}
/* use inverse mask for subtraction */
mod_mask_64 = ~mod_mask_64;
for (i=0; i<N-4; i+=4)
*((uint64_t_alias*)&c->coeffs[i]) |= mod_mask_64;
for (; i<N; i++)
c->coeffs[i] |= mod_mask_64;
/* subtract coefficients that are multiplied by -1 */
overflow_ctr_rem = overflow_ctr_start;
for (i=0; i<b->num_neg_ones; i++) {
int16_t j;
int16_t k = b->neg_ones[i];
uint16_t j_end = N-4<b->neg_ones[i] ? 0 : N-4-b->neg_ones[i];
for (j=0; j<j_end; j+=4,k+=4)
*((uint64_t_alias*)&c->coeffs[k]) -= *((uint64_t_alias*)&a->coeffs[j]);
for (; k<N; k++,j++)
c->coeffs[k] -= a->coeffs[j];
for (k=0; j<N-4; j+=4,k+=4)
*((uint64_t_alias*)&c->coeffs[k]) -= *((uint64_t_alias*)&a->coeffs[j]);
for (; j<N; j++,k++)
c->coeffs[k] -= a->coeffs[j];
overflow_ctr_rem--;
if (!overflow_ctr_rem) {
for (j=0; j<N-4; j+=4)
*((uint64_t_alias*)&c->coeffs[j]) |= mod_mask_64;
for (; j<N; j++)
c->coeffs[j] |= mod_mask_64;
overflow_ctr_rem = overflow_ctr_start;
}
}
ntru_mod_mask(c, mod_mask);
return 1;
}
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
uint8_t ntru_mult_prod(NtruIntPoly *a, NtruProdPoly *b, NtruIntPoly *c, uint16_t mod_mask) {
uint16_t N = a->N;
if (N != b->N)
return 0;
c->N = N;
memset(&c->coeffs, 0, N * sizeof c->coeffs[0]);
NtruIntPoly temp;
ntru_mult_tern(a, &b->f1, &temp, mod_mask);
ntru_mult_tern(&temp, &b->f2, c, mod_mask);
NtruIntPoly f3a;
ntru_mult_tern(a, &b->f3, &f3a, mod_mask);
ntru_add(c, &f3a);
ntru_mod_mask(c, mod_mask);
return 1;
}
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
uint8_t ntru_mult_priv(NtruPrivPoly *a, NtruIntPoly *b, NtruIntPoly *c, uint16_t mod_mask) {
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
if (a->prod_flag)
return ntru_mult_prod(b, &a->poly.prod, c, mod_mask);
else
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
return ntru_mult_tern(b, &a->poly.tern, c, mod_mask);
}
/** NtruPrivPoly to binary (coefficients reduced mod 2), 64 bit version */
void ntru_priv_to_mod2_64(NtruPrivPoly *a, uint64_t *b_coeffs64) {
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
if (a->prod_flag) {
NtruProdPoly *prod = &a->poly.prod;
uint16_t N = prod->N;
uint16_t N64 = (prod->N+63) / 64;
memset(b_coeffs64, 0, N64*8);
uint16_t i, j, bidx;
for (i=0; i<prod->f1.num_ones; i++) {
for (j=0; j<prod->f2.num_ones; j++) {
bidx = prod->f1.ones[i] + prod->f2.ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs64[bidx/64] ^= ((uint64_t)1) << (bidx%64);
}
for (j=0; j<prod->f2.num_neg_ones; j++) {
bidx = prod->f1.ones[i] + prod->f2.neg_ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs64[bidx/64] ^= ((uint64_t)1) << (bidx%64);
}
}
for (i=0; i<prod->f1.num_neg_ones; i++) {
for (j=0; j<prod->f2.num_ones; j++) {
bidx = prod->f1.neg_ones[i] + prod->f2.ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs64[bidx/64] ^= ((uint64_t)1) << (bidx%64);
}
for (j=0; j<prod->f2.num_neg_ones; j++) {
bidx = prod->f1.neg_ones[i] + prod->f2.neg_ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs64[bidx/64] ^= ((uint64_t)1) << (bidx%64);
}
}
for (i=0; i<prod->f3.num_ones; i++) {
uint16_t ai = prod->f3.ones[i];
b_coeffs64[ai/64] ^= ((uint64_t)1) << (ai%64);
}
for (i=0; i<prod->f3.num_neg_ones; i++) {
uint16_t ai = prod->f3.neg_ones[i];
b_coeffs64[ai/64] ^= ((uint64_t)1) << (ai%64);
}
}
else
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
{
NtruTernPoly *tern = &a->poly.tern;
uint16_t N64 = (tern->N+63) / 64;
memset(b_coeffs64, 0, N64*8);
uint16_t i;
for (i=0; i<tern->num_ones; i++) {
uint16_t ai = tern->ones[i];
b_coeffs64[ai/64] ^= ((uint64_t)1) << (ai%64);
}
for (i=0; i<tern->num_neg_ones; i++) {
uint16_t ai = tern->neg_ones[i];
b_coeffs64[ai/64] ^= ((uint64_t)1) << (ai%64);
}
}
}
/** NtruPrivPoly to binary (coefficients reduced mod 2), 32 bit version */
void ntru_priv_to_mod2_32(NtruPrivPoly *a, uint32_t *b_coeffs32) {
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
if (a->prod_flag) {
NtruProdPoly *prod = &a->poly.prod;
uint16_t N = prod->N;
uint16_t N32 = (prod->N+31) / 32;
memset(b_coeffs32, 0, N32*4);
uint16_t i, j, bidx;
for (i=0; i<prod->f1.num_ones; i++) {
for (j=0; j<prod->f2.num_ones; j++) {
bidx = prod->f1.ones[i] + prod->f2.ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs32[bidx/32] ^= ((uint32_t)1) << (bidx%32);
}
for (j=0; j<prod->f2.num_neg_ones; j++) {
bidx = prod->f1.ones[i] + prod->f2.neg_ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs32[bidx/32] ^= ((uint32_t)1) << (bidx%32);
}
}
for (i=0; i<prod->f1.num_neg_ones; i++) {
for (j=0; j<prod->f2.num_ones; j++) {
bidx = prod->f1.neg_ones[i] + prod->f2.ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs32[bidx/32] ^= ((uint32_t)1) << (bidx%32);
}
for (j=0; j<prod->f2.num_neg_ones; j++) {
bidx = prod->f1.neg_ones[i] + prod->f2.neg_ones[j];
if (bidx >= N)
bidx -= N;
b_coeffs32[bidx/32] ^= ((uint32_t)1) << (bidx%32);
}
}
for (i=0; i<prod->f3.num_ones; i++) {
uint16_t ai = prod->f3.ones[i];
b_coeffs32[ai/32] ^= ((uint32_t)1) << (ai%32);
}
for (i=0; i<prod->f3.num_neg_ones; i++) {
uint16_t ai = prod->f3.neg_ones[i];
b_coeffs32[ai/32] ^= ((uint32_t)1) << (ai%32);
}
}
else
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
{
NtruTernPoly *tern = &a->poly.tern;
uint16_t N32 = (tern->N+31) / 32;
memset(b_coeffs32, 0, N32*4);
uint16_t i;
for (i=0; i<tern->num_ones; i++) {
uint16_t ai = tern->ones[i];
b_coeffs32[ai/32] ^= ((uint32_t)1) << (ai%32);
}
for (i=0; i<tern->num_neg_ones; i++) {
uint16_t ai = tern->neg_ones[i];
b_coeffs32[ai/32] ^= ((uint32_t)1) << (ai%32);
}
}
}
void ntru_to_arr_64(NtruIntPoly *p, uint16_t q, uint8_t *a) {
uint16_t N = p->N;
uint8_t log_q = ntru_log2(q);
uint16_t enc_bytes = ntru_enc_len_Nq(N, q);
uint16_t rem = enc_bytes % sizeof(uint64_t);
uint16_t quo = enc_bytes / sizeof(uint64_t);
uint16_t enc_last_int = rem ? quo : quo - 1;
uint16_t enc_last_int_valid = rem ? rem : sizeof(uint64_t);
uint64_t last = 0;
typedef uint64_t __attribute__((__may_alias__)) *uint64_t_alias;
uint64_t *a64 = (uint64_t_alias)a;
uint16_t a_idx = 0; /* index into a64 */
uint8_t bit_idx = 0; /* next unused bit of a64[a_idx] */
a64[0] = 0;
uint16_t p_idx;
uint64_t mod_mask = q - 1;
for (p_idx=0; p_idx<N; p_idx++) {
uint64_t coeff = p->coeffs[p_idx] & mod_mask;
if (bit_idx < 64-log_q) {
if (a_idx == enc_last_int)
last |= coeff << bit_idx;
else
a64[a_idx] |= coeff << bit_idx;
bit_idx += log_q;
}
else {
a64[a_idx] |= coeff << bit_idx;
a_idx++;
bit_idx += log_q - 64;
if (a_idx == enc_last_int)
last = coeff >> (log_q - bit_idx);
else
a64[a_idx] = coeff >> (log_q-bit_idx);
}
}
/* reverse byte order on big-endian machines */
uint16_t i;
for (i = 0; i <= a_idx; i++)
{
if (i == enc_last_int) {
last = htole64(last);
memcpy(&a64[i], &last, enc_last_int_valid);
}
else
a64[i] = htole64(a64[i]);
}
}
void ntru_to_arr_32(NtruIntPoly *p, uint16_t q, uint8_t *a) {
uint16_t N = p->N;
uint8_t log_q = ntru_log2(q);
uint16_t enc_bytes = ntru_enc_len_Nq(N, q);
uint16_t rem = enc_bytes % sizeof(uint32_t);
uint16_t quo = enc_bytes / sizeof(uint32_t);
uint16_t enc_last_int = rem ? quo : quo - 1;
uint16_t enc_last_int_valid = rem ? rem : sizeof(uint32_t);
uint32_t last = 0;
typedef uint32_t __attribute__((__may_alias__)) *uint32_t_alias;
uint32_t *a32 = (uint32_t_alias)a;
uint16_t a_idx = 0; /* index into a32 */
uint8_t bit_idx = 0; /* next unused bit of a32[a_idx] */
a32[0] = 0;
uint16_t p_idx;
uint32_t mod_mask = q - 1;
for (p_idx=0; p_idx<N; p_idx++) {
uint32_t coeff = p->coeffs[p_idx] & mod_mask;
if (bit_idx < 32-log_q) {
if (a_idx == enc_last_int)
last |= coeff << bit_idx;
else
a32[a_idx] |= coeff << bit_idx;
bit_idx += log_q;
}
else {
a32[a_idx] |= coeff << bit_idx;
a_idx++;
bit_idx += log_q - 32;
if (a_idx == enc_last_int)
last = coeff >> (log_q - bit_idx);
else
a32[a_idx] = coeff >> (log_q-bit_idx);
}
}
/* reverse byte order on big-endian machines */
uint16_t i;
for (i = 0; i <= a_idx; i++)
{
if (i == enc_last_int) {
last = htole32(last);
memcpy(&a32[i], &last, enc_last_int_valid);
}
else
a32[i] = htole32(a32[i]);
}
}
void ntru_to_arr4(NtruIntPoly *p, uint8_t *arr) {
uint16_t i = 0;
while (i < p->N-3) {
int8_t c0 = p->coeffs[i] & 3;
int8_t c1 = p->coeffs[i+1] & 3;
int8_t c2 = p->coeffs[i+2] & 3;
int8_t c3 = p->coeffs[i+3] & 3;
int16_t d = c0 + (c1<<2) + (c2<<4) + (c3<<6);
arr[i/4] = d;
i += 4;
}
/* handle the last 0 to 3 coefficients */
if (i >= p->N)
return;
uint16_t last = i / 4;
arr[last] = p->coeffs[i] & 3;
i++;
if (i >= p->N)
return;
arr[last] |= (p->coeffs[i]&3) << 2;
i++;
if (i >= p->N)
return;
arr[last] |= (p->coeffs[i]&3) << 4;
i++;
if (i >= p->N)
return;
arr[last] |= (p->coeffs[i]&3) << 6;
}
void ntru_from_arr(uint8_t *arr, uint16_t N, uint16_t q, NtruIntPoly *p) {
p->N = N;
memset(&p->coeffs, 0, N * sizeof p->coeffs[0]);
uint8_t bits_per_coeff = ntru_log2(q);
uint32_t mask = 0xFFFFFFFF >> (32-bits_per_coeff); /* for truncating values to bitsPerCoeff bits */
uint16_t byte_idx = 0;
uint8_t bit_idx = 0; /* next bit in arr[byte_idx] */
uint32_t coeff_buf = 0; /* contains (bit_idx) bits */
uint8_t coeff_bits = 0; /* length of coeffBuf */
uint16_t coeff_idx = 0; /* index into coeffs */
while (coeff_idx < N) {
/* copy bits_per_coeff or more into coeff_buf */
while (coeff_bits < bits_per_coeff) {
coeff_buf += (arr[byte_idx]&0xFF) << coeff_bits;
coeff_bits += 8 - bit_idx;
byte_idx++;
bit_idx = 0;
}
/* low bits_per_coeff bits = next coefficient */
p->coeffs[coeff_idx] = coeff_buf & mask;
coeff_idx++;
coeff_buf >>= bits_per_coeff;
coeff_bits -= bits_per_coeff;
}
}
void ntru_mult_fac(NtruIntPoly *a, int16_t factor) {
uint16_t i;
for (i=0; i<a->N; i++)
a->coeffs[i] *= factor;
}
void ntru_mod_64(NtruIntPoly *p, uint16_t mod_mask) {
typedef uint64_t __attribute__((__may_alias__)) uint64_t_alias;
uint64_t mod_mask_64 = mod_mask;
mod_mask_64 += mod_mask_64 << 16;
mod_mask_64 += mod_mask_64 << 32;
uint16_t i;
for (i=0; i<p->N; i+=4)
*((uint64_t_alias*)&p->coeffs[i]) &= mod_mask_64;
}
void ntru_mod_32(NtruIntPoly *p, uint16_t mod_mask) {
typedef uint32_t __attribute__((__may_alias__)) uint32_t_alias;
uint32_t mod_mask_32 = mod_mask;
mod_mask_32 += mod_mask_32 << 16;
uint16_t i;
for (i=0; i<p->N; i+=2)
*((uint32_t_alias*)&p->coeffs[i]) &= mod_mask_32;
}
void ntru_mod3_standard(NtruIntPoly *p) {
uint16_t i;
for (i=0; i<p->N; i++) {
int8_t c = p->coeffs[i] % 3;
if (c == -2)
c = 1;
if (c == -1)
c = 2;
p->coeffs[i] = c;
}
}
void ntru_mod3(NtruIntPoly *p) {
#ifdef NTRU_DETECT_SIMD
if (__builtin_cpu_supports("avx2"))
return ntru_mod3_avx2(p);
else if (__builtin_cpu_supports("ssse3"))
return ntru_mod3_sse(p);
else
return ntru_mod3_standard(p);
#else
#ifdef __AVX2__
ntru_mod3_avx2(p);
#elif __SSSE3__
ntru_mod3_sse(p);
#else
ntru_mod3_standard(p);
#endif /* __SSSE3__ */
#endif /* NTRU_DETECT_SIMD */
}
void ntru_mod_center(NtruIntPoly *p, uint16_t modulus) {
uint16_t m2 = modulus / 2;
uint16_t mod_mask = modulus - 1;
uint16_t i;
for (i=0; i<p->N; i++) {
uint16_t c = p->coeffs[i] & mod_mask; // note that c is unsigned
if (c > m2)
c -= modulus;
p->coeffs[i] = c;
}
}
uint8_t ntru_equals_int(NtruIntPoly *a, NtruIntPoly *b) {
if (a->N != b->N)
return 0;
uint16_t i;
for (i=0; i<a->N; i++)
if (a->coeffs[i] != b->coeffs[i])
return 0;
return 1;
}
uint16_t ntru_deg_64(uint64_t *coeffs, uint16_t len) {
uint16_t deg = 64*len - 1;
len--;
while (len>0 && coeffs[len]==0) {
len--;
deg -= 64;
}
while (coeffs[len]>>(deg%64)==0 && deg>0)
deg--;
return deg;
}
uint16_t ntru_deg_32(uint32_t *coeffs, uint16_t len) {
uint16_t deg = 32*len - 1;
len--;
while (len>0 && coeffs[len]==0) {
len--;
deg -= 32;
}
while (coeffs[len]>>(deg%32)==0 && deg>0)
deg--;
return deg;
}
void ntru_clear_tern(NtruTernPoly *p) {
memset(&p->ones, 0, p->num_ones * sizeof p->ones[0]);
memset(&p->neg_ones, 0, p->num_neg_ones * sizeof p->neg_ones[0]);
}
void ntru_clear_priv(NtruPrivPoly *p) {
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
if (p->prod_flag) {
ntru_clear_tern(&p->poly.prod.f1);
ntru_clear_tern(&p->poly.prod.f2);
ntru_clear_tern(&p->poly.prod.f3);
}
else
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
ntru_clear_tern(&p->poly.tern);
}
void ntru_clear_int(NtruIntPoly *p) {
uint16_t i;
for (i=0; i<p->N; i++)
p->coeffs[i] = 0;
}
/**
* @brief Lift inverse
*
* Given a polynomial a and the inverse of (1+3a) mod 2, this function
* calculates the inverse of (1+3a) mod q.
*
* @param a a polynomial such that Fq = (1+3a)^(-1) (mod 2)
* @param Fq the inverse of 1+3a modulo 2
* @param q the modulus
*/
void ntru_lift_inverse(NtruPrivPoly *a, NtruIntPoly *Fq, uint16_t q) {
NtruIntPoly temp1, temp2;
uint32_t v = 2;
while (v < q) {
v *= v;
/* temp1 = (1+3a)*Fq */
ntru_mult_priv(a, Fq, &temp1, q-1);
ntru_mult_fac(&temp1, 3);
ntru_add(&temp1, Fq);
ntru_neg_mod(&temp1, q);
temp1.coeffs[0] += 2;
memcpy(&temp2, Fq, sizeof *Fq);
ntru_mult_int(&temp1, &temp2, Fq, q-1);
}
}
uint8_t ntru_invert_32(NtruPrivPoly *a, uint16_t mod_mask, NtruIntPoly *Fq) {
int16_t i;
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
uint16_t N = a->prod_flag ? a->poly.prod.N : a->poly.tern.N;
#else
uint16_t N = a->poly.tern.N;
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
uint16_t k = 0;
uint16_t N32 = (N+1+31) / 32; /* #uint32_t's needed for N+1 coeffs */
/* b = 1 */
uint32_t b_coeffs32_arr[N32];
uint32_t *b_coeffs32 = b_coeffs32_arr;
memset(b_coeffs32+1, 0, (N32-1)*4);
b_coeffs32[0] = 1;
/* c = 0 */
uint32_t c_coeffs32_arr[N32];
uint32_t *c_coeffs32 = c_coeffs32_arr;
memset(c_coeffs32, 0, N32*4);
/* f=3a+1; skip multiplication by 3 because f=3f (mod 2) */
uint32_t f_coeffs32_arr[N32];
uint32_t *f_coeffs32 = f_coeffs32_arr;
ntru_priv_to_mod2_32(a, f_coeffs32);
f_coeffs32[0] ^= 1;
/* g(x) = x^N − 1 */
uint32_t g_coeffs32_arr[N32];
uint32_t *g_coeffs32 = g_coeffs32_arr;
memset(g_coeffs32, 0, N32*4);
g_coeffs32[0] = 1;
g_coeffs32[N/32] |= ((uint32_t)1) << (N%32);
uint16_t deg_f = ntru_deg_32(f_coeffs32, N32);
uint16_t deg_g = N;
for (;;) {
uint16_t num_zeros = 0;
/* while f[0]==0 */
while ((f_coeffs32[num_zeros/32]&(((uint32_t)1)<<(num_zeros%32)))==0 && num_zeros<=N)
num_zeros++;
if (num_zeros >= N) /* not invertible */
return 0;
k += num_zeros;
/* right-shift f, left-shift c num_zeros coefficients each */
if (num_zeros >= 32) {
memmove(c_coeffs32+num_zeros/32, c_coeffs32, N32*4-num_zeros/32*4);
memset(c_coeffs32, 0, num_zeros/32*4);
memmove(f_coeffs32, f_coeffs32+num_zeros/32, N32*4-num_zeros/32*4);
memset(f_coeffs32+N32-num_zeros/32, 0, num_zeros/32*4);
deg_f -= num_zeros / 32 * 32;
num_zeros %= 32;
}
if (num_zeros > 0) {
/* c(x) = c(x)*(x^num_zeros) */
for (i=N32-1; i>0; i--) {
c_coeffs32[i] <<= num_zeros;
c_coeffs32[i] |= c_coeffs32[i-1] >> (32-num_zeros);
}
c_coeffs32[0] <<= num_zeros;
/* f(x) = f(x)/(x^num_zeros) */
for (i=1; i<N32; i++) {
f_coeffs32[i-1] >>= num_zeros;
f_coeffs32[i-1] |= f_coeffs32[i] << (32-num_zeros);
}
f_coeffs32[i-1] >>= num_zeros;
}
deg_f -= num_zeros;
if (deg_f==0 && f_coeffs32[0]==1) /* if f==1 */
break;
if (deg_f < deg_g) {
/* exchange f and g */
uint32_t *temp_coeffs = f_coeffs32;
f_coeffs32 = g_coeffs32;
g_coeffs32 = temp_coeffs;
uint16_t temp = deg_f;
deg_f = deg_g;
deg_g = temp;
/* exchange b and c */
temp_coeffs = b_coeffs32;
b_coeffs32 = c_coeffs32;
c_coeffs32 = temp_coeffs;
}
ntru_add_mod2_32(f_coeffs32, g_coeffs32, N32);
/* adding f+g may have lowered the degree of f */
while (deg_f>0 && (f_coeffs32[deg_f/32]&(((uint32_t)1)<<(deg_f%32)))==0)
deg_f--;
ntru_add_mod2_32(b_coeffs32, c_coeffs32, N32);
}
if ((b_coeffs32[(N+1-1)/32]&(((uint32_t)1)<<((N+1-1)%32))) != 0) /* if (b[N]!=0) */
return 0;
/* Fq(x) = x^(N-k) * b(x) */
memset(&Fq->coeffs, 0, N * sizeof Fq->coeffs[0]);
Fq->N = N;
int16_t j = 0;
while (k >= N)
k -= N;
for (i=N-1; i>=0; i--) {
j = i - k;
if (j < 0)
j += N;
Fq->coeffs[j] = (b_coeffs32[i/32]>>(i%32)) & 1; /* Fq->coeffs[j]=b[i] */
}
ntru_lift_inverse(a, Fq, mod_mask+1);
return 1;
}
uint8_t ntru_invert_64(NtruPrivPoly *a, uint16_t mod_mask, NtruIntPoly *Fq) {
#ifndef NTRU_AVOID_HAMMING_WT_PATENT
uint16_t N = a->prod_flag ? a->poly.prod.N : a->poly.tern.N;
#else
uint16_t N = a->poly.tern.N;
#endif /* NTRU_AVOID_HAMMING_WT_PATENT */
uint16_t k = 0;
uint16_t N64 = (N+1+63) / 64; /* #uint64_t's needed for N+1 coeffs */
/* b = 1 */
uint64_t b_coeffs64_arr[N64];
uint64_t *b_coeffs64 = b_coeffs64_arr;
memset(b_coeffs64+1, 0, (N64-1)*8);
b_coeffs64[0] = 1;
/* c = 0 */
uint64_t c_coeffs64_arr[N64];
uint64_t *c_coeffs64 = c_coeffs64_arr;
memset(c_coeffs64, 0, N64*8);
/* f=3a+1; skip multiplication by 3 because f=3f (mod 2) */
uint64_t f_coeffs64_arr[N64];
uint64_t *f_coeffs64 = f_coeffs64_arr;
ntru_priv_to_mod2_64(a, f_coeffs64);
f_coeffs64[0] ^= 1;
/* g(x) = x^N − 1 */
uint64_t g_coeffs64_arr[N64];
uint64_t *g_coeffs64 = g_coeffs64_arr;
memset(g_coeffs64, 0, N64*8);
g_coeffs64[0] = 1;
g_coeffs64[N/64] |= ((uint64_t)1) << (N%64);
uint16_t deg_f = ntru_deg_64(f_coeffs64, N64);
uint16_t deg_g = N;
for (;;) {
uint16_t num_zeros = 0;
/* while f[0]==0 */
while ((f_coeffs64[num_zeros/64]&(((uint64_t)1)<<(num_zeros%64)))==0 && num_zeros<=N)
num_zeros++;
if (num_zeros >= N) /* not invertible */
return 0;
k += num_zeros;
/* right-shift f, left-shift c num_zeros coefficients each */
if (num_zeros >= 64) {
memmove(c_coeffs64+num_zeros/64, c_coeffs64, N64*8-num_zeros/64*8);
memset(c_coeffs64, 0, num_zeros/64*8);
memmove(f_coeffs64, f_coeffs64+num_zeros/64, N64*8-num_zeros/64*8);
memset(f_coeffs64+N64-num_zeros/64, 0, num_zeros/64*8);
deg_f -= num_zeros / 64 * 64;
num_zeros %= 64;
}
if (num_zeros > 0) {
int16_t i;
/* c(x) = c(x)*(x^num_zeros) */
for (i=N64-1; i>0; i--) {
c_coeffs64[i] <<= num_zeros;
c_coeffs64[i] |= c_coeffs64[i-1] >> (64-num_zeros);
}
c_coeffs64[0] <<= num_zeros;
/* f(x) = f(x)/(x^num_zeros) */
for (i=1; i<N64; i++) {
f_coeffs64[i-1] >>= num_zeros;
f_coeffs64[i-1] |= f_coeffs64[i] << (64-num_zeros);
}
f_coeffs64[i-1] >>= num_zeros;
}
deg_f -= num_zeros;
if (deg_f==0 && f_coeffs64[0]==1) /* if f==1 */
break;
if (deg_f < deg_g) {
/* exchange f and g */
uint64_t *temp_coeffs = f_coeffs64;
f_coeffs64 = g_coeffs64;
g_coeffs64 = temp_coeffs;
uint16_t temp = deg_f;
deg_f = deg_g;
deg_g = temp;
/* exchange b and c */
temp_coeffs = b_coeffs64;
b_coeffs64 = c_coeffs64;
c_coeffs64 = temp_coeffs;
}
ntru_add_mod2_64(f_coeffs64, g_coeffs64, N64);
/* adding f+g may have lowered the degree of f */
while (deg_f>0 && (f_coeffs64[deg_f/64]&(((uint64_t)1)<<(deg_f%64)))==0)
deg_f--;
ntru_add_mod2_64(b_coeffs64, c_coeffs64, N64);
}
if ((b_coeffs64[(N+1-1)/64]&(((uint64_t)1)<<((N+1-1)%64))) != 0) /* if (b[N]!=0) */
return 0;
/* Fq(x) = x^(N-k) * b(x) */
memset(&Fq->coeffs, 0, N * sizeof Fq->coeffs[0]);
Fq->N = N;
int16_t j = 0;
while (k >= N)
k -= N;
int16_t i;
for (i=N-1; i>=0; i--) {
j = i - k;
if (j < 0)
j += N;
Fq->coeffs[j] = (b_coeffs64[i/64]>>(i%64)) & 1; /* Fq->coeffs[j]=b[i] */
}
ntru_lift_inverse(a, Fq, mod_mask+1);
return 1;
}
uint8_t (*ntru_invert)(NtruPrivPoly *a, uint16_t mod_mask, NtruIntPoly *Fq);
void ntru_set_optimized_impl_poly() {
#ifdef NTRU_DETECT_SIMD
if (__builtin_cpu_supports("avx2")) {
ntru_mult_int = ntru_mult_int_avx2;
ntru_mult_tern = ntru_mult_tern_avx2;
ntru_to_arr = ntru_to_arr_sse;
ntru_mod_mask = ntru_mod_avx2;
}
else if (__builtin_cpu_supports("ssse3")) {
ntru_mult_int = ntru_mult_int_sse;
ntru_mult_tern = ntru_mult_tern_sse;
ntru_to_arr = ntru_to_arr_sse;
ntru_mod_mask = ntru_mod_sse;
}
else if (sizeof(void*) >= 8) { /* 64-bit arch */
ntru_mult_int = ntru_mult_int_64;
ntru_mult_tern = ntru_mult_tern_64;
ntru_to_arr = ntru_to_arr_64;
ntru_mod_mask = ntru_mod_64;
}
else {
ntru_mult_int = ntru_mult_int_16;
ntru_mult_tern = ntru_mult_tern_32;
ntru_to_arr = ntru_to_arr_32;
ntru_mod_mask = ntru_mod_32;
}
if (sizeof(void*) >= 8) /* 64-bit arch */
ntru_invert = ntru_invert_64;
else
ntru_invert = ntru_invert_32;
#else
#ifdef __AVX2__
ntru_mult_int = ntru_mult_int_avx2;
ntru_mult_tern = ntru_mult_tern_avx2;
ntru_to_arr = ntru_to_arr_sse;
ntru_mod_mask = ntru_mod_avx2;
#elif __SSSE3__
ntru_mult_int = ntru_mult_int_sse;
ntru_mult_tern = ntru_mult_tern_sse;
ntru_to_arr = ntru_to_arr_sse;
ntru_mod_mask = ntru_mod_sse;
#elif _LP64
ntru_mult_int = ntru_mult_int_64;
ntru_mult_tern = ntru_mult_tern_64;
ntru_to_arr = ntru_to_arr_64;
ntru_mod_mask = ntru_mod_64;
#else
ntru_mult_int = ntru_mult_int_16;
ntru_mult_tern = ntru_mult_tern_32;
ntru_to_arr = ntru_to_arr_32;
ntru_mod_mask = ntru_mod_32;
#endif
#if _LP64
ntru_invert = ntru_invert_64;
#else
ntru_invert = ntru_invert_32;
#endif
#endif /* NTRU_DETECT_SIMD */
}