/* * Copyright (C) 2021 Denys Vlasenko * * Licensed under GPLv2, see file LICENSE in this source tree. */ #include "tls.h" #define SP_DEBUG 0 #define FIXED_SECRET 0 #define FIXED_PEER_PUBKEY 0 #define ALLOW_ASM 1 #if SP_DEBUG # define dbg(...) fprintf(stderr, __VA_ARGS__) static void dump_hex(const char *fmt, const void *vp, int len) { char hexbuf[32 * 1024 + 4]; const uint8_t *p = vp; bin2hex(hexbuf, (void*)p, len)[0] = '\0'; dbg(fmt, hexbuf); } #else # define dbg(...) ((void)0) # define dump_hex(...) ((void)0) #endif typedef uint32_t sp_digit; typedef int32_t signed_sp_digit; /* The code below is taken from parts of * wolfssl-3.15.3/wolfcrypt/src/sp_c32.c * and heavily modified. */ typedef struct sp_point { sp_digit x[2 * 8]; sp_digit y[2 * 8]; sp_digit z[2 * 8]; int infinity; } sp_point; /* The modulus (prime) of the curve P256. */ static const sp_digit p256_mod[8] = { 0xffffffff,0xffffffff,0xffffffff,0x00000000, 0x00000000,0x00000000,0x00000001,0xffffffff, }; #define p256_mp_mod ((sp_digit)0x000001) /* Normalize the values in each word to 32 bits - NOP */ #define sp_256_norm_8(a) ((void)0) /* Write r as big endian to byte array. * Fixed length number of bytes written: 32 * * r A single precision integer. * a Byte array. */ static void sp_256_to_bin_8(const sp_digit* r, uint8_t* a) { int i; sp_256_norm_8(r); r += 8; for (i = 0; i < 8; i++) { r--; move_to_unaligned32(a, SWAP_BE32(*r)); a += 4; } } /* Read big endian unsigned byte array into r. * * r A single precision integer. * a Byte array. * n Number of bytes in array to read. */ static void sp_256_from_bin_8(sp_digit* r, const uint8_t* a) { int i; r += 8; for (i = 0; i < 8; i++) { sp_digit v; move_from_unaligned32(v, a); *--r = SWAP_BE32(v); a += 4; } } #if SP_DEBUG static void dump_256(const char *fmt, const sp_digit* r) { uint8_t b32[32]; sp_256_to_bin_8(r, b32); dump_hex(fmt, b32, 32); } static void dump_512(const char *fmt, const sp_digit* r) { uint8_t b64[64]; sp_256_to_bin_8(r, b64 + 32); sp_256_to_bin_8(r+8, b64); dump_hex(fmt, b64, 64); } #else # define dump_256(...) ((void)0) # define dump_512(...) ((void)0) #endif /* Convert a point of big-endian 32-byte x,y pair to type sp_point. */ static void sp_256_point_from_bin2x32(sp_point* p, const uint8_t *bin2x32) { memset(p, 0, sizeof(*p)); /*p->infinity = 0;*/ sp_256_from_bin_8(p->x, bin2x32); sp_256_from_bin_8(p->y, bin2x32 + 32); p->z[0] = 1; /* p->z = 1 */ } /* Compare a with b. * * return -ve, 0 or +ve if a is less than, equal to or greater than b * respectively. */ static signed_sp_digit sp_256_cmp_8(const sp_digit* a, const sp_digit* b) { int i; for (i = 7; i >= 0; i--) { /* signed_sp_digit r = a[i] - b[i]; * if (r != 0) * return r; * does not work: think about a[i]=0, b[i]=0xffffffff */ if (a[i] == b[i]) continue; return (a[i] > b[i]) * 2 - 1; } return 0; } /* Compare two numbers to determine if they are equal. * * return 1 when equal and 0 otherwise. */ static int sp_256_cmp_equal_8(const sp_digit* a, const sp_digit* b) { return sp_256_cmp_8(a, b) == 0; } /* Add b to a into r. (r = a + b). Return !0 on overflow */ static int sp_256_add_8(sp_digit* r, const sp_digit* a, const sp_digit* b) { #if ALLOW_ASM && defined(__GNUC__) && defined(__i386__) sp_digit reg; asm volatile ( "\n movl (%0), %3" "\n addl (%1), %3" "\n movl %3, (%2)" "\n" "\n movl 1*4(%0), %3" "\n adcl 1*4(%1), %3" "\n movl %3, 1*4(%2)" "\n" "\n movl 2*4(%0), %3" "\n adcl 2*4(%1), %3" "\n movl %3, 2*4(%2)" "\n" "\n movl 3*4(%0), %3" "\n adcl 3*4(%1), %3" "\n movl %3, 3*4(%2)" "\n" "\n movl 4*4(%0), %3" "\n adcl 4*4(%1), %3" "\n movl %3, 4*4(%2)" "\n" "\n movl 5*4(%0), %3" "\n adcl 5*4(%1), %3" "\n movl %3, 5*4(%2)" "\n" "\n movl 6*4(%0), %3" "\n adcl 6*4(%1), %3" "\n movl %3, 6*4(%2)" "\n" "\n movl 7*4(%0), %3" "\n adcl 7*4(%1), %3" "\n movl %3, 7*4(%2)" "\n" "\n sbbl %3, %3" "\n" : "=r" (a), "=r" (b), "=r" (r), "=r" (reg) : "0" (a), "1" (b), "2" (r) : "memory" ); return reg; #else int i; sp_digit carry; carry = 0; for (i = 0; i < 8; i++) { sp_digit w, v; w = b[i] + carry; v = a[i]; if (w != 0) { v = a[i] + w; carry = (v < a[i]); /* hope compiler detects above as "carry flag set" */ } /* else: b + carry == 0, two cases: * b:ffffffff, carry:1 * b:00000000, carry:0 * in either case, r[i] = a[i] and carry remains unchanged */ r[i] = v; } return carry; #endif } /* Sub b from a into r. (r = a - b). Return !0 on underflow */ static int sp_256_sub_8(sp_digit* r, const sp_digit* a, const sp_digit* b) { #if ALLOW_ASM && defined(__GNUC__) && defined(__i386__) sp_digit reg; asm volatile ( "\n movl (%0), %3" "\n subl (%1), %3" "\n movl %3, (%2)" "\n" "\n movl 1*4(%0), %3" "\n sbbl 1*4(%1), %3" "\n movl %3, 1*4(%2)" "\n" "\n movl 2*4(%0), %3" "\n sbbl 2*4(%1), %3" "\n movl %3, 2*4(%2)" "\n" "\n movl 3*4(%0), %3" "\n sbbl 3*4(%1), %3" "\n movl %3, 3*4(%2)" "\n" "\n movl 4*4(%0), %3" "\n sbbl 4*4(%1), %3" "\n movl %3, 4*4(%2)" "\n" "\n movl 5*4(%0), %3" "\n sbbl 5*4(%1), %3" "\n movl %3, 5*4(%2)" "\n" "\n movl 6*4(%0), %3" "\n sbbl 6*4(%1), %3" "\n movl %3, 6*4(%2)" "\n" "\n movl 7*4(%0), %3" "\n sbbl 7*4(%1), %3" "\n movl %3, 7*4(%2)" "\n" "\n sbbl %3, %3" "\n" : "=r" (a), "=r" (b), "=r" (r), "=r" (reg) : "0" (a), "1" (b), "2" (r) : "memory" ); return reg; #else int i; sp_digit borrow; borrow = 0; for (i = 0; i < 8; i++) { sp_digit w, v; w = b[i] + borrow; v = a[i]; if (w != 0) { v = a[i] - w; borrow = (v > a[i]); /* hope compiler detects above as "carry flag set" */ } /* else: b + borrow == 0, two cases: * b:ffffffff, borrow:1 * b:00000000, borrow:0 * in either case, r[i] = a[i] and borrow remains unchanged */ r[i] = v; } return borrow; #endif } #if ALLOW_ASM && defined(__GNUC__) && defined(__i386__) /* Sub p256_mod from a into r. (r = a - p256_mod). */ static void sp_256_sub_8_p256_mod(sp_digit* r, const sp_digit* a) { sp_digit reg; //p256_mod[7..0] = ffffffff 00000001 00000000 00000000 00000000 ffffffff ffffffff ffffffff asm volatile ( "\n movl (%0), %2" "\n subl $0xffffffff, %2" "\n movl %2, (%1)" "\n" "\n movl 1*4(%0), %2" "\n sbbl $0xffffffff, %2" "\n movl %2, 1*4(%1)" "\n" "\n movl 2*4(%0), %2" "\n sbbl $0xffffffff, %2" "\n movl %2, 2*4(%1)" "\n" "\n movl 3*4(%0), %2" "\n sbbl $0, %2" "\n movl %2, 3*4(%1)" "\n" "\n movl 4*4(%0), %2" "\n sbbl $0, %2" "\n movl %2, 4*4(%1)" "\n" "\n movl 5*4(%0), %2" "\n sbbl $0, %2" "\n movl %2, 5*4(%1)" "\n" "\n movl 6*4(%0), %2" "\n sbbl $1, %2" "\n movl %2, 6*4(%1)" "\n" "\n movl 7*4(%0), %2" "\n sbbl $0xffffffff, %2" "\n movl %2, 7*4(%1)" "\n" : "=r" (a), "=r" (r), "=r" (reg) : "0" (a), "1" (r) : "memory" ); } #else # define sp_256_sub_8_p256_mod(r, a) sp_256_sub_8((r), (a), p256_mod) #endif /* Multiply a and b into r. (r = a * b) */ static void sp_256_mul_8(sp_digit* r, const sp_digit* a, const sp_digit* b) { #if ALLOW_ASM && defined(__GNUC__) && defined(__i386__) sp_digit rr[15]; /* in case r coincides with a or b */ int k; uint32_t accl; uint32_t acch; acch = accl = 0; for (k = 0; k < 15; k++) { int i, j; uint32_t acc_hi; i = k - 7; if (i < 0) i = 0; j = k - i; acc_hi = 0; do { //////////////////////// // uint64_t m = ((uint64_t)a[i]) * b[j]; // acc_hi:acch:accl += m; asm volatile ( // a[i] is already loaded in %%eax "\n mull %7" "\n addl %%eax, %0" "\n adcl %%edx, %1" "\n adcl $0, %2" : "=rm" (accl), "=rm" (acch), "=rm" (acc_hi) : "0" (accl), "1" (acch), "2" (acc_hi), "a" (a[i]), "m" (b[j]) : "cc", "dx" ); //////////////////////// j--; i++; } while (i != 8 && i <= k); rr[k] = accl; accl = acch; acch = acc_hi; } r[15] = accl; memcpy(r, rr, sizeof(rr)); #elif 0 //TODO: arm assembly (untested) sp_digit tmp[16]; asm volatile ( "\n mov r5, #0" "\n mov r6, #0" "\n mov r7, #0" "\n mov r8, #0" "\n 1:" "\n subs r3, r5, #28" "\n movcc r3, #0" "\n sub r4, r5, r3" "\n 2:" "\n ldr r14, [%[a], r3]" "\n ldr r12, [%[b], r4]" "\n umull r9, r10, r14, r12" "\n adds r6, r6, r9" "\n adcs r7, r7, r10" "\n adc r8, r8, #0" "\n add r3, r3, #4" "\n sub r4, r4, #4" "\n cmp r3, #32" "\n beq 3f" "\n cmp r3, r5" "\n ble 2b" "\n 3:" "\n str r6, [%[r], r5]" "\n mov r6, r7" "\n mov r7, r8" "\n mov r8, #0" "\n add r5, r5, #4" "\n cmp r5, #56" "\n ble 1b" "\n str r6, [%[r], r5]" : [r] "r" (tmp), [a] "r" (a), [b] "r" (b) : "memory", "r3", "r4", "r5", "r6", "r7", "r8", "r9", "r10", "r14", "r12" ); memcpy(r, tmp, sizeof(tmp)); #else sp_digit rr[15]; /* in case r coincides with a or b */ int i, j, k; uint64_t acc; acc = 0; for (k = 0; k < 15; k++) { uint32_t acc_hi; i = k - 7; if (i < 0) i = 0; j = k - i; acc_hi = 0; do { uint64_t m = ((uint64_t)a[i]) * b[j]; acc += m; if (acc < m) acc_hi++; j--; i++; } while (i != 8 && i <= k); rr[k] = acc; acc = (acc >> 32) | ((uint64_t)acc_hi << 32); } r[15] = acc; memcpy(r, rr, sizeof(rr)); #endif } /* Shift number right one bit. Bottom bit is lost. */ static void sp_256_rshift1_8(sp_digit* r, sp_digit* a, sp_digit carry) { int i; carry = (!!carry << 31); for (i = 7; i >= 0; i--) { sp_digit c = a[i] << 31; r[i] = (a[i] >> 1) | carry; carry = c; } } /* Divide the number by 2 mod the modulus (prime). (r = a / 2 % m) */ static void sp_256_div2_8(sp_digit* r, const sp_digit* a, const sp_digit* m) { int carry = 0; if (a[0] & 1) carry = sp_256_add_8(r, a, m); sp_256_norm_8(r); sp_256_rshift1_8(r, r, carry); } /* Add two Montgomery form numbers (r = a + b % m) */ static void sp_256_mont_add_8(sp_digit* r, const sp_digit* a, const sp_digit* b /*, const sp_digit* m*/) { // const sp_digit* m = p256_mod; int carry = sp_256_add_8(r, a, b); sp_256_norm_8(r); if (carry) { sp_256_sub_8_p256_mod(r, r /*, m*/); sp_256_norm_8(r); } } /* Subtract two Montgomery form numbers (r = a - b % m) */ static void sp_256_mont_sub_8(sp_digit* r, const sp_digit* a, const sp_digit* b /*, const sp_digit* m*/) { const sp_digit* m = p256_mod; int borrow; borrow = sp_256_sub_8(r, a, b); sp_256_norm_8(r); if (borrow) { sp_256_add_8(r, r, m); sp_256_norm_8(r); } } /* Double a Montgomery form number (r = a + a % m) */ static void sp_256_mont_dbl_8(sp_digit* r, const sp_digit* a /*, const sp_digit* m*/) { // const sp_digit* m = p256_mod; int carry = sp_256_add_8(r, a, a); sp_256_norm_8(r); if (carry) sp_256_sub_8_p256_mod(r, r /*, m*/); sp_256_norm_8(r); } /* Triple a Montgomery form number (r = a + a + a % m) */ static void sp_256_mont_tpl_8(sp_digit* r, const sp_digit* a /*, const sp_digit* m*/) { // const sp_digit* m = p256_mod; int carry = sp_256_add_8(r, a, a); sp_256_norm_8(r); if (carry) { sp_256_sub_8_p256_mod(r, r /*, m*/); sp_256_norm_8(r); } carry = sp_256_add_8(r, r, a); sp_256_norm_8(r); if (carry) { sp_256_sub_8_p256_mod(r, r /*, m*/); sp_256_norm_8(r); } } /* Shift the result in the high 256 bits down to the bottom. */ static void sp_256_mont_shift_8(sp_digit* r, const sp_digit* a) { int i; for (i = 0; i < 8; i++) { r[i] = a[i+8]; r[i+8] = 0; } } /* Mul a by scalar b and add into r. (r += a * b) */ static int sp_256_mul_add_8(sp_digit* r /*, const sp_digit* a, sp_digit b*/) { // const sp_digit* a = p256_mod; //a[7..0] = ffffffff 00000001 00000000 00000000 00000000 ffffffff ffffffff ffffffff sp_digit b = r[0]; uint64_t t; // t = 0; // for (i = 0; i < 8; i++) { // uint32_t t_hi; // uint64_t m = ((uint64_t)b * a[i]) + r[i]; // t += m; // t_hi = (t < m); // r[i] = (sp_digit)t; // t = (t >> 32) | ((uint64_t)t_hi << 32); // } // r[8] += (sp_digit)t; // Unroll, then optimize the above loop: //uint32_t t_hi; uint64_t m; uint32_t t32; //m = ((uint64_t)b * a[0]) + r[0]; // Since b is r[0] and a[0] is ffffffff, the above optimizes to: // m = r[0] * ffffffff + r[0] = (r[0] * 100000000 - r[0]) + r[0] = r[0] << 32; //t += m; // t = r[0] << 32 = b << 32; //t_hi = (t < m); // t_hi = 0; //r[0] = (sp_digit)t; r[0] = 0; //t = (t >> 32) | ((uint64_t)t_hi << 32); // t = b; //m = ((uint64_t)b * a[1]) + r[1]; // Since a[1] is ffffffff, the above optimizes to: // m = b * ffffffff + r[1] = (b * 100000000 - b) + r[1] = (b << 32) - b + r[1]; //t += m; // t = b + (b << 32) - b + r[1] = (b << 32) + r[1]; //t_hi = (t < m); // t_hi = 0; //r[1] = (sp_digit)t; // r[1] = r[1]; //t = (t >> 32) | ((uint64_t)t_hi << 32); // t = b; //m = ((uint64_t)b * a[2]) + r[2]; // Since a[2] is ffffffff, the above optimizes to: // m = b * ffffffff + r[2] = (b * 100000000 - b) + r[2] = (b << 32) - b + r[2]; //t += m; // t = b + (b << 32) - b + r[2] = (b << 32) + r[2] //t_hi = (t < m); // t_hi = 0; //r[2] = (sp_digit)t; // r[2] = r[2]; //t = (t >> 32) | ((uint64_t)t_hi << 32); // t = b; //m = ((uint64_t)b * a[3]) + r[3]; // Since a[3] is 00000000, the above optimizes to: // m = b * 0 + r[3] = r[3]; //t += m; // t = b + r[3]; //t_hi = (t < m); // t_hi = 0; //r[3] = (sp_digit)t; r[3] = r[3] + b; //t = (t >> 32) | ((uint64_t)t_hi << 32); t32 = (r[3] < b); // 0 or 1 //m = ((uint64_t)b * a[4]) + r[4]; // Since a[4] is 00000000, the above optimizes to: // m = b * 0 + r[4] = r[4]; //t += m; // t = t32 + r[4]; //t_hi = (t < m); // t_hi = 0; //r[4] = (sp_digit)t; //t = (t >> 32) | ((uint64_t)t_hi << 32); if (t32 != 0) { r[4]++; t32 = (r[4] == 0); // 0 or 1 //m = ((uint64_t)b * a[5]) + r[5]; // Since a[5] is 00000000, the above optimizes to: // m = b * 0 + r[5] = r[5]; //t += m; // t = t32 + r[5]; (t32 is 0 or 1) //t_hi = (t < m); // t_hi = 0; //r[5] = (sp_digit)t; //t = (t >> 32) | ((uint64_t)t_hi << 32); if (t32 != 0) { r[5]++; t32 = (r[5] == 0); // 0 or 1 } } //m = ((uint64_t)b * a[6]) + r[6]; // Since a[6] is 00000001, the above optimizes to: // m = (uint64_t)b + r[6]; // 33 bits at most //t += m; t = t32 + (uint64_t)b + r[6]; //t_hi = (t < m); // t_hi = 0; r[6] = (sp_digit)t; //t = (t >> 32) | ((uint64_t)t_hi << 32); t = (t >> 32); //m = ((uint64_t)b * a[7]) + r[7]; // Since a[7] is ffffffff, the above optimizes to: // m = b * ffffffff + r[7] = (b * 100000000 - b) + r[7] m = ((uint64_t)b << 32) - b + r[7]; t += m; //t_hi = (t < m); // t_hi in fact is always 0 here (256bit * 32bit can't have more than 32 bits of overflow) r[7] = (sp_digit)t; //t = (t >> 32) | ((uint64_t)t_hi << 32); t = (t >> 32); r[8] += (sp_digit)t; return (r[8] < (sp_digit)t); /* 1 if addition overflowed */ } /* Reduce the number back to 256 bits using Montgomery reduction. * * a A single precision number to reduce in place. * m The single precision number representing the modulus. * mp The digit representing the negative inverse of m mod 2^n. */ static void sp_256_mont_reduce_8(sp_digit* a/*, const sp_digit* m, sp_digit mp*/) { // const sp_digit* m = p256_mod; sp_digit mp = p256_mp_mod; int i; // sp_digit mu; if (mp != 1) { sp_digit word16th = 0; for (i = 0; i < 8; i++) { // mu = (sp_digit)(a[i] * mp); if (sp_256_mul_add_8(a+i /*, m, mu*/)) { int j = i + 8; inc_next_word0: if (++j > 15) { /* a[16] array has no more words? */ word16th++; continue; } if (++a[j] == 0) /* did this overflow too? */ goto inc_next_word0; } } sp_256_mont_shift_8(a, a); if (word16th != 0) sp_256_sub_8_p256_mod(a, a /*, m*/); sp_256_norm_8(a); } else { /* Same code for explicit mp == 1 (which is always the case for P256) */ sp_digit word16th = 0; for (i = 0; i < 8; i++) { /*mu = a[i];*/ if (sp_256_mul_add_8(a+i /*, m, mu*/)) { int j = i + 8; inc_next_word: if (++j > 15) { /* a[16] array has no more words? */ word16th++; continue; } if (++a[j] == 0) /* did this overflow too? */ goto inc_next_word; } } sp_256_mont_shift_8(a, a); if (word16th != 0) sp_256_sub_8_p256_mod(a, a /*, m*/); sp_256_norm_8(a); } } #if 0 //TODO: arm32 asm (also adapt for x86?) static void sp_256_mont_reduce_8(sp_digit* a, sp_digit* m, sp_digit mp) { sp_digit ca = 0; asm volatile ( # i = 0 mov r12, #0 # i = 0 ldr r10, [%[a], #0] # r10 = a[0] ldr r14, [%[a], #4] # r14 = a[1] 1: # mu = a[i] * mp # mul r8, %[mp], r10 # mu = a[i] * mp # a[i+0] += m[0] * mu # ldr r7, [%[m], #0] # a[i+0] += m[0] * mu ldr r9, [%[a], #0] # umull r6, r7, r8, r7 # r7:r6 = mu * m[0] adds r10, r10, r6 # r5:r10 += r7:r6 adc r5, r7, #0 # # a[i+1] += m[1] * mu # ldr r7, [%[m], #4] # a[i+1] += m[1] * mu ldr r9, [%[a], #4] # umull r6, r7, r8, r7 # r7:r6 = mu * m[1] adds r10, r14, r6 # r4:r10 = r7:r14 + r7:r6 adc r4, r7, #0 # adds r10, r10, r5 # r4:r10 += r5 adc r4, r4, #0 # # a[i+2] += m[2] * mu # ldr r7, [%[m], #8] # a[i+2] += m[2] * mu ldr r14, [%[a], #8] # umull r6, r7, r8, r7 # adds r14, r14, r6 # adc r5, r7, #0 # adds r14, r14, r4 # adc r5, r5, #0 # # a[i+3] += m[3] * mu # ldr r7, [%[m], #12] # a[i+3] += m[3] * mu ldr r9, [%[a], #12] # umull r6, r7, r8, r7 # adds r9, r9, r6 # adc r4, r7, #0 # adds r9, r9, r5 # str r9, [%[a], #12] # a[3] = r9 adc r4, r4, #0 # # a[i+4] += m[4] * mu # ldr r7, [%[m], #16] # a[i+4] += m[4] * mu ldr r9, [%[a], #16] # umull r6, r7, r8, r7 # adds r9, r9, r6 # adc r5, r7, #0 # adds r9, r9, r4 # str r9, [%[a], #16] # a[4] = r9 adc r5, r5, #0 # # a[i+5] += m[5] * mu # ldr r7, [%[m], #20] # a[i+5] += m[5] * mu ldr r9, [%[a], #20] # umull r6, r7, r8, r7 # adds r9, r9, r6 # adc r4, r7, #0 # adds r9, r9, r5 # str r9, [%[a], #20] # a[5] = r9 adc r4, r4, #0 # # a[i+6] += m[6] * mu # ldr r7, [%[m], #24] # a[i+6] += m[6] * mu ldr r9, [%[a], #24] # umull r6, r7, r8, r7 # adds r9, r9, r6 # adc r5, r7, #0 # adds r9, r9, r4 # str r9, [%[a], #24] # a[6] = r9 adc r5, r5, #0 # # a[i+7] += m[7] * mu # ldr r7, [%[m], #28] # a[i+7] += m[7] * mu ldr r9, [%[a], #28] # umull r6, r7, r8, r7 # adds r5, r5, r6 # adcs r7, r7, %[ca] # mov %[ca], #0 # adc %[ca], %[ca], %[ca] # ca = CF adds r9, r9, r5 # str r9, [%[a], #28] # a[7] = r9 ldr r9, [%[a], #32] # r9 = a[8] adcs r9, r9, r7 # str r9, [%[a], #32] # a[8] = r9 adc %[ca], %[ca], #0 # ca += CF # i += 1 # i++ add %[a], %[a], #4 # a++ add r12, r12, #4 # i += 4 cmp r12, #32 # if (i < 32) blt 1b # goto 1 str r10, [%[a], #0] # a[0] = r10 str r14, [%[a], #4] # a[1] = r14 : [ca] "+r" (ca), [a] "+r" (a) : [m] "r" (m), [mp] "r" (mp) : "memory", "r4", "r5", "r6", "r7", "r8", "r9", "r10", "r14", "r12" ); if (ca) a -= m; } #endif /* Multiply two Montogmery form numbers mod the modulus (prime). * (r = a * b mod m) * * r Result of multiplication. * a First number to multiply in Montogmery form. * b Second number to multiply in Montogmery form. * m Modulus (prime). * mp Montogmery mulitplier. */ static void sp_256_mont_mul_8(sp_digit* r, const sp_digit* a, const sp_digit* b /*, const sp_digit* m, sp_digit mp*/) { //const sp_digit* m = p256_mod; //sp_digit mp = p256_mp_mod; sp_256_mul_8(r, a, b); sp_256_mont_reduce_8(r /*, m, mp*/); } /* Square the Montgomery form number. (r = a * a mod m) * * r Result of squaring. * a Number to square in Montogmery form. * m Modulus (prime). * mp Montogmery mulitplier. */ static void sp_256_mont_sqr_8(sp_digit* r, const sp_digit* a /*, const sp_digit* m, sp_digit mp*/) { //const sp_digit* m = p256_mod; //sp_digit mp = p256_mp_mod; sp_256_mont_mul_8(r, a, a /*, m, mp*/); } /* Invert the number, in Montgomery form, modulo the modulus (prime) of the * P256 curve. (r = 1 / a mod m) * * r Inverse result. * a Number to invert. */ #if 0 /* Mod-2 for the P256 curve. */ static const uint32_t p256_mod_2[8] = { 0xfffffffd,0xffffffff,0xffffffff,0x00000000, 0x00000000,0x00000000,0x00000001,0xffffffff, }; //Bit pattern: //2 2 2 2 2 2 2 1...1 //5 5 4 3 2 1 0 9...0 9...1 //543210987654321098765432109876543210987654321098765432109876543210...09876543210...09876543210 //111111111111111111111111111111110000000000000000000000000000000100...00000111111...11111111101 #endif static void sp_256_mont_inv_8(sp_digit* r, sp_digit* a) { sp_digit t[2*8]; //can be just [8]? int i; memcpy(t, a, sizeof(sp_digit) * 8); for (i = 254; i >= 0; i--) { sp_256_mont_sqr_8(t, t /*, p256_mod, p256_mp_mod*/); /*if (p256_mod_2[i / 32] & ((sp_digit)1 << (i % 32)))*/ if (i >= 224 || i == 192 || (i <= 95 && i != 1)) sp_256_mont_mul_8(t, t, a /*, p256_mod, p256_mp_mod*/); } memcpy(r, t, sizeof(sp_digit) * 8); } /* Multiply a number by Montogmery normalizer mod modulus (prime). * * r The resulting Montgomery form number. * a The number to convert. */ static void sp_256_mod_mul_norm_8(sp_digit* r, const sp_digit* a) { int64_t t[8]; int32_t o; #define A(n) ((uint64_t)a[n]) /* 1 1 0 -1 -1 -1 -1 0 */ t[0] = 0 + A(0) + A(1) - A(3) - A(4) - A(5) - A(6); /* 0 1 1 0 -1 -1 -1 -1 */ t[1] = 0 + A(1) + A(2) - A(4) - A(5) - A(6) - A(7); /* 0 0 1 1 0 -1 -1 -1 */ t[2] = 0 + A(2) + A(3) - A(5) - A(6) - A(7); /* -1 -1 0 2 2 1 0 -1 */ t[3] = 0 - A(0) - A(1) + 2 * A(3) + 2 * A(4) + A(5) - A(7); /* 0 -1 -1 0 2 2 1 0 */ t[4] = 0 - A(1) - A(2) + 2 * A(4) + 2 * A(5) + A(6); /* 0 0 -1 -1 0 2 2 1 */ t[5] = 0 - A(2) - A(3) + 2 * A(5) + 2 * A(6) + A(7); /* -1 -1 0 0 0 1 3 2 */ t[6] = 0 - A(0) - A(1) + A(5) + 3 * A(6) + 2 * A(7); /* 1 0 -1 -1 -1 -1 0 3 */ t[7] = 0 + A(0) - A(2) - A(3) - A(4) - A(5) + 3 * A(7); #undef A t[1] += t[0] >> 32; t[0] &= 0xffffffff; t[2] += t[1] >> 32; t[1] &= 0xffffffff; t[3] += t[2] >> 32; t[2] &= 0xffffffff; t[4] += t[3] >> 32; t[3] &= 0xffffffff; t[5] += t[4] >> 32; t[4] &= 0xffffffff; t[6] += t[5] >> 32; t[5] &= 0xffffffff; t[7] += t[6] >> 32; t[6] &= 0xffffffff; o = t[7] >> 32; //t[7] &= 0xffffffff; t[0] += o; t[3] -= o; t[6] -= o; t[7] += o; r[0] = (sp_digit)t[0]; t[1] += t[0] >> 32; r[1] = (sp_digit)t[1]; t[2] += t[1] >> 32; r[2] = (sp_digit)t[2]; t[3] += t[2] >> 32; r[3] = (sp_digit)t[3]; t[4] += t[3] >> 32; r[4] = (sp_digit)t[4]; t[5] += t[4] >> 32; r[5] = (sp_digit)t[5]; t[6] += t[5] >> 32; r[6] = (sp_digit)t[6]; // t[7] += t[6] >> 32; // r[7] = (sp_digit)t[7]; r[7] = (sp_digit)t[7] + (sp_digit)(t[6] >> 32); } /* Map the Montgomery form projective co-ordinate point to an affine point. * * r Resulting affine co-ordinate point. * p Montgomery form projective co-ordinate point. */ static void sp_256_map_8(sp_point* r, sp_point* p) { sp_digit t1[2*8]; sp_digit t2[2*8]; sp_256_mont_inv_8(t1, p->z); sp_256_mont_sqr_8(t2, t1 /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t1, t2, t1 /*, p256_mod, p256_mp_mod*/); /* x /= z^2 */ sp_256_mont_mul_8(r->x, p->x, t2 /*, p256_mod, p256_mp_mod*/); memset(r->x + 8, 0, sizeof(r->x) / 2); sp_256_mont_reduce_8(r->x /*, p256_mod, p256_mp_mod*/); /* Reduce x to less than modulus */ if (sp_256_cmp_8(r->x, p256_mod) >= 0) sp_256_sub_8_p256_mod(r->x, r->x /*, p256_mod*/); sp_256_norm_8(r->x); /* y /= z^3 */ sp_256_mont_mul_8(r->y, p->y, t1 /*, p256_mod, p256_mp_mod*/); memset(r->y + 8, 0, sizeof(r->y) / 2); sp_256_mont_reduce_8(r->y /*, p256_mod, p256_mp_mod*/); /* Reduce y to less than modulus */ if (sp_256_cmp_8(r->y, p256_mod) >= 0) sp_256_sub_8_p256_mod(r->y, r->y /*, p256_mod*/); sp_256_norm_8(r->y); memset(r->z, 0, sizeof(r->z)); r->z[0] = 1; } /* Double the Montgomery form projective point p. * * r Result of doubling point. * p Point to double. */ static void sp_256_proj_point_dbl_8(sp_point* r, sp_point* p) { sp_digit t1[2*8]; sp_digit t2[2*8]; /* Put point to double into result */ if (r != p) *r = *p; /* struct copy */ if (r->infinity) return; if (SP_DEBUG) { /* unused part of t2, may result in spurios * differences in debug output. Clear it. */ memset(t2, 0, sizeof(t2)); } /* T1 = Z * Z */ sp_256_mont_sqr_8(t1, r->z /*, p256_mod, p256_mp_mod*/); /* Z = Y * Z */ sp_256_mont_mul_8(r->z, r->y, r->z /*, p256_mod, p256_mp_mod*/); /* Z = 2Z */ sp_256_mont_dbl_8(r->z, r->z /*, p256_mod*/); /* T2 = X - T1 */ sp_256_mont_sub_8(t2, r->x, t1 /*, p256_mod*/); /* T1 = X + T1 */ sp_256_mont_add_8(t1, r->x, t1 /*, p256_mod*/); /* T2 = T1 * T2 */ sp_256_mont_mul_8(t2, t1, t2 /*, p256_mod, p256_mp_mod*/); /* T1 = 3T2 */ sp_256_mont_tpl_8(t1, t2 /*, p256_mod*/); /* Y = 2Y */ sp_256_mont_dbl_8(r->y, r->y /*, p256_mod*/); /* Y = Y * Y */ sp_256_mont_sqr_8(r->y, r->y /*, p256_mod, p256_mp_mod*/); /* T2 = Y * Y */ sp_256_mont_sqr_8(t2, r->y /*, p256_mod, p256_mp_mod*/); /* T2 = T2/2 */ sp_256_div2_8(t2, t2, p256_mod); /* Y = Y * X */ sp_256_mont_mul_8(r->y, r->y, r->x /*, p256_mod, p256_mp_mod*/); /* X = T1 * T1 */ sp_256_mont_mul_8(r->x, t1, t1 /*, p256_mod, p256_mp_mod*/); /* X = X - Y */ sp_256_mont_sub_8(r->x, r->x, r->y /*, p256_mod*/); /* X = X - Y */ sp_256_mont_sub_8(r->x, r->x, r->y /*, p256_mod*/); /* Y = Y - X */ sp_256_mont_sub_8(r->y, r->y, r->x /*, p256_mod*/); /* Y = Y * T1 */ sp_256_mont_mul_8(r->y, r->y, t1 /*, p256_mod, p256_mp_mod*/); /* Y = Y - T2 */ sp_256_mont_sub_8(r->y, r->y, t2 /*, p256_mod*/); dump_512("y2 %s\n", r->y); } /* Add two Montgomery form projective points. * * r Result of addition. * p Frist point to add. * q Second point to add. */ static void sp_256_proj_point_add_8(sp_point* r, sp_point* p, sp_point* q) { sp_digit t1[2*8]; sp_digit t2[2*8]; sp_digit t3[2*8]; sp_digit t4[2*8]; sp_digit t5[2*8]; /* Ensure only the first point is the same as the result. */ if (q == r) { sp_point* a = p; p = q; q = a; } /* Check double */ sp_256_sub_8(t1, p256_mod, q->y); sp_256_norm_8(t1); if (sp_256_cmp_equal_8(p->x, q->x) && sp_256_cmp_equal_8(p->z, q->z) && (sp_256_cmp_equal_8(p->y, q->y) || sp_256_cmp_equal_8(p->y, t1)) ) { sp_256_proj_point_dbl_8(r, p); } else { sp_point tp; sp_point *v; v = r; if (p->infinity | q->infinity) { memset(&tp, 0, sizeof(tp)); v = &tp; } *r = p->infinity ? *q : *p; /* struct copy */ /* U1 = X1*Z2^2 */ sp_256_mont_sqr_8(t1, q->z /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t3, t1, q->z /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t1, t1, v->x /*, p256_mod, p256_mp_mod*/); /* U2 = X2*Z1^2 */ sp_256_mont_sqr_8(t2, v->z /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t4, t2, v->z /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t2, t2, q->x /*, p256_mod, p256_mp_mod*/); /* S1 = Y1*Z2^3 */ sp_256_mont_mul_8(t3, t3, v->y /*, p256_mod, p256_mp_mod*/); /* S2 = Y2*Z1^3 */ sp_256_mont_mul_8(t4, t4, q->y /*, p256_mod, p256_mp_mod*/); /* H = U2 - U1 */ sp_256_mont_sub_8(t2, t2, t1 /*, p256_mod*/); /* R = S2 - S1 */ sp_256_mont_sub_8(t4, t4, t3 /*, p256_mod*/); /* Z3 = H*Z1*Z2 */ sp_256_mont_mul_8(v->z, v->z, q->z /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(v->z, v->z, t2 /*, p256_mod, p256_mp_mod*/); /* X3 = R^2 - H^3 - 2*U1*H^2 */ sp_256_mont_sqr_8(v->x, t4 /*, p256_mod, p256_mp_mod*/); sp_256_mont_sqr_8(t5, t2 /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(v->y, t1, t5 /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t5, t5, t2 /*, p256_mod, p256_mp_mod*/); sp_256_mont_sub_8(v->x, v->x, t5 /*, p256_mod*/); sp_256_mont_dbl_8(t1, v->y /*, p256_mod*/); sp_256_mont_sub_8(v->x, v->x, t1 /*, p256_mod*/); /* Y3 = R*(U1*H^2 - X3) - S1*H^3 */ sp_256_mont_sub_8(v->y, v->y, v->x /*, p256_mod*/); sp_256_mont_mul_8(v->y, v->y, t4 /*, p256_mod, p256_mp_mod*/); sp_256_mont_mul_8(t5, t5, t3 /*, p256_mod, p256_mp_mod*/); sp_256_mont_sub_8(v->y, v->y, t5 /*, p256_mod*/); } } /* Multiply the point by the scalar and return the result. * If map is true then convert result to affine co-ordinates. * * r Resulting point. * g Point to multiply. * k Scalar to multiply by. * map Indicates whether to convert result to affine. */ static void sp_256_ecc_mulmod_8(sp_point* r, const sp_point* g, const sp_digit* k /*, int map*/) { enum { map = 1 }; /* we always convert result to affine coordinates */ sp_point t[3]; sp_digit n = n; /* for compiler */ int c, y; memset(t, 0, sizeof(t)); /* t[0] = {0, 0, 1} * norm */ t[0].infinity = 1; /* t[1] = {g->x, g->y, g->z} * norm */ sp_256_mod_mul_norm_8(t[1].x, g->x); sp_256_mod_mul_norm_8(t[1].y, g->y); sp_256_mod_mul_norm_8(t[1].z, g->z); /* For every bit, starting from most significant... */ k += 7; c = 256; for (;;) { if ((c & 0x1f) == 0) { if (c == 0) break; n = *k--; } y = (n >> 31); dbg("y:%d t[%d] = t[0]+t[1]\n", y, y^1); sp_256_proj_point_add_8(&t[y^1], &t[0], &t[1]); dump_512("t[0].x %s\n", t[0].x); dump_512("t[0].y %s\n", t[0].y); dump_512("t[0].z %s\n", t[0].z); dump_512("t[1].x %s\n", t[1].x); dump_512("t[1].y %s\n", t[1].y); dump_512("t[1].z %s\n", t[1].z); dbg("t[2] = t[%d]\n", y); memcpy(&t[2], &t[y], sizeof(sp_point)); dbg("t[2] *= 2\n"); sp_256_proj_point_dbl_8(&t[2], &t[2]); dump_512("t[2].x %s\n", t[2].x); dump_512("t[2].y %s\n", t[2].y); dump_512("t[2].z %s\n", t[2].z); memcpy(&t[y], &t[2], sizeof(sp_point)); n <<= 1; c--; } if (map) sp_256_map_8(r, &t[0]); else memcpy(r, &t[0], sizeof(sp_point)); memset(t, 0, sizeof(t)); //paranoia } /* Multiply the base point of P256 by the scalar and return the result. * If map is true then convert result to affine co-ordinates. * * r Resulting point. * k Scalar to multiply by. * map Indicates whether to convert result to affine. */ static void sp_256_ecc_mulmod_base_8(sp_point* r, sp_digit* k /*, int map*/) { /* Since this function is called only once, save space: * don't have "static const sp_point p256_base = {...}", * it would have more zeros than data. */ static const uint8_t p256_base_bin[] = { /* x (big-endian) */ 0x6b,0x17,0xd1,0xf2,0xe1,0x2c,0x42,0x47,0xf8,0xbc,0xe6,0xe5,0x63,0xa4,0x40,0xf2,0x77,0x03,0x7d,0x81,0x2d,0xeb,0x33,0xa0,0xf4,0xa1,0x39,0x45,0xd8,0x98,0xc2,0x96, /* y */ 0x4f,0xe3,0x42,0xe2,0xfe,0x1a,0x7f,0x9b,0x8e,0xe7,0xeb,0x4a,0x7c,0x0f,0x9e,0x16,0x2b,0xce,0x33,0x57,0x6b,0x31,0x5e,0xce,0xcb,0xb6,0x40,0x68,0x37,0xbf,0x51,0xf5, /* z will be set to 1, infinity flag to "false" */ }; sp_point p256_base; sp_256_point_from_bin2x32(&p256_base, p256_base_bin); sp_256_ecc_mulmod_8(r, &p256_base, k /*, map*/); } /* Multiply the point by the scalar and serialize the X ordinate. * The number is 0 padded to maximum size on output. * * priv Scalar to multiply the point by. * pub2x32 Point to multiply. * out32 Buffer to hold X ordinate. */ static void sp_ecc_secret_gen_256(const sp_digit priv[8], const uint8_t *pub2x32, uint8_t* out32) { sp_point point[1]; #if FIXED_PEER_PUBKEY memset((void*)pub2x32, 0x55, 64); #endif dump_hex("peerkey %s\n", pub2x32, 32); /* in TLS, this is peer's public key */ dump_hex(" %s\n", pub2x32 + 32, 32); sp_256_point_from_bin2x32(point, pub2x32); dump_512("point->x %s\n", point->x); dump_512("point->y %s\n", point->y); sp_256_ecc_mulmod_8(point, point, priv); sp_256_to_bin_8(point->x, out32); dump_hex("out32: %s\n", out32, 32); } /* Generates a random scalar in [1..order-1] range. */ static void sp_256_ecc_gen_k_8(sp_digit k[8]) { /* Since 32-bit words are "dense", no need to use * sp_256_from_bin_8(k, buf) to convert random stream * to sp_digit array - just store random bits there directly. */ tls_get_random(k, 8 * sizeof(k[0])); #if FIXED_SECRET memset(k, 0x77, 8 * sizeof(k[0])); #endif // If scalar is too large, try again (pseudo-code) // if (k >= 0xffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551 - 1) // order of P256 // goto pick_another_random; // k++; // ensure non-zero /* Simpler alternative, at the cost of not choosing some valid * random values, and slightly non-uniform distribution */ if (k[0] == 0) k[0] = 1; if (k[7] >= 0xffffffff) k[7] = 0xfffffffe; } /* Makes a random EC key pair. */ static void sp_ecc_make_key_256(sp_digit privkey[8], uint8_t *pubkey) { sp_point point[1]; sp_256_ecc_gen_k_8(privkey); dump_256("privkey %s\n", privkey); sp_256_ecc_mulmod_base_8(point, privkey); dump_512("point->x %s\n", point->x); dump_512("point->y %s\n", point->y); sp_256_to_bin_8(point->x, pubkey); sp_256_to_bin_8(point->y, pubkey + 32); memset(point, 0, sizeof(point)); //paranoia } void FAST_FUNC curve_P256_compute_pubkey_and_premaster( uint8_t *pubkey2x32, uint8_t *premaster32, const uint8_t *peerkey2x32) { sp_digit privkey[8]; dump_hex("peerkey2x32: %s\n", peerkey2x32, 64); sp_ecc_make_key_256(privkey, pubkey2x32); dump_hex("pubkey: %s\n", pubkey2x32, 32); dump_hex(" %s\n", pubkey2x32 + 32, 32); /* Combine our privkey and peer's public key to generate premaster */ sp_ecc_secret_gen_256(privkey, /*x,y:*/peerkey2x32, premaster32); dump_hex("premaster: %s\n", premaster32, 32); }