284 lines
10 KiB
C
284 lines
10 KiB
C
/********************************************************************************************
|
|
* SIDH: an efficient supersingular isogeny cryptography library
|
|
*
|
|
* Abstract: core functions over GF(p) and GF(p^2)
|
|
*********************************************************************************************/
|
|
#include <openssl/base.h>
|
|
|
|
#include "utils.h"
|
|
#include "fpx.h"
|
|
|
|
extern const struct params_t sike_params;
|
|
|
|
// Multiprecision squaring, c = a^2 mod p.
|
|
static void fpsqr_mont(const felm_t ma, felm_t mc)
|
|
{
|
|
dfelm_t temp = {0};
|
|
sike_mpmul(ma, ma, temp);
|
|
sike_fprdc(temp, mc);
|
|
}
|
|
|
|
// Chain to compute a^(p-3)/4 using Montgomery arithmetic.
|
|
static void fpinv_chain_mont(felm_t a)
|
|
{
|
|
unsigned int i, j;
|
|
felm_t t[31], tt;
|
|
|
|
// Precomputed table
|
|
fpsqr_mont(a, tt);
|
|
sike_fpmul_mont(a, tt, t[0]);
|
|
for (i = 0; i <= 29; i++) sike_fpmul_mont(t[i], tt, t[i+1]);
|
|
|
|
sike_fpcopy(a, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[5], tt, tt);
|
|
for (i = 0; i < 10; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[14], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[3], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[23], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[13], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[24], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[7], tt, tt);
|
|
for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[12], tt, tt);
|
|
for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[30], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[1], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[30], tt, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[21], tt, tt);
|
|
for (i = 0; i < 9; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[2], tt, tt);
|
|
for (i = 0; i < 9; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[19], tt, tt);
|
|
for (i = 0; i < 9; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[1], tt, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[24], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[26], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[16], tt, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[10], tt, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[6], tt, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[0], tt, tt);
|
|
for (i = 0; i < 9; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[20], tt, tt);
|
|
for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[9], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[25], tt, tt);
|
|
for (i = 0; i < 9; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[30], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[26], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(a, tt, tt);
|
|
for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[28], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[6], tt, tt);
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[10], tt, tt);
|
|
for (i = 0; i < 9; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[22], tt, tt);
|
|
for (j = 0; j < 35; j++) {
|
|
for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(t[30], tt, tt);
|
|
}
|
|
sike_fpcopy(tt, a);
|
|
}
|
|
|
|
// Field inversion using Montgomery arithmetic, a = a^(-1)*R mod p.
|
|
static void fpinv_mont(felm_t a)
|
|
{
|
|
felm_t tt = {0};
|
|
sike_fpcopy(a, tt);
|
|
fpinv_chain_mont(tt);
|
|
fpsqr_mont(tt, tt);
|
|
fpsqr_mont(tt, tt);
|
|
sike_fpmul_mont(a, tt, a);
|
|
}
|
|
|
|
// Multiprecision addition, c = a+b, where lng(a) = lng(b) = nwords. Returns the carry bit.
|
|
#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
|
|
inline static unsigned int mp_add(const felm_t a, const felm_t b, felm_t c, const unsigned int nwords) {
|
|
uint8_t carry = 0;
|
|
for (size_t i = 0; i < nwords; i++) {
|
|
ADDC(carry, a[i], b[i], carry, c[i]);
|
|
}
|
|
return carry;
|
|
}
|
|
|
|
// Multiprecision subtraction, c = a-b, where lng(a) = lng(b) = nwords. Returns the borrow bit.
|
|
inline static unsigned int mp_sub(const felm_t a, const felm_t b, felm_t c, const unsigned int nwords) {
|
|
uint32_t borrow = 0;
|
|
for (size_t i = 0; i < nwords; i++) {
|
|
SUBC(borrow, a[i], b[i], borrow, c[i]);
|
|
}
|
|
return borrow;
|
|
}
|
|
#endif
|
|
|
|
// Multiprecision addition, c = a+b.
|
|
inline static void mp_addfast(const felm_t a, const felm_t b, felm_t c)
|
|
{
|
|
#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
|
|
mp_add(a, b, c, NWORDS_FIELD);
|
|
#else
|
|
sike_mpadd_asm(a, b, c);
|
|
#endif
|
|
}
|
|
|
|
// Multiprecision subtraction, c = a-b, where lng(a) = lng(b) = 2*NWORDS_FIELD.
|
|
// If c < 0 then returns mask = 0xFF..F, else mask = 0x00..0
|
|
inline static crypto_word_t mp_subfast(const dfelm_t a, const dfelm_t b, dfelm_t c) {
|
|
#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
|
|
return (0 - (crypto_word_t)mp_sub(a, b, c, 2*NWORDS_FIELD));
|
|
#else
|
|
return sike_mpsubx2_asm(a, b, c);
|
|
#endif
|
|
}
|
|
|
|
// Multiprecision subtraction, c = c-a-b, where lng(a) = lng(b) = 2*NWORDS_FIELD.
|
|
// Inputs should be s.t. c > a and c > b
|
|
inline static void mp_dblsubfast(const dfelm_t a, const dfelm_t b, dfelm_t c) {
|
|
#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
|
|
mp_sub(c, a, c, 2*NWORDS_FIELD);
|
|
mp_sub(c, b, c, 2*NWORDS_FIELD);
|
|
#else
|
|
sike_mpdblsubx2_asm(a, b, c);
|
|
#endif
|
|
}
|
|
|
|
// Copy a field element, c = a.
|
|
void sike_fpcopy(const felm_t a, felm_t c) {
|
|
for (size_t i = 0; i < NWORDS_FIELD; i++) {
|
|
c[i] = a[i];
|
|
}
|
|
}
|
|
|
|
// Field multiplication using Montgomery arithmetic, c = a*b*R^-1 mod prime, where R=2^768
|
|
void sike_fpmul_mont(const felm_t ma, const felm_t mb, felm_t mc)
|
|
{
|
|
dfelm_t temp = {0};
|
|
sike_mpmul(ma, mb, temp);
|
|
sike_fprdc(temp, mc);
|
|
}
|
|
|
|
// Conversion from Montgomery representation to standard representation,
|
|
// c = ma*R^(-1) mod p = a mod p, where ma in [0, p-1].
|
|
void sike_from_mont(const felm_t ma, felm_t c)
|
|
{
|
|
felm_t one = {0};
|
|
one[0] = 1;
|
|
|
|
sike_fpmul_mont(ma, one, c);
|
|
sike_fpcorrection(c);
|
|
}
|
|
|
|
// GF(p^2) squaring using Montgomery arithmetic, c = a^2 in GF(p^2).
|
|
// Inputs: a = a0+a1*i, where a0, a1 are in [0, 2*p-1]
|
|
// Output: c = c0+c1*i, where c0, c1 are in [0, 2*p-1]
|
|
void sike_fp2sqr_mont(const f2elm_t a, f2elm_t c) {
|
|
felm_t t1, t2, t3;
|
|
|
|
mp_addfast(a->c0, a->c1, t1); // t1 = a0+a1
|
|
sike_fpsub(a->c0, a->c1, t2); // t2 = a0-a1
|
|
mp_addfast(a->c0, a->c0, t3); // t3 = 2a0
|
|
sike_fpmul_mont(t1, t2, c->c0); // c0 = (a0+a1)(a0-a1)
|
|
sike_fpmul_mont(t3, a->c1, c->c1); // c1 = 2a0*a1
|
|
}
|
|
|
|
// Modular negation, a = -a mod p503.
|
|
// Input/output: a in [0, 2*p503-1]
|
|
void sike_fpneg(felm_t a) {
|
|
uint32_t borrow = 0;
|
|
for (size_t i = 0; i < NWORDS_FIELD; i++) {
|
|
SUBC(borrow, sike_params.prime_x2[i], a[i], borrow, a[i]);
|
|
}
|
|
}
|
|
|
|
// Modular division by two, c = a/2 mod p503.
|
|
// Input : a in [0, 2*p503-1]
|
|
// Output: c in [0, 2*p503-1]
|
|
void sike_fpdiv2(const felm_t a, felm_t c) {
|
|
uint32_t carry = 0;
|
|
crypto_word_t mask;
|
|
|
|
mask = 0 - (crypto_word_t)(a[0] & 1); // If a is odd compute a+p503
|
|
for (size_t i = 0; i < NWORDS_FIELD; i++) {
|
|
ADDC(carry, a[i], sike_params.prime[i] & mask, carry, c[i]);
|
|
}
|
|
|
|
// Multiprecision right shift by one.
|
|
for (size_t i = 0; i < NWORDS_FIELD-1; i++) {
|
|
c[i] = (c[i] >> 1) ^ (c[i+1] << (RADIX - 1));
|
|
}
|
|
c[NWORDS_FIELD-1] >>= 1;
|
|
}
|
|
|
|
// Modular correction to reduce field element a in [0, 2*p503-1] to [0, p503-1].
|
|
void sike_fpcorrection(felm_t a) {
|
|
uint32_t borrow = 0;
|
|
crypto_word_t mask;
|
|
|
|
for (size_t i = 0; i < NWORDS_FIELD; i++) {
|
|
SUBC(borrow, a[i], sike_params.prime[i], borrow, a[i]);
|
|
}
|
|
mask = 0 - (crypto_word_t)borrow;
|
|
|
|
borrow = 0;
|
|
for (size_t i = 0; i < NWORDS_FIELD; i++) {
|
|
ADDC(borrow, a[i], sike_params.prime[i] & mask, borrow, a[i]);
|
|
}
|
|
}
|
|
|
|
// GF(p^2) multiplication using Montgomery arithmetic, c = a*b in GF(p^2).
|
|
// Inputs: a = a0+a1*i and b = b0+b1*i, where a0, a1, b0, b1 are in [0, 2*p-1]
|
|
// Output: c = c0+c1*i, where c0, c1 are in [0, 2*p-1]
|
|
void sike_fp2mul_mont(const f2elm_t a, const f2elm_t b, f2elm_t c) {
|
|
felm_t t1, t2;
|
|
dfelm_t tt1, tt2, tt3;
|
|
crypto_word_t mask;
|
|
|
|
mp_addfast(a->c0, a->c1, t1); // t1 = a0+a1
|
|
mp_addfast(b->c0, b->c1, t2); // t2 = b0+b1
|
|
sike_mpmul(a->c0, b->c0, tt1); // tt1 = a0*b0
|
|
sike_mpmul(a->c1, b->c1, tt2); // tt2 = a1*b1
|
|
sike_mpmul(t1, t2, tt3); // tt3 = (a0+a1)*(b0+b1)
|
|
mp_dblsubfast(tt1, tt2, tt3); // tt3 = (a0+a1)*(b0+b1) - a0*b0 - a1*b1
|
|
mask = mp_subfast(tt1, tt2, tt1); // tt1 = a0*b0 - a1*b1. If tt1 < 0 then mask = 0xFF..F, else if tt1 >= 0 then mask = 0x00..0
|
|
|
|
for (size_t i = 0; i < NWORDS_FIELD; i++) {
|
|
t1[i] = sike_params.prime[i] & mask;
|
|
}
|
|
|
|
sike_fprdc(tt3, c->c1); // c[1] = (a0+a1)*(b0+b1) - a0*b0 - a1*b1
|
|
mp_addfast(&tt1[NWORDS_FIELD], t1, &tt1[NWORDS_FIELD]);
|
|
sike_fprdc(tt1, c->c0); // c[0] = a0*b0 - a1*b1
|
|
}
|
|
|
|
// GF(p^2) inversion using Montgomery arithmetic, a = (a0-i*a1)/(a0^2+a1^2).
|
|
void sike_fp2inv_mont(f2elm_t a) {
|
|
f2elm_t t1;
|
|
|
|
fpsqr_mont(a->c0, t1->c0); // t10 = a0^2
|
|
fpsqr_mont(a->c1, t1->c1); // t11 = a1^2
|
|
sike_fpadd(t1->c0, t1->c1, t1->c0); // t10 = a0^2+a1^2
|
|
fpinv_mont(t1->c0); // t10 = (a0^2+a1^2)^-1
|
|
sike_fpneg(a->c1); // a = a0-i*a1
|
|
sike_fpmul_mont(a->c0, t1->c0, a->c0);
|
|
sike_fpmul_mont(a->c1, t1->c0, a->c1); // a = (a0-i*a1)*(a0^2+a1^2)^-1
|
|
}
|