cleanup

P503

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "csidh/faster-csidh"]
+	path = csidh/faster-csidh
+	url = https://gitlab.cs.hs-rm.de/pqcrypto/faster-csidh.git
+[submodule "csidh/ref/csidh"]
+	path = csidh/ref/csidh
+	url = https://git.amongbytes.com/kris/csidh.git

README.md

@@ -0,0 +1,3 @@
+Works only against SIKE from MSR and CF and uses P751 only
+
+**work ongoing**

csidh/Makefile

@@ -0,0 +1,59 @@
+objify = $(patsubst %.c,$(BUILD_DIR)/%.$2,$(patsubst %.s,$(BUILD_DIR)/%.$2,$1))
+
+BUILD_DIR=build
+BIN_DIR=bin
+MK_FILE_PATH = $(lastword $(MAKEFILE_LIST))
+PRJ_DIR      = $(abspath $(dir $(MK_FILE_PATH)))
+GOPATH_LOCAL = $(PRJ_DIR)/$(BUILD_DIR)/
+
+CFLAGS=-Wall -Wextra -Wpedantic -O3 -funroll-loops
+LDFLAGS=
+AR=ar rcs
+RANLIB=ranlib
+GOPATH=$(BUILD_DIR)
+GO=go
+IMPORT_PATH:=github.com/henrydcase/sidh_torture/csidh
+
+CODE_SRC_C = \
+	ref/csidh/rng.c \
+	ref/csidh/mont.c \
+	ref/csidh/csidh.c
+CODE_SRC_S = \
+	ref/csidh/u512.s \
+	ref/csidh/fp.s
+CODE_OBJ = \
+ 	$(call objify,$(CODE_SRC_C),o) \
+ 	$(call objify,$(CODE_SRC_S),o)
+
+$(BUILD_DIR)/%.o: %.c
+	case $@ in */*) f=$@; mkdir -p $${f%/*} ;; esac
+	$(CC) -c -o $@ $< $(CFLAGS)
+
+$(BUILD_DIR)/%.o: %.s
+	case $@ in */*) f=$@; mkdir -p $${f%/*} ;; esac
+	$(CC) -c -o $@ $< $(CFLAGS)
+
+all: libcsidh torturer vector-creator
+
+libcsidh: $(CODE_OBJ)
+	mkdir -p $(PRJ_DIR)/$(BIN_DIR)
+	$(AR) $(BIN_DIR)/libcsidh.a $^
+	$(RANLIB) $(BIN_DIR)/libcsidh.a
+	$(CC) -o $(BIN_DIR)/test_ref ref/test/main.c -L$(BIN_DIR) -lcsidh
+
+clean:
+	rm -rf build
+	rm -rf bin
+
+torturer: $(BUILD_DIR)/.ok
+	GOPATH=$(GOPATH_LOCAL) $(GO) install $(IMPORT_PATH)/cmd/$@
+
+vector-creator: $(BUILD_DIR)/.ok
+	GOPATH=$(GOPATH_LOCAL) $(GO) install $(IMPORT_PATH)/cmd/$@
+
+$(BUILD_DIR)/.ok:
+	mkdir -p "$(dir $(BUILD_DIR)/src/$(IMPORT_PATH))"
+	ln -s `pwd` "$(BUILD_DIR)/src/$(IMPORT_PATH)"
+	mkdir -p bin
+	ln -s `pwd`/bin $(BUILD_DIR)/bin
+	touch $@

csidh/cmd/torturer/main.go

@@ -0,0 +1,86 @@
+package main
+
+import (
+    "bytes"
+    "fmt"
+    "time"
+
+    "github.com/henrydcase/nobs/dh/csidh"
+    "github.com/henrydcase/nobs/drbg"
+    "github.com/henrydcase/sidh_torture/csidh/ref/go-wrapper"
+)
+
+var rng *drbg.CtrDrbg
+
+func TestSS() {
+    var r wrapper.Ref
+    var ssRef [csidh.SharedSecretSize]byte
+    var ss [csidh.SharedSecretSize]byte
+    var prB csidh.PrivateKey
+    var pkB csidh.PublicKey
+
+    // Go types
+    prA := r.KeygenPrv()
+    pkA := r.KeygenPub(&prA)
+
+    csidh.GeneratePrivateKey(&prB, rng)
+    csidh.GeneratePublicKey(&pkB, &prB, rng)
+
+    s1 := time.Now()
+    // Generate with Go
+    csidh.DeriveSecret(&ss, &pkA, &prB, rng)
+    s2 := time.Since(s1)
+
+    // Generate with C
+    s1 = time.Now()
+    r.Derive(ssRef[:], &pkB, &prA)
+    s3 := time.Since(s1)
+    fmt.Printf("|TestSharedSecret  | %.8s | %.8s |\n", s2, s3)
+
+    if !bytes.Equal(ss[:], ssRef[:]) {
+        panic("TestSS")
+    }
+}
+
+func TestKeyGeneration() {
+    var r wrapper.Ref
+    var ssRef [csidh.SharedSecretSize]byte
+    var ss [csidh.SharedSecretSize]byte
+
+    // Go types
+    var prA csidh.PrivateKey
+    var pkA csidh.PublicKey
+    var prB csidh.PrivateKey
+    var pkB csidh.PublicKey
+
+    csidh.GeneratePrivateKey(&prA, rng)
+    csidh.GeneratePublicKey(&pkA, &prA, rng)
+    csidh.GeneratePrivateKey(&prB, rng)
+    csidh.GeneratePublicKey(&pkB, &prB, rng)
+
+    s1 := time.Now()
+    // Generate with Go
+    r.Derive(ss[:], &pkA, &prB)
+    s2 := time.Since(s1)
+
+    // Generate with C
+    s1 = time.Now()
+    r.Derive(ssRef[:], &pkB, &prA)
+    s3 := time.Since(s1)
+    fmt.Printf("|TestKeyGeneration | %.8s | %.8s |\n", s2, s3)
+
+    if !bytes.Equal(ss[:], ssRef[:]) {
+        panic("TestKeyGeneration")
+    }
+}
+
+func main() {
+    fmt.Printf("| TestName         |Go        | C        |\n")
+    fmt.Printf("|------------------|----------|----------|\n")
+
+    rng = drbg.NewCtrDrbg()
+    for {
+        TestSS()
+        TestKeyGeneration()
+    }
+}

csidh/cmd/vector-creator/main.go

@@ -0,0 +1,188 @@
+package main
+
+import (
+	"bytes"
+	"encoding/hex"
+	"encoding/json"
+	"io/ioutil"
+
+	"github.com/henrydcase/nobs/dh/csidh"
+	"github.com/henrydcase/sidh_torture/csidh/ref/go-wrapper"
+)
+
+// Possible values for "Status"
+const (
+	Valid               = iota // Indicates that shared secret must be agreed correctly
+	ValidPublicKey2            // Public key 2 must succeed validation
+	InvalidSharedSecret        // Calculated shared secret must be different than test vector
+	InvalidPublicKey1          // Public key 1 generated from private key must be different than test vector
+	InvalidPublicKey2          // Public key 2 must fail validation
+)
+
+var StatusValues = map[int]string{
+	Valid:               "valid",
+	ValidPublicKey2:     "valid_public_key2",
+	InvalidSharedSecret: "invalid_shared_secret",
+	InvalidPublicKey1:   "invalid_public_key1",
+	InvalidPublicKey2:   "invalid_public_key2",
+}
+
+type TestVector struct {
+	Id      int    `json:"Id"`
+	Pk1     string `json:"Pk1"`
+	Pr1     string `json:"Pr1"`
+	Pk2     string `json:"Pk2"`
+	Ss      string `json:"Ss"`
+	Status  string `json:"status"`
+	Comment string `json:"comment"`
+}
+
+type TestVectors struct {
+	Vectors []TestVector `json:"Vectors"`
+}
+
+// R is a reference to C implementation.
+var R wrapper.Ref
+
+// createValid creates 'num' of TestVector's. Each vector contains
+// valid data.
+func createValid(num int) TestVector {
+	var tv TestVector
+	var ss [csidh.SharedSecretSize]byte
+	var pk [csidh.PublicKeySize]byte
+	var pr [csidh.PrivateKeySize]byte
+
+	prA := R.KeygenPrv()
+	pkA := R.KeygenPub(&prA)
+	prB := R.KeygenPrv()
+	pkB := R.KeygenPub(&prB)
+
+	R.Derive(ss[:], &pkB, &prA)
+
+	tv.Id = num
+	tv.Status = StatusValues[Valid]
+
+	tv.Ss = hex.EncodeToString(ss[:])
+
+	prA.Export(pr[:])
+	tv.Pr1 = hex.EncodeToString(pr[:])
+
+	pkA.Export(pk[:])
+	tv.Pk1 = hex.EncodeToString(pk[:])
+
+	for i, _ := range pk {
+		pk[i] = 0
+	}
+
+	pkB.Export(pk[:])
+	tv.Pk2 = hex.EncodeToString(pk[:])
+
+	return tv
+}
+
+func createNegativeSharedSecret(vectors *TestVectors) {
+	n := len(vectors.Vectors)
+	tv := createValid(n)
+	ss, err := hex.DecodeString(tv.Ss)
+	if err != nil {
+		panic("Can't decode shared secret")
+	}
+
+	for i := 0; i < csidh.SharedSecretSize; i++ {
+		var newSs [csidh.SharedSecretSize]byte
+		copy(newSs[:], ss[:])
+		newSs[i] = ss[i] ^ ss[(i+1)%csidh.SharedSecretSize]
+
+		if bytes.Equal(newSs[:], ss) {
+			tv.Status = StatusValues[Valid]
+		} else {
+			tv.Status = StatusValues[InvalidSharedSecret]
+		}
+		tv.Ss = hex.EncodeToString(newSs[:])
+		tv.Id = n + i
+		vectors.Vectors = append(vectors.Vectors, tv)
+	}
+}
+
+// Public key validation fails
+func createNegativePk2(vectors *TestVectors) {
+	n := len(vectors.Vectors)
+	tv := createValid(n)
+	pk, err := hex.DecodeString(tv.Pk2)
+	if err != nil {
+		panic("Can't decode public key 2")
+	}
+
+	for i := 0; i < csidh.PublicKeySize; i++ {
+		var newPk [csidh.PublicKeySize]byte
+
+		// Modify good public key so it's probably no longer valid
+		copy(newPk[:], pk[:])
+		newPk[i] = pk[i] ^ pk[(i+1)%csidh.PublicKeySize]
+
+		// Try to validate and set the status
+		if R.Validate(newPk[:]) {
+			tv.Status = StatusValues[Valid]
+		} else {
+			tv.Status = StatusValues[InvalidPublicKey2]
+		}
+		tv.Pk2 = hex.EncodeToString(newPk[:])
+		tv.Id = n + i
+		vectors.Vectors = append(vectors.Vectors, tv)
+	}
+}
+
+// Private key doesn't correspond to public key
+func createNegativePrivateKey(vectors *TestVectors, num int) {
+	n := len(vectors.Vectors)
+	for i := 0; i < num; i++ {
+		var tv TestVector
+		var pkBytes1 [csidh.PublicKeySize]byte
+		var pkBytes2 [csidh.PublicKeySize]byte
+		var prBytes [csidh.PrivateKeySize]byte
+
+		pr1 := R.KeygenPrv()
+		pk1 := R.KeygenPub(&pr1)
+
+		// Store private key 1
+		pr1.Export(prBytes[:])
+		tv.Pr1 = hex.EncodeToString(prBytes[:])
+
+		// Generate public key which doesn't correspond to pr1
+		pr2 := R.KeygenPrv()
+		pk2 := R.KeygenPub(&pr2)
+
+		pk1.Export(pkBytes1[:])
+		pk2.Export(pkBytes2[:])
+		if bytes.Equal(pkBytes1[:], pkBytes2[:]) {
+			tv.Status = StatusValues[Valid]
+		} else {
+			tv.Status = StatusValues[InvalidPublicKey1]
+		}
+
+		tv.Id = n + i
+		tv.Pk1 = hex.EncodeToString(pkBytes2[:])
+
+		vectors.Vectors = append(vectors.Vectors, tv)
+	}
+}
+
+// TODO: Produce test case for a Pk2 with all zeros and status valid_public_key2.
+// comment should be "Zero key should validate correctly"
+
+func main() {
+	var vectors TestVectors
+
+	for i := 0; i < 100; i++ {
+		vectors.Vectors = append(vectors.Vectors, createValid(i))
+	}
+	createNegativeSharedSecret(&vectors)
+	createNegativePrivateKey(&vectors, 100)
+	createNegativePk2(&vectors)
+
+	marshalled, err := json.MarshalIndent(vectors, "", "  ")
+	if err != nil {
+		panic("Error occured while Marshalling")
+	}
+	ioutil.WriteFile("testvectors.dat", marshalled, 0644)
+}

csidh/ref/csidh/bench.c

@@ -0,0 +1,54 @@
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <time.h>
+#include <assert.h>
+
+#include "u512.h"
+#include "fp.h"
+#include "mont.h"
+#include "csidh.h"
+
+#include <inttypes.h>
+
+static __inline__ uint64_t rdtsc(void)
+{
+    uint32_t hi, lo;
+    __asm__ __volatile__ ("rdtsc" : "=a"(lo), "=d"(hi));
+    return lo | (uint64_t) hi << 32;
+}
+
+unsigned long its = 10000;
+
+int main()
+{
+    clock_t t0, t1, time = 0;
+    uint64_t c0, c1, cycles = 0;
+
+    private_key priv;
+    public_key pub = base;
+
+    for (unsigned long i = 0; i < its; ++i) {
+
+        csidh_private(&priv);
+
+        t0 = clock();
+        c0 = rdtsc();
+
+        /**************************************/
+        assert(validate(&pub));
+        action(&pub, &pub, &priv);
+        /**************************************/
+
+        c1 = rdtsc();
+        t1 = clock();
+        cycles += c1 - c0;
+        time += t1 - t0;
+    }
+
+    printf("iterations: %lu\n", its);
+    printf("clock cycles: %" PRIu64 "\n", (uint64_t) cycles / its);
+    printf("wall-clock time: %.3lf ms\n", 1000. * time / CLOCKS_PER_SEC / its);
+}
+

