kris
/
pqc
1
1
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
Browse Source

pull SIKE

sike
Henry Case 3 years ago
parent
188e672764
commit
ef4895b556
9 changed files with 2855 additions and 0 deletions
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +1095
    -0
      src/kem/sike/sike-p434-sha256/asm/fp-x86_64.S
  2. +179
    -0
      src/kem/sike/sike-p434-sha256/asm/fp_generic.c
  3. +282
    -0
      src/kem/sike/sike-p434-sha256/fpx.c
  4. +112
    -0
      src/kem/sike/sike-p434-sha256/fpx.h
  5. +262
    -0
      src/kem/sike/sike-p434-sha256/isogeny.c
  6. +49
    -0
      src/kem/sike/sike-p434-sha256/isogeny.h
  7. +128
    -0
      src/kem/sike/sike-p434-sha256/params.c
  8. +517
    -0
      src/kem/sike/sike-p434-sha256/sike.c
  9. +231
    -0
      src/kem/sike/sike-p434-sha256/utils.h

+ 1095
- 0
src/kem/sike/sike-p434-sha256/asm/fp-x86_64.S
File diff suppressed because it is too large
View File


+ 179
- 0
src/kem/sike/sike-p434-sha256/asm/fp_generic.c View File

@@ -0,0 +1,179 @@
/********************************************************************************************
* SIDH: an efficient supersingular isogeny cryptography library
*
* Abstract: portable modular arithmetic for P503
*********************************************************************************************/

#if defined(ARCH_GENERIC) || \
(!defined(ARCH_X86_64) && !defined(ARCH_AARCH64))

#include "../utils.h"
#include "../fpx.h"

// Global constants
extern const struct params_t params;

static void digit_x_digit(const crypto_word_t a, const crypto_word_t b, crypto_word_t* c)
{ // Digit multiplication, digit * digit -> 2-digit result
crypto_word_t al, ah, bl, bh, temp;
crypto_word_t albl, albh, ahbl, ahbh, res1, res2, res3, carry;
crypto_word_t mask_low = (crypto_word_t)(-1) >> (sizeof(crypto_word_t)*4);
crypto_word_t mask_high = (crypto_word_t)(-1) << (sizeof(crypto_word_t)*4);

al = a & mask_low; // Low part
ah = a >> (sizeof(crypto_word_t) * 4); // High part
bl = b & mask_low;
bh = b >> (sizeof(crypto_word_t) * 4);

albl = al*bl;
albh = al*bh;
ahbl = ah*bl;
ahbh = ah*bh;
c[0] = albl & mask_low; // C00

res1 = albl >> (sizeof(crypto_word_t) * 4);
res2 = ahbl & mask_low;
res3 = albh & mask_low;
temp = res1 + res2 + res3;
carry = temp >> (sizeof(crypto_word_t) * 4);
c[0] ^= temp << (sizeof(crypto_word_t) * 4); // C01

res1 = ahbl >> (sizeof(crypto_word_t) * 4);
res2 = albh >> (sizeof(crypto_word_t) * 4);
res3 = ahbh & mask_low;
temp = res1 + res2 + res3 + carry;
c[1] = temp & mask_low; // C10
carry = temp & mask_high;
c[1] ^= (ahbh & mask_high) + carry; // C11
}

void sike_fpadd(const felm_t a, const felm_t b, felm_t c)
{ // Modular addition, c = a+b mod p434.
// Inputs: a, b in [0, 2*p434-1]
// Output: c in [0, 2*p434-1]
unsigned int i, carry = 0;
crypto_word_t mask;

for (i = 0; i < NWORDS_FIELD; i++) {
ADDC(carry, a[i], b[i], carry, c[i]);
}

carry = 0;
for (i = 0; i < NWORDS_FIELD; i++) {
SUBC(carry, c[i], params.prime_x2[i], carry, c[i]);
}
mask = 0 - (crypto_word_t)carry;

carry = 0;
for (i = 0; i < NWORDS_FIELD; i++) {
ADDC(carry, c[i], params.prime_x2[i] & mask, carry, c[i]);
}
}

void sike_fpsub(const felm_t a, const felm_t b, felm_t c)
{ // Modular subtraction, c = a-b mod p434.
// Inputs: a, b in [0, 2*p434-1]
// Output: c in [0, 2*p434-1]
unsigned int i, borrow = 0;
crypto_word_t mask;

for (i = 0; i < NWORDS_FIELD; i++) {
SUBC(borrow, a[i], b[i], borrow, c[i]);
}
mask = 0 - (crypto_word_t)borrow;

borrow = 0;
for (i = 0; i < NWORDS_FIELD; i++) {
ADDC(borrow, c[i], params.prime_x2[i] & mask, borrow, c[i]);
}
}

void sike_mpmul(const felm_t a, const felm_t b, dfelm_t c)
{ // Multiprecision comba multiply, c = a*b, where lng(a) = lng(b) = NWORDS_FIELD.
unsigned int i, j;
crypto_word_t t = 0, u = 0, v = 0, UV[2];
unsigned int carry = 0;

for (i = 0; i < NWORDS_FIELD; i++) {
for (j = 0; j <= i; j++) {
MUL(a[j], b[i-j], UV+1, UV[0]);
ADDC(0, UV[0], v, carry, v);
ADDC(carry, UV[1], u, carry, u);
t += carry;
}
c[i] = v;
v = u;
u = t;
t = 0;
}

for (i = NWORDS_FIELD; i < 2*NWORDS_FIELD-1; i++) {
for (j = i-NWORDS_FIELD+1; j < NWORDS_FIELD; j++) {
MUL(a[j], b[i-j], UV+1, UV[0]);
ADDC(0, UV[0], v, carry, v);
ADDC(carry, UV[1], u, carry, u);
t += carry;
}
c[i] = v;
v = u;
u = t;
t = 0;
}
c[2*NWORDS_FIELD-1] = v;
}

void sike_fprdc(const felm_t ma, felm_t mc)
{ // Efficient Montgomery reduction using comba and exploiting the special form of the prime p434.
// mc = ma*R^-1 mod p434x2, where R = 2^448.
// If ma < 2^448*p434, the output mc is in the range [0, 2*p434-1].
// ma is assumed to be in Montgomery representation.
unsigned int i, j, carry, count = ZERO_WORDS;
crypto_word_t UV[2], t = 0, u = 0, v = 0;

for (i = 0; i < NWORDS_FIELD; i++) {
mc[i] = 0;
}

for (i = 0; i < NWORDS_FIELD; i++) {
for (j = 0; j < i; j++) {
if (j < (i-ZERO_WORDS+1)) {
MUL(mc[j], params.prime_p1[i-j], UV+1, UV[0]);
ADDC(0, UV[0], v, carry, v);
ADDC(carry, UV[1], u, carry, u);
t += carry;
}
}
ADDC(0, v, ma[i], carry, v);
ADDC(carry, u, 0, carry, u);
t += carry;
mc[i] = v;
v = u;
u = t;
t = 0;
}

for (i = NWORDS_FIELD; i < 2*NWORDS_FIELD-1; i++) {
if (count > 0) {
count -= 1;
}
for (j = i-NWORDS_FIELD+1; j < NWORDS_FIELD; j++) {
if (j < (NWORDS_FIELD-count)) {
MUL(mc[j], params.prime_p1[i-j], UV+1, UV[0]);
ADDC(0, UV[0], v, carry, v);
ADDC(carry, UV[1], u, carry, u);
t += carry;
}
}
ADDC(0, v, ma[i], carry, v);
ADDC(carry, u, 0, carry, u);
t += carry;
mc[i-NWORDS_FIELD] = v;
v = u;
u = t;
t = 0;
}
ADDC(0, v, ma[2*NWORDS_FIELD-1], carry, v);
mc[NWORDS_FIELD-1] = v;
}

#endif // NO_ASM || (!X86_64 && !AARCH64)

+ 282
- 0
src/kem/sike/sike-p434-sha256/fpx.c View File

@@ -0,0 +1,282 @@
/********************************************************************************************
* SIDH: an efficient supersingular isogeny cryptography library
*
* Abstract: core functions over GF(p) and GF(p^2)
*********************************************************************************************/
#include <stddef.h>
#include "utils.h"
#include "fpx.h"

