useful files for testing HTTP and SMTP protocols

tests/Makefile

@@ -0,0 +1,21 @@
+TARGETS=client smtp
+
+all: $(TARGETS)
+
+.c.o:  ptwist.h slitheen.h crypto.h
+	gcc -g -ggdb -c $< -o $@ -I../sslout/include
+
+slitheen.o crypto.o ptwist168.o :: ptwist.h slitheen.h crypto.h
+
+client: slitheen.o ptwist168.o crypto.o crypto.h ptwist.h slitheen.h
+	gcc -g -ggdb -o $@ $^ client.c -I../sslout/include -L../sslout/lib -lssl -lcrypto -ldl -lpthread
+
+smtp: slitheen.o ptwist168.o crypto.o crypto.h ptwist.h slitheen.h
+	gcc -g -ggdb -o $@ $^ smtpClient.c -I../sslout/include -L../sslout/lib -lssl -lcrypto -ldl -lpthread
+
+clean:
+	-rm *.o
+
+veryclean: clean
+	-rm $(TARGETS)
+

tests/client.c

@@ -0,0 +1,285 @@
+/*
+ * A "simple" SSLeay 0.8.0 demo program
+ *
+ * This program implements a simple SSL v2 or v3 client
+ * which connects to a web server and issues a "HEAD / HTTP/1.0"
+ * command and prints the output.
+ *
+ * Written by Emil Sit <sit@mit.edu>
+ */
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <string.h>
+#include <fcntl.h>
+
+#include <netdb.h>
+#include <netinet/in.h>
+
+#include <sys/types.h>
+#include <sys/socket.h>
+#include <unistd.h>
+
+#include <openssl/ssl.h>
+#include <openssl/err.h>
+
+#include "slitheen.h"
+
+int my_connect( char *, int );
+void apps_ssl_info_callback( SSL *s, int where, int ret );
+
+int my_dumb_callback( int ok, X509_STORE_CTX *ctx ) {
+    return 1;
+}
+
+int main( int argc, char **argv ) {
+    SSL_CTX *ctx = NULL;
+    SSL *session = NULL;
+
+    char *command = "HEAD / HTTP/1.0\r\n\r\n";
+    
+    int s;
+    int status;
+
+    /* We first need to establish what sort of
+     * connection we know how to make. We can use one of
+     * SSLv23_client_method(), SSLv2_client_method() and
+     * SSLv3_client_method().
+     */
+    const SSL_METHOD *meth = TLSv1_2_client_method();
+    if (meth == NULL) {
+	fprintf( stderr, "no method. :(\n" ); exit(1);
+    }
+
+    /* This enables all ciphers in SSLeay, these include:
+     *   DES, RC4, IDEA, RC2, Blowfish,
+     *   MD2, SHA, DSA.
+     * See crypto/c_all.c
+     */
+    SSL_load_error_strings();	
+    OpenSSL_add_ssl_algorithms();
+
+    /* Initialize the context. This is shared between SSL sessions
+     * and can do FH caching.
+     */
+    ctx = SSL_CTX_new( meth );
+    if ( ctx == NULL ) { fprintf( stderr, "no context. :(\n" ); exit(1);
+		ERR_print_errors_fp(stderr);
+	}
+
+    SSL_CTX_set_cipher_list(ctx, "ECDHE-RSA-AES256-GCM-SHA384:DHE-RSA-AES256-GCM-SHA384");
+
+    /* Set slitheen callbacks */
+    slitheen_init();
+
+    SSL_CTX_set_client_hello_callback(ctx, slitheen_tag_hello);
+    SSL_CTX_set_generate_ec_key_callback(ctx, slitheen_ec_seed_from_tag);
+    SSL_CTX_set_generate_key_callback(ctx, slitheen_seed_from_tag);
+    SSL_CTX_set_finished_mac_callback(ctx, slitheen_finished_mac);
+
+    /* Set up a callback for each state change so we can see what's
+     * going on */
+    SSL_CTX_set_info_callback(ctx,apps_ssl_info_callback);
+
+    /* Set it up so tha we will connect to *any* site, regardless
+     * of their certificate. */
+    SSL_CTX_set_verify( ctx, SSL_VERIFY_NONE, my_dumb_callback );
+
+    /* MACRO. Set's CTX options. Not sure. I think this enables bug
+     * support hacks. */
+    SSL_CTX_set_options(ctx,SSL_OP_ALL);
+
+    /* Finally, we're all set so we can set up the session holder */
+    session = SSL_new( ctx );
+    if ( session == NULL ) { fprintf( stderr, "no session. :(\n" ); exit(1);}
+    
+    /* Make connection s.t. s is the appropriate fd */
+    s = my_connect( (argc == 2) ? argv[1] : "scspc430.cs.uwaterloo.ca" , 8888 );
+
+    /* Set up the SSL side of the connection */
+    SSL_set_fd( session, s );
+    status = SSL_connect( session );
+    /* Check the results. */
+    switch (SSL_get_error(session,status)) {
+		case SSL_ERROR_NONE:
+		/* Everything worked :-) */
+		break;
+		case SSL_ERROR_SSL:
+		fprintf( stderr, "ssl handshake failure\n" );
+		ERR_print_errors_fp(stderr);
+		goto byebye;
+		break;
+
+		/* These are for NON-BLOCKING I/O only! */
+		case SSL_ERROR_WANT_READ:
+		case SSL_ERROR_WANT_WRITE:
+		fprintf( stderr, "want read/write. Use blocking?\n" );
+		goto byebye;	break;
+		case SSL_ERROR_WANT_CONNECT:
+		fprintf( stderr, "want connect. sleep a while, try again." );
+		goto byebye;    break;
+		
+		case SSL_ERROR_SYSCALL:
+		perror("SSL_connect");
+		goto byebye;    break;
+		case SSL_ERROR_WANT_X509_LOOKUP:
+		/* not used! */
+		fprintf( stderr, "shouldn't be getting this.\n" );
+		break;
+		case SSL_ERROR_ZERO_RETURN:
+		fprintf( stderr, "connection closed.\n" );
+		goto byebye;
+    }
+
+	/*Resume session*/
+	printf("Resuming session\n");
+	
+	SSL_SESSION *sess = SSL_get1_session(session);
+	SSL_shutdown(session);
+	SSL_free(session);
+	close(s);
+
+    s = my_connect( (argc == 2) ? argv[1] : "scspc430.cs.uwaterloo.ca" , 8888 );
+	session = SSL_new(ctx);
+
+    SSL_set_fd( session, s );
+
+	SSL_set_session(session, sess);
+
+    status = SSL_connect( session );
+
+    switch (SSL_get_error(session,status)) {
+		case SSL_ERROR_NONE:
+		/* Everything worked :-) */
+		break;
+		case SSL_ERROR_SSL:
+		fprintf( stderr, "ssl handshake failure\n" );
+		ERR_print_errors_fp(stderr);
+		goto byebye;
+		break;
+
+		/* These are for NON-BLOCKING I/O only! */
+		case SSL_ERROR_WANT_READ:
+		case SSL_ERROR_WANT_WRITE:
+		fprintf( stderr, "want read/write. Use blocking?\n" );
+		goto byebye;	break;
+		case SSL_ERROR_WANT_CONNECT:
+		fprintf( stderr, "want connect. sleep a while, try again." );
+		goto byebye;    break;
+		
+		case SSL_ERROR_SYSCALL:
+		perror("SSL_connect");
+		goto byebye;    break;
+		case SSL_ERROR_WANT_X509_LOOKUP:
+		/* not used! */
+		fprintf( stderr, "shouldn't be getting this.\n" );
+		break;
+		case SSL_ERROR_ZERO_RETURN:
+		fprintf( stderr, "connection closed.\n" );
+		goto byebye;
+    }
+
+byebye:
+
+    /* close everything down */
+    SSL_shutdown(session);
+    close(s);
+    
+    SSL_free( session ); session = NULL;
+    SSL_CTX_free(ctx);
+    return 0;
+}
+
+/* returns a socket connected to host on port port */
+int my_connect( char *hostname, int port ) {
+    struct hostent *remote_host;
+    char local_hostname[256];
+    struct sockaddr_in address;
+    int s;
+    
+    /* Get a socket to work with.  This socket will be in the Internet domain, and */
+    /* will be a stream socket. */
+    if ((s = socket(AF_INET, SOCK_STREAM, 0)) < 0) {
+		perror( "connect: cannot create socket" );
+		exit(1);
+    }
+
+    /* If hostname is NULL or of zero length, look up the local hostname */
+    if ((hostname==NULL)||(hostname[0]=='\0'))
+    {
+	/* Get the local hostname */
+	if (gethostname(local_hostname, sizeof(local_hostname))==-1)
+	{
+	    perror("connect");
+	    exit(1);
+	}
+
+	/* Look up the remote host (local host) to get its network number */
+	if ((remote_host=gethostbyname(local_hostname)) == NULL)
+	{
+	    perror("connect");
+	    exit(1);
+	}
+    }
+    else
+	/* Look up the remote host to get its network number. */
+	if ((remote_host=gethostbyname(hostname)) == NULL)
+	{
+	    perror("connect");
+	    exit(1);
+	}
+
+    /* Initialize the address varaible, which specifies where
+       connect() should attempt to connect. */
+    bcopy(remote_host->h_addr, &address.sin_addr, remote_host->h_length);
+    address.sin_family = AF_INET;
+    address.sin_port = htons(port);
+
+    if (!(connect(s, (struct sockaddr *)(&address), sizeof(address)) >= 0))
+    {
+	perror("connect");
+	exit(1);
+    }
+
+    /* Set the socket to non-blocking mode */
+    /* fcntl(s, F_SETFL, O_NONBLOCK); */
+
+    return(s);
+}
+
+void apps_ssl_info_callback( SSL *s, int where, int ret )
+{
+    char *str;
+    int w;
+    
+    w=where& ~SSL_ST_MASK;
+    
+    if (w & SSL_ST_CONNECT) str="SSL_connect";
+    else if (w & SSL_ST_ACCEPT) str="SSL_accept";
+    else str="undefined";
+    
+    if (where & SSL_CB_LOOP)
+    {
+	fprintf(stderr,"%s: %s\n",str,SSL_state_string_long(s));
+    }
+    else if (where & SSL_CB_ALERT)
+    {
+	str=(where & SSL_CB_READ)?"read":"write";
+	fprintf(stderr,"SSL3 alert %s:%s:%s\n",
+		   str,
+		   SSL_alert_type_string_long(ret),
+		   SSL_alert_desc_string_long(ret));
+    }
+    else if (where & SSL_CB_EXIT)
+    {
+	if (ret == 0)
+	    fprintf(stderr,"%s:failed in %s\n",
+		       str,SSL_state_string_long(s));
+	else if (ret < 0)
+	{
+	    fprintf(stderr,"%s:error in %s\n",
+		       str,SSL_state_string_long(s));
+	}
+    }
+}
+

tests/crypto.c

@@ -0,0 +1,216 @@
+
+/**
+ * Author: Cecylia Bocovich <cbocovic@uwaterloo.ca>
+ *
+ * This file contains cryptographic helper functions to
+ * tag flows for use with the Slitheen decoy routing system
+ * Some code in this document is based on the OpenSSL source files:
+ *      crypto/ec/ec_key.c
+ *      crypto/dh/dh_key.c
+ *
+ */
+/*
+ * Written by Nils Larsch for the OpenSSL project.
+ */
+/* ====================================================================
+ * Copyright (c) 1998-2005 The OpenSSL Project.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ *
+ * 2. 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.
+ *
+ * 3. All advertising materials mentioning features or use of this
+ *    software must display the following acknowledgment:
+ *    "This product includes software developed by the OpenSSL Project
+ *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
+ *
+ * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
+ *    endorse or promote products derived from this software without
+ *    prior written permission. For written permission, please contact
+ *    openssl-core@openssl.org.
+ *
+ * 5. Products derived from this software may not be called "OpenSSL"
+ *    nor may "OpenSSL" appear in their names without prior written
+ *    permission of the OpenSSL Project.
+ *
+ * 6. Redistributions of any form whatsoever must retain the following
+ *    acknowledgment:
+ *    "This product includes software developed by the OpenSSL Project
+ *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
+ * EXPRESSED 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 OpenSSL PROJECT OR
+ * ITS 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.
+ * ====================================================================
+ *
+ * This product includes cryptographic software written by Eric Young
+ * (eay@cryptsoft.com).  This product includes software written by Tim
+ * Hudson (tjh@cryptsoft.com).
+ *
+ */
+/* ====================================================================
+ * Copyright 2002 Sun Microsystems, Inc. ALL RIGHTS RESERVED.
+ * Portions originally developed by SUN MICROSYSTEMS, INC., and
+ * contributed to the OpenSSL project.
+ */
+
+/* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
+ * All rights reserved.
+ *
+ * This package is an SSL implementation written
+ * by Eric Young (eay@cryptsoft.com).
+ * The implementation was written so as to conform with Netscapes SSL.
+ *
+ * This library is free for commercial and non-commercial use as long as
+ * the following conditions are aheared to.  The following conditions
+ * apply to all code found in this distribution, be it the RC4, RSA,
+ * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
+ * included with this distribution is covered by the same copyright terms
+ * except that the holder is Tim Hudson (tjh@cryptsoft.com).
+ *
+ * Copyright remains Eric Young's, and as such any Copyright notices in
+ * the code are not to be removed.
+ * If this package is used in a product, Eric Young should be given attribution
+ * as the author of the parts of the library used.
+ * This can be in the form of a textual message at program startup or
+ * in documentation (online or textual) provided with the package.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. 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.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *    "This product includes cryptographic software written by
+ *     Eric Young (eay@cryptsoft.com)"
+ *    The word 'cryptographic' can be left out if the rouines from the library
+ *    being used are not cryptographic related :-).
+ * 4. If you include any Windows specific code (or a derivative thereof) from
+ *    the apps directory (application code) you must include an acknowledgement:
+ *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
+ *
+ * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``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 AUTHOR 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.
+ *
+ * The licence and distribution terms for any publically available version or
+ * derivative of this code cannot be changed.  i.e. this code cannot simply be
+ * copied and put under another distribution licence
+ * [including the GNU Public Licence.]
+ */
+
+
+#include <openssl/evp.h>
+#include <openssl/dh.h>
+#include <openssl/bn.h>
+#include <openssl/err.h>
+#include <openssl/rand.h>
+#include <openssl/ssl.h>
+#include "crypto.h"
+
+/* PRF using sha384, as defined in RFC 5246 */
+int PRF(uint8_t *secret, int32_t secret_len,
+		uint8_t *seed1, int32_t seed1_len,
+		uint8_t *seed2, int32_t seed2_len,
+		uint8_t *seed3, int32_t seed3_len,
+		uint8_t *seed4, int32_t seed4_len,
+		uint8_t *output, int32_t output_len){
+
+	EVP_MD_CTX ctx, ctx_tmp, ctx_init;
+	EVP_PKEY *mac_key;
+	const EVP_MD *md = EVP_sha384();
+
+	uint8_t A[EVP_MAX_MD_SIZE];
+	size_t len, A_len;
+	int chunk = EVP_MD_size(md);
+	int remaining = output_len;
+
+	uint8_t *out = output;
+
+	EVP_MD_CTX_init(&ctx);
+	EVP_MD_CTX_init(&ctx_tmp);
+	EVP_MD_CTX_init(&ctx_init);
+	EVP_MD_CTX_set_flags(&ctx_init, EVP_MD_CTX_FLAG_NON_FIPS_ALLOW);
+
+	mac_key = EVP_PKEY_new_mac_key(EVP_PKEY_HMAC, NULL, secret, secret_len);
+
+	/* Calculate first A value */
+	EVP_DigestSignInit(&ctx_init, NULL, md, NULL, mac_key);
+	EVP_MD_CTX_copy_ex(&ctx, &ctx_init);
+	if(seed1 != NULL && seed1_len > 0){
+		EVP_DigestSignUpdate(&ctx, seed1, seed1_len);
+	}
+	if(seed2 != NULL && seed2_len > 0){
+		EVP_DigestSignUpdate(&ctx, seed2, seed2_len);
+	}
+	if(seed3 != NULL && seed3_len > 0){
+		EVP_DigestSignUpdate(&ctx, seed3, seed3_len);
+	}
+	if(seed4 != NULL && seed4_len > 0){
+		EVP_DigestSignUpdate(&ctx, seed4, seed4_len);
+	}
+	EVP_DigestSignFinal(&ctx, A, &A_len);
+
+	//iterate until desired length is achieved
+	while(remaining > 0){
+		/* Now compute SHA384(secret, A+seed) */
+		EVP_MD_CTX_copy_ex(&ctx, &ctx_init);
+		EVP_DigestSignUpdate(&ctx, A, A_len);
+		EVP_MD_CTX_copy_ex(&ctx_tmp, &ctx);
+		if(seed1 != NULL && seed1_len > 0){
+			EVP_DigestSignUpdate(&ctx, seed1, seed1_len);
+		}
+		if(seed2 != NULL && seed2_len > 0){
+			EVP_DigestSignUpdate(&ctx, seed2, seed2_len);
+		}
+		if(seed3 != NULL && seed3_len > 0){
+			EVP_DigestSignUpdate(&ctx, seed3, seed3_len);
+		}
+		if(seed4 != NULL && seed4_len > 0){
+			EVP_DigestSignUpdate(&ctx, seed4, seed4_len);
+		}
+		
+		if(remaining > chunk){
+			EVP_DigestSignFinal(&ctx, out, &len);
+			out += len;
+			remaining -= len;
+
+			/* Next A value */
+			EVP_DigestSignFinal(&ctx_tmp, A, &A_len);
+		} else {
+			EVP_DigestSignFinal(&ctx, A, &A_len);
+			memcpy(out, A, remaining);
+			remaining -= remaining;
+		}
+	}
+	return 1;
+}