csidh/ref/csidh/csidh.c

@@ -0,0 +1,285 @@
+
+#include <string.h>
+#include <assert.h>
+
+#include "csidh.h"
+#include "rng.h"
+
+#include <stdio.h>
+
+/* specific to p, should perhaps be somewhere else */
+const unsigned primes[num_primes] = {
+      3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,
+     61,  67,  71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137,
+    139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
+    229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
+    317, 331, 337, 347, 349, 353, 359, 367, 373, 587,
+};
+
+const u512 four_sqrt_p = {{
+    0x85e2579c786882cf, 0x4e3433657e18da95, 0x850ae5507965a0b3, 0xa15bc4e676475964,
+}};
+
+
+const public_key base = {0}; /* A = 0 */
+
+void csidh_private(private_key *priv)
+{
+    memset(&priv->e, 0, sizeof(priv->e));
+    for (size_t i = 0; i < num_primes; ) {
+        int8_t buf[64];
+        randombytes(buf, sizeof(buf));
+        for (size_t j = 0; j < sizeof(buf); ++j) {
+            if (buf[j] <= max_exponent && buf[j] >= -max_exponent) {
+                priv->e[i / 2] |= (buf[j] & 0xf) << i % 2 * 4;
+                if (++i >= num_primes)
+                    break;
+            }
+        }
+    }
+}
+
+/* compute [(p+1)/l] P for all l in our list of primes. */
+/* divide and conquer is much faster than doing it naively,
+ * but uses more memory. */
+static void cofactor_multiples(proj *P, const proj *A, size_t lower, size_t upper)
+{
+    assert(lower < upper);
+
+    if (upper - lower == 1)
+        return;
+
+    size_t mid = lower + (upper - lower + 1) / 2;
+
+    u512 cl = u512_1, cu = u512_1;
+    for (size_t i = lower; i < mid; ++i)
+        u512_mul3_64(&cu, &cu, primes[i]);
+    for (size_t i = mid; i < upper; ++i)
+        u512_mul3_64(&cl, &cl, primes[i]);
+
+    xMUL(&P[mid], A, &P[lower], &cu);
+    xMUL(&P[lower], A, &P[lower], &cl);
+
+    cofactor_multiples(P, A, lower, mid);
+    cofactor_multiples(P, A, mid, upper);
+}
+
+/* never accepts invalid keys. */
+bool validate(public_key const *in)
+{
+    const proj A = {in->A, fp_1};
+
+    do {
+
+        proj P[num_primes];
+        fp_random(&P->x);
+        P->z = fp_1;
+
+        /* maximal 2-power in p+1 */
+        xDBL(P, &A, P);
+        xDBL(P, &A, P);
+
+        cofactor_multiples(P, &A, 0, num_primes);
+
+        u512 order = u512_1;
+
+        for (size_t i = num_primes - 1; i < num_primes; --i) {
+
+            /* we only gain information if [(p+1)/l] P is non-zero */
+            if (memcmp(&P[i].z, &fp_0, sizeof(fp))) {
+
+                u512 tmp;
+                u512_set(&tmp, primes[i]);
+                xMUL(&P[i], &A, &P[i], &tmp);
+
+                if (memcmp(&P[i].z, &fp_0, sizeof(fp)))
+                    /* P does not have order dividing p+1. */
+                    return false;
+
+                u512_mul3_64(&order, &order, primes[i]);
+
+                if (u512_sub3(&tmp, &four_sqrt_p, &order)) /* returns borrow */
+                    /* order > 4 sqrt(p), hence definitely supersingular */
+                    return true;
+            }
+        }
+
+    /* P didn't have big enough order to prove supersingularity. */
+    } while (1);
+}
+
+/* compute x^3 + Ax^2 + x */
+static void montgomery_rhs(fp *rhs, fp const *A, fp const *x)
+{
+    fp tmp;
+    *rhs = *x;
+    fp_sq1(rhs);
+    fp_mul3(&tmp, A, x);
+    fp_add2(rhs, &tmp);
+    fp_add2(rhs, &fp_1);
+    fp_mul2(rhs, x);
+}
+
+/* totally not constant-time. */
+void action(public_key *out, public_key const *in, private_key const *priv)
+{
+    u512 k[2];
+    u512_set(&k[0], 4); /* maximal 2-power in p+1 */
+    u512_set(&k[1], 4); /* maximal 2-power in p+1 */
+
+    uint8_t e[2][num_primes];
+
+    for (size_t i = 0; i < num_primes; ++i) {
+
+        int8_t t = (int8_t) (priv->e[i / 2] << i % 2 * 4) >> 4;
+
+        if (t > 0) {
+            e[0][i] = t;
+            e[1][i] = 0;
+            u512_mul3_64(&k[1], &k[1], primes[i]);
+        }
+        else if (t < 0) {
+            e[1][i] = -t;
+            e[0][i] = 0;
+            u512_mul3_64(&k[0], &k[0], primes[i]);
+        }
+        else {
+            e[0][i] = 0;
+            e[1][i] = 0;
+            u512_mul3_64(&k[0], &k[0], primes[i]);
+            u512_mul3_64(&k[1], &k[1], primes[i]);
+        }
+    }
+
+    proj A = {in->A, fp_1};
+
+    bool done[2] = {false, false};
+
+    do {
+
+        assert(!memcmp(&A.z, &fp_1, sizeof(fp)));
+
+        proj P;
+        fp_random(&P.x);
+        P.z = fp_1;
+
+        fp rhs;
+        montgomery_rhs(&rhs, &A.x, &P.x);
+        bool sign = !fp_issquare(&rhs);
+
+        if (done[sign])
+            continue;
+
+        xMUL(&P, &A, &P, &k[sign]);
+
+        done[sign] = true;
+
+        for (size_t i = 0; i < num_primes; ++i) {
+
+            if (e[sign][i]) {
+
+                u512 cof = u512_1;
+                for (size_t j = i + 1; j < num_primes; ++j)
+                    if (e[sign][j])
+                        u512_mul3_64(&cof, &cof, primes[j]);
+
+                proj K;
+                xMUL(&K, &A, &P, &cof);
+
+                if (memcmp(&K.z, &fp_0, sizeof(fp))) {
+
+                    xISOG(&A, &P, &K, primes[i]);
+
+                    if (!--e[sign][i])
+                        u512_mul3_64(&k[sign], &k[sign], primes[i]);
+
+                }
+
+            }
+
+            done[sign] &= !e[sign][i];
+        }
+
+        fp_inv(&A.z);
+        fp_mul2(&A.x, &A.z);
+        A.z = fp_1;
+
+    } while (!(done[0] && done[1]));
+
+    out->A = A.x;
+}
+
+/* includes public-key validation. */
+bool csidh(public_key *out, public_key const *in, private_key const *priv)
+{
+    if (!validate(in)) {
+        fp_random(&out->A);
+        return false;
+    }
+    action(out, in, priv);
+    return true;
+}
+
+/// Public API
+
+void export_public(uint8_t *out, const public_key *pub) {
+    for(size_t i=0; i<64; i++) {
+        size_t j = i / 8;
+        size_t k = i % 8;
+        out[i] = ((uint64_t)pub->A.x.c[j] >> (8*k)) & 0xFF;
+    }
+}
+
+bool import_public(public_key *pub, const uint8_t *out) {
+    for (size_t i=0; i<64; i++) {
+        size_t j = i / 8;
+        size_t k = i%8;
+        uint64_t tmp = out[i];
+        pub->A.x.c[j] |= tmp << (8*k);
+    }
+    return validate(pub);
+}
+
+void export_private(uint8_t *out, const private_key *prv) {
+    for(size_t i=0; i<37; i++) {
+        out[i] = (uint8_t) prv->e[i];
+    }
+}
+
+void import_private(private_key *prv, const uint8_t *out) {
+    for (size_t i=0; i<37; i++) {
+        prv->e[i] = (int8_t) out[i];
+    }
+}
+
+void keygen_prv(uint8_t private[37]) {
+    private_key priv;
+    csidh_private(&priv);
+    export_private(&private[0], &priv);
+}
+
+void keygen_pub(uint8_t public[64], const uint8_t private[37]) {
+    private_key prv1 = {0};
+    public_key pub1 = {0};
+    import_private(&prv1, &private[0]);
+    csidh(&pub1, &base, &prv1);
+    export_public(&public[0], &pub1);
+}
+
+void derive(uint8_t result[64], const uint8_t public[64], const uint8_t private[37]) {
+    public_key pub = {0};
+    public_key ss = {0};
+    private_key prv;
+
+    import_private(&prv, &private[0]);
+    import_public(&pub, &public[0]);
+
+    csidh(&ss, &pub, &prv);
+    export_public(&result[0], &ss);
+}
+
+// returns 1 if valid 0 otherwise
+int is_valid(const uint8_t *a) {
+    public_key pub = {0};
+    return import_public(&pub, a);
+}