extern const struct params_t 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(ARCH_GENERIC) || (!defined(ARCH_X86_64) && !defined(ARCH_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(ARCH_GENERIC) || (!defined(ARCH_X86_64) && !defined(ARCH_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(ARCH_GENERIC) || (!defined(ARCH_X86_64) && !defined(ARCH_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(ARCH_GENERIC) || (!defined(ARCH_X86_64) && !defined(ARCH_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, 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], 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], 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], 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] = 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
}

+ 112
- 0
src/kem/sike/sike-p434-sha256/fpx.h View File

@@ -0,0 +1,112 @@
#ifndef FPX_H_
#define FPX_H_

#include "utils.h"

#if defined(__cplusplus)
extern "C" {
#endif

// Modular addition, c = a+b mod p.
void sike_fpadd(const felm_t a, const felm_t b, felm_t c);
// Modular subtraction, c = a-b mod p.
void sike_fpsub(const felm_t a, const felm_t b, felm_t c);
// Modular division by two, c = a/2 mod p.
void sike_fpdiv2(const felm_t a, felm_t c);
// Modular correction to reduce field element a in [0, 2*p-1] to [0, p-1].
void sike_fpcorrection(felm_t a);
// Multiprecision multiply, c = a*b, where lng(a) = lng(b) = nwords.
void sike_mpmul(const felm_t a, const felm_t b, dfelm_t c);
// 443-bit Montgomery reduction, c = a mod p
void sike_fprdc(const dfelm_t a, felm_t c);
// Double 2x443-bit multiprecision subtraction, c = c-a-b
void sike_mpdblsubx2_asm(const felm_t a, const felm_t b, felm_t c);
// Multiprecision subtraction, c = a-b
crypto_word_t sike_mpsubx2_asm(const dfelm_t a, const dfelm_t b, dfelm_t c);
// 443-bit multiprecision addition, c = a+b
void sike_mpadd_asm(const felm_t a, const felm_t b, felm_t c);
// Modular negation, a = -a mod p.
void sike_fpneg(felm_t a);
// Copy of a field element, c = a
void sike_fpcopy(const felm_t a, felm_t c);
// Copy a field element, c = a.
void sike_fpzero(felm_t a);
// If option = 0xFF...FF x=y; y=x, otherwise swap doesn't happen. Constant time.
void sike_cswap_asm(point_proj_t x, point_proj_t y, const crypto_word_t option);
// 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);
// Field multiplication using Montgomery arithmetic, c = a*b*R^-1 mod p443, where R=2^768
void sike_fpmul_mont(const felm_t ma, const felm_t mb, felm_t mc);
// GF(p443^2) multiplication using Montgomery arithmetic, c = a*b in GF(p443^2)
void sike_fp2mul_mont(const f2elm_t a, const f2elm_t b, f2elm_t c);
// GF(p443^2) inversion using Montgomery arithmetic, a = (a0-i*a1)/(a0^2+a1^2)
void sike_fp2inv_mont(f2elm_t a);
// GF(p^2) squaring using Montgomery arithmetic, c = a^2 in GF(p^2).
void sike_fp2sqr_mont(const f2elm_t a, f2elm_t c);
// Modular correction, a = a in GF(p^2).
void sike_fp2correction(f2elm_t a);

#if defined(__cplusplus)
} // extern C
#endif

// GF(p^2) addition, c = a+b in GF(p^2).
#define sike_fp2add(a, b, c) \
do { \
sike_fpadd(a->c0, b->c0, c->c0); \
sike_fpadd(a->c1, b->c1, c->c1); \
} while(0)

// GF(p^2) subtraction, c = a-b in GF(p^2).
#define sike_fp2sub(a,b,c) \
do { \
sike_fpsub(a->c0, b->c0, c->c0); \
sike_fpsub(a->c1, b->c1, c->c1); \
} while(0)

// Copy a GF(p^2) element, c = a.
#define sike_fp2copy(a, c) \
do { \
sike_fpcopy(a->c0, c->c0); \
sike_fpcopy(a->c1, c->c1); \
} while(0)

// GF(p^2) negation, a = -a in GF(p^2).
#define sike_fp2neg(a) \
do { \
sike_fpneg(a->c0); \
sike_fpneg(a->c1); \
} while(0)

// GF(p^2) division by two, c = a/2 in GF(p^2).
#define sike_fp2div2(a, c) \
do { \
sike_fpdiv2(a->c0, c->c0); \
sike_fpdiv2(a->c1, c->c1); \
} while(0)

// Modular correction, a = a in GF(p^2).
#define sike_fp2correction(a) \
do { \
sike_fpcorrection(a->c0); \
sike_fpcorrection(a->c1); \
} while(0)

// Conversion of a GF(p^2) element to Montgomery representation,
// mc_i = a_i*R^2*R^(-1) = a_i*R in GF(p^2).
#define sike_to_fp2mont(a, mc) \
do { \
sike_fpmul_mont(a->c0, params.mont_R2, mc->c0); \
sike_fpmul_mont(a->c1, params.mont_R2, mc->c1); \
} while(0)

// Conversion of a GF(p^2) element from Montgomery representation to standard representation,
// c_i = ma_i*R^(-1) = a_i in GF(p^2).
#define sike_from_fp2mont(ma, c) \
do { \
sike_from_mont(ma->c0, c->c0); \
sike_from_mont(ma->c1, c->c1); \
} while(0)

#endif // FPX_H_

+ 262
- 0
src/kem/sike/sike-p434-sha256/isogeny.c View File

@@ -0,0 +1,262 @@
/********************************************************************************************
* SIDH: an efficient supersingular isogeny cryptography library
*
* Abstract: elliptic curve and isogeny functions
*********************************************************************************************/
#include <stddef.h>
#include <string.h>
#include "utils.h"
#include "isogeny.h"
#include "fpx.h"

static void xDBL(const point_proj_t P, point_proj_t Q, const f2elm_t A24plus, const f2elm_t C24)
{ // Doubling of a Montgomery point in projective coordinates (X:Z).
// Input: projective Montgomery x-coordinates P = (X1:Z1), where x1=X1/Z1 and Montgomery curve constants A+2C and 4C.
// Output: projective Montgomery x-coordinates Q = 2*P = (X2:Z2).
f2elm_t t0, t1;

sike_fp2sub(P->X, P->Z, t0); // t0 = X1-Z1
sike_fp2add(P->X, P->Z, t1); // t1 = X1+Z1
sike_fp2sqr_mont(t0, t0); // t0 = (X1-Z1)^2
sike_fp2sqr_mont(t1, t1); // t1 = (X1+Z1)^2
sike_fp2mul_mont(C24, t0, Q->Z); // Z2 = C24*(X1-Z1)^2
sike_fp2mul_mont(t1, Q->Z, Q->X); // X2 = C24*(X1-Z1)^2*(X1+Z1)^2
sike_fp2sub(t1, t0, t1); // t1 = (X1+Z1)^2-(X1-Z1)^2
sike_fp2mul_mont(A24plus, t1, t0); // t0 = A24plus*[(X1+Z1)^2-(X1-Z1)^2]
sike_fp2add(Q->Z, t0, Q->Z); // Z2 = A24plus*[(X1+Z1)^2-(X1-Z1)^2] + C24*(X1-Z1)^2
sike_fp2mul_mont(Q->Z, t1, Q->Z); // Z2 = [A24plus*[(X1+Z1)^2-(X1-Z1)^2] + C24*(X1-Z1)^2]*[(X1+Z1)^2-(X1-Z1)^2]
}

void xDBLe(const point_proj_t P, point_proj_t Q, const f2elm_t A24plus, const f2elm_t C24, size_t e)
{ // Computes [2^e](X:Z) on Montgomery curve with projective constant via e repeated doublings.
// Input: projective Montgomery x-coordinates P = (XP:ZP), such that xP=XP/ZP and Montgomery curve constants A+2C and 4C.
// Output: projective Montgomery x-coordinates Q <- (2^e)*P.

memmove(Q, P, sizeof(*P));
for (size_t i = 0; i < e; i++) {
xDBL(Q, Q, A24plus, C24);
}
}

void get_4_isog(const point_proj_t P, f2elm_t A24plus, f2elm_t C24, f2elm_t* coeff)
{ // Computes the corresponding 4-isogeny of a projective Montgomery point (X4:Z4) of order 4.
// Input: projective point of order four P = (X4:Z4).
// Output: the 4-isogenous Montgomery curve with projective coefficients A+2C/4C and the 3 coefficients
// that are used to evaluate the isogeny at a point in eval_4_isog().

sike_fp2sub(P->X, P->Z, coeff[1]); // coeff[1] = X4-Z4
sike_fp2add(P->X, P->Z, coeff[2]); // coeff[2] = X4+Z4
sike_fp2sqr_mont(P->Z, coeff[0]); // coeff[0] = Z4^2
sike_fp2add(coeff[0], coeff[0], coeff[0]); // coeff[0] = 2*Z4^2
sike_fp2sqr_mont(coeff[0], C24); // C24 = 4*Z4^4
sike_fp2add(coeff[0], coeff[0], coeff[0]); // coeff[0] = 4*Z4^2
sike_fp2sqr_mont(P->X, A24plus); // A24plus = X4^2
sike_fp2add(A24plus, A24plus, A24plus); // A24plus = 2*X4^2
sike_fp2sqr_mont(A24plus, A24plus); // A24plus = 4*X4^4
}

void eval_4_isog(point_proj_t P, f2elm_t* coeff)
{ // Evaluates the isogeny at the point (X:Z) in the domain of the isogeny, given a 4-isogeny phi defined
// by the 3 coefficients in coeff (computed in the function get_4_isog()).
// Inputs: the coefficients defining the isogeny, and the projective point P = (X:Z).
// Output: the projective point P = phi(P) = (X:Z) in the codomain.
f2elm_t t0, t1;

sike_fp2add(P->X, P->Z, t0); // t0 = X+Z
sike_fp2sub(P->X, P->Z, t1); // t1 = X-Z
sike_fp2mul_mont(t0, coeff[1], P->X); // X = (X+Z)*coeff[1]
sike_fp2mul_mont(t1, coeff[2], P->Z); // Z = (X-Z)*coeff[2]
sike_fp2mul_mont(t0, t1, t0); // t0 = (X+Z)*(X-Z)
sike_fp2mul_mont(t0, coeff[0], t0); // t0 = coeff[0]*(X+Z)*(X-Z)
sike_fp2add(P->X, P->Z, t1); // t1 = (X-Z)*coeff[2] + (X+Z)*coeff[1]
sike_fp2sub(P->X, P->Z, P->Z); // Z = (X-Z)*coeff[2] - (X+Z)*coeff[1]
sike_fp2sqr_mont(t1, t1); // t1 = [(X-Z)*coeff[2] + (X+Z)*coeff[1]]^2
sike_fp2sqr_mont(P->Z, P->Z); // Z = [(X-Z)*coeff[2] - (X+Z)*coeff[1]]^2
sike_fp2add(t1, t0, P->X); // X = coeff[0]*(X+Z)*(X-Z) + [(X-Z)*coeff[2] + (X+Z)*coeff[1]]^2
sike_fp2sub(P->Z, t0, t0); // t0 = [(X-Z)*coeff[2] - (X+Z)*coeff[1]]^2 - coeff[0]*(X+Z)*(X-Z)
sike_fp2mul_mont(P->X, t1, P->X); // Xfinal
sike_fp2mul_mont(P->Z, t0, P->Z); // Zfinal
}


void xTPL(const point_proj_t P, point_proj_t Q, const f2elm_t A24minus, const f2elm_t A24plus)
{ // Tripling of a Montgomery point in projective coordinates (X:Z).
// Input: projective Montgomery x-coordinates P = (X:Z), where x=X/Z and Montgomery curve constants A24plus = A+2C and A24minus = A-2C.
// Output: projective Montgomery x-coordinates Q = 3*P = (X3:Z3).
f2elm_t t0, t1, t2, t3, t4, t5, t6;

sike_fp2sub(P->X, P->Z, t0); // t0 = X-Z
sike_fp2sqr_mont(t0, t2); // t2 = (X-Z)^2
sike_fp2add(P->X, P->Z, t1); // t1 = X+Z
sike_fp2sqr_mont(t1, t3); // t3 = (X+Z)^2
sike_fp2add(t0, t1, t4); // t4 = 2*X
sike_fp2sub(t1, t0, t0); // t0 = 2*Z
sike_fp2sqr_mont(t4, t1); // t1 = 4*X^2
sike_fp2sub(t1, t3, t1); // t1 = 4*X^2 - (X+Z)^2
sike_fp2sub(t1, t2, t1); // t1 = 4*X^2 - (X+Z)^2 - (X-Z)^2
sike_fp2mul_mont(t3, A24plus, t5); // t5 = A24plus*(X+Z)^2
sike_fp2mul_mont(t3, t5, t3); // t3 = A24plus*(X+Z)^3
sike_fp2mul_mont(A24minus, t2, t6); // t6 = A24minus*(X-Z)^2
sike_fp2mul_mont(t2, t6, t2); // t2 = A24minus*(X-Z)^3
sike_fp2sub(t2, t3, t3); // t3 = A24minus*(X-Z)^3 - coeff*(X+Z)^3
sike_fp2sub(t5, t6, t2); // t2 = A24plus*(X+Z)^2 - A24minus*(X-Z)^2
sike_fp2mul_mont(t1, t2, t1); // t1 = [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2]
sike_fp2add(t3, t1, t2); // t2 = [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2] + A24minus*(X-Z)^3 - coeff*(X+Z)^3
sike_fp2sqr_mont(t2, t2); // t2 = t2^2
sike_fp2mul_mont(t4, t2, Q->X); // X3 = 2*X*t2
sike_fp2sub(t3, t1, t1); // t1 = A24minus*(X-Z)^3 - A24plus*(X+Z)^3 - [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2]
sike_fp2sqr_mont(t1, t1); // t1 = t1^2
sike_fp2mul_mont(t0, t1, Q->Z); // Z3 = 2*Z*t1
}

void xTPLe(const point_proj_t P, point_proj_t Q, const f2elm_t A24minus, const f2elm_t A24plus, size_t e)
{ // Computes [3^e](X:Z) on Montgomery curve with projective constant via e repeated triplings.
// Input: projective Montgomery x-coordinates P = (XP:ZP), such that xP=XP/ZP and Montgomery curve constants A24plus = A+2C and A24minus = A-2C.
// Output: projective Montgomery x-coordinates Q <- (3^e)*P.
memmove(Q, P, sizeof(*P));
for (size_t i = 0; i < e; i++) {
xTPL(Q, Q, A24minus, A24plus);
}
}

void get_3_isog(const point_proj_t P, f2elm_t A24minus, f2elm_t A24plus, f2elm_t* coeff)
{ // Computes the corresponding 3-isogeny of a projective Montgomery point (X3:Z3) of order 3.
// Input: projective point of order three P = (X3:Z3).
// Output: the 3-isogenous Montgomery curve with projective coefficient A/C.
f2elm_t t0, t1, t2, t3, t4;

sike_fp2sub(P->X, P->Z, coeff[0]); // coeff0 = X-Z
sike_fp2sqr_mont(coeff[0], t0); // t0 = (X-Z)^2
sike_fp2add(P->X, P->Z, coeff[1]); // coeff1 = X+Z
sike_fp2sqr_mont(coeff[1], t1); // t1 = (X+Z)^2
sike_fp2add(t0, t1, t2); // t2 = (X+Z)^2 + (X-Z)^2
sike_fp2add(coeff[0], coeff[1], t3); // t3 = 2*X
sike_fp2sqr_mont(t3, t3); // t3 = 4*X^2
sike_fp2sub(t3, t2, t3); // t3 = 4*X^2 - (X+Z)^2 - (X-Z)^2
sike_fp2add(t1, t3, t2); // t2 = 4*X^2 - (X-Z)^2
sike_fp2add(t3, t0, t3); // t3 = 4*X^2 - (X+Z)^2
sike_fp2add(t0, t3, t4); // t4 = 4*X^2 - (X+Z)^2 + (X-Z)^2
sike_fp2add(t4, t4, t4); // t4 = 2(4*X^2 - (X+Z)^2 + (X-Z)^2)
sike_fp2add(t1, t4, t4); // t4 = 8*X^2 - (X+Z)^2 + 2*(X-Z)^2
sike_fp2mul_mont(t2, t4, A24minus); // A24minus = [4*X^2 - (X-Z)^2]*[8*X^2 - (X+Z)^2 + 2*(X-Z)^2]
sike_fp2add(t1, t2, t4); // t4 = 4*X^2 + (X+Z)^2 - (X-Z)^2
sike_fp2add(t4, t4, t4); // t4 = 2(4*X^2 + (X+Z)^2 - (X-Z)^2)
sike_fp2add(t0, t4, t4); // t4 = 8*X^2 + 2*(X+Z)^2 - (X-Z)^2
sike_fp2mul_mont(t3, t4, t4); // t4 = [4*X^2 - (X+Z)^2]*[8*X^2 + 2*(X+Z)^2 - (X-Z)^2]
sike_fp2sub(t4, A24minus, t0); // t0 = [4*X^2 - (X+Z)^2]*[8*X^2 + 2*(X+Z)^2 - (X-Z)^2] - [4*X^2 - (X-Z)^2]*[8*X^2 - (X+Z)^2 + 2*(X-Z)^2]
sike_fp2add(A24minus, t0, A24plus); // A24plus = 8*X^2 - (X+Z)^2 + 2*(X-Z)^2
}


void eval_3_isog(point_proj_t Q, f2elm_t* coeff)
{ // Computes the 3-isogeny R=phi(X:Z), given projective point (X3:Z3) of order 3 on a Montgomery curve and
// a point P with 2 coefficients in coeff (computed in the function get_3_isog()).
// Inputs: projective points P = (X3:Z3) and Q = (X:Z).
// Output: the projective point Q <- phi(Q) = (X3:Z3).
f2elm_t t0, t1, t2;

sike_fp2add(Q->X, Q->Z, t0); // t0 = X+Z
sike_fp2sub(Q->X, Q->Z, t1); // t1 = X-Z
sike_fp2mul_mont(t0, coeff[0], t0); // t0 = coeff0*(X+Z)
sike_fp2mul_mont(t1, coeff[1], t1); // t1 = coeff1*(X-Z)
sike_fp2add(t0, t1, t2); // t2 = coeff0*(X+Z) + coeff1*(X-Z)
sike_fp2sub(t1, t0, t0); // t0 = coeff1*(X-Z) - coeff0*(X+Z)
sike_fp2sqr_mont(t2, t2); // t2 = [coeff0*(X+Z) + coeff1*(X-Z)]^2
sike_fp2sqr_mont(t0, t0); // t0 = [coeff1*(X-Z) - coeff0*(X+Z)]^2
sike_fp2mul_mont(Q->X, t2, Q->X); // X3final = X*[coeff0*(X+Z) + coeff1*(X-Z)]^2
sike_fp2mul_mont(Q->Z, t0, Q->Z); // Z3final = Z*[coeff1*(X-Z) - coeff0*(X+Z)]^2
}


void inv_3_way(f2elm_t z1, f2elm_t z2, f2elm_t z3)
{ // 3-way simultaneous inversion
// Input: z1,z2,z3
// Output: 1/z1,1/z2,1/z3 (override inputs).
f2elm_t t0, t1, t2, t3;

sike_fp2mul_mont(z1, z2, t0); // t0 = z1*z2
sike_fp2mul_mont(z3, t0, t1); // t1 = z1*z2*z3
sike_fp2inv_mont(t1); // t1 = 1/(z1*z2*z3)
sike_fp2mul_mont(z3, t1, t2); // t2 = 1/(z1*z2)
sike_fp2mul_mont(t2, z2, t3); // t3 = 1/z1
sike_fp2mul_mont(t2, z1, z2); // z2 = 1/z2
sike_fp2mul_mont(t0, t1, z3); // z3 = 1/z3
sike_fp2copy(t3, z1); // z1 = 1/z1
}


void get_A(const f2elm_t xP, const f2elm_t xQ, const f2elm_t xR, f2elm_t A)
{ // Given the x-coordinates of P, Q, and R, returns the value A corresponding to the Montgomery curve E_A: y^2=x^3+A*x^2+x such that R=Q-P on E_A.
// Input: the x-coordinates xP, xQ, and xR of the points P, Q and R.
// Output: the coefficient A corresponding to the curve E_A: y^2=x^3+A*x^2+x.
f2elm_t t0, t1, one = F2ELM_INIT;

extern const struct params_t params;
sike_fpcopy(params.mont_one, one->c0);
sike_fp2add(xP, xQ, t1); // t1 = xP+xQ
sike_fp2mul_mont(xP, xQ, t0); // t0 = xP*xQ
sike_fp2mul_mont(xR, t1, A); // A = xR*t1
sike_fp2add(t0, A, A); // A = A+t0
sike_fp2mul_mont(t0, xR, t0); // t0 = t0*xR
sike_fp2sub(A, one, A); // A = A-1
sike_fp2add(t0, t0, t0); // t0 = t0+t0
sike_fp2add(t1, xR, t1); // t1 = t1+xR
sike_fp2add(t0, t0, t0); // t0 = t0+t0
sike_fp2sqr_mont(A, A); // A = A^2
sike_fp2inv_mont(t0); // t0 = 1/t0
sike_fp2mul_mont(A, t0, A); // A = A*t0
sike_fp2sub(A, t1, A); // Afinal = A-t1
}


void j_inv(const f2elm_t A, const f2elm_t C, f2elm_t jinv)
{ // Computes the j-invariant of a Montgomery curve with projective constant.
// Input: A,C in GF(p^2).
// Output: j=256*(A^2-3*C^2)^3/(C^4*(A^2-4*C^2)), which is the j-invariant of the Montgomery curve B*y^2=x^3+(A/C)*x^2+x or (equivalently) j-invariant of B'*y^2=C*x^3+A*x^2+C*x.
f2elm_t t0, t1;

sike_fp2sqr_mont(A, jinv); // jinv = A^2
sike_fp2sqr_mont(C, t1); // t1 = C^2
sike_fp2add(t1, t1, t0); // t0 = t1+t1
sike_fp2sub(jinv, t0, t0); // t0 = jinv-t0
sike_fp2sub(t0, t1, t0); // t0 = t0-t1
sike_fp2sub(t0, t1, jinv); // jinv = t0-t1
sike_fp2sqr_mont(t1, t1); // t1 = t1^2
sike_fp2mul_mont(jinv, t1, jinv); // jinv = jinv*t1
sike_fp2add(t0, t0, t0); // t0 = t0+t0
sike_fp2add(t0, t0, t0); // t0 = t0+t0
sike_fp2sqr_mont(t0, t1); // t1 = t0^2
sike_fp2mul_mont(t0, t1, t0); // t0 = t0*t1
sike_fp2add(t0, t0, t0); // t0 = t0+t0
sike_fp2add(t0, t0, t0); // t0 = t0+t0
sike_fp2inv_mont(jinv); // jinv = 1/jinv
sike_fp2mul_mont(jinv, t0, jinv); // jinv = t0*jinv
}


void xDBLADD(point_proj_t P, point_proj_t Q, const f2elm_t xPQ, const f2elm_t A24)
{ // Simultaneous doubling and differential addition.
// Input: projective Montgomery points P=(XP:ZP) and Q=(XQ:ZQ) such that xP=XP/ZP and xQ=XQ/ZQ, affine difference xPQ=x(P-Q) and Montgomery curve constant A24=(A+2)/4.
// Output: projective Montgomery points P <- 2*P = (X2P:Z2P) such that x(2P)=X2P/Z2P, and Q <- P+Q = (XQP:ZQP) such that = x(Q+P)=XQP/ZQP.
f2elm_t t0, t1, t2;

sike_fp2add(P->X, P->Z, t0); // t0 = XP+ZP
sike_fp2sub(P->X, P->Z, t1); // t1 = XP-ZP
sike_fp2sqr_mont(t0, P->X); // XP = (XP+ZP)^2
sike_fp2sub(Q->X, Q->Z, t2); // t2 = XQ-ZQ
sike_fp2correction(t2);
sike_fp2add(Q->X, Q->Z, Q->X); // XQ = XQ+ZQ
sike_fp2mul_mont(t0, t2, t0); // t0 = (XP+ZP)*(XQ-ZQ)
sike_fp2sqr_mont(t1, P->Z); // ZP = (XP-ZP)^2
sike_fp2mul_mont(t1, Q->X, t1); // t1 = (XP-ZP)*(XQ+ZQ)
sike_fp2sub(P->X, P->Z, t2); // t2 = (XP+ZP)^2-(XP-ZP)^2
sike_fp2mul_mont(P->X, P->Z, P->X); // XP = (XP+ZP)^2*(XP-ZP)^2
sike_fp2mul_mont(t2, A24, Q->X); // XQ = A24*[(XP+ZP)^2-(XP-ZP)^2]
sike_fp2sub(t0, t1, Q->Z); // ZQ = (XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)
sike_fp2add(Q->X, P->Z, P->Z); // ZP = A24*[(XP+ZP)^2-(XP-ZP)^2]+(XP-ZP)^2
sike_fp2add(t0, t1, Q->X); // XQ = (XP+ZP)*(XQ-ZQ)+(XP-ZP)*(XQ+ZQ)
sike_fp2mul_mont(P->Z, t2, P->Z); // ZP = [A24*[(XP+ZP)^2-(XP-ZP)^2]+(XP-ZP)^2]*[(XP+ZP)^2-(XP-ZP)^2]
sike_fp2sqr_mont(Q->Z, Q->Z); // ZQ = [(XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)]^2
sike_fp2sqr_mont(Q->X, Q->X); // XQ = [(XP+ZP)*(XQ-ZQ)+(XP-ZP)*(XQ+ZQ)]^2
sike_fp2mul_mont(Q->Z, xPQ, Q->Z); // ZQ = xPQ*[(XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)]^2
}

+ 49
- 0
src/kem/sike/sike-p434-sha256/isogeny.h View File

@@ -0,0 +1,49 @@
#ifndef ISOGENY_H_
#define ISOGENY_H_

// Computes [2^e](X:Z) on Montgomery curve with projective
// constant via e repeated doublings.
void xDBLe(
const point_proj_t P, point_proj_t Q, const f2elm_t A24plus,
const f2elm_t C24, size_t e);
// Simultaneous doubling and differential addition.
void xDBLADD(
point_proj_t P, point_proj_t Q, const f2elm_t xPQ,
const f2elm_t A24);
// Tripling of a Montgomery point in projective coordinates (X:Z).
void xTPL(
const point_proj_t P, point_proj_t Q, const f2elm_t A24minus,
const f2elm_t A24plus);
// Computes [3^e](X:Z) on Montgomery curve with projective constant
// via e repeated triplings.
void xTPLe(
const point_proj_t P, point_proj_t Q, const f2elm_t A24minus,
const f2elm_t A24plus, size_t e);
// Given the x-coordinates of P, Q, and R, returns the value A
// corresponding to the Montgomery curve E_A: y^2=x^3+A*x^2+x such that R=Q-P on E_A.
void get_A(
const f2elm_t xP, const f2elm_t xQ, const f2elm_t xR, f2elm_t A);
// Computes the j-invariant of a Montgomery curve with projective constant.
void j_inv(
const f2elm_t A, const f2elm_t C, f2elm_t jinv);
// Computes the corresponding 4-isogeny of a projective Montgomery
// point (X4:Z4) of order 4.
void get_4_isog(
const point_proj_t P, f2elm_t A24plus, f2elm_t C24, f2elm_t* coeff);
// Computes the corresponding 3-isogeny of a projective Montgomery
// point (X3:Z3) of order 3.
void get_3_isog(
const point_proj_t P, f2elm_t A24minus, f2elm_t A24plus,
f2elm_t* coeff);
// Computes the 3-isogeny R=phi(X:Z), given projective point (X3:Z3)
// of order 3 on a Montgomery curve and a point P with coefficients given in coeff.
void eval_3_isog(
point_proj_t Q, f2elm_t* coeff);
// Evaluates the isogeny at the point (X:Z) in the domain of the isogeny.
void eval_4_isog(
point_proj_t P, f2elm_t* coeff);
// 3-way simultaneous inversion
void inv_3_way(
f2elm_t z1, f2elm_t z2, f2elm_t z3);

#endif // ISOGENY_H_

+ 128
- 0
src/kem/sike/sike-p434-sha256/params.c View File

@@ -0,0 +1,128 @@
/********************************************************************************************
* SIDH: an efficient supersingular isogeny cryptography library
*
* Abstract: supersingular isogeny parameters and generation of functions for P434
*********************************************************************************************/

#include "utils.h"

// Parameters for isogeny system "SIKE"
const struct params_t params = {
.prime = {
U64_TO_WORDS(0xFFFFFFFFFFFFFFFF), U64_TO_WORDS(0xFFFFFFFFFFFFFFFF),
U64_TO_WORDS(0xFFFFFFFFFFFFFFFF), U64_TO_WORDS(0xFDC1767AE2FFFFFF),
U64_TO_WORDS(0x7BC65C783158AEA3), U64_TO_WORDS(0x6CFC5FD681C52056),
U64_TO_WORDS(0x0002341F27177344)
},
.prime_p1 = {
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0xFDC1767AE3000000),
U64_TO_WORDS(0x7BC65C783158AEA3), U64_TO_WORDS(0x6CFC5FD681C52056),
U64_TO_WORDS(0x0002341F27177344)
},
.prime_x2 = {
U64_TO_WORDS(0xFFFFFFFFFFFFFFFE), U64_TO_WORDS(0xFFFFFFFFFFFFFFFF),
U64_TO_WORDS(0xFFFFFFFFFFFFFFFF), U64_TO_WORDS(0xFB82ECF5C5FFFFFF),
U64_TO_WORDS(0xF78CB8F062B15D47), U64_TO_WORDS(0xD9F8BFAD038A40AC),
U64_TO_WORDS(0x0004683E4E2EE688)
},
.A_gen = {
U64_TO_WORDS(0x05ADF455C5C345BF), U64_TO_WORDS(0x91935C5CC767AC2B),
U64_TO_WORDS(0xAFE4E879951F0257), U64_TO_WORDS(0x70E792DC89FA27B1),
U64_TO_WORDS(0xF797F526BB48C8CD), U64_TO_WORDS(0x2181DB6131AF621F),
U64_TO_WORDS(0x00000A1C08B1ECC4), // XPA0
U64_TO_WORDS(0x74840EB87CDA7788), U64_TO_WORDS(0x2971AA0ECF9F9D0B),
U64_TO_WORDS(0xCB5732BDF41715D5), U64_TO_WORDS(0x8CD8E51F7AACFFAA),
U64_TO_WORDS(0xA7F424730D7E419F), U64_TO_WORDS(0xD671EB919A179E8C),
U64_TO_WORDS(0x0000FFA26C5A924A), // XPA1
U64_TO_WORDS(0xFEC6E64588B7273B), U64_TO_WORDS(0xD2A626D74CBBF1C6),
U64_TO_WORDS(0xF8F58F07A78098C7), U64_TO_WORDS(0xE23941F470841B03),
U64_TO_WORDS(0x1B63EDA2045538DD), U64_TO_WORDS(0x735CFEB0FFD49215),
U64_TO_WORDS(0x0001C4CB77542876), // XQA0
U64_TO_WORDS(0xADB0F733C17FFDD6), U64_TO_WORDS(0x6AFFBD037DA0A050),
U64_TO_WORDS(0x680EC43DB144E02F), U64_TO_WORDS(0x1E2E5D5FF524E374),
U64_TO_WORDS(0xE2DDA115260E2995), U64_TO_WORDS(0xA6E4B552E2EDE508),
U64_TO_WORDS(0x00018ECCDDF4B53E), // XQA1
U64_TO_WORDS(0x01BA4DB518CD6C7D), U64_TO_WORDS(0x2CB0251FE3CC0611),
U64_TO_WORDS(0x259B0C6949A9121B), U64_TO_WORDS(0x60E17AC16D2F82AD),
U64_TO_WORDS(0x3AA41F1CE175D92D), U64_TO_WORDS(0x413FBE6A9B9BC4F3),
U64_TO_WORDS(0x00022A81D8D55643), // XRA0
U64_TO_WORDS(0xB8ADBC70FC82E54A), U64_TO_WORDS(0xEF9CDDB0D5FADDED),
U64_TO_WORDS(0x5820C734C80096A0), U64_TO_WORDS(0x7799994BAA96E0E4),
U64_TO_WORDS(0x044961599E379AF8), U64_TO_WORDS(0xDB2B94FBF09F27E2),
U64_TO_WORDS(0x0000B87FC716C0C6) // XRA1
},
.B_gen = {
U64_TO_WORDS(0x6E5497556EDD48A3), U64_TO_WORDS(0x2A61B501546F1C05),
U64_TO_WORDS(0xEB919446D049887D), U64_TO_WORDS(0x5864A4A69D450C4F),
U64_TO_WORDS(0xB883F276A6490D2B), U64_TO_WORDS(0x22CC287022D5F5B9),
U64_TO_WORDS(0x0001BED4772E551F), // XPB0
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), // XPB1
U64_TO_WORDS(0xFAE2A3F93D8B6B8E), U64_TO_WORDS(0x494871F51700FE1C),
U64_TO_WORDS(0xEF1A94228413C27C), U64_TO_WORDS(0x498FF4A4AF60BD62),
U64_TO_WORDS(0xB00AD2A708267E8A), U64_TO_WORDS(0xF4328294E017837F),
U64_TO_WORDS(0x000034080181D8AE), // XQB0
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), // XQB1
U64_TO_WORDS(0x283B34FAFEFDC8E4), U64_TO_WORDS(0x9208F44977C3E647),
U64_TO_WORDS(0x7DEAE962816F4E9A), U64_TO_WORDS(0x68A2BA8AA262EC9D),
U64_TO_WORDS(0x8176F112EA43F45B), U64_TO_WORDS(0x02106D022634F504),
U64_TO_WORDS(0x00007E8A50F02E37), // XRB0
U64_TO_WORDS(0xB378B7C1DA22CCB1), U64_TO_WORDS(0x6D089C99AD1D9230),
U64_TO_WORDS(0xEBE15711813E2369), U64_TO_WORDS(0x2B35A68239D48A53),
U64_TO_WORDS(0x445F6FD138407C93), U64_TO_WORDS(0xBEF93B29A3F6B54B),
U64_TO_WORDS(0x000173FA910377D3) // XRB1
},
.mont_R2 = {
U64_TO_WORDS(0x28E55B65DCD69B30), U64_TO_WORDS(0xACEC7367768798C2),
U64_TO_WORDS(0xAB27973F8311688D), U64_TO_WORDS(0x175CC6AF8D6C7C0B),
U64_TO_WORDS(0xABCD92BF2DDE347E), U64_TO_WORDS(0x69E16A61C7686D9A),
U64_TO_WORDS(0x000025A89BCDD12A)
},
.mont_one = {
U64_TO_WORDS(0x000000000000742C), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0xB90FF404FC000000),
U64_TO_WORDS(0xD801A4FB559FACD4), U64_TO_WORDS(0xE93254545F77410C),
U64_TO_WORDS(0x0000ECEEA7BD2EDA)
},
.mont_six = {
U64_TO_WORDS(0x000000000002B90A), U64_TO_WORDS(0x0000000000000000),
U64_TO_WORDS(0x0000000000000000), U64_TO_WORDS(0x5ADCCB2822000000),
U64_TO_WORDS(0x187D24F39F0CAFB4), U64_TO_WORDS(0x9D353A4D394145A0),
U64_TO_WORDS(0x00012559A0403298)
},
.A_strat = {
0x30, 0x1C, 0x10, 0x08, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01,
0x01, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x08, 0x04,
0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x04, 0x02, 0x01, 0x01,
0x02, 0x01, 0x01, 0x0D, 0x07, 0x04, 0x02, 0x01, 0x01, 0x02,
0x01, 0x01, 0x03, 0x02, 0x01, 0x01, 0x01, 0x01, 0x05, 0x04,
0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x01,
0x15, 0x0C, 0x07, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01,
0x03, 0x02, 0x01, 0x01, 0x01, 0x01, 0x05, 0x03, 0x02, 0x01,
0x01, 0x01, 0x01, 0x02, 0x01, 0x01, 0x01, 0x09, 0x05, 0x03,
0x02, 0x01, 0x01, 0x01, 0x01, 0x02, 0x01, 0x01, 0x01, 0x04,
0x02, 0x01, 0x01, 0x01, 0x02, 0x01, 0x01
},
.B_strat = {
0x42, 0x21, 0x11, 0x09, 0x05, 0x03, 0x02, 0x01, 0x01, 0x01,
0x01, 0x02, 0x01, 0x01, 0x01, 0x04, 0x02, 0x01, 0x01, 0x01,
0x02, 0x01, 0x01, 0x08, 0x04, 0x02, 0x01, 0x01, 0x01, 0x02,
0x01, 0x01, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x10,
0x08, 0x04, 0x02, 0x01, 0x01, 0x01, 0x02, 0x01, 0x01, 0x04,
0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x08, 0x04, 0x02, 0x01,
0x01, 0x02, 0x01, 0x01, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01,
0x01, 0x20, 0x10, 0x08, 0x04, 0x03, 0x01, 0x01, 0x01, 0x01,
0x02, 0x01, 0x01, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01,
0x08, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x04, 0x02,
0x01, 0x01, 0x02, 0x01, 0x01, 0x10, 0x08, 0x04, 0x02, 0x01,
0x01, 0x02, 0x01, 0x01, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01,
0x01, 0x08, 0x04, 0x02, 0x01, 0x01, 0x02, 0x01, 0x01, 0x04,
0x02, 0x01, 0x01, 0x02, 0x01, 0x01
}
};

