Infrastructure/ftu
1
0
Fork 0
Code Issues Pull requests Actions Packages Projects Releases Wiki Activity
ftu/blst/no_asm.h

1288 lines
36 KiB
C
Raw Normal View History

initial stuff
2022-09-09 06:47:49 +00:00
/*
* Copyright Supranational LLC
* Licensed under the Apache License, Version 2.0, see LICENSE for details.
* SPDX-License-Identifier: Apache-2.0
*/
#if LIMB_T_BITS==32
typedef unsigned long long llimb_t;
#endif
#if defined(__clang__)
# pragma GCC diagnostic ignored "-Wstatic-in-inline"
#endif
static void mul_mont_n(limb_t ret[], const limb_t a[], const limb_t b[],
const limb_t p[], limb_t n0, size_t n)
{
llimb_t limbx;
limb_t mask, borrow, mx, hi, tmp[n+1], carry;
size_t i, j;
for (mx=b[0], hi=0, i=0; i<n; i++) {
limbx = (mx * (llimb_t)a[i]) + hi;
tmp[i] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
mx = n0*tmp[0];
tmp[i] = hi;
for (carry=0, j=0; ; ) {
limbx = (mx * (llimb_t)p[0]) + tmp[0];
hi = (limb_t)(limbx >> LIMB_T_BITS);
for (i=1; i<n; i++) {
limbx = (mx * (llimb_t)p[i] + hi) + tmp[i];
tmp[i-1] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
limbx = tmp[i] + (hi + (llimb_t)carry);
tmp[i-1] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
if (++j==n)
break;
for (mx=b[j], hi=0, i=0; i<n; i++) {
limbx = (mx * (llimb_t)a[i] + hi) + tmp[i];
tmp[i] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
mx = n0*tmp[0];
limbx = hi + (llimb_t)carry;
tmp[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = carry - borrow;
for(i=0; i<n; i++)
ret[i] = (ret[i] & ~mask) | (tmp[i] & mask);
}
#define MUL_MONT_IMPL(bits) \
inline void mul_mont_##bits(vec##bits ret, const vec##bits a, \
const vec##bits b, const vec##bits p, limb_t n0) \
{ mul_mont_n(ret, a, b, p, n0, NLIMBS(bits)); } \
\
inline void sqr_mont_##bits(vec##bits ret, const vec##bits a, \
const vec##bits p, limb_t n0) \
{ mul_mont_n(ret, a, a, p, n0, NLIMBS(bits)); }
/*
* 256-bit subroutines can handle arbitrary modulus, even non-"sparse",
* but we have to harmonize the naming with assembly.
*/
#define mul_mont_256 mul_mont_sparse_256
#define sqr_mont_256 sqr_mont_sparse_256
MUL_MONT_IMPL(256)
#undef mul_mont_256
#undef sqr_mont_256
MUL_MONT_IMPL(384)
static void add_mod_n(limb_t ret[], const limb_t a[], const limb_t b[],
const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t mask, carry, borrow, tmp[n];
size_t i;
for (carry=0, i=0; i<n; i++) {
limbx = a[i] + (b[i] + (llimb_t)carry);
tmp[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = carry - borrow;
for(i=0; i<n; i++)
ret[i] = (ret[i] & ~mask) | (tmp[i] & mask);
}
#define ADD_MOD_IMPL(bits) \
inline void add_mod_##bits(vec##bits ret, const vec##bits a, \
const vec##bits b, const vec##bits p) \
{ add_mod_n(ret, a, b, p, NLIMBS(bits)); }
ADD_MOD_IMPL(256)
ADD_MOD_IMPL(384)
static void sub_mod_n(limb_t ret[], const limb_t a[], const limb_t b[],
const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t mask, carry, borrow;
size_t i;
for (borrow=0, i=0; i<n; i++) {
limbx = a[i] - (b[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = 0 - borrow;
for (carry=0, i=0; i<n; i++) {
limbx = ret[i] + ((p[i] & mask) + (llimb_t)carry);
ret[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
}
#define SUB_MOD_IMPL(bits) \
inline void sub_mod_##bits(vec##bits ret, const vec##bits a, \
const vec##bits b, const vec##bits p) \
{ sub_mod_n(ret, a, b, p, NLIMBS(bits)); }
SUB_MOD_IMPL(256)
SUB_MOD_IMPL(384)
static void mul_by_3_mod_n(limb_t ret[], const limb_t a[], const limb_t p[],
size_t n)
{
llimb_t limbx;
limb_t mask, carry, borrow, tmp[n], two_a[n];
size_t i;
for (carry=0, i=0; i<n; i++) {
limb_t a_i = a[i];
tmp[i] = a_i<<1 | carry;
carry = a_i>>(LIMB_T_BITS-1);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
two_a[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = carry - borrow;
for(i=0; i<n; i++)
two_a[i] = (two_a[i] & ~mask) | (tmp[i] & mask);
for (carry=0, i=0; i<n; i++) {
limbx = a[i] + (two_a[i] + (llimb_t)carry);
tmp[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = carry - borrow;
for(i=0; i<n; i++)
ret[i] = (ret[i] & ~mask) | (tmp[i] & mask);
}
#define MUL_BY_3_MOD_IMPL(bits) \
inline void mul_by_3_mod_##bits(vec##bits ret, const vec##bits a, \
const vec##bits p) \
{ mul_by_3_mod_n(ret, a, p, NLIMBS(bits)); }
MUL_BY_3_MOD_IMPL(256)
MUL_BY_3_MOD_IMPL(384)
static void lshift_mod_n(limb_t ret[], const limb_t a[], size_t count,
const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t mask, carry, borrow, tmp[n];
size_t i;
while (count--) {
for (carry=0, i=0; i<n; i++) {
limb_t a_i = a[i];
tmp[i] = a_i<<1 | carry;
carry = a_i>>(LIMB_T_BITS-1);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = carry - borrow;
for(i=0; i<n; i++)
ret[i] = (ret[i] & ~mask) | (tmp[i] & mask);
a = ret;
}
}
#define LSHIFT_MOD_IMPL(bits) \
inline void lshift_mod_##bits(vec##bits ret, const vec##bits a, size_t count, \
const vec##bits p) \
{ lshift_mod_n(ret, a, count, p, NLIMBS(bits)); }
LSHIFT_MOD_IMPL(256)
LSHIFT_MOD_IMPL(384)
static void cneg_mod_n(limb_t ret[], const limb_t a[], bool_t flag,
const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t borrow, mask, tmp[n];
size_t i;
for (borrow=0, i=0; i<n; i++) {
limbx = p[i] - (a[i] + (llimb_t)borrow);
tmp[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
flag &= vec_is_zero(a, sizeof(tmp)) ^ 1;
mask = (limb_t)0 - flag;
for(i=0; i<n; i++)
ret[i] = (a[i] & ~mask) | (tmp[i] & mask);
}
#define CNEG_MOD_IMPL(bits) \
inline void cneg_mod_##bits(vec##bits ret, const vec##bits a, bool_t flag, \
const vec##bits p) \
{ cneg_mod_n(ret, a, flag, p, NLIMBS(bits)); }
CNEG_MOD_IMPL(256)
CNEG_MOD_IMPL(384)
static limb_t check_mod_n(const byte a[], const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t borrow, ai, acc;
size_t i, j;
for (acc=borrow=0, i=0; i<n; i++) {
for (ai=0, j=0; j<8*sizeof(limb_t); j+=8)
ai |= (limb_t)(*a++) << j;
acc |= ai;
limbx = ai - (p[i] + (llimb_t)borrow);
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
return borrow & (is_zero(acc) ^ 1);
}
#define CHECK_MOD_IMPL(bits) \
inline limb_t check_mod_##bits(const pow##bits a, const vec##bits p) \
{ return check_mod_n(a, p, NLIMBS(bits)); }
CHECK_MOD_IMPL(256)
static limb_t add_n_check_mod_n(byte ret[], const byte a[], const byte b[],
const limb_t p[], size_t n)
{
limb_t ret_[n], a_[n], b_[n], zero;
limbs_from_le_bytes(a_, a, sizeof(a_));
limbs_from_le_bytes(b_, b, sizeof(b_));
add_mod_n(ret_, a_, b_, p, n);
zero = vec_is_zero(ret_, sizeof(ret_));
le_bytes_from_limbs(ret, ret_, sizeof(ret_));
return zero^1;
}
#define ADD_N_CHECK_MOD_IMPL(bits) \
inline limb_t add_n_check_mod_##bits(pow##bits ret, const pow##bits a, \
const pow##bits b, const vec##bits p) \
{ return add_n_check_mod_n(ret, a, b, p, NLIMBS(bits)); }
ADD_N_CHECK_MOD_IMPL(256)
static limb_t sub_n_check_mod_n(byte ret[], const byte a[], const byte b[],
const limb_t p[], size_t n)
{
limb_t ret_[n], a_[n], b_[n], zero;
limbs_from_le_bytes(a_, a, sizeof(a_));
limbs_from_le_bytes(b_, b, sizeof(b_));
sub_mod_n(ret_, a_, b_, p, n);
zero = vec_is_zero(ret_, sizeof(ret_));
le_bytes_from_limbs(ret, ret_, sizeof(ret_));
return zero^1;
}
#define SUB_N_CHECK_MOD_IMPL(bits) \
inline limb_t sub_n_check_mod_##bits(pow##bits ret, const pow##bits a, \
const pow##bits b, const vec##bits p) \
{ return sub_n_check_mod_n(ret, a, b, p, NLIMBS(bits)); }
SUB_N_CHECK_MOD_IMPL(256)
static void from_mont_n(limb_t ret[], const limb_t a[],
const limb_t p[], limb_t n0, size_t n)
{
llimb_t limbx;
limb_t mask, borrow, mx, hi, tmp[n];
size_t i, j;
for (j=0; j<n; j++) {
mx = n0*a[0];
limbx = (mx * (llimb_t)p[0]) + a[0];
hi = (limb_t)(limbx >> LIMB_T_BITS);
for (i=1; i<n; i++) {
limbx = (mx * (llimb_t)p[i] + hi) + a[i];
tmp[i-1] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
tmp[i-1] = hi;
a = tmp;
}
/* this is needed only if input can be non-fully-reduced */
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = 0 - borrow;
for(i=0; i<n; i++)
ret[i] = (ret[i] & ~mask) | (tmp[i] & mask);
}
#define FROM_MONT_IMPL(bits) \
inline void from_mont_##bits(vec##bits ret, const vec##bits a, \
const vec##bits p, limb_t n0) \
{ from_mont_n(ret, a, p, n0, NLIMBS(bits)); }
FROM_MONT_IMPL(256)
FROM_MONT_IMPL(384)
static void redc_mont_n(limb_t ret[], const limb_t a[],
const limb_t p[], limb_t n0, size_t n)
{
llimb_t limbx;
limb_t mask, carry, borrow, mx, hi, tmp[n];
const limb_t *b = a;
size_t i, j;
for (j=0; j<n; j++) {
mx = n0*b[0];
limbx = (mx * (llimb_t)p[0]) + b[0];
hi = (limb_t)(limbx >> LIMB_T_BITS);
for (i=1; i<n; i++) {
limbx = (mx * (llimb_t)p[i] + hi) + b[i];
tmp[i-1] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
tmp[i-1] = hi;
b = tmp;
}
for (carry=0, i=0; i<n; i++) {
limbx = a[n+i] + (tmp[i] + (llimb_t)carry);
tmp[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
ret[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = carry - borrow;
for(i=0; i<n; i++)
ret[i] = (ret[i] & ~mask) | (tmp[i] & mask);
}
#define REDC_MONT_IMPL(bits, bits2) \
inline void redc_mont_##bits(vec##bits ret, const vec##bits2 a, \
const vec##bits p, limb_t n0) \
{ redc_mont_n(ret, a, p, n0, NLIMBS(bits)); }
REDC_MONT_IMPL(256, 512)
REDC_MONT_IMPL(384, 768)
static void rshift_mod_n(limb_t ret[], const limb_t a[], size_t count,
const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t mask, carry, limb, next;
size_t i;
while (count--) {
mask = 0 - (a[0] & 1);
for (carry=0, i=0; i<n; i++) {
limbx = a[i] + ((p[i]&mask) + (llimb_t)carry);
ret[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
for (next=ret[0], i=0; i<n-1; i++) {
limb = next >> 1;
next = ret[i+1];
ret[i] = limb | next << (LIMB_T_BITS-1);
}
ret[i] = next >> 1 | carry << (LIMB_T_BITS-1);
a = ret;
}
}
#define RSHIFT_MOD_IMPL(bits) \
inline void rshift_mod_##bits(vec##bits ret, const vec##bits a, size_t count, \
const vec##bits p) \
{ rshift_mod_n(ret, a, count, p, NLIMBS(bits)); }
RSHIFT_MOD_IMPL(256)
RSHIFT_MOD_IMPL(384)
#define DIV_BY_2_MOD_IMPL(bits) \
inline void div_by_2_mod_##bits(vec##bits ret, const vec##bits a, \
const vec##bits p) \
{ rshift_mod_n(ret, a, 1, p, NLIMBS(bits)); }
DIV_BY_2_MOD_IMPL(384)
static limb_t sgn0_pty_mod_n(const limb_t a[], const limb_t p[], size_t n)
{
llimb_t limbx;
limb_t carry, borrow, ret, tmp[n];
size_t i;
ret = a[0] & 1; /* parity */
for (carry=0, i=0; i<n; i++) {
limb_t a_i = a[i];
tmp[i] = a_i<<1 | carry;
carry = a_i>>(LIMB_T_BITS-1);
}
for (borrow=0, i=0; i<n; i++) {
limbx = tmp[i] - (p[i] + (llimb_t)borrow);
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
ret |= ((carry - borrow) & 2) ^ 2;
return ret;
}
inline limb_t sgn0_pty_mod_384(const vec384 a, const vec384 p)
{ return sgn0_pty_mod_n(a, p, NLIMBS(384)); }
inline limb_t sgn0_pty_mont_384(const vec384 a, const vec384 p, limb_t n0)
{
vec384 tmp;
from_mont_n(tmp, a, p, n0, NLIMBS(384));
return sgn0_pty_mod_n(tmp, p, NLIMBS(384));
}
inline limb_t sgn0_pty_mod_384x(const vec384x a, const vec384 p)
{
limb_t re, im, sign, prty;
re = sgn0_pty_mod_n(a[0], p, NLIMBS(384));
im = sgn0_pty_mod_n(a[1], p, NLIMBS(384));
/* a->im!=0 ? sgn0(a->im) : sgn0(a->re) */
sign = (limb_t)0 - vec_is_zero(a[1], sizeof(vec384));
sign = (re & sign) | (im & ~sign);
/* a->re==0 ? prty(a->im) : prty(a->re) */
prty = (limb_t)0 - vec_is_zero(a[0], sizeof(vec384));
prty = (im & prty) | (re & ~prty);
return (sign & 2) | (prty & 1);
}
inline limb_t sgn0_pty_mont_384x(const vec384x a, const vec384 p, limb_t n0)
{
vec384x tmp;
from_mont_n(tmp[0], a[0], p, n0, NLIMBS(384));
from_mont_n(tmp[1], a[1], p, n0, NLIMBS(384));
return sgn0_pty_mod_384x(tmp, p);
}
void mul_mont_384x(vec384x ret, const vec384x a, const vec384x b,
const vec384 p, limb_t n0)
{
vec384 aa, bb, cc;
add_mod_n(aa, a[0], a[1], p, NLIMBS(384));
add_mod_n(bb, b[0], b[1], p, NLIMBS(384));
mul_mont_n(bb, bb, aa, p, n0, NLIMBS(384));
mul_mont_n(aa, a[0], b[0], p, n0, NLIMBS(384));
mul_mont_n(cc, a[1], b[1], p, n0, NLIMBS(384));
sub_mod_n(ret[0], aa, cc, p, NLIMBS(384));
sub_mod_n(ret[1], bb, aa, p, NLIMBS(384));
sub_mod_n(ret[1], ret[1], cc, p, NLIMBS(384));
}
/*
* mul_mont_n without final conditional subtraction, which implies
* that modulus is one bit short, which in turn means that there are
* no carries to handle between iterations...
*/
static void mul_mont_nonred_n(limb_t ret[], const limb_t a[], const limb_t b[],
const limb_t p[], limb_t n0, size_t n)
{
llimb_t limbx;
limb_t mx, hi, tmp[n+1];
size_t i, j;
for (mx=b[0], hi=0, i=0; i<n; i++) {
limbx = (mx * (llimb_t)a[i]) + hi;
tmp[i] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
mx = n0*tmp[0];
tmp[i] = hi;
for (j=0; ; ) {
limbx = (mx * (llimb_t)p[0]) + tmp[0];
hi = (limb_t)(limbx >> LIMB_T_BITS);
for (i=1; i<n; i++) {
limbx = (mx * (llimb_t)p[i] + hi) + tmp[i];
tmp[i-1] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
tmp[i-1] = tmp[i] + hi;
if (++j==n)
break;
for (mx=b[j], hi=0, i=0; i<n; i++) {
limbx = (mx * (llimb_t)a[i] + hi) + tmp[i];
tmp[i] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
mx = n0*tmp[0];
tmp[i] = hi;
}
vec_copy(ret, tmp, sizeof(tmp)-sizeof(limb_t));
}
void sqr_n_mul_mont_383(vec384 ret, const vec384 a, size_t count,
const vec384 p, limb_t n0, const vec384 b)
{
while(count--) {
mul_mont_nonred_n(ret, a, a, p, n0, NLIMBS(384));
a = ret;
}
mul_mont_n(ret, ret, b, p, n0, NLIMBS(384));
}
void sqr_mont_382x(vec384x ret, const vec384x a,
const vec384 p, limb_t n0)
{
llimb_t limbx;
limb_t mask, carry, borrow;
size_t i;
vec384 t0, t1;
/* "add_mod_n(t0, a[0], a[1], p, NLIMBS(384));" */
for (carry=0, i=0; i<NLIMBS(384); i++) {
limbx = a[0][i] + (a[1][i] + (llimb_t)carry);
t0[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
/* "sub_mod_n(t1, a[0], a[1], p, NLIMBS(384));" */
for (borrow=0, i=0; i<NLIMBS(384); i++) {
limbx = a[0][i] - (a[1][i] + (llimb_t)borrow);
t1[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = 0 - borrow;
/* "mul_mont_n(ret[1], a[0], a[1], p, n0, NLIMBS(384));" */
mul_mont_nonred_n(ret[1], a[0], a[1], p, n0, NLIMBS(384));
/* "add_mod_n(ret[1], ret[1], ret[1], p, NLIMBS(384));" */
for (carry=0, i=0; i<NLIMBS(384); i++) {
limb_t a_i = ret[1][i];
ret[1][i] = a_i<<1 | carry;
carry = a_i>>(LIMB_T_BITS-1);
}
/* "mul_mont_n(ret[0], t0, t1, p, n0, NLIMBS(384));" */
mul_mont_nonred_n(ret[0], t0, t1, p, n0, NLIMBS(384));
/* account for t1's sign... */
for (borrow=0, i=0; i<NLIMBS(384); i++) {
limbx = ret[0][i] - ((t0[i] & mask) + (llimb_t)borrow);
ret[0][i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = 0 - borrow;
for (carry=0, i=0; i<NLIMBS(384); i++) {
limbx = ret[0][i] + ((p[i] & mask) + (llimb_t)carry);
ret[0][i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
}
#define MSB(x) ((x) >> (LIMB_T_BITS-1))
static size_t num_bits(limb_t l)
{
limb_t x, mask;
size_t bits = is_zero(l) ^ 1;
if (sizeof(limb_t) == 8) {
x = l >> (32 & (8*sizeof(limb_t)-1));
mask = 0 - MSB(0 - x);
bits += 32 & mask;
l ^= (x ^ l) & mask;
}
x = l >> 16;
mask = 0 - MSB(0 - x);
bits += 16 & mask;
l ^= (x ^ l) & mask;
x = l >> 8;
mask = 0 - MSB(0 - x);
bits += 8 & mask;
l ^= (x ^ l) & mask;
x = l >> 4;
mask = 0 - MSB(0 - x);
bits += 4 & mask;
l ^= (x ^ l) & mask;
x = l >> 2;
mask = 0 - MSB(0 - x);
bits += 2 & mask;
l ^= (x ^ l) & mask;
bits += l >> 1;
return bits;
}
#if defined(__clang_major__) && __clang_major__>7
__attribute__((optnone))
#endif
static limb_t lshift_2(limb_t hi, limb_t lo, size_t l)
{
size_t r = LIMB_T_BITS - l;
limb_t mask = 0 - (is_zero(l)^1);
return (hi << (l&(LIMB_T_BITS-1))) | ((lo & mask) >> (r&(LIMB_T_BITS-1)));
}
/*
* https://eprint.iacr.org/2020/972 with 'k' being LIMB_T_BITS-1.
*/
static void ab_approximation_n(limb_t a_[2], const limb_t a[],
limb_t b_[2], const limb_t b[], size_t n)
{
limb_t a_hi, a_lo, b_hi, b_lo, mask;
size_t i;
i = n-1;
a_hi = a[i], a_lo = a[i-1];
b_hi = b[i], b_lo = b[i-1];
for (i--; --i;) {
mask = 0 - is_zero(a_hi | b_hi);
a_hi = ((a_lo ^ a_hi) & mask) ^ a_hi;
b_hi = ((b_lo ^ b_hi) & mask) ^ b_hi;
a_lo = ((a[i] ^ a_lo) & mask) ^ a_lo;
b_lo = ((b[i] ^ b_lo) & mask) ^ b_lo;
}
i = LIMB_T_BITS - num_bits(a_hi | b_hi);
/* |i| can be LIMB_T_BITS if all a[2..]|b[2..] were zeros */
a_[0] = a[0], a_[1] = lshift_2(a_hi, a_lo, i);
b_[0] = b[0], b_[1] = lshift_2(b_hi, b_lo, i);
}
typedef struct { limb_t f0, g0, f1, g1; } factors;
static void inner_loop_n(factors *fg, const limb_t a_[2], const limb_t b_[2],
size_t n)
{
llimb_t limbx;
limb_t f0 = 1, g0 = 0, f1 = 0, g1 = 1;
limb_t a_lo, a_hi, b_lo, b_hi, t_lo, t_hi, odd, borrow, xorm;
a_lo = a_[0], a_hi = a_[1];
b_lo = b_[0], b_hi = b_[1];
while(n--) {
odd = 0 - (a_lo&1);
/* a_ -= b_ if a_ is odd */
t_lo = a_lo, t_hi = a_hi;
limbx = a_lo - (llimb_t)(b_lo & odd);
a_lo = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
limbx = a_hi - ((llimb_t)(b_hi & odd) + borrow);
a_hi = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS);
/* negate a_-b_ if it borrowed */
a_lo ^= borrow;
a_hi ^= borrow;
limbx = a_lo + (llimb_t)(borrow & 1);
a_lo = (limb_t)limbx;
a_hi += (limb_t)(limbx >> LIMB_T_BITS) & 1;
/* b_=a_ if a_-b_ borrowed */
b_lo = ((t_lo ^ b_lo) & borrow) ^ b_lo;
b_hi = ((t_hi ^ b_hi) & borrow) ^ b_hi;
/* exchange f0 and f1 if a_-b_ borrowed */
xorm = (f0 ^ f1) & borrow;
f0 ^= xorm;
f1 ^= xorm;
/* exchange g0 and g1 if a_-b_ borrowed */
xorm = (g0 ^ g1) & borrow;
g0 ^= xorm;
g1 ^= xorm;
/* subtract if a_ was odd */
f0 -= f1 & odd;
g0 -= g1 & odd;
f1 <<= 1;
g1 <<= 1;
a_lo >>= 1; a_lo |= a_hi << (LIMB_T_BITS-1);
a_hi >>= 1;
}
fg->f0 = f0, fg->g0 = g0, fg->f1 = f1, fg->g1= g1;
}
static limb_t cneg_n(limb_t ret[], const limb_t a[], limb_t neg, size_t n)
{
llimb_t limbx = 0;
limb_t carry;
size_t i;
for (carry=neg&1, i=0; i<n; i++) {
limbx = (llimb_t)(a[i] ^ neg) + carry;
ret[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
return 0 - MSB((limb_t)limbx);
}
static limb_t add_n(limb_t ret[], const limb_t a[], limb_t b[], size_t n)
{
llimb_t limbx;
limb_t carry;
size_t i;
for (carry=0, i=0; i<n; i++) {
limbx = a[i] + (b[i] + (llimb_t)carry);
ret[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
return carry;
}
static limb_t umul_n(limb_t ret[], const limb_t a[], limb_t b, size_t n)
{
llimb_t limbx;
limb_t hi;
size_t i;
for (hi=0, i=0; i<n; i++) {
limbx = (b * (llimb_t)a[i]) + hi;
ret[i] = (limb_t)limbx;
hi = (limb_t)(limbx >> LIMB_T_BITS);
}
return hi;
}
static limb_t smul_n_shift_n(limb_t ret[], const limb_t a[], limb_t *f_,
const limb_t b[], limb_t *g_,
size_t n)
{
limb_t a_[n+1], b_[n+1], f, g, neg, carry, hi;
size_t i;
/* |a|*|f_| */
f = *f_;
neg = 0 - MSB(f);
f = (f ^ neg) - neg; /* ensure |f| is positive */
(void)cneg_n(a_, a, neg, n);
hi = umul_n(a_, a_, f, n);
a_[n] = hi - (f & neg);
/* |b|*|g_| */
g = *g_;
neg = 0 - MSB(g);
g = (g ^ neg) - neg; /* ensure |g| is positive */
(void)cneg_n(b_, b, neg, n);
hi = umul_n(b_, b_, g, n);
b_[n] = hi - (g & neg);
/* |a|*|f_| + |b|*|g_| */
(void)add_n(a_, a_, b_, n+1);
/* (|a|*|f_| + |b|*|g_|) >> k */
for (carry=a_[0], i=0; i<n; i++) {
hi = carry >> (LIMB_T_BITS-2);
carry = a_[i+1];
ret[i] = hi | (carry << 2);
}
/* ensure result is non-negative, fix up |f_| and |g_| accordingly */
neg = 0 - MSB(carry);
*f_ = (*f_ ^ neg) - neg;
*g_ = (*g_ ^ neg) - neg;
(void)cneg_n(ret, ret, neg, n);
return neg;
}
static limb_t smul_2n(limb_t ret[], const limb_t u[], limb_t f,
const limb_t v[], limb_t g, size_t n)
{
limb_t u_[n], v_[n], neg, hi;
/* |u|*|f_| */
neg = 0 - MSB(f);
f = (f ^ neg) - neg; /* ensure |f| is positive */
neg = cneg_n(u_, u, neg, n);
hi = umul_n(u_, u_, f, n) - (f&neg);
/* |v|*|g_| */
neg = 0 - MSB(g);
g = (g ^ neg) - neg; /* ensure |g| is positive */
neg = cneg_n(v_, v, neg, n);
hi += umul_n(v_, v_, g, n) - (g&neg);
/* |u|*|f_| + |v|*|g_| */
hi += add_n(ret, u_, v_, n);
return hi;
}
static void ct_inverse_mod_n(limb_t ret[], const limb_t inp[],
const limb_t mod[], const limb_t modx[], size_t n)
{
llimb_t limbx;
limb_t a[n], b[n], u[2*n], v[2*n], t[2*n];
limb_t a_[2], b_[2], sign, carry, top;
factors fg;
size_t i;
vec_copy(a, inp, sizeof(a));
vec_copy(b, mod, sizeof(b));
vec_zero(u, sizeof(u)); u[0] = 1;
vec_zero(v, sizeof(v));
for (i=0; i<(2*n*LIMB_T_BITS)/(LIMB_T_BITS-2); i++) {
ab_approximation_n(a_, a, b_, b, n);
inner_loop_n(&fg, a_, b_, LIMB_T_BITS-2);
(void)smul_n_shift_n(t, a, &fg.f0, b, &fg.g0, n);
(void)smul_n_shift_n(b, a, &fg.f1, b, &fg.g1, n);
vec_copy(a, t, sizeof(a));
smul_2n(t, u, fg.f0, v, fg.g0, 2*n);
smul_2n(v, u, fg.f1, v, fg.g1, 2*n);
vec_copy(u, t, sizeof(u));
}
inner_loop_n(&fg, a, b, (2*n*LIMB_T_BITS)%(LIMB_T_BITS-2));
top = smul_2n(ret, u, fg.f1, v, fg.g1, 2*n);
sign = 0 - MSB(top); /* top is 1, 0 or -1 */
for (carry=0, i=0; i<n; i++) {
limbx = ret[n+i] + ((modx[i] & sign) + (llimb_t)carry);
ret[n+i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
top += carry;
sign = 0 - top; /* top is 1, 0 or -1 */
top |= sign;
for (i=0; i<n; i++)
a[i] = modx[i] & top;
(void)cneg_n(a, a, 0 - MSB(sign), n);
add_n(ret+n, ret+n, a, n);
}
#define CT_INVERSE_MOD_IMPL(bits) \
inline void ct_inverse_mod_##bits(vec##bits ret, const vec##bits inp, \
const vec##bits mod, const vec##bits modx) \
{ ct_inverse_mod_n(ret, inp, mod, modx, NLIMBS(bits)); }
CT_INVERSE_MOD_IMPL(256)
CT_INVERSE_MOD_IMPL(384)
/*
* Copy of inner_loop_n above, but with |L| updates.
*/
static limb_t legendre_loop_n(limb_t L, factors *fg, const limb_t a_[2],
const limb_t b_[2], size_t n)
{
llimb_t limbx;
limb_t f0 = 1, g0 = 0, f1 = 0, g1 = 1;
limb_t a_lo, a_hi, b_lo, b_hi, t_lo, t_hi, odd, borrow, xorm;
a_lo = a_[0], a_hi = a_[1];
b_lo = b_[0], b_hi = b_[1];
while(n--) {
odd = 0 - (a_lo&1);
/* a_ -= b_ if a_ is odd */
t_lo = a_lo, t_hi = a_hi;
limbx = a_lo - (llimb_t)(b_lo & odd);
a_lo = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
limbx = a_hi - ((llimb_t)(b_hi & odd) + borrow);
a_hi = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS);
L += ((t_lo & b_lo) >> 1) & borrow;
/* negate a_-b_ if it borrowed */
a_lo ^= borrow;
a_hi ^= borrow;
limbx = a_lo + (llimb_t)(borrow & 1);
a_lo = (limb_t)limbx;
a_hi += (limb_t)(limbx >> LIMB_T_BITS) & 1;
/* b_=a_ if a_-b_ borrowed */
b_lo = ((t_lo ^ b_lo) & borrow) ^ b_lo;
b_hi = ((t_hi ^ b_hi) & borrow) ^ b_hi;
/* exchange f0 and f1 if a_-b_ borrowed */
xorm = (f0 ^ f1) & borrow;
f0 ^= xorm;
f1 ^= xorm;
/* exchange g0 and g1 if a_-b_ borrowed */
xorm = (g0 ^ g1) & borrow;
g0 ^= xorm;
g1 ^= xorm;
/* subtract if a_ was odd */
f0 -= f1 & odd;
g0 -= g1 & odd;
f1 <<= 1;
g1 <<= 1;
a_lo >>= 1; a_lo |= a_hi << (LIMB_T_BITS-1);
a_hi >>= 1;
L += (b_lo + 2) >> 2;
}
fg->f0 = f0, fg->g0 = g0, fg->f1 = f1, fg->g1 = g1;
return L;
}
static bool_t ct_is_sqr_mod_n(const limb_t inp[], const limb_t mod[], size_t n)
{
limb_t a[n], b[n], t[n];
limb_t a_[2], b_[2], neg, L = 0;
factors fg;
size_t i;
vec_copy(a, inp, sizeof(a));
vec_copy(b, mod, sizeof(b));
for (i=0; i<(2*n*LIMB_T_BITS)/(LIMB_T_BITS-2); i++) {
ab_approximation_n(a_, a, b_, b, n);
L = legendre_loop_n(L, &fg, a_, b_, LIMB_T_BITS-2);
neg = smul_n_shift_n(t, a, &fg.f0, b, &fg.g0, n);
(void)smul_n_shift_n(b, a, &fg.f1, b, &fg.g1, n);
vec_copy(a, t, sizeof(a));
L += (b[0] >> 1) & neg;
}
L = legendre_loop_n(L, &fg, a, b, (2*n*LIMB_T_BITS)%(LIMB_T_BITS-2));
return (L & 1) ^ 1;
}
#define CT_IS_SQR_MOD_IMPL(bits) \
inline bool_t ct_is_square_mod_##bits(const vec##bits inp, \
const vec##bits mod) \
{ return ct_is_sqr_mod_n(inp, mod, NLIMBS(bits)); }
CT_IS_SQR_MOD_IMPL(384)
/*
* |div_top| points at two most significant limbs of the dividend, |d_hi|
* and |d_lo| are two most significant limbs of the divisor. If divisor
* is only one limb, it is to be passed in |d_hi| with zero in |d_lo|.
* The divisor is required to be "bitwise left-aligned," and dividend's
* top limbs to be not larger than the divisor's. The latter limitation
* can be problematic in the first iteration of multi-precision division,
* where in most general case the condition would have to be "smaller."
* The subroutine considers four limbs, two of which are "overlapping,"
* hence the name... Another way to look at it is to think of the pair
* of the dividend's limbs being suffixed with a zero:
* +-------+-------+-------+
* R | | | 0 |
* +-------+-------+-------+
* +-------+-------+
* D | | |
* +-------+-------+
*/
limb_t div_3_limbs(const limb_t div_top[2], limb_t d_lo, limb_t d_hi)
{
llimb_t Rx;
limb_t r_lo = div_top[0], r_hi = div_top[1];
limb_t Q = 0, mask, borrow, rx;
size_t i;
for (i = 0; i < LIMB_T_BITS; i++) {
/* "borrow, Rx = R - D" */
Rx = (llimb_t)r_lo - d_lo;
rx = (limb_t)Rx;
borrow = (limb_t)(Rx >> LIMB_T_BITS) & 1;
Rx = r_hi - (d_hi + (llimb_t)borrow);
borrow = (limb_t)(Rx >> LIMB_T_BITS);
/* "if (R >= D) R -= D" */
r_lo = ((r_lo ^ rx) & borrow) ^ rx;
rx = (limb_t)Rx;
r_hi = ((r_hi ^ rx) & borrow) ^ rx;
Q <<= 1;
Q |= ~borrow & 1;
/* "D >>= 1" */
d_lo >>= 1; d_lo |= d_hi << (LIMB_T_BITS - 1);
d_hi >>= 1;
}
mask = 0 - MSB(Q); /* does it overflow? */
/* "borrow, Rx = R - D" */
Rx = (llimb_t)r_lo - d_lo;
rx = (limb_t)Rx;
borrow = (limb_t)(Rx >> LIMB_T_BITS) & 1;
Rx = r_hi - (d_hi + (llimb_t)borrow);
borrow = (limb_t)(Rx >> LIMB_T_BITS) & 1;
Q <<= 1;
Q |= borrow ^ 1;
return (Q | mask);
}
static limb_t quot_rem_n(limb_t *div_rem, const limb_t *divisor,
limb_t quotient, size_t n)
{
llimb_t limbx;
limb_t tmp[n+1], carry, mask, borrow;
size_t i;
/* divisor*quotient */
for (carry=0, i=0; i<n; i++) {
limbx = (quotient * (llimb_t)divisor[i]) + carry;
tmp[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS);
}
tmp[i] = carry;
/* remainder = dividend - divisor*quotient */
for (borrow=0, i=0; i<=n; i++) {
limbx = div_rem[i] - (tmp[i] + (llimb_t)borrow);
tmp[i] = (limb_t)limbx;
borrow = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
mask = 0 - borrow;
/* if quotient was off by one, add divisor to the remainder */
for (carry=0, i=0; i<n; i++) {
limbx = tmp[i] + ((divisor[i] & mask) + (llimb_t)carry);
div_rem[i] = (limb_t)limbx;
carry = (limb_t)(limbx >> LIMB_T_BITS) & 1;
}
return (div_rem[i] = quotient + mask);
}
inline limb_t quot_rem_128(limb_t *div_rem, const limb_t *divisor,
limb_t quotient)
{ return quot_rem_n(div_rem, divisor, quotient, NLIMBS(128)); }
inline limb_t quot_rem_64(limb_t *div_rem, const limb_t *divisor,
limb_t quotient)
{ return quot_rem_n(div_rem, divisor, quotient, NLIMBS(64)); }
/*
* Unlock reference implementations in vect.c
*/
#define mul_by_8_mod_384 mul_by_8_mod_384
#define mul_by_8_mod_384x mul_by_8_mod_384x
#define mul_by_3_mod_384x mul_by_3_mod_384x
#define mul_by_1_plus_i_mod_384x mul_by_1_plus_i_mod_384x
#define add_mod_384x add_mod_384x
#define sub_mod_384x sub_mod_384x
#define lshift_mod_384x lshift_mod_384x
#define sqr_mont_384x sqr_mont_384x
inline void vec_prefetch(const void *ptr, size_t len)
{ (void)ptr; (void)len; }
/*
* SHA-256
*/
#define ROTR(x,n) ((x)>>n | (x)<<(32-n))
#define Sigma0(x) (ROTR((x),2) ^ ROTR((x),13) ^ ROTR((x),22))
#define Sigma1(x) (ROTR((x),6) ^ ROTR((x),11) ^ ROTR((x),25))
#define sigma0(x) (ROTR((x),7) ^ ROTR((x),18) ^ ((x)>>3))
#define sigma1(x) (ROTR((x),17) ^ ROTR((x),19) ^ ((x)>>10))
#define Ch(x,y,z) (((x) & (y)) ^ ((~(x)) & (z)))
#define Maj(x,y,z) (((x) & (y)) ^ ((x) & (z)) ^ ((y) & (z)))
void blst_sha256_block_data_order(unsigned int *v, const void *inp,
size_t blocks)
{
static const unsigned int K256[64] = {
0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5,
0x3956c25b, 0x59f111f1, 0x923f82a4, 0xab1c5ed5,
0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3,
0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174,
0xe49b69c1, 0xefbe4786, 0x0fc19dc6, 0x240ca1cc,
0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7,
0xc6e00bf3, 0xd5a79147, 0x06ca6351, 0x14292967,
0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13,
0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85,
0xa2bfe8a1, 0xa81a664b, 0xc24b8b70, 0xc76c51a3,
0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5,
0x391c0cb3, 0x4ed8aa4a, 0x5b9cca4f, 0x682e6ff3,
0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208,
0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2
};
unsigned int X[16], l, a, b, c, d, e, f, g, h, s0, s1, T1, T2;
const unsigned char *data = inp;
size_t round;
a = v[0];
b = v[1];
c = v[2];
d = v[3];
e = v[4];
f = v[5];
g = v[6];
h = v[7];
while (blocks--) {
for (round = 0; round < 16; round++) {
l = (unsigned int)data[0] << 24;
l |= (unsigned int)data[1] << 16;
l |= (unsigned int)data[2] << 8;
l |= (unsigned int)data[3];
data += 4;
T1 = X[round] = l;
T1 += h + Sigma1(e) + Ch(e, f, g) + K256[round];
T2 = Sigma0(a) + Maj(a, b, c);
h = g;
g = f;
f = e;
e = d + T1;
d = c;
c = b;
b = a;
a = T1 + T2;
}
for (; round < 64; round++) {
s0 = X[(round + 1) & 0x0f];
s0 = sigma0(s0);
s1 = X[(round + 14) & 0x0f];
s1 = sigma1(s1);
T1 = X[round & 0xf] += s0 + s1 + X[(round + 9) & 0xf];
T1 += h + Sigma1(e) + Ch(e, f, g) + K256[round];
T2 = Sigma0(a) + Maj(a, b, c);
h = g;
g = f;
f = e;
e = d + T1;
d = c;
c = b;
b = a;
a = T1 + T2;
}
a += v[0]; v[0] = a;
b += v[1]; v[1] = b;
c += v[2]; v[2] = c;
d += v[3]; v[3] = d;
e += v[4]; v[4] = e;
f += v[5]; v[5] = f;
g += v[6]; v[6] = g;
h += v[7]; v[7] = h;
}
}
#undef ROTR
#undef Sigma0
#undef Sigma1
#undef sigma0
#undef sigma1
#undef Ch
#undef Maj
void blst_sha256_hcopy(unsigned int dst[8], const unsigned int src[8])
{
size_t i;
for (i=0; i<8; i++)
dst[i] = src[i];
}
void blst_sha256_emit(unsigned char md[32], const unsigned int h[8])
{
size_t i;
for (i=0; i<8; i++, md+=4) {
unsigned int h_i = h[i];
md[0] = (unsigned char)(h_i >> 24);
md[1] = (unsigned char)(h_i >> 16);
md[2] = (unsigned char)(h_i >> 8);
md[3] = (unsigned char)h_i;
}
}
void blst_sha256_bcopy(void *dst_, const void *src_, size_t len)
{
unsigned char *dst = dst_;
const unsigned char *src = src_;
size_t i;
for (i=0; i<len; i++)
dst[i] = src[i];
}
Reference in a new issue Copy permalink