csidh/ref/csidh/csidh.h

@@ -0,0 +1,35 @@
+#ifndef CSIDH_H
+#define CSIDH_H
+
+#include "u512.h"
+#include "fp.h"
+#include "mont.h"
+
+/* specific to p, should perhaps be somewhere else */
+#define num_primes 74
+#define max_exponent 5 /* (2*5+1)^74 is roughly 2^256 */
+
+
+typedef struct private_key {
+    int8_t e[(num_primes + 1) / 2]; /* packed int4_t */
+} private_key;
+
+typedef struct public_key {
+    fp A; /* Montgomery coefficient: represents y^2 = x^3 + Ax^2 + x */
+} public_key;
+
+extern const public_key base;
+
+void csidh_private(private_key *priv);
+bool csidh(public_key *out, public_key const *in, private_key const *priv);
+
+void keygen_prv(uint8_t private[37]);
+void keygen_pub(uint8_t public[64], const uint8_t private[37]);
+void derive(uint8_t result[64], const uint8_t public[64], const uint8_t private_key[37]);
+
+void export_public(uint8_t *out, const public_key *pub);
+bool import_public(public_key *pub, const uint8_t *out);
+void export_private(uint8_t *out, const private_key *prv);
+void import_private(private_key *prv, const uint8_t *out);
+int is_valid(const uint8_t *a);
+#endif