tests/crypto.h

@@ -0,0 +1,21 @@
+#ifndef _CRYPTO_H_
+#define _CRYPTO_H_
+
+#include <stdint.h>
+
+# define n2s(c,s)        ((s=(((unsigned int)(c[0]))<< 8)| \
+							(((unsigned int)(c[1]))    )),c+=2)
+
+
+int PRF(uint8_t *secret, int32_t secret_len,
+		uint8_t *seed1, int32_t seed1_len,
+		uint8_t *seed2, int32_t seed2_len,
+		uint8_t *seed3, int32_t seed3_len,
+		uint8_t *seed4, int32_t seed4_len,
+		uint8_t *output, int32_t output_len);
+
+#define PRE_MASTER_LEN 256
+
+#define SSL_MAX_DIGEST 6 //this is from ssl_locl.h
+
+#endif /* _CRYPTO_H_ */

tests/domain.crt

@@ -0,0 +1,22 @@
+-----BEGIN CERTIFICATE-----
+MIIDlTCCAn2gAwIBAgIJAOegHMxTnxsfMA0GCSqGSIb3DQEBCwUAMGExCzAJBgNV
+BAYTAkNBMQswCQYDVQQIDAJPTjELMAkGA1UEBwwCTkExCzAJBgNVBAoMAk5BMQsw
+CQYDVQQLDAJOQTELMAkGA1UEAwwCTkExETAPBgkqhkiG9w0BCQEWAk5BMB4XDTE2
+MTEwNDE0NDEwOFoXDTE3MTEwNDE0NDEwOFowYTELMAkGA1UEBhMCQ0ExCzAJBgNV
+BAgMAk9OMQswCQYDVQQHDAJOQTELMAkGA1UECgwCTkExCzAJBgNVBAsMAk5BMQsw
+CQYDVQQDDAJOQTERMA8GCSqGSIb3DQEJARYCTkEwggEiMA0GCSqGSIb3DQEBAQUA
+A4IBDwAwggEKAoIBAQDSBOGs0Ckg6G/1VRHKIpQrLUfvsrUKv8SpW4RxIpmop67P
+oy/dfRAU+6YfVE1BdXZnyjGHmt3pJ57PilN6Vm6CQbbRB5cfLlp/GexkNVmptRV1
+LgLGzCqTpkemFtGnV55eVMrEPvy291lUID6As/evYO75AgfutoOPFBBqsPUYOW19
+fyntG/5LedqkVSJa2nNnVDd+uRtX2E0f11dAha0+ld/wdLMxgzTfQM1xQwcIMmYx
+TsavqAUQTOU2e7zB/9XNrlIqVvP07vDdBe8tySoW2+pHhpJ+rZY3BngAjstvSc0X
+YD/XdKbdkTHbLtPmzlrz43zujEsUOyQeNGgyslRXAgMBAAGjUDBOMB0GA1UdDgQW
+BBRrtIC0DvaaS/oKodIXFBBV5MdjMTAfBgNVHSMEGDAWgBRrtIC0DvaaS/oKodIX
+FBBV5MdjMTAMBgNVHRMEBTADAQH/MA0GCSqGSIb3DQEBCwUAA4IBAQCNzW7HdGd7
+swAquhaTpvYylK8/wVhz2F/rHZLHIsIXs87GReMtt+tT826xsBQOucKWD8EfZw3A
+WffPBAFv4mbRssX2OMTFtvrjZeTSJ468NQ+L0j8sYtFiWLQFe9R+LHZimD1iaJuH
+JiNddCBO0hsnqnBjSvc+VaNPdGUnC3bxK54ecw/msjVCQMWSARfPKLug0ZaDlyXN
+BCXCH3CvvNdeGJojJAqEqhuGn/ZtORhIYqhm5DovnvteQzb+K+7r6oE3mb7QWfne
+yrnUZ3gsfAiYX4iQrJjyNuuzzf3Kuw33E6yNcU0BnKv4XKZSslYJG0cfler70iuZ
+iPbz5D+Y+mxB
+-----END CERTIFICATE-----

tests/domain.key

@@ -0,0 +1,28 @@
+-----BEGIN PRIVATE KEY-----
+MIIEvwIBADANBgkqhkiG9w0BAQEFAASCBKkwggSlAgEAAoIBAQDSBOGs0Ckg6G/1
+VRHKIpQrLUfvsrUKv8SpW4RxIpmop67Poy/dfRAU+6YfVE1BdXZnyjGHmt3pJ57P
+ilN6Vm6CQbbRB5cfLlp/GexkNVmptRV1LgLGzCqTpkemFtGnV55eVMrEPvy291lU
+ID6As/evYO75AgfutoOPFBBqsPUYOW19fyntG/5LedqkVSJa2nNnVDd+uRtX2E0f
+11dAha0+ld/wdLMxgzTfQM1xQwcIMmYxTsavqAUQTOU2e7zB/9XNrlIqVvP07vDd
+Be8tySoW2+pHhpJ+rZY3BngAjstvSc0XYD/XdKbdkTHbLtPmzlrz43zujEsUOyQe
+NGgyslRXAgMBAAECggEAFnYBMLs+es+BEcrcUTLrh8aa9+6z4AdsX/RysKrTpi7w
+yjIYr2ZtruDO/nDEpmiKlkUQkQEmFIZxPJ79mIrf3LxsceSO6/3IHenG11zBTx4X
+wmV+VIaHUy0oYoH0QCJlxI7q9JwL6QGBIzcjmB4nMjF6vVTTRltsrIOn0dFFFSb1
+hRo9N+gpKlo3sxIpsT/AbVNgp8v4ERRSaek/NdOxZa2JrdjtOFd6fyIapdVcTNlq
+BnflAIDI800waY0bI61jvbkK28/s9VZX+LTNgoM8s3WUg8gSqTnnM3P+ExjAAYlz
+oaLnJLHzHa4+Ge0Umc/bq81IWGv4OC2Gfz2sVppVKQKBgQDxLTXklNAvt961/RQ1
+Z+TEyZSuyFzLBsGCLWmO4BKIW/5NJwXDs4IlvNVDcSKhx0JngytcJFzsmrI/3PYx
+PZB4ArhwoFllNbUwReejbH6zfbwMu9cll2a9iAN4qFmZ1vdhZyrO54Gs+7rzPal3
+48+ot46UZmwYC4EVPlKbmZHb6wKBgQDe7Wxb8r6HFFc+Bv1BH3eAe6Q25CO0rylw
+Uj6ywGcxQN02ZwOpIuL4B1VgcBw6Evau85xeOIljxQoqYZx1BrLP0W7m4Hz/9pXt
+xO2jLt+fFA45nFGmoz9hPFFcE7G9CPJa651TJeRyPmnPuIra+YSv1mTBK047uLYE
+CQn0eQiqRQKBgQCeuO1d+R4mOXFF/j/dgfmKFUxX1x9SsqiJvs9Yq1knCrwLFKJz
+kGEX6vXDw+JR7CxkqWvQ4VGSmaCz7PzNnUHF6axu2ujadgU/ttHiuNWN1TL8gmbu
+o3yHQctfDhSh7VqzbM08jBnIuYbo7GzqrTAOTnd0bQ+4p2GVlIWTR/c5bwKBgQCd
+MUVrQKvgjBNIhZ0mwo19Z7HPkucqjSjnLtdPLCSXsmlJ5zHJnfP1eRZWvjbpLF4N
+13G3bp5SX6pfYBz22BZWPBO7cSgzRSlLU8VNf8vnXDj7nqJkolAavHGSP8DgzuI2
+AlYZP2rqs9/gvPwMlE1fwAmjvGmZ6xFNnyIJ5R6bEQKBgQDRwYS+tOwkHI1TiWD+
+VHOo7sUa+M5VMNTvmJxN3nDlRtCmOgxxhcJbYxiC1lEGBwh4GLwnL70nQHz+41xx
+BfLfko9t7g/OIgxDMd/7CcEuby7Aisn3RS3juvLup5duQNqfoRDlZa8Ch8JddAh0
+Wzp9qZYdAVANZnyEXtjoKodoDQ==
+-----END PRIVATE KEY-----

tests/key

@@ -0,0 +1,13 @@
+97 57 83 1f 83 df 99 d6 a3 ab d8 b7 09 6c d0 77 00 00 00 00 00 00 00 00 00 00 00 00 
+cc c3 a5 85 b7 b3 4e bc f2 ee a3 a9 e3 ba 88 84 00 00 00 00 00 00 00 00 00 00 00 00 
+78 fb 72 98 ee 65 0f 70 8f fe e1 60 fe 9b a6 b2 00 00 00 00 00 00 00 00 00 00 00 00 
+dd e6 8e db 63 80 86 5c 70 32 0b 7b 31 cd 4c 1d 00 00 00 00 00 00 00 00 00 00 00 00 
+88 43 58 7b 45 c0 bf d4 78 c1 fe dc d2 77 dd 15 00 00 00 00 00 00 00 00 00 00 00 00 
+68 6a e2 d3 b2 ff 1a 47 e3 93 86 ad 7f 49 52 ff 00 00 00 00 00 00 00 00 00 00 00 00 
+74 f4 82 93 8c d9 ac 58 23 79 b4 2b bc 25 77 0f 00 00 00 00 00 00 00 00 00 00 00 00 
+fe da 27 1d 56 a4 13 7f 95 61 12 b0 48 05 40 1f 00 00 00 00 00 00 00 00 00 00 00 00 
+8d f7 44 1e 86 ac 42 7a cf 60 b0 fb fe 54 ac 39 00 00 00 00 00 00 00 00 00 00 00 00 
+7d 63 ae b2 2b e6 8f e6 ea f5 ed 40 0c 58 0e b8 00 00 00 00 00 00 00 00 00 00 00 00 
+c4 70 15 99 27 36 73 37 b9 91 f6 99 fa 2b 66 02 00 00 00 00 00 00 00 00 00 00 00 00 
+e7 8a 9a a6 15 ba 6a c7 8c 75 98 81 04 68 8c 8c 00 00 00 00 00 00 00 00 00 00 00 00 
+3f 64 c2 b2 6e 1d 14 db 52 c3 69 83 19 e5 47 e0 00 00 00 00 00 00 00 00 00 00 00 00 

tests/ptwist.h

@@ -0,0 +1,37 @@
+#ifndef __PTWIST_H__
+#define __PTWIST_H__
+
+#define PTWIST_BITS 168  /* must be a multiple of 8 */
+#define PTWIST_BYTES (PTWIST_BITS/8)
+
+#define PTWIST_TAG_BITS 224  /* must be a multiple of 8 */
+#define PTWIST_TAG_BYTES (PTWIST_TAG_BITS/8)
+
+#define PTWIST_PUZZLE_STRENGTH 0 /*21*/  /* set to 0 to disable client puzzle */
+#define PTWIST_PUZZLE_MASK ((1<<PTWIST_PUZZLE_STRENGTH)-1)
+
+#if PTWIST_PUZZLE_STRENGTH == 0
+#define PTWIST_RESP_BITS 0
+#else
+#define PTWIST_RESP_BITS (PTWIST_PUZZLE_STRENGTH+6)
+#endif
+
+#define PTWIST_RESP_BYTES ((PTWIST_RESP_BITS+7)/8)
+#define PTWIST_RESP_MASK ((1<<(((PTWIST_RESP_BITS&7)==0)?8:(PTWIST_RESP_BITS&7)))-1)
+
+#define PTWIST_HASH_SHOWBITS (PTWIST_TAG_BITS-PTWIST_BITS-PTWIST_RESP_BITS)
+#define PTWIST_HASH_TOTBITS (PTWIST_HASH_SHOWBITS+PTWIST_PUZZLE_STRENGTH)
+#define PTWIST_HASH_TOTBYTES ((PTWIST_HASH_TOTBITS+7)/8)
+#define PTWIST_HASH_MASK ((1<<(((PTWIST_HASH_TOTBITS&7)==0)?8:(PTWIST_HASH_TOTBITS&7)))-1)
+
+typedef unsigned char byte;
+
+/* Figure out whether there's a point with x-coordinate x on the main
+ * curve.  If not, then there's one on the twist curve.  (There are
+ * actually two, which are negatives of each other; that doesn't
+ * matter.)  Multiply that point by seckey and set out to the
+ * x-coordinate of the result. */
+void ptwist_pointmul(byte out[PTWIST_BYTES], const byte x[PTWIST_BYTES],
+	const byte seckey[PTWIST_BYTES]);
+
+#endif

tests/ptwist168.c