+ 517
- 0
src/kem/sike/sike-p434-sha256/sike.c View File

@@ -0,0 +1,517 @@
/********************************************************************************************
* SIDH: an efficient supersingular isogeny cryptography library
*
* Abstract: supersingular isogeny key encapsulation (SIKE) protocol
*********************************************************************************************/

#include <assert.h>
#include <stddef.h>
#include <stdint.h>
#include <string.h>
#include <sha2/sha256.h>
#include <random/randombytes.h>

#include "utils.h"
#include "isogeny.h"
#include "fpx.h"

extern const struct params_t params;

// SIDH_JINV_BYTESZ is a number of bytes used for encoding j-invariant.
#define SIDH_JINV_BYTESZ 110U
// SIDH_PRV_A_BITSZ is a number of bits of SIDH private key (2-isogeny)
#define SIDH_PRV_A_BITSZ 216U
// SIDH_PRV_A_BITSZ is a number of bits of SIDH private key (3-isogeny)
#define SIDH_PRV_B_BITSZ 217U
// MAX_INT_POINTS_ALICE is a number of points used in 2-isogeny tree computation
#define MAX_INT_POINTS_ALICE 7U
// MAX_INT_POINTS_ALICE is a number of points used in 3-isogeny tree computation
#define MAX_INT_POINTS_BOB 8U