csidh/ref/csidh/fp.h

@@ -0,0 +1,37 @@
+#ifndef FP_H
+#define FP_H
+
+#include "u512.h"
+
+/* fp is in the Montgomery domain, so interpreting that
+   as an integer should never make sense.
+   enable compiler warnings when mixing up u512 and fp. */
+typedef struct fp {
+    u512 x;
+} fp;
+
+extern const fp fp_0;
+extern const fp fp_1;
+
+void fp_set(fp *x, uint64_t y);
+void fp_cswap(fp *x, fp *y, bool c);
+
+void fp_enc(fp *x, u512 const *y); /* encode to Montgomery representation */
+void fp_dec(u512 *x, fp const *y); /* decode from Montgomery representation */
+
+void fp_add2(fp *x, fp const *y);
+void fp_sub2(fp *x, fp const *y);
+void fp_mul2(fp *x, fp const *y);
+
+void fp_add3(fp *x, fp const *y, fp const *z);
+void fp_sub3(fp *x, fp const *y, fp const *z);
+void fp_mul3(fp *x, fp const *y, fp const *z);
+
+void fp_sq1(fp *x);
+void fp_sq2(fp *x, fp const *y);
+void fp_inv(fp *x);
+bool fp_issquare(fp const *x);
+
+void fp_random(fp *x);
+
+#endif

csidh/ref/csidh/fp.s