@@ -0,0 +1,1651 @@
+#include "ptwist.h"
+
+/* ptwist168.c by Ian Goldberg. Based on: */
+
+/* crypto/ec/ecp_nistp224.c */
+/*
+ * Written by Emilia Kasper (Google) for the OpenSSL project.
+ */
+/* ====================================================================
+ * Copyright (c) 2000-2010 The OpenSSL Project.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ *
+ * 2. 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.
+ *
+ * 3. All advertising materials mentioning features or use of this
+ *    software must display the following acknowledgment:
+ *    "This product includes software developed by the OpenSSL Project
+ *    for use in the OpenSSL Toolkit. (http://www.OpenSSL.org/)"
+ *
+ * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
+ *    endorse or promote products derived from this software without
+ *    prior written permission. For written permission, please contact
+ *    licensing@OpenSSL.org.
+ *
+ * 5. Products derived from this software may not be called "OpenSSL"
+ *    nor may "OpenSSL" appear in their names without prior written
+ *    permission of the OpenSSL Project.
+ *
+ * 6. Redistributions of any form whatsoever must retain the following
+ *    acknowledgment:
+ *    "This product includes software developed by the OpenSSL Project
+ *    for use in the OpenSSL Toolkit (http://www.OpenSSL.org/)"
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
+ * EXPRESSED 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 OpenSSL PROJECT OR
+ * ITS 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.
+ * ====================================================================
+ *
+ * This product includes cryptographic software written by Eric Young
+ * (eay@cryptsoft.com).  This product includes software written by Tim
+ * Hudson (tjh@cryptsoft.com).
+ *
+ */
+
+/*
+ * A 64-bit implementation of the NIST P-224 elliptic curve point multiplication
+ *
+ * Inspired by Daniel J. Bernstein's public domain nistp224 implementation
+ * and Adam Langley's public domain 64-bit C implementation of curve25519
+ */
+#include <stdint.h>
+#include <string.h>
+
+#if defined(__GNUC__) && (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 1))
+  /* even with gcc, the typedef won't work for 32-bit platforms */
+  typedef __uint128_t uint128_t; /* nonstandard; implemented by gcc on 64-bit platforms */
+#else
+  #error "Need GCC 3.1 or later to define type uint128_t"
+#endif
+
+typedef uint8_t u8;
+
+
+/******************************************************************************/
+/*		    INTERNAL REPRESENTATION OF FIELD ELEMENTS
+ *
+ * Field elements are represented as a_0 + 2^56*a_1 + 2^112*a_2
+ * where each slice a_i is a 64-bit word, i.e., a field element is an fslice
+ * array a with 3 elements, where a[i] = a_i.
+ * Outputs from multiplications are represented as unreduced polynomials
+ * b_0 + 2^56*b_1 + 2^112*b_2 + 2^168*b_3 + 2^224*b_4
+ * where each b_i is a 128-bit word. We ensure that inputs to each field
+ * multiplication satisfy a_i < 2^60, so outputs satisfy b_i < 4*2^60*2^60,
+ * and fit into a 128-bit word without overflow. The coefficients are then
+ * again partially reduced to a_i < 2^57. We only reduce to the unique minimal
+ * representation at the end of the computation.
+ *
+ */
+
+typedef uint64_t fslice;
+
+typedef fslice coord[3];
+typedef coord point[3];
+
+
+#include <stdio.h>
+#include <stdlib.h>
+
+/*static void dump_coord(const char *label, const coord c)
+{
+    if (label) fprintf(stderr, "%s: ", label);
+    printf("%016lx %016lx %016lx\n", c[2], c[1], c[0]);
+}*/
+
+/*static void dump_point(const char *label, point p)
+{
+    if (label) fprintf(stderr, "%s:\n", label);
+    dump_coord(" x", p[0]);
+    dump_coord(" y", p[1]);
+    dump_coord(" z", p[2]);
+}*/
+
+/* Field element represented as a byte arrary.
+ * 21*8 = 168 bits is also the group order size for the elliptic curve.  */
+typedef u8 felem_bytearray[21];
+
+static const felem_bytearray ptwist168_curve_params[5] = {
+	{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,    /* p */
+	 0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,
+	 0xFF},
+	{0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,    /* a */
+	 0xFF,0xFF,0xFF,0xFF,0xFF,0xFE,0xFF,0xFF,0xFF,0xFE,
+	 0xFC},
+	{0x4E,0x35,0x5E,0x95,0xCA,0xFE,0xDD,0x48,0x6E,0xBC,    /* b */
+	 0x69,0xBA,0xD3,0x16,0x46,0xD3,0x20,0xE0,0x1D,0xC7,
+	 0xD6},
+	{0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,    /* x */
+	 0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
+	 0x02},
+	{0xEA,0x67,0x47,0xB7,0x5A,0xF8,0xC7,0xF9,0x3C,0x1F,    /* y */
+	 0x5E,0x6D,0x32,0x0F,0x88,0xB9,0xBE,0x15,0x66,0xD2,
+	 0xF2}
+};
+
+/* Helper functions to convert field elements to/from internal representation */
+static void bin21_to_felem(fslice out[3], const u8 in[21])
+	{
+	out[0] = *((const uint64_t *)(in)) & 0x00ffffffffffffff;
+	out[1] = (*((const uint64_t *)(in+7))) & 0x00ffffffffffffff;
+	out[2] = (*((const uint64_t *)(in+14))) & 0x00ffffffffffffff;
+	}
+
+static void felem_to_bin21(u8 out[21], const fslice in[3])
+	{
+	unsigned i;
+	for (i = 0; i < 7; ++i)
+		{
+		out[i]	  = in[0]>>(8*i);
+		out[i+7]  = in[1]>>(8*i);
+		out[i+14] = in[2]>>(8*i);
+		}
+	}
+
+#if 0
+/* To preserve endianness when using BN_bn2bin and BN_bin2bn */
+static void flip_endian(u8 *out, const u8 *in, unsigned len)
+	{
+	unsigned i;
+	for (i = 0; i < len; ++i)
+		out[i] = in[len-1-i];
+	}
+#endif
+
+/******************************************************************************/
+/*				FIELD OPERATIONS
+ *
+ * Field operations, using the internal representation of field elements.
+ * NB! These operations are specific to our point multiplication and cannot be
+ * expected to be correct in general - e.g., multiplication with a large scalar
+ * will cause an overflow.
+ *
+ */
+
+/* Sum two field elements: out += in */
+static void felem_sum64(fslice out[3], const fslice in[3])
+	{
+	out[0] += in[0];
+	out[1] += in[1];
+	out[2] += in[2];
+	}
+
+/* Subtract field elements: out -= in */
+/* Assumes in[i] < 2^57 */
+static void felem_diff64(fslice out[3], const fslice in[3])
+	{
+	/* a = 3*2^56 - 3 */
+	/* b = 3*2^56 - 3*257 */
+	static const uint64_t a = (((uint64_t) 3) << 56) - ((uint64_t) 3);
+	static const uint64_t b = (((uint64_t) 3) << 56) - ((uint64_t) 771);
+
+	/* Add 0 mod 2^168-2^8-1 to ensure out > in at each element */
+	/* a*2^112 + a*2^56 + b = 3*p */
+	out[0] += b;
+	out[1] += a;
+	out[2] += a;
+
+	out[0] -= in[0];
+	out[1] -= in[1];
+	out[2] -= in[2];
+	}
+
+/* Subtract in unreduced 128-bit mode: out128 -= in128 */
+/* Assumes in[i] < 2^119 */
+static void felem_diff128(uint128_t out[5], const uint128_t in[5])
+	{
+	/* a = 3*2^118 - 192
+	   b = 3*2^118 - 49536
+	   c = 3*2^118
+	   d = 3*2^118 - 12681408
+
+	   a*2^224 + a*2^168 + b*2^112 + c*2^56 + d
+	    = (3*2^174 + 3*2^118 + 49344)*p
+	*/
+	static const uint128_t a = (((uint128_t)3) << 118) - ((uint128_t) 192);
+	static const uint128_t b = (((uint128_t)3) << 118) - ((uint128_t) 49536);
+	static const uint128_t c = (((uint128_t)3) << 118);
+	static const uint128_t d = (((uint128_t)3) << 118) - ((uint128_t) 12681408);;
+
+	/* Add 0 mod 2^168-2^8-1 to ensure out > in */
+	out[0] += d;
+	out[1] += c;
+	out[2] += b;
+	out[3] += a;
+	out[4] += a;
+
+	out[0] -= in[0];
+	out[1] -= in[1];
+	out[2] -= in[2];
+	out[3] -= in[3];
+	out[4] -= in[4];
+	}
+
+/* Subtract in mixed mode: out128 -= in64 */
+/* in[i] < 2^63 */
+static void felem_diff_128_64(uint128_t out[5], const fslice in[3])
+	{
+	/* a = 3*2^62 - 192
+	   b = 3*2^62 - 49344
+	   a*2^112 + a*2^56 + b = 192*p
+	*/
+	static const uint128_t a = (((uint128_t) 3) << 62) - ((uint128_t) 192);
+	static const uint128_t b = (((uint128_t) 3) << 62) - ((uint128_t) 49344);
+
+	/* Add 0 mod 2^168-2^8-1 to ensure out > in */
+	out[0] += b;
+	out[1] += a;
+	out[2] += a;
+
+	out[0] -= in[0];
+	out[1] -= in[1];
+	out[2] -= in[2];
+	}
+
+/* Multiply a field element by a scalar: out64 = out64 * scalar
+ * The scalars we actually use are small, so results fit without overflow */
+static void felem_scalar64(fslice out[3], const fslice scalar)
+	{
+	out[0] *= scalar;
+	out[1] *= scalar;
+	out[2] *= scalar;
+	}
+
+/* Multiply an unreduced field element by a scalar: out128 = out128 * scalar
+ * The scalars we actually use are small, so results fit without overflow */
+static void felem_scalar128(uint128_t out[5], const uint128_t scalar)
+	{
+	out[0] *= scalar;
+	out[1] *= scalar;
+	out[2] *= scalar;
+	out[3] *= scalar;
+	out[4] *= scalar;
+	}
+
+/* Square a field element: out = in^2 */
+static void felem_square(uint128_t out[5], const fslice in[3])
+	{
+	out[0] = ((uint128_t) in[0]) * in[0];
+	out[1] = ((uint128_t) in[0]) * in[1] * 2;
+	out[2] = ((uint128_t) in[0]) * in[2] * 2 + ((uint128_t) in[1]) * in[1];
+	out[3] = ((uint128_t) in[1]) * in[2] * 2;
+	out[4] = ((uint128_t) in[2]) * in[2];
+	}
+
+/* Multiply two field elements: out = in1 * in2 */
+static void felem_mul(uint128_t out[5], const fslice in1[3], const fslice in2[3])
+	{
+	out[0] = ((uint128_t) in1[0]) * in2[0];
+	out[1] = ((uint128_t) in1[0]) * in2[1] + ((uint128_t) in1[1]) * in2[0];
+	out[2] = ((uint128_t) in1[0]) * in2[2] + ((uint128_t) in1[1]) * in2[1] +
+		((uint128_t) in1[2]) * in2[0];
+	out[3] = ((uint128_t) in1[1]) * in2[2] +
+		((uint128_t) in1[2]) * in2[1];
+	out[4] = ((uint128_t) in1[2]) * in2[2];
+	}
+
+#define M257(x) (((x)<<8)+(x))
+
+/* Reduce 128-bit coefficients to 64-bit coefficients. Requires in[i] < 2^126,
+ * ensures out[0] < 2^56, out[1] < 2^56, out[2] < 2^57 */
+static void felem_reduce(fslice out[3], const uint128_t in[5])
+	{
+	static const uint128_t two56m1 = (((uint128_t) 1)<<56) -
+		((uint128_t)1);
+	uint128_t output[3];
+
+	output[0] = in[0];  /* < 2^126 */
+	output[1] = in[1];  /* < 2^126 */
+	output[2] = in[2];  /* < 2^126 */
+
+	/* Eliminate in[3], in[4] */
+	output[2] += M257(in[4] >> 56);       /* < 2^126 + 2^79 */
+	output[1] += M257(in[4] & two56m1);   /* < 2^126 + 2^65 */
+
+	output[1] += M257(in[3] >> 56);       /* < 2^126 + 2^65 + 2^79 */
+	output[0] += M257(in[3] & two56m1);   /* < 2^126 + 2^65 */
+
+	/* Eliminate the top part of output[2] */
+	output[0] += M257(output[2] >> 56);   /* < 2^126 + 2^65 + 2^79 */
+	output[2] &= two56m1;                 /* < 2^56 */
+
+	/* Carry 0 -> 1 -> 2 */
+	output[1] += output[0] >> 56;         /* < 2^126 + 2^71 */
+	output[0] &= two56m1;                 /* < 2^56 */
+
+	output[2] += output[1] >> 56;         /* < 2^71 */
+	output[1] &= two56m1;                 /* < 2^56 */
+
+	/* Eliminate the top part of output[2] */
+	output[0] += M257(output[2] >> 56);   /* < 2^57 */
+	output[2] &= two56m1;                 /* < 2^56 */
+
+	/* Carry 0 -> 1 -> 2 */
+	output[1] += output[0] >> 56;         /* <= 2^56 */
+	out[0] = output[0] & two56m1;         /* < 2^56 */
+
+	out[2] = output[2] + (output[1] >> 56);  /* <= 2^56 */
+	out[1] = output[1] & two56m1;         /* < 2^56 */
+
+	}
+
+/* Reduce to unique minimal representation */
+static void felem_contract(fslice out[3], const fslice in[3])
+	{
+	static const uint64_t two56m1 = (((uint64_t) 1)<<56) -
+		((uint64_t)1);
+	static const uint64_t two56m257 = (((uint64_t) 1)<<56) -
+		((uint64_t)257);
+	uint64_t a;
+
+	/* in[0] < 2^56, in[1] < 2^56, in[2] <= 2^56 */
+	/* so in < 2*p for sure */
+
+	/* Eliminate the top part of in[2] */
+	out[0] = in[0] + M257(in[2] >> 56);   /* < 2^57 */
+	out[2] = in[2] & two56m1;             /* < 2^56, but if out[0] >= 2^56
+	                                         then out[2] now = 0 */
+
+	/* Carry 0 -> 1 -> 2 */
+	out[1] = in[1] + (out[0] >> 56);      /* < 2^56 + 2, but if
+	                                         out[1] >= 2^56 then
+						 out[2] = 0 */
+	out[0] &= two56m1;                    /* < 2^56 */
+
+	out[2] += out[1] >> 56;               /* < 2^56 due to the above */
+	out[1] &= two56m1;                    /* < 2^56 */
+
+	/* Now out < 2^168, but it could still be > p */
+	a = ((out[2] == two56m1) & (out[1] == two56m1) & (out[0] >= two56m257));
+	out[2] -= two56m1*a;
+	out[1] -= two56m1*a;
+	out[0] -= two56m257*a;
+	}
+
+/* Negate a field element: out = -in */
+/* Assumes in[i] < 2^57 */
+static void felem_neg(fslice out[3], const fslice in[3])
+	{
+	/* a = 3*2^56 - 3 */
+	/* b = 3*2^56 - 3*257 */
+	static const uint64_t a = (((uint64_t) 3) << 56) - ((uint64_t) 3);
+	static const uint64_t b = (((uint64_t) 3) << 56) - ((uint64_t) 771);
+	static const uint64_t two56m1 = (((uint64_t) 1) << 56) - ((uint64_t) 1);
+	fslice tmp[3];
+
+	/* Add 0 mod 2^168-2^8-1 to ensure out > in at each element */
+	/* a*2^112 + a*2^56 + b = 3*p */
+	tmp[0] = b - in[0];
+	tmp[1] = a - in[1];
+	tmp[2] = a - in[2];
+
+	/* Carry 0 -> 1 -> 2 */
+	tmp[1] += tmp[0] >> 56;
+	tmp[0] &= two56m1;                 /* < 2^56 */
+
+	tmp[2] += tmp[1] >> 56;         /* < 2^71 */
+	tmp[1] &= two56m1;                 /* < 2^56 */
+
+	felem_contract(out, tmp);
+
+	}
+
+/* Zero-check: returns 1 if input is 0, and 0 otherwise.
+ * We know that field elements are reduced to in < 2^169,
+ * so we only need to check three cases: 0, 2^168 - 2^8 - 1,
+ * and 2^169 - 2^9 - 2 */
+static fslice felem_is_zero(const fslice in[3])
+	{
+	fslice zero, two168m8m1, two169m9m2;
+	static const uint64_t two56m1 = (((uint64_t) 1)<<56) -
+		((uint64_t)1);
+	static const uint64_t two56m257 = (((uint64_t) 1)<<56) -
+		((uint64_t)257);
+	static const uint64_t two57m1 = (((uint64_t) 1)<<57) -
+		((uint64_t)1);
+	static const uint64_t two56m514 = (((uint64_t) 1)<<56) -
+		((uint64_t)514);
+
+	zero = (in[0] == 0) & (in[1] == 0) & (in[2] == 0);
+	two168m8m1 = (in[2] == two56m1) & (in[1] == two56m1) &
+			(in[0] == two56m257);
+	two169m9m2 = (in[2] == two57m1) & (in[1] == two56m1) &
+			(in[0] == two56m514);
+
+	return (zero | two168m8m1 | two169m9m2);
+	}
+
+/* Invert a field element */
+static void felem_inv(fslice out[3], const fslice in[3])
+	{
+	fslice ftmp[3], ftmp2[3], ftmp3[3], ftmp4[3];
+	uint128_t tmp[5];
+	unsigned i;
+
+	felem_square(tmp, in); felem_reduce(ftmp, tmp);		/* 2 */
+	felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp);	/* 2^2 - 1 */
+								/* = ftmp */
+
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^3 - 2 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^4 - 2^2 */
+	felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp2, tmp);	/* 2^4 - 1 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^5 - 2 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^6 - 2^2 */
+	felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp);	/* 2^6 - 1 */
+								/* = ftmp */
+
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^7 - 2 */
+	for (i = 0; i < 5; ++i)					/* 2^12 - 2^6 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp, ftmp2); felem_reduce(ftmp3, tmp);	/* 2^12 - 1 */
+								/* = ftmp3 */
+
+	felem_square(tmp, ftmp3); felem_reduce(ftmp2, tmp);	/* 2^13 - 2 */
+	for (i = 0; i < 11; ++i)				/* 2^24 - 2^12 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp2, ftmp3); felem_reduce(ftmp3, tmp);	/* 2^24 - 1 */
+								/* = ftmp3 */
+	felem_square(tmp, ftmp3); felem_reduce(ftmp2, tmp);	/* 2^25 - 2 */
+	for (i = 0; i < 23; ++i)				/* 2^48 - 2^24 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp2, ftmp3); felem_reduce(ftmp4, tmp);	/* 2^48 - 1 */
+								/* = ftmp4 */
+	felem_square(tmp, ftmp4); felem_reduce(ftmp2, tmp);	/* 2^49 - 2 */
+	for (i = 0; i < 23; ++i)				/* 2^72 - 2^24 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp2, ftmp3); felem_reduce(ftmp4, tmp);	/* 2^72 - 1 */
+								/* = ftmp4 */
+
+	felem_square(tmp, ftmp4); felem_reduce(ftmp2, tmp);	/* 2^73 - 2 */
+	for (i = 0; i < 5; ++i)					/* 2^78 - 2^6 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^78 - 1 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^79 - 2 */
+	felem_mul(tmp, in, ftmp2); felem_reduce(ftmp4, tmp);	/* 2^79 - 1 */
+								/* = ftmp4 */
+	felem_square(tmp, ftmp4); felem_reduce(ftmp2, tmp);	/* 2^80 - 2 */
+	for (i = 0; i < 78; ++i)				/* 2^158 - 2^79 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp4, ftmp2); felem_reduce(ftmp2, tmp); /* 2^158 - 1 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^159 - 2 */
+	felem_mul(tmp, in, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^159 - 1 */
+	for (i = 0; i < 7; ++i)					/* 2^166 - 2^7 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^166 - 2^6 - 1 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^167 - 2^7 - 2 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^168 - 2^8 - 4 */
+	felem_mul(tmp, in, ftmp2); felem_reduce(out, tmp);	/* 2^168 - 2^8 - 3 */
+								/* = out */
+	}
+
+/* Take the square root of a field element */
+static void felem_sqrt(fslice out[3], const fslice in[3])
+	{
+	fslice ftmp[3], ftmp2[3];
+	uint128_t tmp[5];
+	unsigned i;
+
+	felem_square(tmp, in); felem_reduce(ftmp, tmp);		/* 2 */
+	felem_mul(tmp, in, ftmp); felem_reduce(ftmp, tmp);	/* 2^2 - 1 */
+								/* = ftmp */
+
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^3 - 2 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^4 - 2^2 */
+	felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp2, tmp);	/* 2^4 - 1 */
+	felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^5 - 2 */
+	felem_mul(tmp, ftmp2, in); felem_reduce(ftmp, tmp);	/* 2^5 - 1 */
+								/* = ftmp */
+
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^6 - 2 */
+	for (i = 0; i < 4; ++i)					/* 2^10 - 2^5 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp, ftmp2); felem_reduce(ftmp, tmp);	/* 2^10 - 1 */
+								/* = ftmp */
+
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^11 - 2 */
+	for (i = 0; i < 9; ++i)					/* 2^20 - 2^10 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp);	/* 2^20 - 1 */
+								/* = ftmp */
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^21 - 2 */
+	for (i = 0; i < 19; ++i)				/* 2^40 - 2^20 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp);	/* 2^40 - 1 */
+								/* = ftmp */
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^41 - 2 */
+	for (i = 0; i < 39; ++i)				/* 2^80 - 2^40 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp2, ftmp); felem_reduce(ftmp, tmp);	/* 2^80 - 1 */
+								/* = ftmp */
+
+	felem_square(tmp, ftmp); felem_reduce(ftmp2, tmp);	/* 2^81 - 2 */
+	for (i = 0; i < 79; ++i)				/* 2^160 - 2^80 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_mul(tmp, ftmp, ftmp2); felem_reduce(ftmp2, tmp);	/* 2^160 - 1 */
+	for (i = 0; i < 5; ++i)					/* 2^165 - 2^5 */
+		{
+		felem_square(tmp, ftmp2); felem_reduce(ftmp2, tmp);
+		}
+	felem_square(tmp, ftmp2); felem_reduce(out, tmp);	/* 2^166 - 2^6 */
+								/* = out */
+	}
+
+/* Copy in constant time:
+ * if icopy == 1, copy in to out,
+ * if icopy == 0, copy out to itself. */
+static void
+copy_conditional(fslice *out, const fslice *in, unsigned len, fslice icopy)
+	{
+	unsigned i;
+	/* icopy is a (64-bit) 0 or 1, so copy is either all-zero or all-one */
+	const fslice copy = -icopy;
+	for (i = 0; i < len; ++i)
+		{
+		const fslice tmp = copy & (in[i] ^ out[i]);
+		out[i] ^= tmp;
+		}
+	}
+
+/* Copy in constant time:
+ * if isel == 1, copy in2 to out,
+ * if isel == 0, copy in1 to out. */
+static void select_conditional(fslice *out, const fslice *in1, const fslice *in2,
+	unsigned len, fslice isel)
+	{
+	unsigned i;
+	/* isel is a (64-bit) 0 or 1, so sel is either all-zero or all-one */
+	const fslice sel = -isel;
+	for (i = 0; i < len; ++i)
+		{
+		const fslice tmp = sel & (in1[i] ^ in2[i]);
+		out[i] = in1[i] ^ tmp;
+		}
+}
+
+/******************************************************************************/
+/*			 ELLIPTIC CURVE POINT OPERATIONS
+ *
+ * Points are represented in Jacobian projective coordinates:
+ * (X, Y, Z) corresponds to the affine point (X/Z^2, Y/Z^3),
+ * or to the point at infinity if Z == 0.
+ *
+ */
+
+/* Double an elliptic curve point:
+ * (X', Y', Z') = 2 * (X, Y, Z), where
+ * X' = (3 * (X - Z^2) * (X + Z^2))^2 - 8 * X * Y^2
+ * Y' = 3 * (X - Z^2) * (X + Z^2) * (4 * X * Y^2 - X') - 8 * Y^2
+ * Z' = (Y + Z)^2 - Y^2 - Z^2 = 2 * Y * Z
+ * Outputs can equal corresponding inputs, i.e., x_out == x_in is allowed,
+ * while x_out == y_in is not (maybe this works, but it's not tested). */
+static void
+point_double(fslice x_out[3], fslice y_out[3], fslice z_out[3],
+	     const fslice x_in[3], const fslice y_in[3], const fslice z_in[3])
+	{
+	uint128_t tmp[5], tmp2[5];
+	fslice delta[3];
+	fslice gamma[3];
+	fslice beta[3];
+	fslice alpha[3];
+	fslice ftmp[3], ftmp2[3];
+	memcpy(ftmp, x_in, 3 * sizeof(fslice));
+	memcpy(ftmp2, x_in, 3 * sizeof(fslice));
+
+	/* delta = z^2 */
+	felem_square(tmp, z_in);
+	felem_reduce(delta, tmp);
+
+	/* gamma = y^2 */
+	felem_square(tmp, y_in);
+	felem_reduce(gamma, tmp);
+
+	/* beta = x*gamma */
+	felem_mul(tmp, x_in, gamma);
+	felem_reduce(beta, tmp);
+
+	/* alpha = 3*(x-delta)*(x+delta) */
+	felem_diff64(ftmp, delta);
+	/* ftmp[i] < 2^57 + 2^58 + 2 < 2^59 */
+	felem_sum64(ftmp2, delta);
+	/* ftmp2[i] < 2^57 + 2^57 = 2^58 */
+	felem_scalar64(ftmp2, 3);
+	/* ftmp2[i] < 3 * 2^58 < 2^60 */
+	felem_mul(tmp, ftmp, ftmp2);
+	/* tmp[i] < 2^60 * 2^59 * 4 = 2^121 */
+	felem_reduce(alpha, tmp);
+
+	/* x' = alpha^2 - 8*beta */
+	felem_square(tmp, alpha);
+	/* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
+	memcpy(ftmp, beta, 3 * sizeof(fslice));
+	felem_scalar64(ftmp, 8);
+	/* ftmp[i] < 8 * 2^57 = 2^60 */
+	felem_diff_128_64(tmp, ftmp);
+	/* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
+	felem_reduce(x_out, tmp);
+
+	/* z' = (y + z)^2 - gamma - delta */
+	felem_sum64(delta, gamma);
+	/* delta[i] < 2^57 + 2^57 = 2^58 */
+	memcpy(ftmp, y_in, 3 * sizeof(fslice));
+	felem_sum64(ftmp, z_in);
+	/* ftmp[i] < 2^57 + 2^57 = 2^58 */
+	felem_square(tmp, ftmp);
+	/* tmp[i] < 4 * 2^58 * 2^58 = 2^118 */
+	felem_diff_128_64(tmp, delta);
+	/* tmp[i] < 2^118 + 2^64 + 8 < 2^119 */
+	felem_reduce(z_out, tmp);
+
+	/* y' = alpha*(4*beta - x') - 8*gamma^2 */
+	felem_scalar64(beta, 4);
+	/* beta[i] < 4 * 2^57 = 2^59 */
+	felem_diff64(beta, x_out);
+	/* beta[i] < 2^59 + 2^58 + 2 < 2^60 */
+	felem_mul(tmp, alpha, beta);
+	/* tmp[i] < 4 * 2^57 * 2^60 = 2^119 */
+	felem_square(tmp2, gamma);
+	/* tmp2[i] < 4 * 2^57 * 2^57 = 2^116 */
+	felem_scalar128(tmp2, 8);
+	/* tmp2[i] < 8 * 2^116 = 2^119 */
+	felem_diff128(tmp, tmp2);
+	/* tmp[i] < 2^119 + 2^120 < 2^121 */
+	felem_reduce(y_out, tmp);
+	}
+
+/* Add two elliptic curve points:
+ * (X_1, Y_1, Z_1) + (X_2, Y_2, Z_2) = (X_3, Y_3, Z_3), where
+ * X_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1)^2 - (Z_1^2 * X_2 - Z_2^2 * X_1)^3 -
+ * 2 * Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2
+ * Y_3 = (Z_1^3 * Y_2 - Z_2^3 * Y_1) * (Z_2^2 * X_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^2 - X_3) -
+ *        Z_2^3 * Y_1 * (Z_1^2 * X_2 - Z_2^2 * X_1)^3
+ * Z_3 = (Z_1^2 * X_2 - Z_2^2 * X_1) * (Z_1 * Z_2) */
+
+/* This function is not entirely constant-time:
+ * it includes a branch for checking whether the two input points are equal,
+ * (while not equal to the point at infinity).
+ * This case never happens during single point multiplication,
+ * so there is no timing leak for ECDH or ECDSA signing. */
+static void point_add(fslice x3[3], fslice y3[3], fslice z3[3],
+	const fslice x1[3], const fslice y1[3], const fslice z1[3],
+	const fslice x2[3], const fslice y2[3], const fslice z2[3])
+	{
+	fslice ftmp[3], ftmp2[3], ftmp3[3], ftmp4[3], ftmp5[3];
+	fslice xout[3], yout[3], zout[3];
+	uint128_t tmp[5], tmp2[5];
+	fslice z1_is_zero, z2_is_zero, x_equal, y_equal;
+
+	/* ftmp = z1^2 */
+	felem_square(tmp, z1);
+	felem_reduce(ftmp, tmp);
+
+	/* ftmp2 = z2^2 */
+	felem_square(tmp, z2);
+	felem_reduce(ftmp2, tmp);
+
+	/* ftmp3 = z1^3 */
+	felem_mul(tmp, ftmp, z1);
+	felem_reduce(ftmp3, tmp);
+
+	/* ftmp4 = z2^3 */
+	felem_mul(tmp, ftmp2, z2);
+	felem_reduce(ftmp4, tmp);
+
+	/* ftmp3 = z1^3*y2 */
+	felem_mul(tmp, ftmp3, y2);
+	/* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
+
+	/* ftmp4 = z2^3*y1 */
+	felem_mul(tmp2, ftmp4, y1);
+	felem_reduce(ftmp4, tmp2);
+
+	/* ftmp3 = z1^3*y2 - z2^3*y1 */
+	felem_diff_128_64(tmp, ftmp4);
+	/* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
+	felem_reduce(ftmp3, tmp);
+
+	/* ftmp = z1^2*x2 */
+	felem_mul(tmp, ftmp, x2);
+	/* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
+
+	/* ftmp2 =z2^2*x1 */
+	felem_mul(tmp2, ftmp2, x1);
+	felem_reduce(ftmp2, tmp2);
+
+	/* ftmp = z1^2*x2 - z2^2*x1 */
+	felem_diff128(tmp, tmp2);
+	/* tmp[i] < 2^116 + 2^64 + 8 < 2^117 */
+	felem_reduce(ftmp, tmp);
+
+	/* the formulae are incorrect if the points are equal
+	 * so we check for this and do doubling if this happens */
+	x_equal = felem_is_zero(ftmp);
+	y_equal = felem_is_zero(ftmp3);
+	z1_is_zero = felem_is_zero(z1);
+	z2_is_zero = felem_is_zero(z2);
+	/* In affine coordinates, (X_1, Y_1) == (X_2, Y_2) */
+	if (x_equal && y_equal && !z1_is_zero && !z2_is_zero)
+		{
+		point_double(x3, y3, z3, x1, y1, z1);
+		return;
+		}
+
+	/* ftmp5 = z1*z2 */
+	felem_mul(tmp, z1, z2);
+	felem_reduce(ftmp5, tmp);
+
+	/* zout = (z1^2*x2 - z2^2*x1)*(z1*z2) */
+	felem_mul(tmp, ftmp, ftmp5);
+	felem_reduce(zout, tmp);
+
+	/* ftmp = (z1^2*x2 - z2^2*x1)^2 */
+	memcpy(ftmp5, ftmp, 3 * sizeof(fslice));
+	felem_square(tmp, ftmp);
+	felem_reduce(ftmp, tmp);
+
+	/* ftmp5 = (z1^2*x2 - z2^2*x1)^3 */
+	felem_mul(tmp, ftmp, ftmp5);
+	felem_reduce(ftmp5, tmp);
+
+	/* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
+	felem_mul(tmp, ftmp2, ftmp);
+	felem_reduce(ftmp2, tmp);
+
+	/* ftmp4 = z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
+	felem_mul(tmp, ftmp4, ftmp5);
+	/* tmp[i] < 4 * 2^57 * 2^57 = 2^116 */
+
+	/* tmp2 = (z1^3*y2 - z2^3*y1)^2 */
+	felem_square(tmp2, ftmp3);
+	/* tmp2[i] < 4 * 2^57 * 2^57 < 2^116 */
+
+	/* tmp2 = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 */
+	felem_diff_128_64(tmp2, ftmp5);
+	/* tmp2[i] < 2^116 + 2^64 + 8 < 2^117 */
+
+	/* ftmp5 = 2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
+	memcpy(ftmp5, ftmp2, 3 * sizeof(fslice));
+	felem_scalar64(ftmp5, 2);
+	/* ftmp5[i] < 2 * 2^57 = 2^58 */
+
+	/* xout = (z1^3*y2 - z2^3*y1)^2 - (z1^2*x2 - z2^2*x1)^3 -
+	   2*z2^2*x1*(z1^2*x2 - z2^2*x1)^2 */
+	felem_diff_128_64(tmp2, ftmp5);
+	/* tmp2[i] < 2^117 + 2^64 + 8 < 2^118 */
+	felem_reduce(xout, tmp2);
+
+	/* ftmp2 = z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - xout */
+	felem_diff64(ftmp2, xout);
+	/* ftmp2[i] < 2^57 + 2^58 + 2 < 2^59 */
+
+	/* tmp2 = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - xout) */
+	felem_mul(tmp2, ftmp3, ftmp2);
+	/* tmp2[i] < 4 * 2^57 * 2^59 = 2^118 */
+
+	/* yout = (z1^3*y2 - z2^3*y1)*(z2^2*x1*(z1^2*x2 - z2^2*x1)^2 - xout) -
+	   z2^3*y1*(z1^2*x2 - z2^2*x1)^3 */
+	felem_diff128(tmp2, tmp);
+	/* tmp2[i] < 2^118 + 2^120 < 2^121 */
+	felem_reduce(yout, tmp2);
+
+	/* the result (xout, yout, zout) is incorrect if one of the
+	 * inputs is the point at infinity, so we need to check for this
+	 * separately */
+
+	/* if point 1 is at infinity, copy point 2 to output, and vice versa */
+	copy_conditional(xout, x2, 3, z1_is_zero);
+	select_conditional(x3, xout, x1, 3, z2_is_zero);
+	copy_conditional(yout, y2, 3, z1_is_zero);
+	select_conditional(y3, yout, y1, 3, z2_is_zero);
+	copy_conditional(zout, z2, 3, z1_is_zero);
+	select_conditional(z3, zout, z1, 3, z2_is_zero);
+	}
+
+/*static void affine(point P)
+{
+    coord z1, z2, xin, yin;
+    uint128_t tmp[7];
+
+    if (felem_is_zero(P[2])) return;
+    felem_inv(z2, P[2]);
+    felem_square(tmp, z2); felem_reduce(z1, tmp);
+    felem_mul(tmp, P[0], z1); felem_reduce(xin, tmp);
+    felem_contract(P[0], xin);
+    felem_mul(tmp, z1, z2); felem_reduce(z1, tmp);
+    felem_mul(tmp, P[1], z1); felem_reduce(yin, tmp);
+    felem_contract(P[1], yin);
+    memset(P[2], 0, sizeof(coord));
+    P[2][0] = 1;
+}*/
+
+static void affine_x(coord out, point P)
+{
+    coord z1, z2, xin;
+    uint128_t tmp[7];
+
+    if (felem_is_zero(P[2])) return;
+    felem_inv(z2, P[2]);
+    felem_square(tmp, z2); felem_reduce(z1, tmp);
+    felem_mul(tmp, P[0], z1); felem_reduce(xin, tmp);
+    felem_contract(out, xin);
+}
+
+/* Multiply the given point by s */
+static void point_mul(point out, point in, const felem_bytearray s)
+{
+    int i;
+    point tmp;
+
+    point table[16];
+    memset(table[0], 0, sizeof(point));
+    memmove(table[1], in, sizeof(point));
+    for(i=2; i<16; i+=2) {
+	point_double(table[i][0], table[i][1], table[i][2],
+		     table[i/2][0], table[i/2][1], table[i/2][2]);
+	point_add(table[i+1][0], table[i+1][1], table[i+1][2],
+		  table[i][0], table[i][1], table[i][2],
+		  in[0], in[1], in[2]);
+    }
+    /*
+    for(i=0;i<16;++i) {
+	fprintf(stderr, "table[%d]:\n", i);
+	affine(table[i]);
+	dump_point(NULL, table[i]);
+    }
+    */
+
+    memset(tmp, 0, sizeof(point));
+    for(i=0;i<21;i++) {
+	u8 oh = s[20-i] >> 4;
+	u8 ol = s[20-i] & 0x0f;
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_add(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2],
+		  table[oh][0], table[oh][1], table[oh][2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_double(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2]);
+	point_add(tmp[0], tmp[1], tmp[2], tmp[0], tmp[1], tmp[2],
+		  table[ol][0], table[ol][1], table[ol][2]);
+    }
+    memmove(out, tmp, sizeof(point));
+}
+
+#if 0
+/* Select a point from an array of 16 precomputed point multiples,
+ * in constant time: for bits = {b_0, b_1, b_2, b_3}, return the point
+ * pre_comp[8*b_3 + 4*b_2 + 2*b_1 + b_0] */
+static void select_point(const fslice bits[4], const fslice pre_comp[16][3][4],
+	fslice out[12])
+	{
+	fslice tmp[5][12];
+	select_conditional(tmp[0], pre_comp[7][0], pre_comp[15][0], 12, bits[3]);
+	select_conditional(tmp[1], pre_comp[3][0], pre_comp[11][0], 12, bits[3]);
+	select_conditional(tmp[2], tmp[1], tmp[0], 12, bits[2]);
+	select_conditional(tmp[0], pre_comp[5][0], pre_comp[13][0], 12, bits[3]);
+	select_conditional(tmp[1], pre_comp[1][0], pre_comp[9][0], 12, bits[3]);
+	select_conditional(tmp[3], tmp[1], tmp[0], 12, bits[2]);
+	select_conditional(tmp[4], tmp[3], tmp[2], 12, bits[1]);
+	select_conditional(tmp[0], pre_comp[6][0], pre_comp[14][0], 12, bits[3]);
+	select_conditional(tmp[1], pre_comp[2][0], pre_comp[10][0], 12, bits[3]);
+	select_conditional(tmp[2], tmp[1], tmp[0], 12, bits[2]);
+	select_conditional(tmp[0], pre_comp[4][0], pre_comp[12][0], 12, bits[3]);
+	select_conditional(tmp[1], pre_comp[0][0], pre_comp[8][0], 12, bits[3]);
+	select_conditional(tmp[3], tmp[1], tmp[0], 12, bits[2]);
+	select_conditional(tmp[1], tmp[3], tmp[2], 12, bits[1]);
+	select_conditional(out, tmp[1], tmp[4], 12, bits[0]);
+	}
+
+/* Interleaved point multiplication using precomputed point multiples:
+ * The small point multiples 0*P, 1*P, ..., 15*P are in pre_comp[],
+ * the scalars in scalars[]. If g_scalar is non-NULL, we also add this multiple
+ * of the generator, using certain (large) precomputed multiples in g_pre_comp.
+ * Output point (X, Y, Z) is stored in x_out, y_out, z_out */
+static void batch_mul(fslice x_out[4], fslice y_out[4], fslice z_out[4],
+	const felem_bytearray scalars[], const unsigned num_points, const u8 *g_scalar,
+	const fslice pre_comp[][16][3][4], const fslice g_pre_comp[16][3][4])
+	{
+	unsigned i, j, num;
+	unsigned gen_mul = (g_scalar != NULL);
+	fslice nq[12], nqt[12], tmp[12];
+	fslice bits[4];
+	u8 byte;
+
+	/* set nq to the point at infinity */
+	memset(nq, 0, 12 * sizeof(fslice));
+
+	/* Loop over all scalars msb-to-lsb, 4 bits at a time: for each nibble,
+	 * double 4 times, then add the precomputed point multiples.
+	 * If we are also adding multiples of the generator, then interleave
+	 * these additions with the last 56 doublings. */
+	for (i = (num_points ? 28 : 7); i > 0; --i)
+		{
+		for (j = 0; j < 8; ++j)
+			{
+			/* double once */
+			point_double(nq, nq+4, nq+8, nq, nq+4, nq+8);
+			/* add multiples of the generator */
+			if ((gen_mul) && (i <= 7))
+				{
+				bits[3] = (g_scalar[i+20] >> (7-j)) & 1;
+				bits[2] = (g_scalar[i+13] >> (7-j)) & 1;
+				bits[1] = (g_scalar[i+6] >> (7-j)) & 1;
+				bits[0] = (g_scalar[i-1] >> (7-j)) & 1;
+				/* select the point to add, in constant time */
+				select_point(bits, g_pre_comp, tmp);
+				memcpy(nqt, nq, 12 * sizeof(fslice));
+				point_add(nq, nq+4, nq+8, nqt, nqt+4, nqt+8,
+					tmp, tmp+4, tmp+8);
+				}
+			/* do an addition after every 4 doublings */
+			if (j % 4 == 3)
+				{
+				/* loop over all scalars */
+				for (num = 0; num < num_points; ++num)
+					{
+					byte = scalars[num][i-1];
+					bits[3] = (byte >> (10-j)) & 1;
+					bits[2] = (byte >> (9-j)) & 1;
+					bits[1] = (byte >> (8-j)) & 1;
+					bits[0] = (byte >> (7-j)) & 1;
+					/* select the point to add */
+					select_point(bits,
+						pre_comp[num], tmp);
+					memcpy(nqt, nq, 12 * sizeof(fslice));
+					point_add(nq, nq+4, nq+8, nqt, nqt+4,
+						nqt+8, tmp, tmp+4, tmp+8);
+					}
+				}
+			}
+		}
+	memcpy(x_out, nq, 4 * sizeof(fslice));
+	memcpy(y_out, nq+4, 4 * sizeof(fslice));
+	memcpy(z_out, nq+8, 4 * sizeof(fslice));
+	}
+
+/******************************************************************************/
+/*		       FUNCTIONS TO MANAGE PRECOMPUTATION
+ */
+
+static NISTP224_PRE_COMP *nistp224_pre_comp_new()
+	{
+	NISTP224_PRE_COMP *ret = NULL;
+	ret = (NISTP224_PRE_COMP *)OPENSSL_malloc(sizeof(NISTP224_PRE_COMP));
+	if (!ret)
+		{
+		ECerr(EC_F_NISTP224_PRE_COMP_NEW, ERR_R_MALLOC_FAILURE);
+		return ret;
+		}
+	memset(ret->g_pre_comp, 0, sizeof(ret->g_pre_comp));
+	ret->references = 1;
+	return ret;
+	}
+
+static void *nistp224_pre_comp_dup(void *src_)
+	{
+	NISTP224_PRE_COMP *src = src_;
+
+	/* no need to actually copy, these objects never change! */
+	CRYPTO_add(&src->references, 1, CRYPTO_LOCK_EC_PRE_COMP);
+
+	return src_;
+	}
+
+static void nistp224_pre_comp_free(void *pre_)
+	{
+	int i;
+	NISTP224_PRE_COMP *pre = pre_;
+
+	if (!pre)
+		return;
+
+	i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
+	if (i > 0)
+		return;
+
+	OPENSSL_free(pre);
+	}
+
+static void nistp224_pre_comp_clear_free(void *pre_)
+	{
+	int i;
+	NISTP224_PRE_COMP *pre = pre_;
+
+	if (!pre)
+		return;
+
+	i = CRYPTO_add(&pre->references, -1, CRYPTO_LOCK_EC_PRE_COMP);
+	if (i > 0)
+		return;
+
+	OPENSSL_cleanse(pre, sizeof *pre);
+	OPENSSL_free(pre);
+	}
+
+/******************************************************************************/
+/*			   OPENSSL EC_METHOD FUNCTIONS
+ */
+
+int ec_GFp_nistp224_group_init(EC_GROUP *group)
+	{
+	int ret;
+	ret = ec_GFp_simple_group_init(group);
+	group->a_is_minus3 = 1;
+	return ret;
+	}
+
+int ec_GFp_nistp224_group_set_curve(EC_GROUP *group, const BIGNUM *p,
+	const BIGNUM *a, const BIGNUM *b, BN_CTX *ctx)
+	{
+	int ret = 0;
+	BN_CTX *new_ctx = NULL;
+	BIGNUM *curve_p, *curve_a, *curve_b;
+
+	if (ctx == NULL)
+		if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
+	BN_CTX_start(ctx);
+	if (((curve_p = BN_CTX_get(ctx)) == NULL) ||
+		((curve_a = BN_CTX_get(ctx)) == NULL) ||
+		((curve_b = BN_CTX_get(ctx)) == NULL)) goto err;
+	BN_bin2bn(nistp224_curve_params[0], sizeof(felem_bytearray), curve_p);
+	BN_bin2bn(nistp224_curve_params[1], sizeof(felem_bytearray), curve_a);
+	BN_bin2bn(nistp224_curve_params[2], sizeof(felem_bytearray), curve_b);
+	if ((BN_cmp(curve_p, p)) || (BN_cmp(curve_a, a)) ||
+		(BN_cmp(curve_b, b)))
+		{
+		ECerr(EC_F_EC_GFP_NISTP224_GROUP_SET_CURVE,
+			EC_R_WRONG_CURVE_PARAMETERS);
+		goto err;
+		}
+	group->field_mod_func = BN_nist_mod_224;
+	ret = ec_GFp_simple_group_set_curve(group, p, a, b, ctx);
+err:
+	BN_CTX_end(ctx);
+	if (new_ctx != NULL)
+		BN_CTX_free(new_ctx);
+	return ret;
+	}
+
+/* Takes the Jacobian coordinates (X, Y, Z) of a point and returns
+ * (X', Y') = (X/Z^2, Y/Z^3) */
+int ec_GFp_nistp224_point_get_affine_coordinates(const EC_GROUP *group,
+	const EC_POINT *point, BIGNUM *x, BIGNUM *y, BN_CTX *ctx)
+	{
+	fslice z1[4], z2[4], x_in[4], y_in[4], x_out[4], y_out[4];
+	uint128_t tmp[7];
+
+	if (EC_POINT_is_at_infinity(group, point))
+		{
+		ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
+			EC_R_POINT_AT_INFINITY);
+		return 0;
+		}
+	if ((!BN_to_felem(x_in, &point->X)) || (!BN_to_felem(y_in, &point->Y)) ||
+		(!BN_to_felem(z1, &point->Z))) return 0;
+	felem_inv(z2, z1);
+	felem_square(tmp, z2); felem_reduce(z1, tmp);
+	felem_mul(tmp, x_in, z1); felem_reduce(x_in, tmp);
+	felem_contract(x_out, x_in);
+	if (x != NULL)
+		{
+		if (!felem_to_BN(x, x_out)) {
+		ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
+			ERR_R_BN_LIB);
+		return 0;
+		}
+		}
+	felem_mul(tmp, z1, z2); felem_reduce(z1, tmp);
+	felem_mul(tmp, y_in, z1); felem_reduce(y_in, tmp);
+	felem_contract(y_out, y_in);
+	if (y != NULL)
+		{
+		if (!felem_to_BN(y, y_out)) {
+		ECerr(EC_F_EC_GFP_NISTP224_POINT_GET_AFFINE_COORDINATES,
+			ERR_R_BN_LIB);
+		return 0;
+		}
+		}
+	return 1;
+	}
+
+/* Computes scalar*generator + \sum scalars[i]*points[i], ignoring NULL values
+ * Result is stored in r (r can equal one of the inputs). */
+int ec_GFp_nistp224_points_mul(const EC_GROUP *group, EC_POINT *r,
+	const BIGNUM *scalar, size_t num, const EC_POINT *points[],
+	const BIGNUM *scalars[], BN_CTX *ctx)
+	{
+	int ret = 0;
+	int i, j;
+	BN_CTX *new_ctx = NULL;
+	BIGNUM *x, *y, *z, *tmp_scalar;
+	felem_bytearray g_secret;
+	felem_bytearray *secrets = NULL;
+	fslice (*pre_comp)[16][3][4] = NULL;
+	felem_bytearray tmp;
+	unsigned num_bytes;
+	int have_pre_comp = 0;
+	size_t num_points = num;
+	fslice x_in[4], y_in[4], z_in[4], x_out[4], y_out[4], z_out[4];
+	NISTP224_PRE_COMP *pre = NULL;
+	fslice (*g_pre_comp)[3][4] = NULL;
+	EC_POINT *generator = NULL;
+	const EC_POINT *p = NULL;
+	const BIGNUM *p_scalar = NULL;
+
+	if (ctx == NULL)
+		if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
+	BN_CTX_start(ctx);
+	if (((x = BN_CTX_get(ctx)) == NULL) ||
+		((y = BN_CTX_get(ctx)) == NULL) ||
+		((z = BN_CTX_get(ctx)) == NULL) ||
+		((tmp_scalar = BN_CTX_get(ctx)) == NULL))
+		goto err;
+
+	if (scalar != NULL)
+		{
+		pre = EC_EX_DATA_get_data(group->extra_data,
+			nistp224_pre_comp_dup, nistp224_pre_comp_free,
+			nistp224_pre_comp_clear_free);
+		if (pre)
+			/* we have precomputation, try to use it */
+			g_pre_comp = pre->g_pre_comp;
+		else
+			/* try to use the standard precomputation */
+			g_pre_comp = (fslice (*)[3][4]) gmul;
+		generator = EC_POINT_new(group);
+		if (generator == NULL)
+			goto err;
+		/* get the generator from precomputation */
+		if (!felem_to_BN(x, g_pre_comp[1][0]) ||
+			!felem_to_BN(y, g_pre_comp[1][1]) ||
+			!felem_to_BN(z, g_pre_comp[1][2]))
+			{
+			ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
+			goto err;
+			}
+		if (!EC_POINT_set_Jprojective_coordinates_GFp(group,
+				generator, x, y, z, ctx))
+			goto err;
+		if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
+			/* precomputation matches generator */
+			have_pre_comp = 1;
+		else
+			/* we don't have valid precomputation:
+			 * treat the generator as a random point */
+			num_points = num_points + 1;
+		}
+	secrets = OPENSSL_malloc(num_points * sizeof(felem_bytearray));
+	pre_comp = OPENSSL_malloc(num_points * 16 * 3 * 4 * sizeof(fslice));
+
+	if ((num_points) && ((secrets == NULL) || (pre_comp == NULL)))
+		{
+		ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_MALLOC_FAILURE);
+		goto err;
+		}
+
+	/* we treat NULL scalars as 0, and NULL points as points at infinity,
+	 * i.e., they contribute nothing to the linear combination */
+	memset(secrets, 0, num_points * sizeof(felem_bytearray));
+	memset(pre_comp, 0, num_points * 16 * 3 * 4 * sizeof(fslice));
+	for (i = 0; i < num_points; ++i)
+		{
+		if (i == num)
+			/* the generator */
+			{
+			p = EC_GROUP_get0_generator(group);
+			p_scalar = scalar;
+			}
+		else
+			/* the i^th point */
+			{
+			p = points[i];
+			p_scalar = scalars[i];
+			}
+		if ((p_scalar != NULL) && (p != NULL))
+			{
+			num_bytes = BN_num_bytes(p_scalar);
+			/* reduce scalar to 0 <= scalar < 2^224 */
+			if ((num_bytes > sizeof(felem_bytearray)) || (BN_is_negative(p_scalar)))
+				{
+				/* this is an unusual input, and we don't guarantee
+				 * constant-timeness */
+				if (!BN_nnmod(tmp_scalar, p_scalar, &group->order, ctx))
+					{
+					ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
+					goto err;
+					}
+				num_bytes = BN_bn2bin(tmp_scalar, tmp);
+				}
+			else
+				BN_bn2bin(p_scalar, tmp);
+			flip_endian(secrets[i], tmp, num_bytes);
+			/* precompute multiples */
+			if ((!BN_to_felem(x_out, &p->X)) ||
+				(!BN_to_felem(y_out, &p->Y)) ||
+				(!BN_to_felem(z_out, &p->Z))) goto err;
+			memcpy(pre_comp[i][1][0], x_out, 4 * sizeof(fslice));
+			memcpy(pre_comp[i][1][1], y_out, 4 * sizeof(fslice));
+			memcpy(pre_comp[i][1][2], z_out, 4 * sizeof(fslice));
+			for (j = 1; j < 8; ++j)
+				{
+				point_double(pre_comp[i][2*j][0],
+					pre_comp[i][2*j][1],
+					pre_comp[i][2*j][2],
+					pre_comp[i][j][0],
+					pre_comp[i][j][1],
+					pre_comp[i][j][2]);
+				point_add(pre_comp[i][2*j+1][0],
+					pre_comp[i][2*j+1][1],
+					pre_comp[i][2*j+1][2],
+					pre_comp[i][1][0],
+					pre_comp[i][1][1],
+					pre_comp[i][1][2],
+					pre_comp[i][2*j][0],
+					pre_comp[i][2*j][1],
+					pre_comp[i][2*j][2]);
+				}
+			}
+		}
+
+	/* the scalar for the generator */
+	if ((scalar != NULL) && (have_pre_comp))
+		{
+		memset(g_secret, 0, sizeof g_secret);
+		num_bytes = BN_num_bytes(scalar);
+		/* reduce scalar to 0 <= scalar < 2^224 */
+		if ((num_bytes > sizeof(felem_bytearray)) || (BN_is_negative(scalar)))
+			{
+			/* this is an unusual input, and we don't guarantee
+			 * constant-timeness */
+			if (!BN_nnmod(tmp_scalar, scalar, &group->order, ctx))
+				{
+				ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
+				goto err;
+				}
+			num_bytes = BN_bn2bin(tmp_scalar, tmp);
+			}
+		else
+			BN_bn2bin(scalar, tmp);
+		flip_endian(g_secret, tmp, num_bytes);
+		/* do the multiplication with generator precomputation*/
+		batch_mul(x_out, y_out, z_out,
+			(const felem_bytearray (*)) secrets, num_points,
+			g_secret, (const fslice (*)[16][3][4]) pre_comp,
+			(const fslice (*)[3][4]) g_pre_comp);
+		}
+	else
+		/* do the multiplication without generator precomputation */
+		batch_mul(x_out, y_out, z_out,
+			(const felem_bytearray (*)) secrets, num_points,
+			NULL, (const fslice (*)[16][3][4]) pre_comp, NULL);
+	/* reduce the output to its unique minimal representation */
+	felem_contract(x_in, x_out);
+	felem_contract(y_in, y_out);
+	felem_contract(z_in, z_out);
+	if ((!felem_to_BN(x, x_in)) || (!felem_to_BN(y, y_in)) ||
+		(!felem_to_BN(z, z_in)))
+		{
+		ECerr(EC_F_EC_GFP_NISTP224_POINTS_MUL, ERR_R_BN_LIB);
+		goto err;
+		}
+	ret = EC_POINT_set_Jprojective_coordinates_GFp(group, r, x, y, z, ctx);
+
+err:
+	BN_CTX_end(ctx);
+	if (generator != NULL)
+		EC_POINT_free(generator);
+	if (new_ctx != NULL)
+		BN_CTX_free(new_ctx);
+	if (secrets != NULL)
+		OPENSSL_free(secrets);
+	if (pre_comp != NULL)
+		OPENSSL_free(pre_comp);
+	return ret;
+	}
+
+int ec_GFp_nistp224_precompute_mult(EC_GROUP *group, BN_CTX *ctx)
+	{
+	int ret = 0;
+	NISTP224_PRE_COMP *pre = NULL;
+	int i, j;
+	BN_CTX *new_ctx = NULL;
+	BIGNUM *x, *y;
+	EC_POINT *generator = NULL;
+
+	/* throw away old precomputation */
+	EC_EX_DATA_free_data(&group->extra_data, nistp224_pre_comp_dup,
+		nistp224_pre_comp_free, nistp224_pre_comp_clear_free);
+	if (ctx == NULL)
+		if ((ctx = new_ctx = BN_CTX_new()) == NULL) return 0;
+	BN_CTX_start(ctx);
+	if (((x = BN_CTX_get(ctx)) == NULL) ||
+		((y = BN_CTX_get(ctx)) == NULL))
+		goto err;
+	/* get the generator */
+	if (group->generator == NULL) goto err;
+	generator = EC_POINT_new(group);
+	if (generator == NULL)
+		goto err;
+	BN_bin2bn(nistp224_curve_params[3], sizeof (felem_bytearray), x);
+	BN_bin2bn(nistp224_curve_params[4], sizeof (felem_bytearray), y);
+	if (!EC_POINT_set_affine_coordinates_GFp(group, generator, x, y, ctx))
+		goto err;
+	if ((pre = nistp224_pre_comp_new()) == NULL)
+		goto err;
+	/* if the generator is the standard one, use built-in precomputation */
+	if (0 == EC_POINT_cmp(group, generator, group->generator, ctx))
+		{
+		memcpy(pre->g_pre_comp, gmul, sizeof(pre->g_pre_comp));
+		ret = 1;
+		goto err;
+		}
+	if ((!BN_to_felem(pre->g_pre_comp[1][0], &group->generator->X)) ||
+		(!BN_to_felem(pre->g_pre_comp[1][1], &group->generator->Y)) ||
+		(!BN_to_felem(pre->g_pre_comp[1][2], &group->generator->Z)))
+		goto err;
+	/* compute 2^56*G, 2^112*G, 2^168*G */
+	for (i = 1; i < 5; ++i)
+		{
+		point_double(pre->g_pre_comp[2*i][0], pre->g_pre_comp[2*i][1],
+			pre->g_pre_comp[2*i][2], pre->g_pre_comp[i][0],
+			pre->g_pre_comp[i][1], pre->g_pre_comp[i][2]);
+		for (j = 0; j < 55; ++j)
+			{
+			point_double(pre->g_pre_comp[2*i][0],
+				pre->g_pre_comp[2*i][1],
+				pre->g_pre_comp[2*i][2],
+				pre->g_pre_comp[2*i][0],
+				pre->g_pre_comp[2*i][1],
+				pre->g_pre_comp[2*i][2]);
+			}
+		}
+	/* g_pre_comp[0] is the point at infinity */
+	memset(pre->g_pre_comp[0], 0, sizeof(pre->g_pre_comp[0]));
+	/* the remaining multiples */
+	/* 2^56*G + 2^112*G */
+	point_add(pre->g_pre_comp[6][0], pre->g_pre_comp[6][1],
+		pre->g_pre_comp[6][2], pre->g_pre_comp[4][0],
+		pre->g_pre_comp[4][1], pre->g_pre_comp[4][2],
+		pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
+		pre->g_pre_comp[2][2]);
+	/* 2^56*G + 2^168*G */
+	point_add(pre->g_pre_comp[10][0], pre->g_pre_comp[10][1],
+		pre->g_pre_comp[10][2], pre->g_pre_comp[8][0],
+		pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
+		pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
+		pre->g_pre_comp[2][2]);
+	/* 2^112*G + 2^168*G */
+	point_add(pre->g_pre_comp[12][0], pre->g_pre_comp[12][1],
+		pre->g_pre_comp[12][2], pre->g_pre_comp[8][0],
+		pre->g_pre_comp[8][1], pre->g_pre_comp[8][2],
+		pre->g_pre_comp[4][0], pre->g_pre_comp[4][1],
+		pre->g_pre_comp[4][2]);
+	/* 2^56*G + 2^112*G + 2^168*G */
+	point_add(pre->g_pre_comp[14][0], pre->g_pre_comp[14][1],
+		pre->g_pre_comp[14][2], pre->g_pre_comp[12][0],
+		pre->g_pre_comp[12][1], pre->g_pre_comp[12][2],
+		pre->g_pre_comp[2][0], pre->g_pre_comp[2][1],
+		pre->g_pre_comp[2][2]);
+	for (i = 1; i < 8; ++i)
+		{
+		/* odd multiples: add G */
+		point_add(pre->g_pre_comp[2*i+1][0], pre->g_pre_comp[2*i+1][1],
+			pre->g_pre_comp[2*i+1][2], pre->g_pre_comp[2*i][0],
+			pre->g_pre_comp[2*i][1], pre->g_pre_comp[2*i][2],
+			pre->g_pre_comp[1][0], pre->g_pre_comp[1][1],
+			pre->g_pre_comp[1][2]);
+		}
+
+	if (!EC_EX_DATA_set_data(&group->extra_data, pre, nistp224_pre_comp_dup,
+			nistp224_pre_comp_free, nistp224_pre_comp_clear_free))
+		goto err;
+	ret = 1;
+	pre = NULL;
+ err:
+	BN_CTX_end(ctx);
+	if (generator != NULL)
+		EC_POINT_free(generator);
+	if (new_ctx != NULL)
+		BN_CTX_free(new_ctx);
+	if (pre)
+		nistp224_pre_comp_free(pre);
+	return ret;
+	}
+
+int ec_GFp_nistp224_have_precompute_mult(const EC_GROUP *group)
+	{
+	if (EC_EX_DATA_get_data(group->extra_data, nistp224_pre_comp_dup,
+			nistp224_pre_comp_free, nistp224_pre_comp_clear_free)
+		!= NULL)
+		return 1;
+	else
+		return 0;
+	}
+#endif
+
+#ifdef TESTING
+
+#include <sys/time.h>
+
+static u8 ctoh(char c)
+{
+    if (c >= '0' && c <= '9') return c-'0';
+    if (c >= 'a' && c <= 'f') return c-'a'+10;
+    if (c >= 'A' && c <= 'F') return c-'A'+10;
+    return 0;
+}
+
+static void arg_to_bytearray(felem_bytearray ba, const char *arg)
+{
+    /* Convert the arg, which is a string like "1a2637c8" to a byte
+     * array like 0xc8 0x37 0x26 0x1a. */
+    int size = sizeof(felem_bytearray);
+    int arglen = strlen(arg);
+    int argsize = (arglen+1)/2;
+    const char *argp = arg + arglen;
+    u8 *bap = ba;
+
+    memset(ba, 0, size);
+    if (size < argsize) {
+	fprintf(stderr, "Arg too long: %s\n", arg);
+	exit(1);
+    }
+
+    while (argp > arg+1) {
+	argp -= 2;
+	*bap = (ctoh(argp[0])<<4)|(ctoh(argp[1]));
+	++bap;
+    }
+    if (arglen & 1) {
+	/* Handle the stray top nybble */
+	argp -= 1;
+	*bap = ctoh(argp[0]);
+    }
+}
+
+static void arg_to_coord(coord c, const char *arg)
+{
+    felem_bytearray ba;
+
+    arg_to_bytearray(ba, arg);
+    /* Now convert it to a coord */
+    bin21_to_felem(c, ba);
+}
+
+int main(int argc, char **argv)
+{
+    point infinity, P, Q, P2, PQ;
+    felem_bytearray s;
+    int i;
+    struct timeval st, et;
+    unsigned long el;
+    int niter = 1000;
+
+    memset(infinity, 0, sizeof(infinity));
+    memset(P, 0, sizeof(P));
+    memset(Q, 0, sizeof(Q));
+
+    if (argc != 6) {
+	fprintf(stderr, "Usage: %s Px Py Qx Qy s\n", argv[0]);
+	exit(1);
+    }
+
+    arg_to_coord(P[0], argv[1]);
+    arg_to_coord(P[1], argv[2]);
+    P[2][0] = 1;
+    dump_point("P", P);
+    arg_to_coord(Q[0], argv[3]);
+    arg_to_coord(Q[1], argv[4]);
+    Q[2][0] = 1;
+    dump_point("Q", Q);
+    arg_to_bytearray(s, argv[5]);
+
+    point_double(P2[0], P2[1], P2[2], P[0], P[1], P[2]);
+    affine(P2);
+    point_add(PQ[0], PQ[1], PQ[2], P[0], P[1], P[2], Q[0], Q[1], Q[2]);
+    affine(PQ);
+    dump_point("P2", P2);
+    dump_point("PQ", PQ);
+
+    gettimeofday(&st, NULL);
+    for (i=0;i<niter;++i) {
+	point_mul(P, P, s);
+	affine(P);
+    }
+    gettimeofday(&et, NULL);
+    el = (et.tv_sec-st.tv_sec)*1000000 + (et.tv_usec-st.tv_usec);
+    fprintf(stderr, "%lu / %d = %lu us\n", el, niter, el/niter);
+
+    dump_point("Ps", P);
+
+    return 0;
+}
+#endif
+
+/* Figure out whether there's a point with x-coordinate x on the main
+ * curve.  If not, then there's one on the twist curve.  (There are
+ * actually two, which are negatives of each other; that doesn't
+ * matter.)  Multiply that point by seckey and set out to the
+ * x-coordinate of the result. */
+void ptwist_pointmul(byte out[PTWIST_BYTES], const byte x[PTWIST_BYTES],
+	const byte seckey[PTWIST_BYTES])
+{
+    /* Compute z = x^3 + a*x + b */
+    point P, Q;
+    coord z, r2, Qx;
+    uint128_t tmp[5];
+    int ontwist;
+    static const coord three = { 3, 0, 0 };
+    static const coord b =
+	    { 0x46d320e01dc7d6, 0x486ebc69bad316, 0x4e355e95cafedd };
+
+
+    /* Convert the byte array to a coord */
+    bin21_to_felem(P[0], x);
+
+    /* Compute z = x^3 - 3*x + b */
+    felem_square(tmp, P[0]); felem_reduce(z, tmp);
+    felem_diff64(z, three);
+    felem_mul(tmp, z, P[0]); felem_reduce(z, tmp);
+    felem_sum64(z, b);
+
+    /*
+    dump_coord("z", z);
+    */
+    /* Compute r = P[1] = z ^ ((p+1)/4).  This will be a square root of
+     * z, if one exists. */
+    felem_sqrt(P[1], z);
+    /*
+    dump_coord("r", P[1]);
+    */
+
+    /* Is P[1] a square root of z? */
+    felem_square(tmp, P[1]); felem_diff_128_64(tmp, z); felem_reduce(r2, tmp);
+
+    if (felem_is_zero(r2)) {
+	/* P(x,r) is on the curve */
+	ontwist = 0;
+    } else {
+	/* (-x, r) is on the twist */
+	ontwist = 1;
+	felem_neg(P[0], P[0]);
+    }
+    /*
+    fprintf(stderr, "ontwist = %d\n", ontwist);
+    */
+    memset(P[2], 0, sizeof(coord));
+    P[2][0] = 1;
+
+    /* All set.  Now do the point multiplication. */
+    /*
+    dump_point("P", P);
+    for(i=0;i<21;++i) {
+	fprintf(stderr, "%02x", seckey[20-i]);
+    }
+    fprintf(stderr, "\n");
+    */
+    point_mul(Q, P, seckey);
+    affine_x(Qx, Q);
+    /*
+    dump_point("Q", Q);
+    */
+
+    /* Get the x-coordinate of the result, and negate it if we're on the
+     * twist. */
+    if (ontwist) {
+	felem_neg(Qx, Qx);
+    }
+
+    /* Convert back to bytes */
+    felem_to_bin21(out, Qx);
+    /*
+    fprintf(stderr, "out: ");
+    for(i=0;i<21;++i) {
+	fprintf(stderr, "%02x", out[i]);
+    }
+    fprintf(stderr, "\n");
+    */
+}