// Swap points.
// If option = 0 then P <- P and Q <- Q, else if option = 0xFF...FF then P <- Q and Q <- P
#if !defined(ARCH_X86_64) || defined(ARCH_GENERIC)
static void sike_cswap(point_proj_t P, point_proj_t Q, const crypto_word_t option)
{
crypto_word_t temp;
for (size_t i = 0; i < NWORDS_FIELD; i++) {
temp = option & (P->X->c0[i] ^ Q->X->c0[i]);
P->X->c0[i] = temp ^ P->X->c0[i];
Q->X->c0[i] = temp ^ Q->X->c0[i];
temp = option & (P->Z->c0[i] ^ Q->Z->c0[i]);
P->Z->c0[i] = temp ^ P->Z->c0[i];
Q->Z->c0[i] = temp ^ Q->Z->c0[i];
temp = option & (P->X->c1[i] ^ Q->X->c1[i]);
P->X->c1[i] = temp ^ P->X->c1[i];
Q->X->c1[i] = temp ^ Q->X->c1[i];
temp = option & (P->Z->c1[i] ^ Q->Z->c1[i]);
P->Z->c1[i] = temp ^ P->Z->c1[i];
Q->Z->c1[i] = temp ^ Q->Z->c1[i];
}
}
#endif

// Swap points.
// If option = 0 then P <- P and Q <- Q, else if option = 0xFF...FF then P <- Q and Q <- P
static inline void sike_fp2cswap(point_proj_t P, point_proj_t Q, const crypto_word_t option)
{
#if defined(ARCH_X86_64) && !defined(ARCH_GENERIC)
sike_cswap_asm(P, Q, option);
#else
sike_cswap(P, Q, option);
#endif
}