@@ -0,0 +1,452 @@
+
+.intel_syntax noprefix
+
+.section .rodata
+
+.set pbits, 511
+p:
+    .quad 0x1b81b90533c6c87b, 0xc2721bf457aca835, 0x516730cc1f0b4f25, 0xa7aac6c567f35507
+    .quad 0x5afbfcc69322c9cd, 0xb42d083aedc88c42, 0xfc8ab0d15e3e4c4a, 0x65b48e8f740f89bf
+
+
+.global fp_0
+fp_0: .quad 0, 0, 0, 0, 0, 0, 0, 0
+
+.global fp_1
+fp_1: /* 2^512 mod p */
+    .quad 0xc8fc8df598726f0a, 0x7b1bc81750a6af95, 0x5d319e67c1e961b4, 0xb0aa7275301955f1
+    .quad 0x4a080672d9ba6c64, 0x97a5ef8a246ee77b, 0x06ea9e5d4383676a, 0x3496e2e117e0ec80
+
+
+/* (2^512)^2 mod p */
+.r_squared_mod_p:
+    .quad 0x36905b572ffc1724, 0x67086f4525f1f27d, 0x4faf3fbfd22370ca, 0x192ea214bcc584b1
+    .quad 0x5dae03ee2f5de3d0, 0x1e9248731776b371, 0xad5f166e20e4f52d, 0x4ed759aea6f3917e
+
+/* -p^-1 mod 2^64 */
+.inv_min_p_mod_r:
+    .quad 0x66c1301f632e294d
+
+
+.section .text
+
+.global fp_copy
+fp_copy:
+    cld
+    mov rcx, 8
+    rep movsq
+    ret
+
+.global fp_set
+fp_set:
+    push rdi
+    call u512_set
+    pop rdi
+    mov rsi, rdi
+    jmp fp_enc
+
+.global fp_cswap
+fp_cswap:
+    movzx rax, dl
+    neg rax
+    .set k, 0
+    .rept 8
+        mov rcx, [rdi + 8*k]
+        mov rdx, [rsi + 8*k]
+
+        mov r8, rcx
+        xor r8, rdx
+        and r8, rax
+
+        xor rcx, r8
+        xor rdx, r8
+
+        mov [rdi + 8*k], rcx
+        mov [rsi + 8*k], rdx
+
+        .set k, k+1
+    .endr
+    ret
+
+.reduce_once:
+    push rbp
+    mov rbp, rdi
+
+    mov rdi, [rbp +  0]
+    sub rdi, [rip + p +  0]
+    mov rsi, [rbp +  8]
+    sbb rsi, [rip + p +  8]
+    mov rdx, [rbp + 16]
+    sbb rdx, [rip + p + 16]
+    mov rcx, [rbp + 24]
+    sbb rcx, [rip + p + 24]
+    mov r8,  [rbp + 32]
+    sbb r8,  [rip + p + 32]
+    mov r9,  [rbp + 40]
+    sbb r9,  [rip + p + 40]
+    mov r10, [rbp + 48]
+    sbb r10, [rip + p + 48]
+    mov r11, [rbp + 56]
+    sbb r11, [rip + p + 56]
+
+    setnc al
+    movzx rax, al
+    neg rax
+
+.macro cswap2, r, m
+    xor \r, \m
+    and \r, rax
+    xor \m, \r
+.endm
+
+    cswap2 rdi, [rbp +  0]
+    cswap2 rsi, [rbp +  8]
+    cswap2 rdx, [rbp + 16]
+    cswap2 rcx, [rbp + 24]
+    cswap2 r8,  [rbp + 32]
+    cswap2 r9,  [rbp + 40]
+    cswap2 r10, [rbp + 48]
+    cswap2 r11, [rbp + 56]
+
+    pop rbp
+    ret
+
+.global fp_add3
+fp_add3:
+    push rdi
+    call u512_add3
+    pop rdi
+    jmp .reduce_once
+
+.global fp_add2
+fp_add2:
+    mov rdx, rdi
+    jmp fp_add3
+
+.global fp_sub3
+fp_sub3:
+    push rdi
+    call u512_sub3
+    pop rdi
+    xor rsi, rsi
+    xor rdx, rdx
+    xor rcx, rcx
+    xor r8, r8
+    xor r9, r9
+    xor r10, r10
+    xor r11, r11
+    test rax, rax
+    cmovnz rax, [rip + p +  0]
+    cmovnz rsi, [rip + p +  8]
+    cmovnz rdx, [rip + p + 16]
+    cmovnz rcx, [rip + p + 24]
+    cmovnz r8,  [rip + p + 32]
+    cmovnz r9,  [rip + p + 40]
+    cmovnz r10, [rip + p + 48]
+    cmovnz r11, [rip + p + 56]
+    add [rdi +  0], rax
+    adc [rdi +  8], rsi
+    adc [rdi + 16], rdx
+    adc [rdi + 24], rcx
+    adc [rdi + 32],  r8
+    adc [rdi + 40],  r9
+    adc [rdi + 48], r10
+    adc [rdi + 56], r11
+    ret
+
+.global fp_sub2
+fp_sub2:
+    mov rdx, rdi
+    xchg rsi, rdx
+    jmp fp_sub3
+
+
+/* Montgomery arithmetic */
+
+.global fp_enc
+fp_enc:
+    lea rdx, [rip + .r_squared_mod_p]
+    jmp fp_mul3
+
+.global fp_dec
+fp_dec:
+    lea rdx, [rip + u512_1]
+    jmp fp_mul3
+
+.global fp_mul3
+fp_mul3:
+    push rbp
+    push rbx
+    push r12
+    push r13
+    push r14
+    push r15
+
+    push rdi
+
+    mov rdi, rsi
+    mov rsi, rdx
+
+    xor r8,  r8
+    xor r9,  r9
+    xor r10, r10
+    xor r11, r11
+    xor r12, r12
+    xor r13, r13
+    xor r14, r14
+    xor r15, r15
+    xor rbp, rbp
+
+    /* flags are already cleared */
+
+.macro MULSTEP, k, r0, r1, r2, r3, r4, r5, r6, r7, r8
+
+    mov rdx, [rsi +  0]
+    mulx rcx, rdx, [rdi + 8*\k]
+    add rdx, \r0
+    mulx rcx, rdx, [rip + .inv_min_p_mod_r]
+
+    xor rax, rax /* clear flags */
+
+    mulx rbx, rax, [rip + p +  0]
+    adox \r0, rax
+
+    mulx rcx, rax, [rip + p +  8]
+    adcx \r1, rbx
+    adox \r1, rax
+
+    mulx rbx, rax, [rip + p + 16]
+    adcx \r2, rcx
+    adox \r2, rax
+
+    mulx rcx, rax, [rip + p + 24]
+    adcx \r3, rbx
+    adox \r3, rax
+
+    mulx rbx, rax, [rip + p + 32]
+    adcx \r4, rcx
+    adox \r4, rax
+
+    mulx rcx, rax, [rip + p + 40]
+    adcx \r5, rbx
+    adox \r5, rax
+
+    mulx rbx, rax, [rip + p + 48]
+    adcx \r6, rcx
+    adox \r6, rax
+
+    mulx rcx, rax, [rip + p + 56]
+    adcx \r7, rbx
+    adox \r7, rax
+
+    mov rax, 0
+    adcx \r8, rcx
+    adox \r8, rax
+
+
+    mov rdx, [rdi + 8*\k]
+
+    xor rax, rax /* clear flags */
+
+    mulx rbx, rax, [rsi +  0]
+    adox \r0, rax
+
+    mulx rcx, rax, [rsi +  8]
+    adcx \r1, rbx
+    adox \r1, rax
+
+    mulx rbx, rax, [rsi + 16]
+    adcx \r2, rcx
+    adox \r2, rax
+
+    mulx rcx, rax, [rsi + 24]
+    adcx \r3, rbx
+    adox \r3, rax
+
+    mulx rbx, rax, [rsi + 32]
+    adcx \r4, rcx
+    adox \r4, rax
+
+    mulx rcx, rax, [rsi + 40]
+    adcx \r5, rbx
+    adox \r5, rax
+
+    mulx rbx, rax, [rsi + 48]
+    adcx \r6, rcx
+    adox \r6, rax
+
+    mulx rcx, rax, [rsi + 56]
+    adcx \r7, rbx
+    adox \r7, rax
+
+    mov rax, 0
+    adcx \r8, rcx
+    adox \r8, rax
+
+.endm
+
+    MULSTEP 0, r8,  r9,  r10, r11, r12, r13, r14, r15, rbp
+    MULSTEP 1, r9,  r10, r11, r12, r13, r14, r15, rbp, r8
+    MULSTEP 2, r10, r11, r12, r13, r14, r15, rbp, r8,  r9
+    MULSTEP 3, r11, r12, r13, r14, r15, rbp, r8,  r9,  r10
+    MULSTEP 4, r12, r13, r14, r15, rbp, r8,  r9,  r10, r11
+    MULSTEP 5, r13, r14, r15, rbp, r8,  r9,  r10, r11, r12
+    MULSTEP 6, r14, r15, rbp, r8,  r9,  r10, r11, r12, r13
+    MULSTEP 7, r15, rbp, r8,  r9,  r10, r11, r12, r13, r14
+
+    pop rdi
+
+    mov [rdi +  0], rbp
+    mov [rdi +  8], r8
+    mov [rdi + 16], r9
+    mov [rdi + 24], r10
+    mov [rdi + 32], r11
+    mov [rdi + 40], r12
+    mov [rdi + 48], r13
+    mov [rdi + 56], r14
+
+    pop r15
+    pop r14
+    pop r13
+    pop r12
+    pop rbx
+    pop rbp
+    jmp .reduce_once
+
+.global fp_mul2
+fp_mul2:
+    mov rdx, rdi
+    jmp fp_mul3
+
+.global fp_sq2
+fp_sq2:
+    /* TODO implement optimized Montgomery squaring */
+    mov rdx, rsi
+    jmp fp_mul3
+
+.global fp_sq1
+fp_sq1:
+    mov rsi, rdi
+    jmp fp_sq2
+
+/* (obviously) not constant time in the exponent! */
+.fp_pow:
+    push rbx
+    mov rbx, rsi
+    push r12
+    push r13
+    push rdi
+    sub rsp, 64
+
+    mov rsi, rdi
+    mov rdi, rsp
+    call fp_copy
+
+    mov rdi, [rsp + 64]
+    lea rsi, [rip + fp_1]
+    call fp_copy
+
+.macro POWSTEP, k
+        mov r13, [rbx + 8*\k]
+        xor r12, r12
+
+        0:
+        test r13, 1
+        jz 1f
+
+        mov rdi, [rsp + 64]
+        mov rsi, rsp
+        call fp_mul2
+
+        1:
+        mov rdi, rsp
+        call fp_sq1
+
+        shr r13
+
+        inc r12
+        test r12, 64
+        jz 0b
+.endm
+
+    POWSTEP 0
+    POWSTEP 1
+    POWSTEP 2
+    POWSTEP 3
+    POWSTEP 4
+    POWSTEP 5
+    POWSTEP 6
+    POWSTEP 7
+
+    add rsp, 64+8
+    pop r13
+    pop r12
+    pop rbx
+    ret
+
+.section .rodata
+.p_minus_2:
+    .quad 0x1b81b90533c6c879, 0xc2721bf457aca835, 0x516730cc1f0b4f25, 0xa7aac6c567f35507
+    .quad 0x5afbfcc69322c9cd, 0xb42d083aedc88c42, 0xfc8ab0d15e3e4c4a, 0x65b48e8f740f89bf
+
+.section .text
+
+/* TODO use a better addition chain? */
+.global fp_inv
+fp_inv:
+    lea rsi, [rip + .p_minus_2]
+    jmp .fp_pow
+
+.section .rodata
+.p_minus_1_halves:
+    .quad 0x8dc0dc8299e3643d, 0xe1390dfa2bd6541a, 0xa8b398660f85a792, 0xd3d56362b3f9aa83
+    .quad 0x2d7dfe63499164e6, 0x5a16841d76e44621, 0xfe455868af1f2625, 0x32da4747ba07c4df
+
+.section .text
+
+/* TODO use a better addition chain? */
+.global fp_issquare
+fp_issquare:
+    push rdi
+    lea rsi, [rip + .p_minus_1_halves]
+    call .fp_pow
+    pop rdi
+
+    xor rax, rax
+    .set k, 0
+    .rept 8
+        mov rsi, [rdi + 8*k]
+        xor rsi, [rip + fp_1 + 8*k]
+        or rax, rsi
+        .set k, k+1
+    .endr
+    test rax, rax
+    setz al
+    movzx rax, al
+    ret
+
+
+/* not constant time (but this shouldn't leak anything of importance) */
+.global fp_random
+fp_random:
+
+    push rdi
+    mov rsi, 64
+    call randombytes
+    pop rdi
+    mov rax, 1
+    shl rax, (pbits % 64)
+    dec rax
+    and [rdi + 56], rax
+
+    .set k, 7
+    .rept 8
+        mov rax, [rip + p + 8*k]
+        cmp [rdi + 8*k], rax
+        jge fp_random
+        jl 0f
+        .set k, k-1
+    .endr
+    0:
+    ret
+

csidh/ref/csidh/main.c