tests/pubkey

@@ -0,0 +1 @@
+ぼ<𣧂PZ鷩l
<BS螃𣶹疄俑)𧨾櫥皿金灡��

tests/server.c

@@ -0,0 +1,124 @@
+#include <stdio.h>
+#include <unistd.h>
+#include <sys/socket.h>
+#include <arpa/inet.h>
+#include <openssl/ssl.h>
+#include <openssl/err.h>
+
+int create_socket(int port)
+{
+    int s;
+    struct sockaddr_in addr;
+
+    addr.sin_family = AF_INET;
+    addr.sin_port = htons(port);
+    addr.sin_addr.s_addr = htonl(INADDR_ANY);
+
+    s = socket(AF_INET, SOCK_STREAM, 0);
+    if (s < 0) {
+	perror("Unable to create socket");
+	exit(EXIT_FAILURE);
+    }
+
+    if (bind(s, (struct sockaddr*)&addr, sizeof(addr)) < 0) {
+	perror("Unable to bind");
+	exit(EXIT_FAILURE);
+    }
+
+    if (listen(s, 1) < 0) {
+	perror("Unable to listen");
+	exit(EXIT_FAILURE);
+    }
+
+    return s;
+}
+
+void init_openssl()
+{ 
+    SSL_load_error_strings();	
+    OpenSSL_add_ssl_algorithms();
+}
+
+void cleanup_openssl()
+{
+    EVP_cleanup();
+}
+
+SSL_CTX *create_context()
+{
+    const SSL_METHOD *method;
+    SSL_CTX *ctx;
+
+    method = TLSv1_2_server_method();
+
+    ctx = SSL_CTX_new(method);
+    if (!ctx) {
+	perror("Unable to create SSL context");
+	ERR_print_errors_fp(stderr);
+	exit(EXIT_FAILURE);
+    }
+
+    return ctx;
+}
+
+void configure_context(SSL_CTX *ctx)
+{
+    SSL_CTX_set_ecdh_auto(ctx, 1);
+
+    /* Set the key and cert */
+    if (SSL_CTX_use_certificate_file(ctx, "domain.crt", SSL_FILETYPE_PEM) < 0) {
+        ERR_print_errors_fp(stderr);
+	exit(EXIT_FAILURE);
+    }
+
+    if (SSL_CTX_use_PrivateKey_file(ctx, "domain.key", SSL_FILETYPE_PEM) < 0 ) {
+        ERR_print_errors_fp(stderr);
+	exit(EXIT_FAILURE);
+    }
+    SSL_CTX_set_cipher_list(ctx, "ECDHE-RSA-AES256-GCM-SHA384:DHE-RSA-AES256-GCM-SHA384");
+}
+
+int main(int argc, char **argv)
+{
+    int sock;
+    SSL_CTX *ctx;
+
+    init_openssl();
+    ctx = create_context();
+
+    configure_context(ctx);
+
+    sock = create_socket(8888);
+
+    /* Handle connections */
+    while(1) {
+        struct sockaddr_in addr;
+        uint len = sizeof(addr);
+        SSL *ssl;
+        const char reply[] = "test\n";
+
+        int client = accept(sock, (struct sockaddr*)&addr, &len);
+        if (client < 0) {
+            perror("Unable to accept");
+            exit(EXIT_FAILURE);
+        }
+
+        ssl = SSL_new(ctx);
+        SSL_set_fd(ssl, client);
+
+        if (SSL_accept(ssl) <= 0) {
+            ERR_print_errors_fp(stderr);
+        }
+        else {
+            SSL_write(ssl, reply, strlen(reply));
+        }
+
+        SSL_free(ssl);
+        close(client);
+    }
+
+    close(sock);
+    SSL_CTX_free(ctx);
+    cleanup_openssl();
+}
+