static void ladder3Pt(
const f2elm_t xP, const f2elm_t xQ, const f2elm_t xPQ, const uint8_t* m,
int is_A, point_proj_t R, const f2elm_t A) {
point_proj_t R0 = POINT_PROJ_INIT, R2 = POINT_PROJ_INIT;
f2elm_t A24 = F2ELM_INIT;
crypto_word_t mask;
int bit, swap, prevbit = 0;

const size_t nbits = is_A?SIDH_PRV_A_BITSZ:SIDH_PRV_B_BITSZ;

// Initializing constant
sike_fpcopy(params.mont_one, A24[0].c0);
sike_fp2add(A24, A24, A24);
sike_fp2add(A, A24, A24);
sike_fp2div2(A24, A24);
sike_fp2div2(A24, A24); // A24 = (A+2)/4

// Initializing points
sike_fp2copy(xQ, R0->X);
sike_fpcopy(params.mont_one, R0->Z[0].c0);
sike_fp2copy(xPQ, R2->X);
sike_fpcopy(params.mont_one, R2->Z[0].c0);
sike_fp2copy(xP, R->X);
sike_fpcopy(params.mont_one, R->Z[0].c0);
memset(R->Z->c1, 0, sizeof(R->Z->c1));

// Main loop
for (size_t i = 0; i < nbits; i++) {
bit = (m[i >> 3] >> (i & 7)) & 1;
swap = bit ^ prevbit;
prevbit = bit;
mask = 0 - (crypto_word_t)swap;

sike_fp2cswap(R, R2, mask);
xDBLADD(R0, R2, R->X, A24);
sike_fp2mul_mont(R2->X, R->Z, R2->X);
}
swap = 0 ^ prevbit;
mask = 0 - (crypto_word_t)swap;
sike_fp2cswap(R, R2, mask);
}

