libsec: curve25519 (X25519) and AES-GCM AEAD primitives
Adds the crypto primitives libsec needs for modern TLS 1.2:
X25519 scalar multiplication for ECDHE, and AES-GCM AEAD for
the record layer. Also usable standalone for code wanting
X25519 key agreement or AES-GCM.
curve25519.c and curve25519_dh.c are ported from 9front
libsec/port/; curve25519.c is Adam Langley's curve25519-donna.
aes_gcm.c is the 9front implementation of NIST SP 800-38D.
The only change from 9front is the include header swap
("os.h" -> <u.h> + <libc.h> + <libsec.h>).
RFC 7748 §5 (Curve25519/X25519 scalar multiplication). AES-GCM
follows NIST SP 800-38D. TLS binding (record layer, cipher
suite IDs) is in the tls-aead-record-layer and
tls-ecdhe-sni-client patches.
--- sys/include/libsec.h
+++ sys/include/libsec.h
@@ -44,6 +44,20 @@
void setupAESXCBCstate(AESstate *s);
uchar* aesXCBCmac(uchar *p, int len, AESstate *s);
+/* AES-GCM AEAD (RFC 5288, NIST SP 800-38D) */
+typedef struct AESGCMstate AESGCMstate;
+struct AESGCMstate
+{
+ AESstate; /* anonymous embedded AES state */
+ ulong H[4];
+ ulong M[16][256][4]; /* precomputed GHASH tables */
+};
+
+void setupAESGCMstate(AESGCMstate *s, uchar *key, int keylen, uchar *iv, int ivlen);
+void aesgcm_setiv(AESGCMstate *s, uchar *iv, int ivlen);
+void aesgcm_encrypt(uchar *dat, ulong ndat, uchar *aad, ulong naad, uchar tag[16], AESGCMstate *s);
+int aesgcm_decrypt(uchar *dat, ulong ndat, uchar *aad, ulong naad, uchar tag[16], AESGCMstate *s);
+
/*
* Blowfish Definitions
*/
@@ -113,6 +127,14 @@
void ccpoly_encrypt(uchar *dat, ulong ndat, uchar *aad, ulong naad, uchar tag[16], Chachastate *cs);
int ccpoly_decrypt(uchar *dat, ulong ndat, uchar *aad, ulong naad, uchar tag[16], Chachastate *cs);
+
+/*
+ * Curve25519 (RFC 7748)
+ */
+void curve25519(uchar mypublic[32], uchar secret[32], uchar basepoint[32]);
+int x25519(uchar out[32], uchar s[32], uchar u[32]);
+void curve25519_dh_new(uchar x[32], uchar y[32]);
+int curve25519_dh_finish(uchar x[32], uchar y[32], uchar z[32]);
/*
* DES definitions
--- sys/src/libsec/port/aes_gcm.c
+++ sys/src/libsec/port/aes_gcm.c
@@ -0,0 +1,200 @@
+#include <u.h>
+#include <libc.h>
+#include <libsec.h>
+
+static void
+load128(uchar b[16], ulong W[4])
+{
+ W[0] = (ulong)b[15] | (ulong)b[14]<<8 | (ulong)b[13]<<16 | (ulong)b[12]<<24;
+ W[1] = (ulong)b[11] | (ulong)b[10]<<8 | (ulong)b[ 9]<<16 | (ulong)b[ 8]<<24;
+ W[2] = (ulong)b[ 7] | (ulong)b[ 6]<<8 | (ulong)b[ 5]<<16 | (ulong)b[ 4]<<24;
+ W[3] = (ulong)b[ 3] | (ulong)b[ 2]<<8 | (ulong)b[ 1]<<16 | (ulong)b[ 0]<<24;
+}
+
+static void
+store128(ulong W[4], uchar b[16])
+{
+ b[15] = W[0], b[14] = W[0]>>8, b[13] = W[0]>>16, b[12] = W[0]>>24;
+ b[11] = W[1], b[10] = W[1]>>8, b[ 9] = W[1]>>16, b[ 8] = W[1]>>24;
+ b[ 7] = W[2], b[ 6] = W[2]>>8, b[ 5] = W[2]>>16, b[ 4] = W[2]>>24;
+ b[ 3] = W[3], b[ 2] = W[3]>>8, b[ 1] = W[3]>>16, b[ 0] = W[3]>>24;
+}
+
+static void
+gfmul(ulong X[4], ulong Y[4], ulong Z[4])
+{
+ long m, i;
+
+ Z[0] = Z[1] = Z[2] = Z[3] = 0;
+ for(i=127; i>=0; i--){
+ m = ((long)Y[i>>5] << 31-(i&31)) >> 31;
+ Z[0] ^= X[0] & m;
+ Z[1] ^= X[1] & m;
+ Z[2] ^= X[2] & m;
+ Z[3] ^= X[3] & m;
+ m = ((long)X[0]<<31) >> 31;
+ X[0] = X[0]>>1 | X[1]<<31;
+ X[1] = X[1]>>1 | X[2]<<31;
+ X[2] = X[2]>>1 | X[3]<<31;
+ X[3] = X[3]>>1 ^ (0xE1000000 & m);
+ }
+}
+
+static void
+prepareM(ulong H[4], ulong M[16][256][4])
+{
+ ulong X[4], i, j;
+
+ for(i=0; i<16; i++){
+ for(j=0; j<256; j++){
+ X[0] = X[1] = X[2] = X[3] = 0;
+ X[i>>2] = j<<((i&3)<<3);
+ gfmul(X, H, M[i][j]);
+ }
+ }
+}
+
+static void
+ghash1(AESGCMstate *s, ulong X[4], ulong Y[4])
+{
+ ulong *Xi, i;
+
+ X[0] ^= Y[0], X[1] ^= Y[1], X[2] ^= Y[2], X[3] ^= Y[3];
+ if(0){
+ gfmul(X, s->H, Y);
+ return;
+ }
+
+ Y[0] = Y[1] = Y[2] = Y[3] = 0;
+ for(i=0; i<16; i++){
+ Xi = s->M[i][(X[i>>2]>>((i&3)<<3))&0xFF];
+ Y[0] ^= Xi[0];
+ Y[1] ^= Xi[1];
+ Y[2] ^= Xi[2];
+ Y[3] ^= Xi[3];
+ }
+}
+
+static void
+ghashn(AESGCMstate *s, uchar *dat, ulong len, ulong Y[4])
+{
+ uchar tmp[16];
+ ulong X[4];
+
+ while(len >= 16){
+ load128(dat, X);
+ ghash1(s, X, Y);
+ dat += 16, len -= 16;
+ }
+ if(len > 0){
+ memmove(tmp, dat, len);
+ memset(tmp+len, 0, 16-len);
+ load128(tmp, X);
+ ghash1(s, X, Y);
+ }
+}
+
+static ulong
+aesxctr1(AESstate *s, uchar ctr[AESbsize], uchar *dat, ulong len)
+{
+ uchar tmp[AESbsize];
+ ulong i;
+
+ aes_encrypt(s->ekey, s->rounds, ctr, tmp);
+ if(len > AESbsize)
+ len = AESbsize;
+ for(i=0; i<len; i++)
+ dat[i] ^= tmp[i];
+ return len;
+}
+
+static void
+aesxctrn(AESstate *s, uchar *dat, ulong len)
+{
+ uchar ctr[AESbsize];
+ ulong i;
+
+ memmove(ctr, s->ivec, AESbsize);
+ while(len > 0){
+ for(i=AESbsize-1; i>=AESbsize-4; i--)
+ if(++ctr[i] != 0)
+ break;
+
+ if(aesxctr1(s, ctr, dat, len) < AESbsize)
+ break;
+ dat += AESbsize;
+ len -= AESbsize;
+ }
+}
+
+void
+aesgcm_setiv(AESGCMstate *s, uchar *iv, int ivlen)
+{
+ if(ivlen == 96/8){
+ memmove(s->ivec, iv, ivlen);
+ memset(s->ivec+ivlen, 0, AESbsize-ivlen);
+ s->ivec[AESbsize-1] = 1;
+ } else {
+ ulong L[4], Y[4] = {0};
+
+ ghashn(s, iv, ivlen, Y);
+ L[0] = ivlen << 3;
+ L[1] = ivlen >> 29;
+ L[2] = L[3] = 0;
+ ghash1(s, L, Y);
+ store128(Y, s->ivec);
+ }
+}
+
+void
+setupAESGCMstate(AESGCMstate *s, uchar *key, int keylen, uchar *iv, int ivlen)
+{
+ setupAESstate(s, key, keylen, nil);
+
+ memset(s->ivec, 0, AESbsize);
+ aes_encrypt(s->ekey, s->rounds, s->ivec, s->ivec);
+ load128(s->ivec, s->H);
+ memset(s->ivec, 0, AESbsize);
+ prepareM(s->H, s->M);
+
+ if(iv != nil && ivlen > 0)
+ aesgcm_setiv(s, iv, ivlen);
+}
+
+void
+aesgcm_encrypt(uchar *dat, ulong ndat, uchar *aad, ulong naad, uchar tag[16], AESGCMstate *s)
+{
+ ulong L[4], Y[4] = {0};
+
+ ghashn(s, aad, naad, Y);
+ aesxctrn(s, dat, ndat);
+ ghashn(s, dat, ndat, Y);
+ L[0] = ndat << 3;
+ L[1] = ndat >> 29;
+ L[2] = naad << 3;
+ L[3] = naad >> 29;
+ ghash1(s, L, Y);
+ store128(Y, tag);
+ aesxctr1(s, s->ivec, tag, 16);
+}
+
+int
+aesgcm_decrypt(uchar *dat, ulong ndat, uchar *aad, ulong naad, uchar tag[16], AESGCMstate *s)
+{
+ ulong L[4], Y[4] = {0};
+ uchar tmp[16];
+
+ ghashn(s, aad, naad, Y);
+ ghashn(s, dat, ndat, Y);
+ L[0] = ndat << 3;
+ L[1] = ndat >> 29;
+ L[2] = naad << 3;
+ L[3] = naad >> 29;
+ ghash1(s, L, Y);
+ store128(Y, tmp);
+ aesxctr1(s, s->ivec, tmp, 16);
+ if(tsmemcmp(tag, tmp, 16) != 0)
+ return -1;
+ aesxctrn(s, dat, ndat);
+ return 0;
+}
--- sys/src/libsec/port/curve25519.c
+++ sys/src/libsec/port/curve25519.c
@@ -0,0 +1,571 @@
+/* Copyright 2008, Google Inc.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following disclaimer
+ * in the documentation and/or other materials provided with the
+ * distribution.
+ * * Neither the name of Google Inc. nor the names of its
+ * contributors may be used to endorse or promote products derived from
+ * this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+ * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * curve25519: Curve25519 elliptic curve, public key function
+ *
+ * http://code.google.com/p/curve25519-donna/
+ *
+ * Adam Langley <agl@imperialviolet.org>
+ *
+ * Derived from public domain C code by Daniel J. Bernstein <djb@cr.yp.to>
+ *
+ * More information about curve25519 can be found here
+ * http://cr.yp.to/ecdh.html
+ *
+ * djb's sample implementation of curve25519 is written in a special assembly
+ * language called qhasm and uses the floating point registers.
+ *
+ * This is, almost, a clean room reimplementation from the curve25519 paper. It
+ * uses many of the tricks described therein. Only the crecip function is taken
+ * from the sample implementation.
+ */
+#include <u.h>
+#include <libc.h>
+#include <libsec.h>
+
+typedef vlong felem;
+
+/* Sum two numbers: output += in */
+static void fsum(felem *output, felem *in) {
+ unsigned i;
+ for (i = 0; i < 10; i += 2) {
+ output[0+i] = (output[0+i] + in[0+i]);
+ output[1+i] = (output[1+i] + in[1+i]);
+ }
+}
+
+/* Find the difference of two numbers: output = in - output
+ * (note the order of the arguments!)
+ */
+static void fdifference(felem *output, felem *in) {
+ unsigned i;
+ for (i = 0; i < 10; ++i) {
+ output[i] = (in[i] - output[i]);
+ }
+}
+
+/* Multiply a number my a scalar: output = in * scalar */
+static void fscalar_product(felem *output, felem *in, felem scalar) {
+ unsigned i;
+ for (i = 0; i < 10; ++i) {
+ output[i] = in[i] * scalar;
+ }
+}
+
+/* Multiply two numbers: output = in2 * in
+ *
+ * output must be distinct to both inputs. The inputs are reduced coefficient
+ * form, the output is not.
+ */
+static void fproduct(felem *output, felem *in2, felem *in) {
+ output[0] = in2[0] * in[0];
+ output[1] = in2[0] * in[1] +
+ in2[1] * in[0];
+ output[2] = 2 * in2[1] * in[1] +
+ in2[0] * in[2] +
+ in2[2] * in[0];
+ output[3] = in2[1] * in[2] +
+ in2[2] * in[1] +
+ in2[0] * in[3] +
+ in2[3] * in[0];
+ output[4] = in2[2] * in[2] +
+ 2 * (in2[1] * in[3] +
+ in2[3] * in[1]) +
+ in2[0] * in[4] +
+ in2[4] * in[0];
+ output[5] = in2[2] * in[3] +
+ in2[3] * in[2] +
+ in2[1] * in[4] +
+ in2[4] * in[1] +
+ in2[0] * in[5] +
+ in2[5] * in[0];
+ output[6] = 2 * (in2[3] * in[3] +
+ in2[1] * in[5] +
+ in2[5] * in[1]) +
+ in2[2] * in[4] +
+ in2[4] * in[2] +
+ in2[0] * in[6] +
+ in2[6] * in[0];
+ output[7] = in2[3] * in[4] +
+ in2[4] * in[3] +
+ in2[2] * in[5] +
+ in2[5] * in[2] +
+ in2[1] * in[6] +
+ in2[6] * in[1] +
+ in2[0] * in[7] +
+ in2[7] * in[0];
+ output[8] = in2[4] * in[4] +
+ 2 * (in2[3] * in[5] +
+ in2[5] * in[3] +
+ in2[1] * in[7] +
+ in2[7] * in[1]) +
+ in2[2] * in[6] +
+ in2[6] * in[2] +
+ in2[0] * in[8] +
+ in2[8] * in[0];
+ output[9] = in2[4] * in[5] +
+ in2[5] * in[4] +
+ in2[3] * in[6] +
+ in2[6] * in[3] +
+ in2[2] * in[7] +
+ in2[7] * in[2] +
+ in2[1] * in[8] +
+ in2[8] * in[1] +
+ in2[0] * in[9] +
+ in2[9] * in[0];
+ output[10] = 2 * (in2[5] * in[5] +
+ in2[3] * in[7] +
+ in2[7] * in[3] +
+ in2[1] * in[9] +
+ in2[9] * in[1]) +
+ in2[4] * in[6] +
+ in2[6] * in[4] +
+ in2[2] * in[8] +
+ in2[8] * in[2];
+ output[11] = in2[5] * in[6] +
+ in2[6] * in[5] +
+ in2[4] * in[7] +
+ in2[7] * in[4] +
+ in2[3] * in[8] +
+ in2[8] * in[3] +
+ in2[2] * in[9] +
+ in2[9] * in[2];
+ output[12] = in2[6] * in[6] +
+ 2 * (in2[5] * in[7] +
+ in2[7] * in[5] +
+ in2[3] * in[9] +
+ in2[9] * in[3]) +
+ in2[4] * in[8] +
+ in2[8] * in[4];
+ output[13] = in2[6] * in[7] +
+ in2[7] * in[6] +
+ in2[5] * in[8] +
+ in2[8] * in[5] +
+ in2[4] * in[9] +
+ in2[9] * in[4];
+ output[14] = 2 * (in2[7] * in[7] +
+ in2[5] * in[9] +
+ in2[9] * in[5]) +
+ in2[6] * in[8] +
+ in2[8] * in[6];
+ output[15] = in2[7] * in[8] +
+ in2[8] * in[7] +
+ in2[6] * in[9] +
+ in2[9] * in[6];
+ output[16] = in2[8] * in[8] +
+ 2 * (in2[7] * in[9] +
+ in2[9] * in[7]);
+ output[17] = in2[8] * in[9] +
+ in2[9] * in[8];
+ output[18] = 2 * in2[9] * in[9];
+}
+
+/* Reduce a long form to a short form by taking the input mod 2^255 - 19. */
+static void freduce_degree(felem *output) {
+ output[8] += 19 * output[18];
+ output[7] += 19 * output[17];
+ output[6] += 19 * output[16];
+ output[5] += 19 * output[15];
+ output[4] += 19 * output[14];
+ output[3] += 19 * output[13];
+ output[2] += 19 * output[12];
+ output[1] += 19 * output[11];
+ output[0] += 19 * output[10];
+}
+
+/* Reduce all coefficients of the short form input to be -2**25 <= x <= 2**25
+ */
+static void freduce_coefficients(felem *output) {
+ unsigned i;
+ do {
+ output[10] = 0;
+
+ for (i = 0; i < 10; i += 2) {
+ felem over = output[i] / 0x2000000l;
+ felem over2 = (over + ((over >> 63) * 2) + 1) / 2;
+ output[i+1] += over2;
+ output[i] -= over2 * 0x4000000l;
+
+ over = output[i+1] / 0x2000000;
+ output[i+2] += over;
+ output[i+1] -= over * 0x2000000;
+ }
+ output[0] += 19 * output[10];
+ } while (output[10]);
+}
+
+/* A helpful wrapper around fproduct: output = in * in2.
+ *
+ * output must be distinct to both inputs. The output is reduced degree and
+ * reduced coefficient.
+ */
+static void
+fmult(felem *output, felem *in, felem *in2) {
+ felem t[19];
+ fproduct(t, in, in2);
+ freduce_degree(t);
+ freduce_coefficients(t);
+ memcpy(output, t, sizeof(felem) * 10);
+}
+
+static void fsquare_inner(felem *output, felem *in) {
+ felem tmp;
+ output[0] = in[0] * in[0];
+ output[1] = 2 * in[0] * in[1];
+ output[2] = 2 * (in[1] * in[1] +
+ in[0] * in[2]);
+ output[3] = 2 * (in[1] * in[2] +
+ in[0] * in[3]);
+ output[4] = in[2] * in[2] +
+ 4 * in[1] * in[3] +
+ 2 * in[0] * in[4];
+ output[5] = 2 * (in[2] * in[3] +
+ in[1] * in[4] +
+ in[0] * in[5]);
+ output[6] = 2 * (in[3] * in[3] +
+ in[2] * in[4] +
+ in[0] * in[6] +
+ 2 * in[1] * in[5]);
+ output[7] = 2 * (in[3] * in[4] +
+ in[2] * in[5] +
+ in[1] * in[6] +
+ in[0] * in[7]);
+ tmp = in[1] * in[7] + in[3] * in[5];
+ output[8] = in[4] * in[4] +
+ 2 * (in[2] * in[6] +
+ in[0] * in[8] +
+ 2 * tmp);
+ output[9] = 2 * (in[4] * in[5] +
+ in[3] * in[6] +
+ in[2] * in[7] +
+ in[1] * in[8] +
+ in[0] * in[9]);
+ tmp = in[3] * in[7] + in[1] * in[9];
+ output[10] = 2 * (in[5] * in[5] +
+ in[4] * in[6] +
+ in[2] * in[8] +
+ 2 * tmp);
+ output[11] = 2 * (in[5] * in[6] +
+ in[4] * in[7] +
+ in[3] * in[8] +
+ in[2] * in[9]);
+ output[12] = in[6] * in[6] +
+ 2 * (in[4] * in[8] +
+ 2 * (in[5] * in[7] +
+ in[3] * in[9]));
+ output[13] = 2 * (in[6] * in[7] +
+ in[5] * in[8] +
+ in[4] * in[9]);
+ output[14] = 2 * (in[7] * in[7] +
+ in[6] * in[8] +
+ 2 * in[5] * in[9]);
+ output[15] = 2 * (in[7] * in[8] +
+ in[6] * in[9]);
+ output[16] = in[8] * in[8] +
+ 4 * in[7] * in[9];
+ output[17] = 2 * in[8] * in[9];
+ output[18] = 2 * in[9] * in[9];
+}
+
+static void
+fsquare(felem *output, felem *in) {
+ felem t[19];
+ fsquare_inner(t, in);
+ freduce_degree(t);
+ freduce_coefficients(t);
+ memcpy(output, t, sizeof(felem) * 10);
+}
+
+/* Take a little-endian, 32-byte number and expand it into polynomial form */
+static void
+fexpand(felem *output, uchar *input) {
+#define F(n,start,shift,mask) \
+ output[n] = ((((felem) input[start + 0]) | \
+ ((felem) input[start + 1]) << 8 | \
+ ((felem) input[start + 2]) << 16 | \
+ ((felem) input[start + 3]) << 24) >> shift) & mask;
+ F(0, 0, 0, 0x3ffffff);
+ F(1, 3, 2, 0x1ffffff);
+ F(2, 6, 3, 0x3ffffff);
+ F(3, 9, 5, 0x1ffffff);
+ F(4, 12, 6, 0x3ffffff);
+ F(5, 16, 0, 0x1ffffff);
+ F(6, 19, 1, 0x3ffffff);
+ F(7, 22, 3, 0x1ffffff);
+ F(8, 25, 4, 0x3ffffff);
+ F(9, 28, 6, 0x1ffffff);
+#undef F
+}
+
+/* Take a fully reduced polynomial form number and contract it into a
+ * little-endian, 32-byte array
+ */
+static void
+fcontract(uchar *output, felem *input) {
+ int i;
+
+ do {
+ for (i = 0; i < 9; ++i) {
+ if ((i & 1) == 1) {
+ while (input[i] < 0) {
+ input[i] += 0x2000000;
+ input[i + 1]--;
+ }
+ } else {
+ while (input[i] < 0) {
+ input[i] += 0x4000000;
+ input[i + 1]--;
+ }
+ }
+ }
+ while (input[9] < 0) {
+ input[9] += 0x2000000;
+ input[0] -= 19;
+ }
+ } while (input[0] < 0);
+
+ input[1] <<= 2;
+ input[2] <<= 3;
+ input[3] <<= 5;
+ input[4] <<= 6;
+ input[6] <<= 1;
+ input[7] <<= 3;
+ input[8] <<= 4;
+ input[9] <<= 6;
+#define F(i, s) \
+ output[s+0] |= input[i] & 0xff; \
+ output[s+1] = (input[i] >> 8) & 0xff; \
+ output[s+2] = (input[i] >> 16) & 0xff; \
+ output[s+3] = (input[i] >> 24) & 0xff;
+ output[0] = 0;
+ output[16] = 0;
+ F(0,0);
+ F(1,3);
+ F(2,6);
+ F(3,9);
+ F(4,12);
+ F(5,16);
+ F(6,19);
+ F(7,22);
+ F(8,25);
+ F(9,28);
+#undef F
+}
+
+/* Input: Q, Q', Q-Q'
+ * Output: 2Q, Q+Q'
+ *
+ * x2 z3: long form
+ * x3 z3: long form
+ * x z: short form, destroyed
+ * xprime zprime: short form, destroyed
+ * qmqp: short form, preserved
+ */
+static void fmonty(felem *x2, felem *z2, /* output 2Q */
+ felem *x3, felem *z3, /* output Q + Q' */
+ felem *x, felem *z, /* input Q */
+ felem *xprime, felem *zprime, /* input Q' */
+ felem *qmqp /* input Q - Q' */) {
+ felem origx[10], origxprime[10], zzz[19], xx[19], zz[19], xxprime[19],
+ zzprime[19], zzzprime[19], xxxprime[19];
+
+ memcpy(origx, x, 10 * sizeof(felem));
+ fsum(x, z);
+ fdifference(z, origx); /* does x - z */
+
+ memcpy(origxprime, xprime, sizeof(felem) * 10);
+ fsum(xprime, zprime);
+ fdifference(zprime, origxprime);
+ fproduct(xxprime, xprime, z);
+ fproduct(zzprime, x, zprime);
+ freduce_degree(xxprime);
+ freduce_coefficients(xxprime);
+ freduce_degree(zzprime);
+ freduce_coefficients(zzprime);
+ memcpy(origxprime, xxprime, sizeof(felem) * 10);
+ fsum(xxprime, zzprime);
+ fdifference(zzprime, origxprime);
+ fsquare(xxxprime, xxprime);
+ fsquare(zzzprime, zzprime);
+ fproduct(zzprime, zzzprime, qmqp);
+ freduce_degree(zzprime);
+ freduce_coefficients(zzprime);
+ memcpy(x3, xxxprime, sizeof(felem) * 10);
+ memcpy(z3, zzprime, sizeof(felem) * 10);
+
+ fsquare(xx, x);
+ fsquare(zz, z);
+ fproduct(x2, xx, zz);
+ freduce_degree(x2);
+ freduce_coefficients(x2);
+ fdifference(zz, xx); /* does zz = xx - zz */
+ memset(zzz + 10, 0, sizeof(felem) * 9);
+ fscalar_product(zzz, zz, 121665);
+ freduce_degree(zzz);
+ freduce_coefficients(zzz);
+ fsum(zzz, xx);
+ fproduct(z2, zz, zzz);
+ freduce_degree(z2);
+ freduce_coefficients(z2);
+}
+
+/* Calculates nQ where Q is the x-coordinate of a point on the curve
+ *
+ * resultx/resultz: the x coordinate of the resulting curve point (short form)
+ * n: a little endian, 32-byte number
+ * q: a point of the curve (short form)
+ */
+static void
+cmult(felem *resultx, felem *resultz, uchar *n, felem *q) {
+ felem a[19] = {0}, b[19] = {1}, c[19] = {1}, d[19] = {0};
+ felem *nqpqx = a, *nqpqz = b, *nqx = c, *nqz = d, *t;
+ felem e[19] = {0}, f[19] = {1}, g[19] = {0}, h[19] = {1};
+ felem *nqpqx2 = e, *nqpqz2 = f, *nqx2 = g, *nqz2 = h;
+
+ unsigned i, j;
+
+ memcpy(nqpqx, q, sizeof(felem) * 10);
+
+ for (i = 0; i < 32; ++i) {
+ uchar byte = n[31 - i];
+ for (j = 0; j < 8; ++j) {
+ if (byte & 0x80) {
+ fmonty(nqpqx2, nqpqz2,
+ nqx2, nqz2,
+ nqpqx, nqpqz,
+ nqx, nqz,
+ q);
+ } else {
+ fmonty(nqx2, nqz2,
+ nqpqx2, nqpqz2,
+ nqx, nqz,
+ nqpqx, nqpqz,
+ q);
+ }
+
+ t = nqx;
+ nqx = nqx2;
+ nqx2 = t;
+ t = nqz;
+ nqz = nqz2;
+ nqz2 = t;
+ t = nqpqx;
+ nqpqx = nqpqx2;
+ nqpqx2 = t;
+ t = nqpqz;
+ nqpqz = nqpqz2;
+ nqpqz2 = t;
+
+ byte <<= 1;
+ }
+ }
+
+ memcpy(resultx, nqx, sizeof(felem) * 10);
+ memcpy(resultz, nqz, sizeof(felem) * 10);
+}
+
+/* ----------------------------------------------------------------------------- */
+/* Shamelessly copied from djb's code */
+/* ----------------------------------------------------------------------------- */
+static void
+crecip(felem *out, felem *z) {
+ felem z2[10];
+ felem z9[10];
+ felem z11[10];
+ felem z2_5_0[10];
+ felem z2_10_0[10];
+ felem z2_20_0[10];
+ felem z2_50_0[10];
+ felem z2_100_0[10];
+ felem t0[10];
+ felem t1[10];
+ int i;
+
+ /* 2 */ fsquare(z2,z);
+ /* 4 */ fsquare(t1,z2);
+ /* 8 */ fsquare(t0,t1);
+ /* 9 */ fmult(z9,t0,z);
+ /* 11 */ fmult(z11,z9,z2);
+ /* 22 */ fsquare(t0,z11);
+ /* 2^5 - 2^0 = 31 */ fmult(z2_5_0,t0,z9);
+
+ /* 2^6 - 2^1 */ fsquare(t0,z2_5_0);
+ /* 2^7 - 2^2 */ fsquare(t1,t0);
+ /* 2^8 - 2^3 */ fsquare(t0,t1);
+ /* 2^9 - 2^4 */ fsquare(t1,t0);
+ /* 2^10 - 2^5 */ fsquare(t0,t1);
+ /* 2^10 - 2^0 */ fmult(z2_10_0,t0,z2_5_0);
+
+ /* 2^11 - 2^1 */ fsquare(t0,z2_10_0);
+ /* 2^12 - 2^2 */ fsquare(t1,t0);
+ /* 2^20 - 2^10 */ for (i = 2;i < 10;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+ /* 2^20 - 2^0 */ fmult(z2_20_0,t1,z2_10_0);
+
+ /* 2^21 - 2^1 */ fsquare(t0,z2_20_0);
+ /* 2^22 - 2^2 */ fsquare(t1,t0);
+ /* 2^40 - 2^20 */ for (i = 2;i < 20;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+ /* 2^40 - 2^0 */ fmult(t0,t1,z2_20_0);
+
+ /* 2^41 - 2^1 */ fsquare(t1,t0);
+ /* 2^42 - 2^2 */ fsquare(t0,t1);
+ /* 2^50 - 2^10 */ for (i = 2;i < 10;i += 2) { fsquare(t1,t0); fsquare(t0,t1); }
+ /* 2^50 - 2^0 */ fmult(z2_50_0,t0,z2_10_0);
+
+ /* 2^51 - 2^1 */ fsquare(t0,z2_50_0);
+ /* 2^52 - 2^2 */ fsquare(t1,t0);
+ /* 2^100 - 2^50 */ for (i = 2;i < 50;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+ /* 2^100 - 2^0 */ fmult(z2_100_0,t1,z2_50_0);
+
+ /* 2^101 - 2^1 */ fsquare(t1,z2_100_0);
+ /* 2^102 - 2^2 */ fsquare(t0,t1);
+ /* 2^200 - 2^100 */ for (i = 2;i < 100;i += 2) { fsquare(t1,t0); fsquare(t0,t1); }
+ /* 2^200 - 2^0 */ fmult(t1,t0,z2_100_0);
+
+ /* 2^201 - 2^1 */ fsquare(t0,t1);
+ /* 2^202 - 2^2 */ fsquare(t1,t0);
+ /* 2^250 - 2^50 */ for (i = 2;i < 50;i += 2) { fsquare(t0,t1); fsquare(t1,t0); }
+ /* 2^250 - 2^0 */ fmult(t0,t1,z2_50_0);
+
+ /* 2^251 - 2^1 */ fsquare(t1,t0);
+ /* 2^252 - 2^2 */ fsquare(t0,t1);
+ /* 2^253 - 2^3 */ fsquare(t1,t0);
+ /* 2^254 - 2^4 */ fsquare(t0,t1);
+ /* 2^255 - 2^5 */ fsquare(t1,t0);
+ /* 2^255 - 21 */ fmult(out,t1,z11);
+}
+
+void
+curve25519(uchar mypublic[32], uchar secret[32], uchar basepoint[32]) {
+ felem bp[10], x[10], z[10], zmone[10];
+ fexpand(bp, basepoint);
+ cmult(x, z, secret, bp);
+ crecip(zmone, z);
+ fmult(z, x, zmone);
+ fcontract(mypublic, z);
+}
--- sys/src/libsec/port/curve25519_dh.c
+++ sys/src/libsec/port/curve25519_dh.c
@@ -0,0 +1,70 @@
+#include <u.h>
+#include <libc.h>
+#include <libsec.h>
+
+static uchar nine[32] = {9};
+static uchar zero[32] = {0};
+
+int
+x25519(uchar out[32], uchar s[32], uchar u[32])
+{
+ uchar sf, sl, ul;
+
+ sf = s[0];
+ sl = s[31];
+ ul = u[31];
+
+ /* clamp */
+ s[0] &= ~7; /* clear bit 0,1,2 */
+ s[31] = 0x40 | (s[31] & 0x7f); /* set bit 254, clear bit 255 */
+
+ /*
+ Implementations MUST accept non-canonical values and process them as
+ if they had been reduced modulo the field prime. The non-canonical
+ values are 2^255 - 19 through 2^255 - 1 for X25519
+ */
+ u[31] &= 0x7f;
+
+ curve25519(out, s, u);
+
+ s[0] = sf;
+ s[31] = sl;
+ u[31] = ul;
+
+ return tsmemcmp(out, zero, 32) != 0;
+}
+
+void
+curve25519_dh_new(uchar x[32], uchar y[32])
+{
+ /* new public/private key pair */
+ uchar b;
+
+ genrandom(x, 32);
+ b = x[31];
+
+ /* don't check for zero: the scalar is never
+ zero because of clamping, and the basepoint is not the identity
+ in the prime-order subgroup(s). */
+ x25519(y, x, nine);
+
+ /* bit 255 is always 0, so make it random */
+ y[31] |= b & 0x80;
+}
+
+int
+curve25519_dh_finish(uchar x[32], uchar y[32], uchar z[32])
+{
+ int r;
+
+ /* remove the random bit */
+ y[31] &= 0x7f;
+
+ /* calculate dhx key */
+ r = x25519(z, x, y);
+
+ memset(x, 0, 32);
+ memset(y, 0, 32);
+
+ return r;
+}
--- sys/src/libsec/port/mkfile
+++ sys/src/libsec/port/mkfile
@@ -3,7 +3,8 @@
LIB=/$objtype/lib/libsec.a
CFILES = des.c desmodes.c desECB.c desCBC.c des3ECB.c des3CBC.c\
- aes.c blowfish.c chacha.c \
+ aes.c aes_gcm.c blowfish.c chacha.c \
+ curve25519.c curve25519_dh.c\
hmac.c md5.c md5block.c md4.c sha1.c sha1block.c\
sha2_64.c sha2_128.c sha2block64.c sha2block128.c\
sha1pickle.c md5pickle.c\
|