tests/slitheen.c

@@ -0,0 +1,607 @@
+/**
+ * Author: Cecylia Bocovich <cbocovic@uwaterloo.ca>
+ *
+ * This file contains callback functions and all necessary helper functions to
+ * tag flows for use with the Slitheen decoy routing system
+ *
+ */
+
+#include "slitheen.h"
+
+#include <stdio.h>
+#include <string.h>
+#include <semaphore.h>
+#include <openssl/rand.h>
+#include <openssl/sha.h>
+#include <openssl/ssl.h>
+#include <openssl/err.h>
+#include <openssl/ec.h>
+#include <openssl/bn.h>
+
+#include "ptwist.h"
+#include "crypto.h"
+
+tag_pair *keys;
+
+//look for key list
+sem_t key_lock;
+
+byte maingen[PTWIST_BYTES];
+byte twistgen[PTWIST_BYTES];
+byte mainpub[PTWIST_BYTES];
+byte twistpub[PTWIST_BYTES];
+
+//called once when phantomjs loads
+void slitheen_init(){
+	int test;
+
+	if(sem_init(&key_lock, 0, 1) == -1){
+		printf("Initialization of semaphore failed\n");
+		exit(1);
+	}
+
+	sem_getvalue(&key_lock, &test);
+}
+
+static void gen_tag(byte *tag, byte stored_key[16],
+                    const byte *context, size_t context_len){
+    byte seckey[PTWIST_BYTES];
+    byte sharedsec[PTWIST_BYTES+context_len];
+    byte usetwist;
+    byte taghashout[32];
+#if PTWIST_PUZZLE_STRENGTH > 0
+    size_t puzzle_len = 16+PTWIST_RESP_BYTES;
+    byte value_to_hash[puzzle_len];
+    byte hashout[32];
+    bn_t Rbn, Hbn;
+    int i, len, sign;
+#endif
+
+    memset(tag, 0xAA, PTWIST_TAG_BYTES);
+    memset(stored_key, 0, 16);
+
+    /* Use the main or the twist curve? */
+    RAND_bytes(&usetwist, 1);
+    usetwist &= 1;
+
+    /* Create seckey*G and seckey*Y */
+    RAND_bytes(seckey, PTWIST_BYTES);
+    ptwist_pointmul(tag, usetwist ? twistgen : maingen, seckey);
+    ptwist_pointmul(sharedsec, usetwist ? twistpub : mainpub, seckey);
+
+    /* Create the tag hash keys */
+    memmove(sharedsec+PTWIST_BYTES, context, context_len);
+    SHA256(sharedsec, PTWIST_BYTES, taghashout);
+
+#if PTWIST_PUZZLE_STRENGTH > 0
+    /* The puzzle is to find a response R such that SHA256(K || R)
+       starts with PTWIST_PUZZLE_STRENGTH bits of 0s.  K is the first
+       128 bits of the above hash tag keys. */
+
+    /* Construct our response to the puzzle.  Start looking for R in a
+     * random place. */
+    memmove(value_to_hash, taghashout, 16);
+    RAND_bytes(value_to_hash+16, PTWIST_RESP_BYTES);
+    value_to_hash[16+PTWIST_RESP_BYTES-1] &= PTWIST_RESP_MASK;
+
+    while(1) {
+	unsigned int firstbits;
+
+	md_map_sh256(hashout, value_to_hash, puzzle_len);
+#if PTWIST_PUZZLE_STRENGTH < 32
+	/* This assumes that you're on an architecture that doesn't care
+	 * about alignment, and is little endian. */
+	firstbits = *(unsigned int*)hashout;
+	if ((firstbits & PTWIST_PUZZLE_MASK) == 0) {
+	    break;
+	}
+	/* Increment R and try again. */
+	for(i=0;i<PTWIST_RESP_BYTES;++i) {
+	    if (++value_to_hash[16+i]) break;
+	}
+	value_to_hash[16+PTWIST_RESP_BYTES-1] &= PTWIST_RESP_MASK;
+#else
+#error "Code assumes PTWIST_PUZZLE_STRENGTH < 32"
+#endif
+    }
+
+	/*
+	for(i=0;i<puzzle_len;++i) {
+	    printf("%02x", value_to_hash[i]);
+	    if ((i%4) == 3) printf(" ");
+	}
+	printf("\n");
+	for(i=0;i<32;++i) {
+	    printf("%02x", hashout[i]);
+	    if ((i%4) == 3) printf(" ");
+	}
+	printf("\n");
+	*/
+    /* When we get here, we have solved the puzzle.  R is in
+     * value_to_hash[16..16+PTWIST_RESP_BYTES-1], the hash output
+     * hashout starts with PTWIST_PUZZLE_STRENGTH bits of 0s, and we'll
+     * want to copy out H (the next PTWIST_HASH_SHOWBITS bits of the
+     * hash output).  The final tag is [seckey*G]_x || R || H . */
+    bn_new(Rbn);
+    bn_new(Hbn);
+
+    bn_read_bin(Rbn, value_to_hash+16, PTWIST_RESP_BYTES, BN_POS);
+    hashout[PTWIST_HASH_TOTBYTES-1] &= PTWIST_HASH_MASK;
+    bn_read_bin(Hbn, hashout, PTWIST_HASH_TOTBYTES, BN_POS);
+    bn_lsh(Hbn, Hbn, PTWIST_RESP_BITS-PTWIST_PUZZLE_STRENGTH);
+    bn_add(Hbn, Hbn, Rbn);
+    len = PTWIST_TAG_BYTES-PTWIST_BYTES;
+    bn_write_bin(tag+PTWIST_BYTES, &len, &sign, Hbn);
+	/*
+	for(i=0;i<PTWIST_TAG_BYTES;++i) {
+	    printf("%02x", tag[i]);
+	    if ((i%4) == 3) printf(" ");
+	}
+	printf("\n");
+	*/
+
+    bn_free(Rbn);
+    bn_free(Hbn);
+#elif PTWIST_HASH_SHOWBITS <= 128
+    /* We're not using a client puzzle, so the tag is [seckey*G]_x || H
+     * where H is the first PTWIST_HASH_SHOWBITS bits of the above hash
+     * output.  The key generated is the last 128 bits of that output.
+     * If there's no client puzzle, PTWIST_HASH_SHOWBITS must be a multiple
+     * of 8. */
+    memmove(tag+PTWIST_BYTES, taghashout, PTWIST_HASH_SHOWBITS/8);
+#else
+#error "No client puzzle used, but PWTIST_HASH_SHOWBITS > 128"
+#endif
+
+    memmove(stored_key, taghashout+16, 16);
+}
+
+int tag_hello(unsigned char *target, byte stored_key[16]){
+    FILE *fp;
+    int res;
+    byte *tag;
+
+    /* Create the generators */
+    memset(maingen, 0, PTWIST_BYTES);
+    maingen[0] = 2;
+    memset(twistgen, 0, PTWIST_BYTES);
+
+
+    /* Read the public keys */
+    fp = fopen("pubkey", "rb");
+    if (fp == NULL) {
+		perror("fopen");
+		exit(1);
+    }
+    res = fread(mainpub, PTWIST_BYTES, 1, fp);
+    if (res < 1) {
+		perror("fread");
+		exit(1);
+    }
+    res = fread(twistpub, PTWIST_BYTES, 1, fp);
+    if (res < 1) {
+		perror("fread");
+		exit(1);
+    }
+    fclose(fp);
+
+    tag = target;
+
+    gen_tag(tag, stored_key, (const byte *)"context", 7);
+
+    return 0;
+}
+
+//Client hello callback
+int slitheen_tag_hello(SSL *s){
+    unsigned char *result;
+    int len;
+
+
+    result = s->s3->client_random;
+    len = sizeof(s->s3->client_random);
+
+    if(len < PTWIST_TAG_BYTES) {
+            printf("Uhoh\n");
+            return 1;
+    }
+    unsigned long Time = (unsigned long)time(NULL);
+    unsigned char *p = result;
+    l2n(Time, p);
+
+    //
+    tag_pair *new_pair = calloc(1, sizeof(tag_pair));
+    tag_hello((byte *) result+4, new_pair->key);
+
+    new_pair->next = NULL;
+    memcpy(new_pair->client_random, s->s3->client_random, SSL3_RANDOM_SIZE);
+
+    int test;
+    sem_getvalue(&key_lock, &test);
+
+    sem_wait(&key_lock);
+    tag_pair *last_pair;
+    if(keys == NULL){
+        keys = new_pair;
+    } else {
+        last_pair = keys;
+        while(last_pair->next != NULL){
+            last_pair = last_pair->next;
+        }
+        last_pair->next = new_pair;
+    }
+    sem_post(&key_lock);
+
+    return 0;
+}
+
+//dh callback
+int slitheen_seed_from_tag(SSL *s, DH *dh)
+{
+    int ok = 0;
+    int generate_new_key = 0;
+    unsigned l;
+    BN_CTX *ctx;
+    BN_MONT_CTX *mont = NULL;
+    BIGNUM *pub_key = NULL, *priv_key= NULL;
+    unsigned char *buf = NULL, *seed = NULL;
+    int bytes = 0;
+    byte key[16];
+
+    //find key from keys list
+    sem_wait(&key_lock);
+    tag_pair *pair = keys;
+    while(pair != NULL){
+        if(!memcmp(pair->client_random, s->s3->client_random, SSL3_RANDOM_SIZE)){
+            memcpy(key, pair->key, 16);
+            break;
+        }
+        pair = pair->next;
+    }
+    if(pair == NULL){
+        printf("ERROR: KEY NOT FOUND\n");
+            sem_post(&key_lock);
+        return 1;
+    }
+    sem_post(&key_lock);
+
+    seed = (unsigned char *) key;
+
+    ctx = BN_CTX_new();
+    if (ctx == NULL)
+        goto err;
+
+    if (dh->priv_key == NULL) {
+        priv_key = BN_new();
+        if (priv_key == NULL)
+            goto err;
+        generate_new_key = 1;
+    } else
+        priv_key = dh->priv_key;
+
+    if (dh->pub_key == NULL) {
+        pub_key = BN_new();
+        if (pub_key == NULL)
+            goto err;
+    } else
+        pub_key = dh->pub_key;
+
+    if (dh->flags & DH_FLAG_CACHE_MONT_P) {
+        mont = BN_MONT_CTX_set_locked(&dh->method_mont_p,
+                                      CRYPTO_LOCK_DH, dh->p, ctx);
+        if (!mont)
+            goto err;
+    }
+
+    if (generate_new_key) {
+	/* secret exponent length */
+	l = dh->length ? dh->length : BN_num_bits(dh->p) - 1;
+	bytes = (l+7) / 8;
+
+	/* set exponent to seeded prg value */
+	buf = (unsigned char *)OPENSSL_malloc(bytes);
+	if (buf == NULL){
+	    BNerr(BN_F_BNRAND, ERR_R_MALLOC_FAILURE);
+	    goto err;
+	}
+
+    PRF(seed, 16,
+        (uint8_t *) SLITHEEN_KEYGEN_CONST, SLITHEEN_KEYGEN_CONST_SIZE,
+        NULL, 0, NULL, 0, NULL, 0,
+        buf, bytes);
+
+#ifdef DEBUG
+    printf("Generated the following rand bytes: ");
+    for(i=0; i< bytes; i++){
+        printf(" %02x ", buf[i]);
+    }
+    printf("\n");
+#endif
+
+	if (!BN_bin2bn(buf, bytes, priv_key))
+	    goto err;
+
+    }
+
+    {
+        BIGNUM local_prk;
+        BIGNUM *prk;
+
+        if ((dh->flags & DH_FLAG_NO_EXP_CONSTTIME) == 0) {
+            BN_init(&local_prk);
+            prk = &local_prk;
+            BN_with_flags(prk, priv_key, BN_FLG_CONSTTIME);
+        } else
+            prk = priv_key;
+
+        if (!dh->meth->bn_mod_exp(dh, pub_key, dh->g, prk, dh->p, ctx, mont))
+            goto err;
+    }
+
+    dh->pub_key = pub_key;
+    dh->priv_key = priv_key;
+    ok = 1;
+ err:
+    if (buf != NULL){
+		OPENSSL_cleanse(buf, bytes);
+		OPENSSL_free(buf);
+    }
+    if (ok != 1)
+        DHerr(DH_F_GENERATE_KEY, ERR_R_BN_LIB);
+
+    if ((pub_key != NULL) && (dh->pub_key == NULL))
+        BN_free(pub_key);
+    if ((priv_key != NULL) && (dh->priv_key == NULL))
+        BN_free(priv_key);
+    BN_CTX_free(ctx);
+    return (ok);
+}
+
+int slitheen_ec_seed_from_tag(SSL *s, EC_KEY *eckey){
+    int ok = 0;
+    BN_CTX *ctx = NULL;
+    BIGNUM *priv_key = NULL, *order = NULL;
+    EC_POINT *pub_key = NULL;
+	unsigned l;
+    unsigned char *buf = NULL, *seed = NULL;
+    int bytes = 0;
+    byte key[16];
+
+    //find key from keys list
+    sem_wait(&key_lock);
+    tag_pair *pair = keys;
+    while(pair != NULL){
+        if(!memcmp(pair->client_random, s->s3->client_random, SSL3_RANDOM_SIZE)){
+            memcpy(key, pair->key, 16);
+            break;
+        }
+        pair = pair->next;
+    }
+    if(pair == NULL){
+        printf("ERROR: KEY NOT FOUND\n");
+		sem_post(&key_lock);
+        return 1;
+    }
+    sem_post(&key_lock);
+
+#ifdef DEBUG
+    printf("IN SLITHEEN EC GENERATE CALLBACK (key =");
+    for(i=0; i< 16; i++){
+        printf(" %02x", key[i]);
+    }
+    printf(")\n");
+#endif
+
+    seed = (unsigned char *) key;
+
+    if (!eckey || !EC_KEY_get0_group(eckey)) {
+        ECerr(EC_F_EC_KEY_GENERATE_KEY, ERR_R_PASSED_NULL_PARAMETER);
+        return 0;
+    }
+
+    if ((order = BN_new()) == NULL)
+        goto err;
+    if ((ctx = BN_CTX_new()) == NULL)
+        goto err;
+
+    if (EC_KEY_get0_private_key(eckey) == NULL) {
+        priv_key = BN_new();
+        if (priv_key == NULL)
+            goto err;
+    } else
+        priv_key = EC_KEY_get0_private_key(eckey);
+
+    if (!EC_GROUP_get_order(EC_KEY_get0_group(eckey), order, ctx))
+        goto err;
+
+	/* secret exponent length */
+	l = BN_num_bits(order) - 1;
+	bytes = (l+7) / 8;
+
+	/* set exponent to seeded prg value */
+	buf = (unsigned char *)OPENSSL_malloc(bytes);
+	if (buf == NULL){
+	    BNerr(BN_F_BNRAND, ERR_R_MALLOC_FAILURE);
+	    goto err;
+	}
+
+    PRF(seed, 16,
+        (uint8_t *) SLITHEEN_KEYGEN_CONST, SLITHEEN_KEYGEN_CONST_SIZE,
+        NULL, 0, NULL, 0, NULL, 0,
+        buf, bytes);
+
+#ifdef DEBUG
+    printf("Generated the following rand bytes: ");
+    for(i=0; i< bytes; i++){
+        printf(" %02x ", buf[i]);
+    }
+    printf("\n");
+#endif
+
+	if (!BN_bin2bn(buf, bytes, priv_key))
+	    goto err;
+
+    if (EC_KEY_get0_public_key(eckey) == NULL) {
+        pub_key = EC_POINT_new(EC_KEY_get0_group(eckey));
+        if (pub_key == NULL)
+            goto err;
+    } else
+        pub_key = EC_KEY_get0_public_key(eckey);
+
+    if (!EC_POINT_mul(EC_KEY_get0_group(eckey), pub_key, priv_key, NULL, NULL, ctx))
+        goto err;
+
+    EC_KEY_set_private_key(eckey, priv_key);
+    EC_KEY_set_public_key(eckey, pub_key);
+
+    ok = 1;
+
+ err:
+    if (buf != NULL){
+		OPENSSL_cleanse(buf, bytes);
+		OPENSSL_free(buf);
+    }
+    if (order)
+        BN_free(order);
+    if (pub_key != NULL && EC_KEY_get0_public_key(eckey) == NULL)
+        EC_POINT_free(pub_key);
+    if (priv_key != NULL && EC_KEY_get0_private_key(eckey) == NULL)
+        BN_free(priv_key);
+    if (ctx != NULL)
+        BN_CTX_free(ctx);
+    return (ok);
+
+}
+
+/* Finished_mac_callback
+ *
+ * This function checks the MAC both against the expected input, and additionally
+ * with an extra shared-secret-based input to the MAC. If it passes the former, a
+ * warning flag is set to indicate non-usage of the decoy routing protocol
+ *
+ * Returns 1 (ordinary) or 2 (modified) on success, 0 on failure
+ */
+int slitheen_finished_mac(SSL *ssl, unsigned char *finished_mac){
+
+    EVP_MD_CTX ctx;
+    uint8_t output[2*EVP_MAX_MD_SIZE];
+    uint32_t output_len;
+    int i;
+
+    uint8_t *modified_mac, *extra_input;
+    uint32_t extra_input_len;
+    byte key[16];
+
+    EVP_MD_CTX_init(&ctx);
+
+    modified_mac = calloc(1, ssl->s3->tmp.peer_finish_md_len);
+
+    //find shared secret from keys list
+    sem_wait(&key_lock);
+    tag_pair *pair = keys;
+    tag_pair *prev = NULL;
+    while(pair != NULL){
+        if(!memcmp(pair->client_random, ssl->s3->client_random, SSL3_RANDOM_SIZE)){
+            memcpy(key, pair->key, 16);
+            break;
+        }
+        prev = pair;
+        pair = pair->next;
+    }
+    if(pair == NULL){
+        printf("ERROR: KEY NOT FOUND\n");
+		sem_post(&key_lock);
+        return 1;
+    }
+    sem_post(&key_lock);
+
+    //remove that key from the keys list: we are now done with it
+    sem_wait(&key_lock);
+    if (prev == NULL){
+        keys = keys->next;
+        free(pair);
+    } else {
+        prev->next = pair->next;
+        free(pair);
+    }
+    sem_post(&key_lock);
+
+
+    //compute extra input to finished from shared secret
+    extra_input_len = SSL3_RANDOM_SIZE;
+    extra_input = calloc(1, extra_input_len);
+
+    PRF(ssl->session->master_key, ssl->session->master_key_length,
+        (uint8_t *) SLITHEEN_FINISHED_INPUT_CONST, SLITHEEN_FINISHED_INPUT_CONST_SIZE,
+        NULL, 0, NULL, 0, NULL, 0,
+        extra_input, extra_input_len);
+
+    //compute modified finish message from relay station
+    EVP_MD_CTX *hdgst = NULL;
+    for(i=0; i< SSL_MAX_DIGEST; i++){
+        if(ssl->s3->handshake_dgst[i]){
+            hdgst = ssl->s3->handshake_dgst[i]; //there are many of these, find right one(s)?
+        }
+    }
+    if(!hdgst){
+        printf("Could not find digest, skipping extra check\n");
+    } else {
+
+        EVP_MD_CTX_copy_ex(&ctx, hdgst);
+        EVP_DigestUpdate(&ctx, extra_input, extra_input_len);
+
+        EVP_DigestFinal_ex(&ctx, output, &output_len);
+
+        PRF(ssl->session->master_key, ssl->session->master_key_length,
+            (uint8_t *) TLS_MD_SERVER_FINISH_CONST, TLS_MD_SERVER_FINISH_CONST_SIZE,
+            output, output_len, NULL, 0, NULL, 0,
+            modified_mac, ssl->s3->tmp.peer_finish_md_len);
+    }
+
+#ifdef DEBUG
+    printf("modified mac:\n");
+    for(i=0; i<ssl->s3->tmp.peer_finish_md_len; i++){
+        printf("%02x ", modified_mac[i]);
+    }
+    printf("\n");
+#endif
+
+    //now add unmodified Finish message to the finish mac computation
+    uint8_t *unmodified_finish = malloc(ssl->s3->tmp.peer_finish_md_len+4);
+    memcpy(unmodified_finish, ssl->init_buf->data, 4);//copy header
+    memcpy(unmodified_finish+4, ssl->s3->tmp.peer_finish_md, ssl->s3->tmp.peer_finish_md_len);
+
+	printf("Unmodified hash:\n");
+	for(i=0; i< ssl->s3->tmp.peer_finish_md_len; i++){
+		printf("%02x ", ssl->s3->tmp.peer_finish_md[i]);
+	}
+	printf("\n");
+
+    ssl3_finish_mac(ssl, unmodified_finish, ssl->s3->tmp.peer_finish_md_len + 4);
+
+    free(unmodified_finish);
+
+    int32_t retval = 0;
+
+    //Compare MAC to what it should be
+    if(CRYPTO_memcmp(finished_mac, ssl->s3->tmp.peer_finish_md, ssl->s3->tmp.peer_finish_md_len) != 0){
+        //check to see if we have the modified MAC instead.
+        if(CRYPTO_memcmp(finished_mac, modified_mac, ssl->s3->tmp.peer_finish_md_len) != 0){
+            printf("MAC unknown\n");
+        } else {
+            retval = 2;
+            printf("MAC was correctly modified!\n");
+        }
+    } else {
+        retval = 1;
+        printf("MAC was ordinary\n");
+    }
+
+    //clean up
+    EVP_MD_CTX_cleanup(&ctx);
+    OPENSSL_cleanse(extra_input, extra_input_len);
+    OPENSSL_cleanse(output, output_len);
+
+    return !retval;
+
+}