@@ -0,0 +1,99 @@
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <time.h>
+#include <assert.h>
+
+#include "u512.h"
+#include "fp.h"
+#include "mont.h"
+#include "csidh.h"
+
+void u512_print(u512 const *x)
+{
+    for (size_t i = 63; i < 64; --i)
+        printf("%02hhx", i[(unsigned char *) x->c]);
+}
+
+void fp_print(fp const *x)
+{
+    u512 y;
+    fp_dec(&y, x);
+    u512_print(&y);
+}
+
+int main()
+{
+    clock_t t0, t1;
+
+    private_key priv_alice, priv_bob;
+    public_key pub_alice, pub_bob;
+    public_key shared_alice, shared_bob;
+
+    printf("\n");
+
+
+    t0 = clock();
+    csidh_private(&priv_alice);
+    t1 = clock();
+
+    printf("Alice's private key   (%7.3lf ms):\n  ", 1000. * (t1 - t0) / CLOCKS_PER_SEC);
+    for (size_t i = 0; i < sizeof(priv_alice); ++i)
+        printf("%02hhx", i[(uint8_t *) &priv_alice]);
+    printf("\n\n");
+
+    t0 = clock();
+    csidh_private(&priv_bob);
+    t1 = clock();
+
+    printf("Bob's private key     (%7.3lf ms):\n  ", 1000. * (t1 - t0) / CLOCKS_PER_SEC);
+    for (size_t i = 0; i < sizeof(priv_bob); ++i)
+        printf("%02hhx", i[(uint8_t *) &priv_bob]);
+    printf("\n\n");
+
+
+    t0 = clock();
+    assert(csidh(&pub_alice, &base, &priv_alice));
+    t1 = clock();
+
+    printf("Alice's public key    (%7.3lf ms):\n  ", 1000. * (t1 - t0) / CLOCKS_PER_SEC);
+    fp_print(&pub_alice.A);
+    printf("\n\n");
+
+    t0 = clock();
+    assert(csidh(&pub_bob, &base, &priv_bob));
+    t1 = clock();
+
+    printf("Bob's public key      (%7.3lf ms):\n  ", 1000. * (t1 - t0) / CLOCKS_PER_SEC);
+    fp_print(&pub_bob.A);
+    printf("\n\n");
+
+
+    t0 = clock();
+    assert(csidh(&shared_alice, &pub_bob, &priv_alice));
+    t1 = clock();
+
+    printf("Alice's shared secret (%7.3lf ms):\n  ", 1000. * (t1 - t0) / CLOCKS_PER_SEC);
+    fp_print(&shared_alice.A);
+    printf("\n\n");
+
+    t0 = clock();
+    assert(csidh(&shared_bob, &pub_alice, &priv_bob));
+    t1 = clock();
+
+    printf("Bob's shared secret   (%7.3lf ms):\n  ", 1000. * (t1 - t0) / CLOCKS_PER_SEC);
+    fp_print(&shared_bob.A);
+    printf("\n\n");
+
+    printf("    ");
+    if (memcmp(&shared_alice, &shared_bob, sizeof(public_key)))
+        printf("\x1b[31mNOT EQUAL!\x1b[0m\n");
+    else
+        printf("\x1b[32mequal.\x1b[0m\n");
+    printf("\n");
+
+    printf("\n");
+}
+

csidh/ref/csidh/mont.c

@@ -0,0 +1,188 @@
+
+#include <assert.h>
+
+#include "mont.h"
+
+void xDBLADD(proj *R, proj *S, proj const *P, proj const *Q, proj const *PQ, proj const *A24)
+{
+    fp tmp0, tmp1, tmp2;
+
+    fp_add3(&tmp0, &P->x, &P->z);
+    fp_sub3(&tmp1, &P->x, &P->z);
+    fp_sq2(&R->x, &tmp0);
+    fp_sub3(&tmp2, &Q->x, &Q->z);
+    fp_add3(&S->x, &Q->x, &Q->z);
+    fp_mul2(&tmp0, &tmp2);
+    fp_sq2(&R->z, &tmp1);
+    fp_mul2(&tmp1, &S->x);
+    fp_sub3(&tmp2, &R->x, &R->z);
+    fp_mul2(&R->z, &A24->z);
+    fp_mul2(&R->x, &R->z);
+    fp_mul3(&S->x, &A24->x, &tmp2);
+    fp_sub3(&S->z, &tmp0, &tmp1);
+    fp_add2(&R->z, &S->x);
+    fp_add3(&S->x, &tmp0, &tmp1);
+    fp_mul2(&R->z, &tmp2);
+    fp_sq1(&S->z);
+    fp_sq1(&S->x);
+    fp_mul2(&S->z, &PQ->x);
+    fp_mul2(&S->x, &PQ->z);
+}
+
+void xDBL(proj *Q, proj const *A, proj const *P)
+{
+    fp a, b, c;
+    fp_add3(&a, &P->x, &P->z);
+    fp_sq1(&a);
+    fp_sub3(&b, &P->x, &P->z);
+    fp_sq1(&b);
+    fp_sub3(&c, &a, &b);
+    fp_add2(&b, &b); fp_add2(&b, &b); /* multiplication by 4 */
+    fp_mul2(&b, &A->z);
+    fp_mul3(&Q->x, &a, &b);
+    fp_add3(&a, &A->z, &A->z); /* multiplication by 2 */
+    fp_add2(&a, &A->x);
+    fp_mul2(&a, &c);
+    fp_add2(&a, &b);
+    fp_mul3(&Q->z, &a, &c);
+}
+
+void xADD(proj *S, proj const *P, proj const *Q, proj const *PQ)
+{
+    fp a, b, c, d;
+    fp_add3(&a, &P->x, &P->z);
+    fp_sub3(&b, &P->x, &P->z);
+    fp_add3(&c, &Q->x, &Q->z);
+    fp_sub3(&d, &Q->x, &Q->z);
+    fp_mul2(&a, &d);
+    fp_mul2(&b, &c);
+    fp_add3(&c, &a, &b);
+    fp_sub3(&d, &a, &b);
+    fp_sq1(&c);
+    fp_sq1(&d);
+    fp_mul3(&S->x, &PQ->z, &c);
+    fp_mul3(&S->z, &PQ->x, &d);
+}
+
+/* Montgomery ladder. */
+/* P must not be the unique point of order 2. */
+/* not constant-time! */
+void xMUL(proj *Q, proj const *A, proj const *P, u512 const *k)
+{
+    proj R = *P;
+    proj A24;
+    const proj Pcopy = *P; /* in case Q = P */
+
+    Q->x = fp_1;
+    Q->z = fp_0;
+
+    fp_add3(&A24.x, &A->z, &A->z);
+    fp_add3(&A24.z, &A24.x, &A24.x);
+    fp_add2(&A24.x, &A->x);
+
+    unsigned long i = 512;
+    while (--i && !u512_bit(k, i));
+
+    do {
+
+        bool bit = u512_bit(k, i);
+
+        if (bit) { proj T = *Q; *Q = R; R = T; } /* not constant-time */
+        //fp_cswap(&Q->x, &R.x, bit);
+        //fp_cswap(&Q->z, &R.z, bit);
+
+        xDBLADD(Q, &R, Q, &R, &Pcopy, &A24);
+
+        if (bit) { proj T = *Q; *Q = R; R = T; } /* not constant-time */
+        //fp_cswap(&Q->x, &R.x, bit);
+        //fp_cswap(&Q->z, &R.z, bit);
+
+    } while (i--);
+}
+
+
+/* computes the isogeny with kernel point K of order k */
+/* returns the new curve coefficient A and the image of P */
+/* (obviously) not constant time in k */
+void xISOG(proj *A, proj *P, proj const *K, uint64_t k)
+{
+    assert (k >= 3);
+    assert (k % 2 == 1);
+
+    fp tmp0, tmp1;
+    fp T[4] = {K->z, K->x, K->x, K->z};
+    proj Q;
+
+    fp_mul3(&Q.x,  &P->x, &K->x);
+    fp_mul3(&tmp0, &P->z, &K->z);
+    fp_sub2(&Q.x,  &tmp0);
+
+    fp_mul3(&Q.z,  &P->x, &K->z);
+    fp_mul3(&tmp0, &P->z, &K->x);
+    fp_sub2(&Q.z,  &tmp0);
+
+    proj M[3] = {*K};
+    xDBL(&M[1], A, K);
+
+    for (uint64_t i = 1; i < k / 2; ++i) {
+
+        if (i >= 2)
+            xADD(&M[i % 3], &M[(i - 1) % 3], K, &M[(i - 2) % 3]);
+
+        fp_mul3(&tmp0, &M[i % 3].x, &T[0]);
+        fp_mul3(&tmp1, &M[i % 3].z, &T[1]);
+        fp_add3(&T[0], &tmp0, &tmp1);
+
+        fp_mul2(&T[1], &M[i % 3].x);
+
+        fp_mul3(&tmp0, &M[i % 3].z, &T[2]);
+        fp_mul3(&tmp1, &M[i % 3].x, &T[3]);
+        fp_add3(&T[2], &tmp0, &tmp1);
+
+        fp_mul2(&T[3], &M[i % 3].z);
+
+
+        fp_mul3(&tmp0, &P->x, &M[i % 3].x);
+        fp_mul3(&tmp1, &P->z, &M[i % 3].z);
+        fp_sub2(&tmp0, &tmp1);
+        fp_mul2(&Q.x,  &tmp0);
+
+        fp_mul3(&tmp0, &P->x, &M[i % 3].z);
+        fp_mul3(&tmp1, &P->z, &M[i % 3].x);
+        fp_sub2(&tmp0, &tmp1);
+        fp_mul2(&Q.z,  &tmp0);
+    }
+
+    fp_mul2(&T[0], &T[1]);
+    fp_add2(&T[0], &T[0]); /* multiplication by 2 */
+
+    fp_sq1(&T[1]);
+
+    fp_mul2(&T[2], &T[3]);
+    fp_add2(&T[2], &T[2]); /* multiplication by 2 */
+
+    fp_sq1(&T[3]);
+
+    /* Ax := T[1] * T[3] * Ax - 3 * Az * (T[1] * T[2] - T[0] * T[3]) */
+    fp_mul3(&tmp0, &T[1], &T[2]);
+    fp_mul3(&tmp1, &T[0], &T[3]);
+    fp_sub2(&tmp0, &tmp1);
+    fp_mul2(&tmp0, &A->z);
+    fp_add3(&tmp1, &tmp0, &tmp0); fp_add2(&tmp0, &tmp1); /* multiplication by 3 */
+
+    fp_mul3(&tmp1, &T[1], &T[3]);
+    fp_mul2(&tmp1, &A->x);
+
+    fp_sub3(&A->x, &tmp1, &tmp0);
+
+    /* Az := Az * T[3]^2 */
+    fp_sq1(&T[3]);
+    fp_mul2(&A->z, &T[3]);
+
+    /* X := X * Xim^2, Z := Z * Zim^2 */
+    fp_sq1(&Q.x);
+    fp_sq1(&Q.z);
+    fp_mul2(&P->x, &Q.x);
+    fp_mul2(&P->z, &Q.z);
+}
+