// Initialization of basis points
static inline void sike_init_basis(const crypto_word_t *gen, f2elm_t XP, f2elm_t XQ, f2elm_t XR) {
sike_fpcopy(gen, XP->c0);
sike_fpcopy(gen + NWORDS_FIELD, XP->c1);
sike_fpcopy(gen + 2*NWORDS_FIELD, XQ->c0);
sike_fpcopy(gen + 3*NWORDS_FIELD, XQ->c1);
sike_fpcopy(gen + 4*NWORDS_FIELD, XR->c0);
sike_fpcopy(gen + 5*NWORDS_FIELD, XR->c1);
}

// Conversion of GF(p^2) element from Montgomery to standard representation.
static inline void sike_fp2_encode(const f2elm_t x, uint8_t *enc) {
f2elm_t t;
sike_from_fp2mont(x, t);

// convert to bytes in little endian form
for (size_t i=0; i<FIELD_BYTESZ; i++) {
enc[i+ 0] = (t[0].c0[i/LSZ] >> (8*(i%LSZ))) & 0xFF;
enc[i+FIELD_BYTESZ] = (t[0].c1[i/LSZ] >> (8*(i%LSZ))) & 0xFF;
}
}

// Parse byte sequence back into GF(p^2) element, and conversion to Montgomery representation.
// Elements over GF(p503) are encoded in 63 octets in little endian format
// (i.e., the least significant octet is located in the lowest memory address).
static inline void fp2_decode(const uint8_t *enc, f2elm_t t) {
memset(t[0].c0, 0, sizeof(t[0].c0));
memset(t[0].c1, 0, sizeof(t[0].c1));
// convert bytes in little endian form to f2elm_t
for (size_t i = 0; i < FIELD_BYTESZ; i++) {
t[0].c0[i/LSZ] |= ((crypto_word_t)enc[i+ 0]) << (8*(i%LSZ));
t[0].c1[i/LSZ] |= ((crypto_word_t)enc[i+FIELD_BYTESZ]) << (8*(i%LSZ));
}
sike_to_fp2mont(t, t);
}

// Alice's ephemeral public key generation
// Input: a private key prA in the range [0, 2^250 - 1], stored in 32 bytes.
// Output: the public key pkA consisting of 3 GF(p503^2) elements encoded in 378 bytes.
static void gen_iso_A(const uint8_t* skA, uint8_t* pkA)
{
point_proj_t R, pts[MAX_INT_POINTS_ALICE];
point_proj_t phiP = POINT_PROJ_INIT;
point_proj_t phiQ = POINT_PROJ_INIT;
point_proj_t phiR = POINT_PROJ_INIT;
f2elm_t XPA, XQA, XRA, coeff[3];
f2elm_t A24plus = F2ELM_INIT;
f2elm_t C24 = F2ELM_INIT;
f2elm_t A = F2ELM_INIT;
unsigned int m, index = 0, pts_index[MAX_INT_POINTS_ALICE], npts = 0, ii = 0;

// Initialize basis points
sike_init_basis(params.A_gen, XPA, XQA, XRA);
sike_init_basis(params.B_gen, phiP->X, phiQ->X, phiR->X);
sike_fpcopy(params.mont_one, (phiP->Z)->c0);
sike_fpcopy(params.mont_one, (phiQ->Z)->c0);
sike_fpcopy(params.mont_one, (phiR->Z)->c0);

// Initialize constants: A24plus = A+2C, C24 = 4C, where A=6, C=1
sike_fpcopy(params.mont_one, A24plus->c0);
sike_fp2add(A24plus, A24plus, A24plus);
sike_fp2add(A24plus, A24plus, C24);
sike_fp2add(A24plus, C24, A);
sike_fp2add(C24, C24, A24plus);

// Retrieve kernel point
ladder3Pt(XPA, XQA, XRA, skA, 1, R, A);

// Traverse tree
index = 0;
for (size_t row = 1; row < A_max; row++) {
while (index < A_max-row) {
sike_fp2copy(R->X, pts[npts]->X);
sike_fp2copy(R->Z, pts[npts]->Z);
pts_index[npts++] = index;
m = params.A_strat[ii++];
xDBLe(R, R, A24plus, C24, (2*m));
index += m;
}
get_4_isog(R, A24plus, C24, coeff);

for (size_t i = 0; i < npts; i++) {
eval_4_isog(pts[i], coeff);
}
eval_4_isog(phiP, coeff);
eval_4_isog(phiQ, coeff);
eval_4_isog(phiR, coeff);

sike_fp2copy(pts[npts-1]->X, R->X);
sike_fp2copy(pts[npts-1]->Z, R->Z);
index = pts_index[npts-1];
npts -= 1;
}

get_4_isog(R, A24plus, C24, coeff);
eval_4_isog(phiP, coeff);
eval_4_isog(phiQ, coeff);
eval_4_isog(phiR, coeff);

inv_3_way(phiP->Z, phiQ->Z, phiR->Z);
sike_fp2mul_mont(phiP->X, phiP->Z, phiP->X);
sike_fp2mul_mont(phiQ->X, phiQ->Z, phiQ->X);
sike_fp2mul_mont(phiR->X, phiR->Z, phiR->X);

// Format public key
sike_fp2_encode(phiP->X, pkA);
sike_fp2_encode(phiQ->X, pkA + SIDH_JINV_BYTESZ);
sike_fp2_encode(phiR->X, pkA + 2*SIDH_JINV_BYTESZ);
}

// Bob's ephemeral key-pair generation
// It produces a private key skB and computes the public key pkB.
// The private key is an integer in the range [0, 2^Floor(Log(2,3^159)) - 1], stored in 32 bytes.
// The public key consists of 3 GF(p503^2) elements encoded in 378 bytes.
static void gen_iso_B(const uint8_t* skB, uint8_t* pkB)
{
point_proj_t R, pts[MAX_INT_POINTS_BOB];
point_proj_t phiP = POINT_PROJ_INIT;
point_proj_t phiQ = POINT_PROJ_INIT;
point_proj_t phiR = POINT_PROJ_INIT;
f2elm_t XPB, XQB, XRB, coeff[3];
f2elm_t A24plus = F2ELM_INIT;
f2elm_t A24minus = F2ELM_INIT;
f2elm_t A = F2ELM_INIT;
unsigned int m, index = 0, pts_index[MAX_INT_POINTS_BOB], npts = 0, ii = 0;

// Initialize basis points
sike_init_basis(params.B_gen, XPB, XQB, XRB);
sike_init_basis(params.A_gen, phiP->X, phiQ->X, phiR->X);
sike_fpcopy(params.mont_one, (phiP->Z)->c0);
sike_fpcopy(params.mont_one, (phiQ->Z)->c0);
sike_fpcopy(params.mont_one, (phiR->Z)->c0);

// Initialize constants: A24minus = A-2C, A24plus = A+2C, where A=6, C=1
sike_fpcopy(params.mont_one, A24plus->c0);
sike_fp2add(A24plus, A24plus, A24plus);
sike_fp2add(A24plus, A24plus, A24minus);
sike_fp2add(A24plus, A24minus, A);
sike_fp2add(A24minus, A24minus, A24plus);

// Retrieve kernel point
ladder3Pt(XPB, XQB, XRB, skB, 0, R, A);

// Traverse tree
index = 0;
for (size_t row = 1; row < B_max; row++) {
while (index < B_max-row) {
sike_fp2copy(R->X, pts[npts]->X);
sike_fp2copy(R->Z, pts[npts]->Z);
pts_index[npts++] = index;
m = params.B_strat[ii++];
xTPLe(R, R, A24minus, A24plus, m);
index += m;
}
get_3_isog(R, A24minus, A24plus, coeff);

for (size_t i = 0; i < npts; i++) {
eval_3_isog(pts[i], coeff);
}
eval_3_isog(phiP, coeff);
eval_3_isog(phiQ, coeff);
eval_3_isog(phiR, coeff);

sike_fp2copy(pts[npts-1]->X, R->X);
sike_fp2copy(pts[npts-1]->Z, R->Z);
index = pts_index[npts-1];
npts -= 1;
}

get_3_isog(R, A24minus, A24plus, coeff);
eval_3_isog(phiP, coeff);
eval_3_isog(phiQ, coeff);
eval_3_isog(phiR, coeff);

inv_3_way(phiP->Z, phiQ->Z, phiR->Z);
sike_fp2mul_mont(phiP->X, phiP->Z, phiP->X);
sike_fp2mul_mont(phiQ->X, phiQ->Z, phiQ->X);
sike_fp2mul_mont(phiR->X, phiR->Z, phiR->X);

// Format public key
sike_fp2_encode(phiP->X, pkB);
sike_fp2_encode(phiQ->X, pkB + SIDH_JINV_BYTESZ);
sike_fp2_encode(phiR->X, pkB + 2*SIDH_JINV_BYTESZ);
}