tests/slitheen.h

@@ -0,0 +1,35 @@
+#ifndef _SLITHEEN_H_
+#define _SLITHEEN_H_
+
+#include <openssl/ssl.h>
+#include "ptwist.h"
+
+# define l2n(l,c)        (*((c)++)=(unsigned char)(((l)>>24)&0xff), \
+                         *((c)++)=(unsigned char)(((l)>>16)&0xff), \
+                         *((c)++)=(unsigned char)(((l)>> 8)&0xff), \
+                         *((c)++)=(unsigned char)(((l)    )&0xff))
+
+void slitheen_init();
+
+int slitheen_tag_hello(SSL *s);
+
+int slitheen_seed_from_tag(SSL *s, DH *dh);
+
+int slitheen_ec_seed_from_tag(SSL *s, EC_KEY *eckey);
+
+int slitheen_finished_mac(SSL *ssl, unsigned char *finished_mac);
+
+typedef struct tag_pair_st {
+    byte key[16];
+    u_char client_random[SSL3_RANDOM_SIZE];
+    struct tag_pair_st *next;
+} tag_pair;
+
+
+#define SLITHEEN_KEYGEN_CONST "SLITHEEN_KEYGEN"
+#define SLITHEEN_KEYGEN_CONST_SIZE 15
+
+#define SLITHEEN_FINISHED_INPUT_CONST "SLITHEEN_FINISH"
+#define SLITHEEN_FINISHED_INPUT_CONST_SIZE 15
+
+#endif /* _SLITHEEN_H_ */