csidh/ref/csidh/mont.h

@@ -0,0 +1,19 @@
+#ifndef MONT_H
+#define MONT_H
+
+#include "u512.h"
+#include "fp.h"
+
+/* P^1 over fp. */
+typedef struct proj {
+    fp x;
+    fp z;
+} proj;
+
+void xDBL(proj *Q, proj const *A, proj const *P);
+void xADD(proj *S, proj const *P, proj const *Q, proj const *PQ);
+void xDBLADD(proj *R, proj *S, proj const *P, proj const *Q, proj const *PQ, proj const *A);
+void xMUL(proj *Q, proj const *A, proj const *P, u512 const *k);
+void xISOG(proj *A, proj *P, proj const *K, uint64_t k);
+
+#endif

csidh/ref/csidh/rng.c

@@ -0,0 +1,18 @@
+
+#include "rng.h"
+
+#include <stdlib.h>
+#include <unistd.h>
+#include <fcntl.h>
+
+void randombytes(void *x, size_t l)
+{
+    static int fd = -1;
+    ssize_t n;
+    if (fd < 0 && 0 > (fd = open("/dev/urandom", O_RDONLY)))
+        exit(1);
+    for (size_t i = 0; i < l; i += n)
+        if (0 >= (n = read(fd, (char *) x + i, l - i)))
+            exit(2);
+}
+

csidh/ref/csidh/rng.h

@@ -0,0 +1,8 @@
+#ifndef RNG_H
+#define RNG_H
+
+#include <stdlib.h>
+
+void randombytes(void *x, size_t l);
+
+#endif

csidh/ref/csidh/u512.h

@@ -0,0 +1,22 @@
+#ifndef UINT_H
+#define UINT_H
+
+#include <stdbool.h>
+#include <stdint.h>
+
+typedef struct u512 {
+    uint64_t c[8];
+} u512;
+
+extern const u512 u512_1;
+
+void u512_set(u512 *x, uint64_t y);
+
+bool u512_bit(u512 const *x, uint64_t k);
+
+bool u512_add3(u512 *x, u512 const *y, u512 const *z); /* returns carry */
+bool u512_sub3(u512 *x, u512 const *y, u512 const *z); /* returns borrow */
+
+void u512_mul3_64(u512 *x, u512 const *y, uint64_t z);
+
+#endif

csidh/ref/csidh/u512.s

@@ -0,0 +1,102 @@
+
+.intel_syntax noprefix
+
+.section .rodata
+
+.global u512_1
+u512_1: .quad 1, 0, 0, 0, 0, 0, 0, 0
+
+
+.section .text
+
+.global u512_set
+u512_set:
+    cld
+    mov rax, rsi
+    stosq
+    xor rax, rax
+    mov rcx, 7
+    rep stosq
+    ret
+
+
+.global u512_bit
+u512_bit:
+    mov rcx, rsi
+    and rcx, 0x3f
+    shr rsi, 6
+    mov rax, [rdi + 8*rsi]
+    shr rax, cl
+    and rax, 1
+    ret
+
+
+.global u512_add3
+u512_add3:
+    mov rax, [rsi +  0]
+    add rax, [rdx +  0]
+    mov [rdi +  0], rax
+    .set k, 1
+    .rept 7
+        mov rax, [rsi + 8*k]
+        adc rax, [rdx + 8*k]
+        mov [rdi + 8*k], rax
+        .set k, k+1
+    .endr
+    setc al
+    movzx rax, al
+    ret
+
+.global u512_sub3
+u512_sub3:
+    mov rax, [rsi +  0]
+    sub rax, [rdx +  0]
+    mov [rdi +  0], rax
+    .set k, 1
+    .rept 7
+        mov rax, [rsi + 8*k]
+        sbb rax, [rdx + 8*k]
+        mov [rdi + 8*k], rax
+        .set k, k+1
+    .endr
+    setc al
+    movzx rax, al
+    ret
+
+
+.global u512_mul3_64
+u512_mul3_64:
+
+    mulx r10, rax, [rsi +  0]
+    mov [rdi +  0], rax
+
+    mulx r11, rax, [rsi +  8]
+    add  rax, r10
+    mov [rdi +  8], rax
+
+    mulx r10, rax, [rsi + 16]
+    adcx rax, r11
+    mov [rdi + 16], rax
+
+    mulx r11, rax, [rsi + 24]
+    adcx rax, r10
+    mov [rdi + 24], rax
+
+    mulx r10, rax, [rsi + 32]
+    adcx rax, r11
+    mov [rdi + 32],rax
+
+    mulx r11, rax, [rsi + 40]
+    adcx rax, r10
+    mov [rdi + 40],rax
+
+    mulx r10, rax, [rsi + 48]
+    adcx rax, r11
+    mov [rdi + 48],rax
+
+    mulx r11, rax, [rsi + 56]
+    adcx rax, r10
+    mov [rdi + 56],rax
+
+    ret
+