// Alice's ephemeral shared secret computation
// It produces a shared secret key ssA using her secret key skA and Bob's public key pkB
// Inputs: Alice's skA is an integer in the range [0, 2^250 - 1], stored in 32 bytes.
// Bob's pkB consists of 3 GF(p503^2) elements encoded in 378 bytes.
// Output: a shared secret ssA that consists of one element in GF(p503^2) encoded in 126 bytes.
static void ex_iso_A(const uint8_t* skA, const uint8_t* pkB, uint8_t* ssA)
{
point_proj_t R, pts[MAX_INT_POINTS_ALICE];
f2elm_t coeff[3], PKB[3], jinv;
f2elm_t A24plus = F2ELM_INIT;
f2elm_t C24 = F2ELM_INIT;
f2elm_t A = F2ELM_INIT;
unsigned int m, index = 0, pts_index[MAX_INT_POINTS_ALICE], npts = 0, ii = 0;

// Initialize images of Bob's basis
fp2_decode(pkB, PKB[0]);
fp2_decode(pkB + SIDH_JINV_BYTESZ, PKB[1]);
fp2_decode(pkB + 2*SIDH_JINV_BYTESZ, PKB[2]);

// Initialize constants
get_A(PKB[0], PKB[1], PKB[2], A);
sike_fpadd(params.mont_one, params.mont_one, C24->c0);
sike_fp2add(A, C24, A24plus);
sike_fpadd(C24->c0, C24->c0, C24->c0);

// Retrieve kernel point
ladder3Pt(PKB[0], PKB[1], PKB[2], skA, 1, R, A);

// Traverse tree
index = 0;
for (size_t row = 1; row < A_max; row++) {
while (index < A_max-row) {
sike_fp2copy(R->X, pts[npts]->X);
sike_fp2copy(R->Z, pts[npts]->Z);
pts_index[npts++] = index;
m = params.A_strat[ii++];
xDBLe(R, R, A24plus, C24, (2*m));
index += m;
}
get_4_isog(R, A24plus, C24, coeff);

for (size_t i = 0; i < npts; i++) {
eval_4_isog(pts[i], coeff);
}

sike_fp2copy(pts[npts-1]->X, R->X);
sike_fp2copy(pts[npts-1]->Z, R->Z);
index = pts_index[npts-1];
npts -= 1;
}

get_4_isog(R, A24plus, C24, coeff);
sike_fp2add(A24plus, A24plus, A24plus);
sike_fp2sub(A24plus, C24, A24plus);
sike_fp2add(A24plus, A24plus, A24plus);
j_inv(A24plus, C24, jinv);
sike_fp2_encode(jinv, ssA);
}

// Bob's ephemeral shared secret computation
// It produces a shared secret key ssB using his secret key skB and Alice's public key pkA
// Inputs: Bob's skB is an integer in the range [0, 2^Floor(Log(2,3^159)) - 1], stored in 32 bytes.
// Alice's pkA consists of 3 GF(p503^2) elements encoded in 378 bytes.
// Output: a shared secret ssB that consists of one element in GF(p503^2) encoded in 126 bytes.
static void ex_iso_B(const uint8_t* skB, const uint8_t* pkA, uint8_t* ssB)
{
point_proj_t R, pts[MAX_INT_POINTS_BOB];
f2elm_t coeff[3], PKB[3], jinv;
f2elm_t A24plus = F2ELM_INIT;
f2elm_t A24minus = F2ELM_INIT;
f2elm_t A = F2ELM_INIT;
unsigned int m, index = 0, pts_index[MAX_INT_POINTS_BOB], npts = 0, ii = 0;

// Initialize images of Alice's basis
fp2_decode(pkA, PKB[0]);
fp2_decode(pkA + SIDH_JINV_BYTESZ, PKB[1]);
fp2_decode(pkA + 2*SIDH_JINV_BYTESZ, PKB[2]);

// Initialize constants
get_A(PKB[0], PKB[1], PKB[2], A);
sike_fpadd(params.mont_one, params.mont_one, A24minus->c0);
sike_fp2add(A, A24minus, A24plus);
sike_fp2sub(A, A24minus, A24minus);

// Retrieve kernel point
ladder3Pt(PKB[0], PKB[1], PKB[2], skB, 0, R, A);

// Traverse tree
index = 0;
for (size_t row = 1; row < B_max; row++) {
while (index < B_max-row) {
sike_fp2copy(R->X, pts[npts]->X);
sike_fp2copy(R->Z, pts[npts]->Z);
pts_index[npts++] = index;
m = params.B_strat[ii++];
xTPLe(R, R, A24minus, A24plus, m);
index += m;
}
get_3_isog(R, A24minus, A24plus, coeff);

for (size_t i = 0; i < npts; i++) {
eval_3_isog(pts[i], coeff);
}

sike_fp2copy(pts[npts-1]->X, R->X);
sike_fp2copy(pts[npts-1]->Z, R->Z);
index = pts_index[npts-1];
npts -= 1;
}

get_3_isog(R, A24minus, A24plus, coeff);
sike_fp2add(A24plus, A24minus, A);
sike_fp2add(A, A, A);
sike_fp2sub(A24plus, A24minus, A24plus);
j_inv(A, A24plus, jinv);
sike_fp2_encode(jinv, ssB);
}

int SIKE_keypair(uint8_t out_priv[SIKE_PRV_BYTESZ],
uint8_t out_pub[SIKE_PUB_BYTESZ]) {
// Calculate private key for Alice. Needs to be in range [0, 2^0xFA - 1] and <
// 253 bits
randombytes(out_priv, SIKE_PRV_BYTESZ);
out_priv[31] = (out_priv[31] | 0x01) & 0x03;

gen_iso_B(out_priv, out_pub);
return 1;
}

void SIKE_encaps(uint8_t out_shared_key[SIKE_SS_BYTESZ],
uint8_t out_ciphertext[SIKE_CT_BYTESZ],
const uint8_t pub_key[SIKE_PUB_BYTESZ]) {
// Secret buffer is reused by the function to store some ephemeral
// secret data. It's size must be maximum of 64,
// SIKE_MSG_BYTESZ and SIDH_PRV_A_BITSZ in bytes.
uint8_t secret[32]; // OZAPTF, why?
uint8_t j[SIDH_JINV_BYTESZ];
uint8_t temp[SIKE_MSG_BYTESZ + SIKE_CT_BYTESZ];
SHA256_CTX ctx;

// Generate secret key for A
// secret key A = SHA256({0,1}^n || pub_key)) mod SIDH_PRV_A_BITSZ
randombytes(temp, SIKE_MSG_BYTESZ);

sha256_init(&ctx);
sha256_update(&ctx, temp, SIKE_MSG_BYTESZ);
sha256_update(&ctx, pub_key, SIKE_PUB_BYTESZ);
sha256_final(&ctx, secret);

// Generate public key for A - first part of the ciphertext
gen_iso_A(secret, out_ciphertext);

// Generate c1:
// h = SHA256(j-invariant)
// c1 = h ^ m
ex_iso_A(secret, pub_key, j);
sha256_init(&ctx);
sha256_update(&ctx, j, sizeof(j));
sha256_final(&ctx, secret);

// c1 = h ^ m
uint8_t *c1 = &out_ciphertext[SIKE_PUB_BYTESZ];
for (size_t i = 0; i < SIKE_MSG_BYTESZ; i++) {
c1[i] = temp[i] ^ secret[i];
}

sha256_init(&ctx);
sha256_update(&ctx, temp, SIKE_MSG_BYTESZ);
sha256_update(&ctx, out_ciphertext, SIKE_CT_BYTESZ);
sha256_final(&ctx, secret);
// Generate shared secret out_shared_key = SHA256(m||out_ciphertext)
memcpy(out_shared_key, secret, SIKE_SS_BYTESZ);
}

void SIKE_decaps(uint8_t out_shared_key[SIKE_SS_BYTESZ],
const uint8_t ciphertext[SIKE_CT_BYTESZ],
const uint8_t pub_key[SIKE_PUB_BYTESZ],
const uint8_t priv_key[SIKE_PRV_BYTESZ]) {
// Secret buffer is reused by the function to store some ephemeral
// secret data. It's size must be maximum of 64,
// SIKE_MSG_BYTESZ and SIDH_PRV_A_BITSZ in bytes.
uint8_t secret[32];
uint8_t j[SIDH_JINV_BYTESZ];
uint8_t c0[SIKE_PUB_BYTESZ];
uint8_t temp[SIKE_MSG_BYTESZ];
uint8_t shared_nok[SIKE_MSG_BYTESZ];
SHA256_CTX ctx;

// This is OK as we are only using ephemeral keys in BoringSSL
randombytes(shared_nok, SIKE_MSG_BYTESZ);

// Recover m
// Let ciphertext = c0 || c1 - both have fixed sizes
// m = F(j-invariant(c0, priv_key)) ^ c1
ex_iso_B(priv_key, ciphertext, j);

sha256_init(&ctx);
sha256_update(&ctx, j, sizeof(j));
sha256_final(&ctx, secret);

const uint8_t *c1 = &ciphertext[sizeof(c0)];
for (size_t i = 0; i < SIKE_MSG_BYTESZ; i++) {
temp[i] = c1[i] ^ secret[i];
}

sha256_init(&ctx);
sha256_update(&ctx, temp, SIKE_MSG_BYTESZ);
sha256_update(&ctx, pub_key, SIKE_PUB_BYTESZ);
sha256_final(&ctx, secret);

// Recover c0 = public key A
gen_iso_A(secret, c0);
crypto_word_t ok = ct_uint_eq(
ct_mem_eq(c0, ciphertext, SIKE_PUB_BYTESZ), 1);
for (size_t i = 0; i < SIKE_MSG_BYTESZ; i++) {
temp[i] = ct_select_8(ok, temp[i], shared_nok[i]);
}

sha256_init(&ctx);
sha256_update(&ctx, temp, SIKE_MSG_BYTESZ);
sha256_update(&ctx, ciphertext, SIKE_CT_BYTESZ);
sha256_final(&ctx, secret);

// Generate shared secret out_shared_key = SHA256(m||ciphertext)
memcpy(out_shared_key, secret, SIKE_SS_BYTESZ);
}