tests/smtpClient.c

@@ -0,0 +1,181 @@
+//
+// smtpClient.c - example SMTP Client shell to send HELO to SMTP server
+//
+//   compile using: gcc -o smtpClient smtpClient.c
+//   run using: ./smtpClient <mail host> [ port ]
+//   example:   ./smtpClient mailhost.ksu.edu 25
+//
+#include <stdio.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/socket.h>
+#include <netinet/in.h>
+#include <arpa/inet.h>
+#include <netdb.h>      // for struct hostent
+
+#include <openssl/ssl.h>
+#include <openssl/err.h>
+
+#include "slitheen.h"
+#define MAX_REQUEST 1024
+#define MAX_REPLY BUFSIZ
+int getLine(int fd, char line[], int max);
+
+void apps_ssl_info_callback( SSL *s, int where, int ret )
+{
+    char *str;
+    int w;
+    
+    w=where& ~SSL_ST_MASK;
+    
+    if (w & SSL_ST_CONNECT) str="SSL_connect";
+    else if (w & SSL_ST_ACCEPT) str="SSL_accept";
+    else str="undefined";
+    
+    if (where & SSL_CB_LOOP)
+    {
+	fprintf(stderr,"%s: %s\n",str,SSL_state_string_long(s));
+    }
+    else if (where & SSL_CB_ALERT)
+    {
+	str=(where & SSL_CB_READ)?"read":"write";
+	fprintf(stderr,"SSL3 alert %s:%s:%s\n",
+		   str,
+		   SSL_alert_type_string_long(ret),
+		   SSL_alert_desc_string_long(ret));
+    }
+    else if (where & SSL_CB_EXIT)
+    {
+	if (ret == 0)
+	    fprintf(stderr,"%s:failed in %s\n",
+		       str,SSL_state_string_long(s));
+	else if (ret < 0)
+	{
+	    fprintf(stderr,"%s:error in %s\n",
+		       str,SSL_state_string_long(s));
+	}
+    }
+}
+
+main(int argc, char **argv)
+{
+    int                 sockfd;
+    struct sockaddr_in  serv_addr;
+    char                request[MAX_REQUEST+1];
+    char                reply[MAX_REPLY+1];
+    unsigned int	server_port = 25;
+    struct hostent	*hostptr;
+    struct in_addr      *ptr;
+    unsigned short	port_number;
+    char userinput[801]; //MAX_IN is initialized to 800
+
+    SSL_CTX *ctx = NULL;
+    SSL *session = NULL;
+
+	//TLS setup
+    const SSL_METHOD *meth = TLSv1_2_client_method();
+    //const SSL_METHOD *meth = SSLv23_client_method();
+    if (meth == NULL) {
+		fprintf( stderr, "no method. :(\n" ); exit(1);
+    }
+
+    SSL_load_error_strings();	
+    OpenSSL_add_ssl_algorithms();
+
+    ctx = SSL_CTX_new( meth );
+    if ( ctx == NULL ) { fprintf( stderr, "no context. :(\n" ); exit(1);
+		ERR_print_errors_fp(stderr);
+	}
+
+    //SSL_CTX_set_cipher_list(ctx, "ECDHE-RSA-AES256-GCM-SHA384:DHE-RSA-AES256-GCM-SHA384");
+
+	//set Slitheen callbacks
+    slitheen_init();
+
+    SSL_CTX_set_client_hello_callback(ctx, slitheen_tag_hello);
+    //SSL_CTX_set_generate_ec_key_callback(ctx, slitheen_ec_seed_from_tag);
+    //SSL_CTX_set_generate_key_callback(ctx, slitheen_seed_from_tag);
+    //SSL_CTX_set_finished_mac_callback(ctx, slitheen_finished_mac);
+
+    /* Set up a callback for each state change so we can see what's
+     * going on */
+    SSL_CTX_set_info_callback(ctx,apps_ssl_info_callback);
+
+
+    // Read hostname and port from the command-line
+    if (argc < 2)
+    {
+        printf("Usage: smtpClient <hostname> [port-number]\n");
+        return(1);
+    }
+    if (argc > 2)
+        port_number = atoi(argv[2]);
+    else
+        port_number = 25;
+
+    if ( (hostptr = (struct hostent *) gethostbyname(argv[1])) == NULL) {
+        perror("gethostbyname error for host");
+	return(1);
+    }
+    ptr = (struct in_addr *) *(hostptr->h_addr_list);
+    printf ("DEBUG: server address: %u %s\n",ptr->s_addr,inet_ntoa(*ptr));
+    memset((char *) &serv_addr, 0, sizeof(serv_addr));
+    serv_addr.sin_family      = AF_INET;
+    serv_addr.sin_addr.s_addr = ptr->s_addr;
+    serv_addr.sin_port        = htons(port_number);
+
+    // Create communication endpoint.
+    if ( (sockfd = socket(AF_INET, SOCK_STREAM, 0)) < 0) {
+        perror("server: can't open stream socket");
+        return(1);
+    }
+
+    // Connect to the server.
+    if (connect(sockfd, (struct sockaddr *) &serv_addr, sizeof(serv_addr)) < 0) {
+        perror("client: can't connect to server");
+        return(1);
+    }
+    getLine(sockfd, reply, MAX_REPLY);
+    puts(reply);
+    printf("SENDING HELO %s\n", argv[1]);
+    sprintf(request,"HELO %s\r\n", argv[1]);
+    write(sockfd, request, strlen(request));
+    read(sockfd, reply, BUFSIZ);
+    puts(reply);
+    printf("SENDING STARTTLS\n");
+    sprintf(request, "STARTTLS\r\n");
+    write(sockfd, request, strlen(request));
+    getLine(sockfd, reply, MAX_REPLY);
+    puts(reply);
+
+
+	//set up ssl socket
+	if (sockfd != -1){
+		session = SSL_new( ctx );
+		SSL_set_fd( session, sockfd );
+		SSL_connect( session );
+	}
+		
+    close(sockfd);
+    return(0);
+}
+
+//
+// get a line of data from fd
+//
+int getLine(int fd, char* line, int lim)
+{
+    int i;
+    char c;
+
+    i =  0;
+    while (--lim > 0 && read(fd,&c,1)>0 && c!='\n' && c!='\0')
+    {
+        line[i++] = c;
+    }
+    if (c=='\n')
+        line[i++] = c;
+    line[i] = '\0';
+    return i;
+}
+