csidh/ref/go-wrapper/wrapper.go

@@ -0,0 +1,80 @@
+package wrapper
+
+/*
+#cgo CFLAGS: -I../../ref/csidh
+#cgo LDFLAGS: -L../../bin -lcsidh
+#include <csidh.h>
+*/
+import "C"
+import (
+    "unsafe"
+
+    csidh "github.com/henrydcase/nobs/dh/csidh"
+)
+
+const (
+    PrvSz = 37
+    PubSz = 64
+    SsSz  = 64
+)
+
+type Ref struct {}
+
+// Converts Golang's byte array to C's unsigned char. 'c' must
+// have same size as byte
+func (Ref) toBytes(c []C.uchar, b []byte) {
+    if len(c) != len(b) {
+        panic("len c!= len b")
+    }
+    for i,v := range(b) {
+        c[i] = C.uchar(v)
+    }
+}
+
+func (Ref) KeygenPrv() (prv csidh.PrivateKey) {
+    var buf [PrvSz]C.uchar
+    C.keygen_prv(&buf[0])
+    prv.Import(C.GoBytes(unsafe.Pointer(&buf[0]), C.int(len(buf))))
+    return
+}
+
+func (c Ref) KeygenPub(prv *csidh.PrivateKey) (pub csidh.PublicKey) {
+    var prv_c_buf [PrvSz]C.uchar
+    var pub_c_buf [PubSz]C.uchar
+    var prv_g_buf [PrvSz]byte
+
+    if !prv.Export(prv_g_buf[:]) {
+        panic("Can't export private key")
+    }
+    c.toBytes(prv_c_buf[:], prv_g_buf[:])
+    C.keygen_pub(&pub_c_buf[0], &prv_c_buf[0])
+    pub.Import(C.GoBytes(unsafe.Pointer(&pub_c_buf[0]), C.int(len(pub_c_buf))))
+    return
+}
+
+func (c Ref) Derive(ss []byte, pub *csidh.PublicKey, prv *csidh.PrivateKey) {
+    var prv_c_buf [PrvSz]C.uchar
+    var pub_c_buf [PubSz]C.uchar
+    var prv_g_buf [PrvSz]byte
+    var pub_g_buf [PubSz]byte
+    var ss_c_buf [SsSz]C.uchar
+
+    if !pub.Export(pub_g_buf[:]) {
+        panic("Can't export public key")
+    }
+
+    if !prv.Export(prv_g_buf[:]) {
+        panic("Can't export private key")
+    }
+
+    c.toBytes(pub_c_buf[:], pub_g_buf[:])
+    c.toBytes(prv_c_buf[:], prv_g_buf[:])
+    C.derive(&ss_c_buf[0], &pub_c_buf[0], &prv_c_buf[0])
+    copy(ss, C.GoBytes(unsafe.Pointer(&ss_c_buf[0]), C.int(len(ss_c_buf))))
+}
+
+func (c Ref) Validate(a []byte) bool {
+    var pub_c_buf [PubSz]C.uchar
+    c.toBytes(pub_c_buf[:], a[:])
+    return C.is_valid(&pub_c_buf[0]) != 0
+}

csidh/ref/test/main.c

@@ -0,0 +1,115 @@
+#include <string.h>
+#include <assert.h>
+#include <stdio.h>
+#include "../csidh/csidh.h"
+#include "../csidh/rng.h"
+
+static void u512_print(u512 const *x)
+{
+    for (size_t i=0; i<8; i++) {
+        printf("0x%016lX,", x->c[i]);
+    }
+    printf("\n");
+}
+
+void print_bytes(uint8_t *out, size_t sz) {
+    for(size_t i=0; i<sz; i++) {
+        printf("%02X", out[i]);
+    }
+    printf("\n");
+}
+
+static void int_print(const int8_t v[37]) {
+    for (size_t i = 0; i<37; i++)
+        printf("0x%X,", (unsigned char)v[i]);
+    printf("\n");
+}
+
+// Print after convertion from Montgomery domain
+static void fp_print(fp const *x)
+{
+    u512 y;
+    fp_dec(&y, x);
+    u512_print(&y);
+}
+
+// Print without converting from Montgomery domain
+static void fp_print_n(fp const *x) {
+    for (size_t i=0; i<8; i++) {
+        printf("0x%016lX, ", x->x.c[i]);
+    }
+    printf("\n");
+}
+
+static bool fp_cmp(fp const *l, fp const *r)
+{
+    u512 lu,ru;
+    fp_dec(&lu, l);
+    fp_dec(&ru, r);
+    return ( sizeof(lu.c) == sizeof(ru.c) ) && (memcmp(&lu.c, &ru.c, sizeof(lu.c)) == 0);
+}
+
+static bool testLoopRef() {
+    for(size_t i=0; i<10; i++) {
+        private_key prA, prB;
+        public_key pkA, pkB;
+        public_key shA, shB;
+
+        // private key
+        csidh_private(&prA);
+        csidh_private(&prB);
+
+        // public key
+        csidh(&pkA, &base, &prA);
+        csidh(&pkB, &base, &prB);
+
+        //csidh
+        csidh(&shA, &pkA, &prB);
+        csidh(&shB, &pkB, &prA);
+        return (memcmp(&shA, &shB, sizeof(shB))==0);
+    }
+}
+
+static bool testImportExport() {
+    uint8_t buf_pub[64]={0};
+    private_key prv1 = {0};
+    private_key prv2 = {0};
+    public_key pub1 = {0};
+    public_key pub2 = {0};
+
+    csidh_private(&prv1);
+    csidh_private(&prv2);
+
+    csidh(&pub1, &base, &prv1);
+
+    export_public(buf_pub, &pub1);
+    import_public(&pub2, buf_pub);
+    return fp_cmp(&pub1.A, &pub2.A);
+}
+
+// special cases:
+//  * diag 0
+//  * only 0
+static void validator() {
+    uint8_t rnd[64] = {0};
+    private_key prA = {0};
+    public_key pubA={0};
+
+    csidh_private(&prA);
+    csidh(&pubA, &base, &prA);
+    export_public(rnd, &pubA);
+
+    for (size_t i=0; i<64; i++) {
+        uint8_t r[64] = {0};
+        memcpy(r, rnd, sizeof(r));
+        r[i]=0x01;
+        printf("%u -> ", import_public(&pubA, r));
+    }
+}
+
+int main() {
+//    validator();
+    return !(
+        testImportExport() &&
+        testLoopRef());
+}

go.mod

@@ -0,0 +1,5 @@
+module github.com/henrydcase/sidh_torture
+
+go 1.13
+
+require github.com/henrydcase/nobs v0.0.0-20191125150329-7efbbf474559

go.sum

@@ -0,0 +1,4 @@
+github.com/henrydcase/nobs v0.0.0-20191125150329-7efbbf474559 h1:tr5PqEaPopB2KYEm/Xji8Eo3j6TR6fuKiN6dOTyQgRY=
+github.com/henrydcase/nobs v0.0.0-20191125150329-7efbbf474559/go.mod h1:+liTPsuK0xSOSyNKhVz4h7Khig8zW4NcvxdVbzS0Jyw=
+golang.org/x/sys v0.0.0-20191120155948-bd437916bb0e h1:N7DeIrjYszNmSW409R3frPPwglRwMkXSBzwVbkOjLLA=
+golang.org/x/sys v0.0.0-20191120155948-bd437916bb0e/go.mod h1:h1NjWce9XRLGQEsW7wpKNCjG9DtNlClVuFLEZdDNbEs=