+ 231
- 0
src/kem/sike/sike-p434-sha256/utils.h View File

@@ -0,0 +1,231 @@
/********************************************************************************************
* SIDH: an efficient supersingular isogeny cryptography library
*
* Abstract: internal header file for P434
*********************************************************************************************/

#ifndef UTILS_H_
#define UTILS_H_

#include <stddef.h>
#include <sike/sike.h>

// Conversion macro from number of bits to number of bytes
#define BITS_TO_BYTES(nbits) (((nbits)+7)/8)

// Bit size of the field
#define BITS_FIELD 434
// Byte size of the field
#define FIELD_BYTESZ BITS_TO_BYTES(BITS_FIELD)
// Number of 64-bit words of a 224-bit element
#define NBITS_ORDER 224
#define NWORDS64_ORDER ((NBITS_ORDER+63)/64)
// Number of elements in Alice's strategy
#define A_max 108
// Number of elements in Bob's strategy
#define B_max 137
// Word size size
#define RADIX sizeof(crypto_word_t)*8
// Byte size of a limb
#define LSZ sizeof(crypto_word_t)

#if defined(CPU_64_BIT)
typedef uint64_t crypto_word_t;
// Number of words of a 434-bit field element
#define NWORDS_FIELD 7
// Number of "0" digits in the least significant part of p434 + 1
#define ZERO_WORDS 3
// U64_TO_WORDS expands |x| for a |crypto_word_t| array literal.
#define U64_TO_WORDS(x) UINT64_C(x)
#else
typedef uint32_t crypto_word_t;
// Number of words of a 434-bit field element
#define NWORDS_FIELD 14
// Number of "0" digits in the least significant part of p434 + 1
#define ZERO_WORDS 6
// U64_TO_WORDS expands |x| for a |crypto_word_t| array literal.
#define U64_TO_WORDS(x) \
(uint32_t)(UINT64_C(x) & 0xffffffff), (uint32_t)(UINT64_C(x) >> 32)
#endif

// Extended datatype support
#if !defined(HAS_UINT128)
typedef uint64_t uint128_t[2];
#endif

// The following functions return 1 (TRUE) if condition is true, 0 (FALSE) otherwise
// Digit multiplication
#define MUL(multiplier, multiplicand, hi, lo) digit_x_digit((multiplier), (multiplicand), &(lo));

// If mask |x|==0xff.ff set |x| to 1, otherwise 0
#define M2B(x) ((x)>>(RADIX-1))

// Digit addition with carry
#define ADDC(carryIn, addend1, addend2, carryOut, sumOut) \
do { \
crypto_word_t tempReg = (addend1) + (crypto_word_t)(carryIn); \
(sumOut) = (addend2) + tempReg; \
(carryOut) = M2B(ct_uint_lt(tempReg, (crypto_word_t)(carryIn)) | \
ct_uint_lt((sumOut), tempReg)); \
} while(0)

// Digit subtraction with borrow
#define SUBC(borrowIn, minuend, subtrahend, borrowOut, differenceOut) \
do { \
crypto_word_t tempReg = (minuend) - (subtrahend); \
crypto_word_t borrowReg = M2B(ct_uint_lt((minuend), (subtrahend))); \
borrowReg |= ((borrowIn) & ct_uint_eq(tempReg, 0)); \
(differenceOut) = tempReg - (crypto_word_t)(borrowIn); \
(borrowOut) = borrowReg; \
} while(0)

/* Old GCC 4.9 (jessie) doesn't implement {0} initialization properly,
which violates C11 as described in 6.7.9, 21 (similarily C99, 6.7.8).
Defines below are used to work around the bug, and provide a way
to initialize f2elem_t and point_proj_t structs.
Bug has been fixed in GCC6 (debian stretch).
*/
#define F2ELM_INIT {{ {0}, {0} }}
#define POINT_PROJ_INIT {{ F2ELM_INIT, F2ELM_INIT }}

// Datatype for representing 434-bit field elements (448-bit max.)
// Elements over GF(p434) are encoded in 63 octets in little endian format
// (i.e., the least significant octet is located in the lowest memory address).
typedef crypto_word_t felm_t[NWORDS_FIELD];

// An element in F_{p^2}, is composed of two coefficients from F_p, * i.e.
// Fp2 element = c0 + c1*i in F_{p^2}
// Datatype for representing double-precision 2x434-bit field elements (448-bit max.)
// Elements (a+b*i) over GF(p434^2), where a and b are defined over GF(p434), are
// encoded as {a, b}, with a in the lowest memory portion.
typedef struct {
felm_t c0;
felm_t c1;
} fp2;

// Our F_{p^2} element type is a pointer to the struct.
typedef fp2 f2elm_t[1];

// Datatype for representing double-precision 2x434-bit
// field elements in contiguous memory.
typedef crypto_word_t dfelm_t[2*NWORDS_FIELD];

// Constants used during SIKE computation.
struct params_t {
// Stores a prime
const crypto_word_t prime[NWORDS_FIELD];
// Stores prime + 1
const crypto_word_t prime_p1[NWORDS_FIELD];
// Stores prime * 2
const crypto_word_t prime_x2[NWORDS_FIELD];
// Alice's generator values {XPA0 + XPA1*i, XQA0 + XQA1*i, XRA0 + XRA1*i}
// in GF(prime^2), expressed in Montgomery representation
const crypto_word_t A_gen[6*NWORDS_FIELD];
// Bob's generator values {XPB0 + XPB1*i, XQB0 + XQB1*i, XRB0 + XRB1*i}
// in GF(prime^2), expressed in Montgomery representation
const crypto_word_t B_gen[6*NWORDS_FIELD];
// Montgomery constant mont_R2 = (2^448)^2 mod prime
const crypto_word_t mont_R2[NWORDS_FIELD];
// Value 'one' in Montgomery representation
const crypto_word_t mont_one[NWORDS_FIELD];
// Value '6' in Montgomery representation
const crypto_word_t mont_six[NWORDS_FIELD];
// Fixed parameters for isogeny tree computation
const unsigned int A_strat[A_max-1];
const unsigned int B_strat[B_max-1];
};

// Point representation in projective XZ Montgomery coordinates.
typedef struct {
f2elm_t X;
f2elm_t Z;
} point_proj;
typedef point_proj point_proj_t[1];

// Checks whether two words are equal. Returns 1 in case it is,
// otherwise 0.
static inline crypto_word_t ct_uint_eq(crypto_word_t x, crypto_word_t y)
{
// if x==y then t = 0
crypto_word_t t = x ^ y;
// if x!=y t will have first bit set
t = (t >> 1) - t;
// return MSB - 1 in case x==y, otherwise 0
return ((~t) >> (RADIX-1));
}
// Constant time select.
// if pick == 1 (out = in1)
// if pick == 0 (out = in2)
// else out is undefined
static inline uint8_t ct_select_8(uint8_t flag, uint8_t in1, uint8_t in2) {
uint8_t mask = ((int8_t)(flag << 7))>>7;
return (in1&mask) | (in2&(~mask));
}

// Constant time memcmp. Returns 1 if p==q, otherwise 0
static inline int ct_mem_eq(const void *p, const void *q, size_t n)
{
const uint8_t *pp = (uint8_t*)p, *qq = (uint8_t*)q;
uint8_t a = 0;

while (n--) a |= *pp++ ^ *qq++;
return (ct_uint_eq(a, 0));
}

/*
// Returns 1 if x<y, otherwise 0
static inline crypto_word_t ct_uint_lt(crypto_word_t x, crypto_word_t y) {
const crypto_word_t t1 = x^y;
const crypto_word_t t2 = x - y;
const crypto_word_t tt = x ^ (t1 | (t2^y));
return (tt >> (RADIX-1));
}
*/

/// OZAPTF: coppied from boringssl
static inline crypto_word_t constant_time_msb_w(crypto_word_t a) {
return 0u - (a >> (sizeof(a) * 8 - 1));
}

// constant_time_lt_w returns 0xff..f if a < b and 0 otherwise.
static inline crypto_word_t ct_uint_lt(crypto_word_t x, crypto_word_t y)
{
/*
const crypto_word_t t1 = x^y;
const crypto_word_t t2 = x - y;
const crypto_word_t tt = x ^ (t1 | (t2^y));
return (tt >> (RADIX-1));
*/
// Consider the two cases of the problem:
// msb(a) == msb(b): a < b iff the MSB of a - b is set.
// msb(a) != msb(b): a < b iff the MSB of b is set.
//
// If msb(a) == msb(b) then the following evaluates as:
// msb(a^((a^b)|((a-b)^a))) ==
// msb(a^((a-b) ^ a)) == (because msb(a^b) == 0)
// msb(a^a^(a-b)) == (rearranging)
// msb(a-b) (because ∀x. x^x == 0)
//
// Else, if msb(a) != msb(b) then the following evaluates as:
// msb(a^((a^b)|((a-b)^a))) ==
// msb(a^(𝟙 | ((a-b)^a))) == (because msb(a^b) == 1 and 𝟙
// represents a value s.t. msb(𝟙) = 1)
// msb(a^𝟙) == (because ORing with 1 results in 1)
// msb(b)
//
//
// Here is an SMT-LIB verification of this formula:
//
// (define-fun lt ((a (_ BitVec 32)) (b (_ BitVec 32))) (_ BitVec 32)
// (bvxor a (bvor (bvxor a b) (bvxor (bvsub a b) a)))
// )
//
// (declare-fun a () (_ BitVec 32))
// (declare-fun b () (_ BitVec 32))
//
// (assert (not (= (= #x00000001 (bvlshr (lt a b) #x0000001f)) (bvult a b))))
// (check-sat)
// (get-model)
return constant_time_msb_w(x^((x^y)|((x-y)^x)));
}
#endif // UTILS_H_

Write Preview
Loading…
Cancel
Save