diff --git a/LICENSE b/LICENSE
index 49c41fa..2f4dfcd 100644
--- a/LICENSE
+++ b/LICENSE
@@ -181,6 +181,29 @@
 SOFTWARE.
 
 
+The code in third_party/sike also carries the MIT license:
+
+Copyright (c) Microsoft Corporation. All rights reserved.
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE
+
+
 Licenses for support code
 -------------------------
 
diff --git a/crypto/CMakeLists.txt b/crypto/CMakeLists.txt
index 5cdfa40..1c505bc 100644
--- a/crypto/CMakeLists.txt
+++ b/crypto/CMakeLists.txt
@@ -115,6 +115,7 @@
 
     chacha/chacha-armv8.${ASM_EXT}
     test/trampoline-armv8.${ASM_EXT}
+    third_party/sike/asm/fp-armv8.${ASM_EXT}
   )
 endif()
 
@@ -136,6 +137,7 @@
     cipher_extra/chacha20_poly1305_x86_64.${ASM_EXT}
     hrss/asm/poly_rq_mul.S
     test/trampoline-x86_64.${ASM_EXT}
+    third_party/sike/asm/fp-x86_64.${ASM_EXT}
   )
 endif()
 
@@ -145,6 +147,8 @@
 perlasm(chacha/chacha-x86_64.${ASM_EXT} chacha/asm/chacha-x86_64.pl)
 perlasm(cipher_extra/aes128gcmsiv-x86_64.${ASM_EXT} cipher_extra/asm/aes128gcmsiv-x86_64.pl)
 perlasm(cipher_extra/chacha20_poly1305_x86_64.${ASM_EXT} cipher_extra/asm/chacha20_poly1305_x86_64.pl)
+perlasm(third_party/sike/asm/fp-x86_64.${ASM_EXT} ../third_party/sike/asm/fp-x86_64.pl)
+perlasm(third_party/sike/asm/fp-armv8.${ASM_EXT} ../third_party/sike/asm/fp-armv8.pl)
 perlasm(test/trampoline-armv4.${ASM_EXT} test/asm/trampoline-armv4.pl)
 perlasm(test/trampoline-armv8.${ASM_EXT} test/asm/trampoline-armv8.pl)
 perlasm(test/trampoline-x86.${ASM_EXT} test/asm/trampoline-x86.pl)
@@ -404,6 +408,11 @@
   x509v3/v3_sxnet.c
   x509v3/v3_utl.c
   ../third_party/fiat/curve25519.c
+  ../third_party/sike/fpx.c
+  ../third_party/sike/isogeny.c
+  ../third_party/sike/P503.c
+  ../third_party/sike/sike.c
+  ../third_party/sike/asm/fp_generic.c
 
   $<TARGET_OBJECTS:fipsmodule>
 
@@ -489,6 +498,7 @@
   x509/x509_time_test.cc
   x509v3/tab_test.cc
   x509v3/v3name_test.cc
+  ../third_party/sike/sike_test.cc
 
   $<TARGET_OBJECTS:crypto_test_data>
   $<TARGET_OBJECTS:boringssl_gtest_main>
diff --git a/third_party/sike/LICENSE b/third_party/sike/LICENSE
new file mode 100644
index 0000000..5cf7c8d
--- /dev/null
+++ b/third_party/sike/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) Microsoft Corporation. All rights reserved.
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE
diff --git a/third_party/sike/P503.c b/third_party/sike/P503.c
new file mode 100644
index 0000000..aecf623
--- /dev/null
+++ b/third_party/sike/P503.c
@@ -0,0 +1,70 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: supersingular isogeny parameters and generation of functions for P503
+*********************************************************************************************/
+
+#include "utils.h"
+
+// Parameters for isogeny system "SIKEp503"
+const struct params_t p503 = {
+    .prime = {
+        0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0xABFFFFFFFFFFFFFF,
+        0x13085BDA2211E7A0, 0x1B9BF6C87B7E7DAF, 0x6045C6BDDA77A4D0, 0x004066F541811E1E
+    },
+    .prime_p1 = {
+        0x0000000000000000, 0x0000000000000000, 0x0000000000000000, 0xAC00000000000000,
+        0x13085BDA2211E7A0, 0x1B9BF6C87B7E7DAF, 0x6045C6BDDA77A4D0, 0x004066F541811E1E
+    },
+    .prime_x2 = {
+        0xFFFFFFFFFFFFFFFE, 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0x57FFFFFFFFFFFFFF,
+        0x2610B7B44423CF41, 0x3737ED90F6FCFB5E, 0xC08B8D7BB4EF49A0, 0x0080CDEA83023C3C
+    },
+    .A_gen = {
+        0xE7EF4AA786D855AF, 0xED5758F03EB34D3B, 0x09AE172535A86AA9, 0x237B9CC07D622723,
+        0xE3A284CBA4E7932D, 0x27481D9176C5E63F, 0x6A323FF55C6E71BF, 0x002ECC31A6FB8773,   // XPA0
+        0x64D02E4E90A620B8, 0xDAB8128537D4B9F1, 0x4BADF77B8A228F98, 0x0F5DBDF9D1FB7D1B,
+        0xBEC4DB288E1A0DCC, 0xE76A8665E80675DB, 0x6D6F252E12929463, 0x003188BD1463FACC,   // XPA1
+        0xB79D41025DE85D56, 0x0B867DA9DF169686, 0x740E5368021C827D, 0x20615D72157BF25C,
+        0xFF1590013C9B9F5B, 0xC884DCADE8C16CEA, 0xEBD05E53BF724E01, 0x0032FEF8FDA5748C,   // XQA0
+        0x12E2E849AA0A8006, 0x41CF47008635A1E8, 0x9CD720A70798AED7, 0x42A820B42FCF04CF,
+        0x7BF9BAD32AAE88B1, 0xF619127A54090BBE, 0x1CB10D8F56408EAA, 0x001D6B54C3C0EDEB,   // XRA0
+        0x34DB54931CBAAC36, 0x420A18CB8DD5F0C4, 0x32008C1A48C0F44D, 0x3B3BA772B1CFD44D,
+        0xA74B058FDAF13515, 0x095FC9CA7EEC17B4, 0x448E829D28F120F8, 0x00261EC3ED16A489    // XRA1
+    },
+    .B_gen = {
+        0x7EDE37F4FA0BC727, 0xF7F8EC5C8598941C, 0xD15519B516B5F5C8, 0xF6D5AC9B87A36282,
+        0x7B19F105B30E952E, 0x13BD8B2025B4EBEE, 0x7B96D27F4EC579A2, 0x00140850CAB7E5DE,   // XPB0
+        0x7764909DAE7B7B2D, 0x578ABB16284911AB, 0x76E2BFD146A6BF4D, 0x4824044B23AA02F0,
+        0x1105048912A321F3, 0xB8A2E482CF0F10C1, 0x42FF7D0BE2152085, 0x0018E599C5223352,   // XPB1
+        0x4256C520FB388820, 0x744FD7C3BAAF0A13, 0x4B6A2DDDB12CBCB8, 0xE46826E27F427DF8,
+        0xFE4A663CD505A61B, 0xD6B3A1BAF025C695, 0x7C3BB62B8FCC00BD, 0x003AFDDE4A35746C,   // XQB0
+        0x75601CD1E6C0DFCB, 0x1A9007239B58F93E, 0xC1F1BE80C62107AC, 0x7F513B898F29FF08,
+        0xEA0BEDFF43E1F7B2, 0x2C6D94018CBAE6D0, 0x3A430D31BCD84672, 0x000D26892ECCFE83,   // XRB0
+        0x1119D62AEA3007A1, 0xE3702AA4E04BAE1B, 0x9AB96F7D59F990E7, 0xF58440E8B43319C0,
+        0xAF8134BEE1489775, 0xE7F7774E905192AA, 0xF54AE09308E98039, 0x001EF7A041A86112    // XRB1
+    },
+    .mont_R2 = {
+        0x5289A0CF641D011F, 0x9B88257189FED2B9, 0xA3B365D58DC8F17A, 0x5BC57AB6EFF168EC,
+        0x9E51998BD84D4423, 0xBF8999CBAC3B5695, 0x46E9127BCE14CDB6, 0x003F6CFCE8B81771
+    },
+    .mont_one = {
+        0x00000000000003F9, 0x0000000000000000, 0x0000000000000000, 0xB400000000000000,
+        0x63CB1A6EA6DED2B4, 0x51689D8D667EB37D, 0x8ACD77C71AB24142, 0x0026FBAEC60F5953
+    },
+    .A_strat = {
+        61, 32, 16, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1, 8, 4, 2, 1, 1, 2, 1, 1,
+        4, 2, 1, 1, 2, 1, 1, 16, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1, 8, 4, 2, 1,
+        1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1, 29, 16, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1,
+        1, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1, 13, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2,
+        1, 1, 2, 1, 1, 5, 4, 2, 1, 1, 2, 1, 1, 2, 1, 1, 1
+    },
+    .B_strat = {
+        71, 38, 21, 13, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1, 5, 4, 2, 1, 1, 2, 1,
+        1, 2, 1, 1, 1, 9, 5, 3, 2, 1, 1, 1, 1, 2, 1, 1, 1, 4, 2, 1, 1, 1, 2, 1, 1, 17, 9,
+        5, 3, 2, 1, 1, 1, 1, 2, 1, 1, 1, 4, 2, 1, 1, 1, 2, 1, 1, 8, 4, 2, 1, 1, 1, 2, 1,
+        1, 4, 2, 1, 1, 2, 1, 1, 33, 17, 9, 5, 3, 2, 1, 1, 1, 1, 2, 1, 1, 1, 4, 2, 1, 1, 1,
+        2, 1, 1, 8, 4, 2, 1, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1, 16, 8, 4, 2, 1, 1, 1, 2,
+        1, 1, 4, 2, 1, 1, 2, 1, 1, 8, 4, 2, 1, 1, 2, 1, 1, 4, 2, 1, 1, 2, 1, 1
+    }
+};
diff --git a/third_party/sike/asm/fp-armv8.pl b/third_party/sike/asm/fp-armv8.pl
new file mode 100644
index 0000000..a1728d1
--- /dev/null
+++ b/third_party/sike/asm/fp-armv8.pl
@@ -0,0 +1,968 @@
+#! /usr/bin/env perl
+#
+# April 2019
+#
+# Abstract: field arithmetic in aarch64 assembly for SIDH/p503
+
+$flavour = shift;
+$output  = shift;
+if ($flavour =~ /\./) { $output = $flavour; undef $flavour; }
+
+$0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
+( $xlate="${dir}arm-xlate.pl" and -f $xlate ) or
+( $xlate="${dir}../../../crypto/perlasm/arm-xlate.pl" and -f $xlate) or
+die "can't locate arm-xlate.pl";
+
+open OUT,"| \"$^X\" \"$xlate\" $flavour \"$output\"";
+*STDOUT=*OUT;
+
+$PREFIX="sike";
+
+$code.=<<___;
+.section  .rodata
+
+.Lp503p1_nz_s8:
+    .quad  0x085BDA2211E7A0AC, 0x9BF6C87B7E7DAF13
+    .quad  0x45C6BDDA77A4D01B, 0x4066F541811E1E60
+
+.Lp503x2:
+    .quad  0xFFFFFFFFFFFFFFFE, 0xFFFFFFFFFFFFFFFF
+    .quad  0x57FFFFFFFFFFFFFF, 0x2610B7B44423CF41
+    .quad  0x3737ED90F6FCFB5E, 0xC08B8D7BB4EF49A0
+    .quad  0x0080CDEA83023C3C
+
+.text
+___
+
+# C[0-2] = A[0] * B[0-1]
+sub mul64x128_comba_cut {
+    my ($A0,$B0,$B1,$C0,$C1,$C2,$T0,$T1)=@_;
+    my $body=<<___;
+        mul     $T1, $A0, $B0
+        umulh   $B0, $A0, $B0
+        adds    $C0, $C0, $C2
+        adc     $C1, $C1, xzr
+
+        mul     $T0, $A0, $B1
+        umulh   $B1, $A0, $B1
+        adds    $C0, $C0, $T1
+        adcs    $C1, $C1, $B0
+        adc     $C2, xzr, xzr
+
+        adds    $C1, $C1, $T0
+        adc     $C2, $C2, $B1
+___
+    return $body;
+}
+
+sub mul256_karatsuba_comba {
+    my ($M,$A0,$A1,$A2,$A3,$B0,$B1,$B2,$B3,$C0,$C1,$C2,$C3,$C4,$C5,$C6,$C7,$T0,$T1)=@_;
+    # (AH+AL) x (BH+BL), low part
+    my $mul_low=&mul64x128_comba_cut($A1, $C6, $T1, $C3, $C4, $C5, $C7, $A0);
+    # AL x BL
+    my $mul_albl=&mul64x128_comba_cut($A1, $B0, $B1, $C1, $T1, $C7, $C6, $A0);
+    # AH x BH
+    my $mul_ahbh=&mul64x128_comba_cut($A3, $B2, $B3, $A1, $C6, $B0, $B1, $A2);
+    my $body=<<___;
+        // A0-A1 <- AH + AL, T0 <- mask
+        adds    $A0, $A0, $A2
+        adcs    $A1, $A1, $A3
+        adc     $T0, xzr, xzr
+
+        // C6, T1 <- BH + BL, C7 <- mask
+        adds    $C6, $B0, $B2
+        adcs    $T1, $B1, $B3
+        adc     $C7, xzr, xzr
+
+        // C0-C1 <- masked (BH + BL)
+        sub     $C2, xzr, $T0
+        sub     $C3, xzr, $C7
+        and     $C0, $C6, $C2
+        and     $C1, $T1, $C2
+
+        // C4-C5 <- masked (AH + AL), T0 <- combined carry
+        and     $C4, $A0, $C3
+        and     $C5, $A1, $C3
+        mul     $C2, $A0, $C6
+        mul     $C3, $A0, $T1
+        and     $T0, $T0, $C7
+
+        // C0-C1, T0 <- (AH+AL) x (BH+BL), part 1
+        adds    $C0, $C4, $C0
+        umulh   $C4, $A0, $T1
+        adcs    $C1, $C5, $C1
+        umulh   $C5, $A0, $C6
+        adc     $T0, $T0, xzr
+
+        // C2-C5 <- (AH+AL) x (BH+BL), low part
+        $mul_low
+        ldp     $A0, $A1, [$M,#0]
+
+        // C2-C5, T0 <- (AH+AL) x (BH+BL), final part
+        adds    $C4, $C0, $C4
+        umulh   $C7, $A0, $B0
+        umulh   $T1, $A0, $B1
+        adcs    $C5, $C1, $C5
+        mul     $C0, $A0, $B0
+        mul     $C1, $A0, $B1
+        adc     $T0, $T0, xzr
+
+        // C0-C1, T1, C7 <- AL x BL
+        $mul_albl
+
+        // C2-C5, T0 <- (AH+AL) x (BH+BL) - ALxBL
+        mul     $A0, $A2, $B2
+        umulh   $B0, $A2, $B2
+        subs    $C2, $C2, $C0
+        sbcs    $C3, $C3, $C1
+        sbcs    $C4, $C4, $T1
+        mul     $A1, $A2, $B3
+        umulh   $C6, $A2, $B3
+        sbcs    $C5, $C5, $C7
+        sbc     $T0, $T0, xzr
+
+        // A0, A1, C6, B0 <- AH x BH
+        $mul_ahbh
+
+        // C2-C5, T0 <- (AH+AL) x (BH+BL) - ALxBL - AHxBH
+        subs    $C2, $C2, $A0
+        sbcs    $C3, $C3, $A1
+        sbcs    $C4, $C4, $C6
+        sbcs    $C5, $C5, $B0
+        sbc     $T0, $T0, xzr
+
+        adds    $C2, $C2, $T1
+        adcs    $C3, $C3, $C7
+        adcs    $C4, $C4, $A0
+        adcs    $C5, $C5, $A1
+        adcs    $C6, $T0, $C6
+        adc     $C7, $B0, xzr
+___
+    return $body;
+}
+
+# 512-bit integer multiplication using Karatsuba (two levels),
+# Comba (lower level).
+# Operation: c [x2] = a [x0] * b [x1]
+sub mul {
+    # (AH+AL) x (BH+BL), low part
+    my $mul_kc_low=&mul256_karatsuba_comba(
+        "x2",                                           # M0
+        "x3","x4","x5","x6",                            # A0-A3
+        "x11","x12","x13","x14",                        # B0-B3
+        "x8","x9","x10","x20","x21","x22","x23","x24",  # C0-C7
+        "x25","x26");                                   # TMP
+    # AL x BL
+    my $mul_albl=&mul256_karatsuba_comba(
+        "x0",                                           # M0
+        "x3","x4","x5","x6",                            # A0-A3
+        "x11","x12","x13","x14",                        # B0-B3
+        "x21","x22","x23","x24","x25","x26","x27","x28",# C0-C7
+        "x8","x9");                                     # TMP
+    # AH x BH
+    my $mul_ahbh=&mul256_karatsuba_comba(
+        "x0",                                           # M0
+        "x3","x4","x5","x6",                            # A0-A3
+        "x11","x12","x13","x14",                        # B0-B3
+        "x21","x22","x23","x24","x25","x26","x27","x28",# C0-C7
+        "x8","x9");                                     # TMP
+
+    my $body=<<___;
+        .global ${PREFIX}_mpmul
+        .align 4
+        ${PREFIX}_mpmul:
+        stp     x29, x30, [sp,#-96]!
+        add     x29, sp, #0
+        stp     x19, x20, [sp,#16]
+        stp     x21, x22, [sp,#32]
+        stp     x23, x24, [sp,#48]
+        stp     x25, x26, [sp,#64]
+        stp     x27, x28, [sp,#80]
+
+        ldp     x3, x4, [x0]
+        ldp     x5, x6, [x0,#16]
+        ldp     x7, x8, [x0,#32]
+        ldp     x9, x10, [x0,#48]
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+        ldp     x15, x16, [x1,#32]
+        ldp     x17, x19, [x1,#48]
+
+        // x3-x7 <- AH + AL, x7 <- carry
+        adds    x3, x3, x7
+        adcs    x4, x4, x8
+        adcs    x5, x5, x9
+        adcs    x6, x6, x10
+        adc     x7, xzr, xzr
+
+        // x11-x14 <- BH + BL, x8 <- carry
+        adds    x11, x11, x15
+        adcs    x12, x12, x16
+        adcs    x13, x13, x17
+        adcs    x14, x14, x19
+        adc     x8, xzr, xzr
+
+        // x9 <- combined carry
+        and      x9, x7, x8
+        // x7-x8 <- mask
+        sub      x7, xzr, x7
+        sub      x8, xzr, x8
+
+
+        // x15-x19 <- masked (BH + BL)
+        and     x15, x11, x7
+        and     x16, x12, x7
+        and     x17, x13, x7
+        and     x19, x14, x7
+
+        // x20-x23 <- masked (AH + AL)
+        and     x20, x3, x8
+        and     x21, x4, x8
+        and     x22, x5, x8
+        and     x23, x6, x8
+
+        // x15-x19, x7 <- masked (AH+AL) + masked (BH+BL), step 1
+        adds    x15, x15, x20
+        adcs    x16, x16, x21
+        adcs    x17, x17, x22
+        adcs    x19, x19, x23
+        adc     x7, x9, xzr
+
+        // x8-x10,x20-x24 <- (AH+AL) x (BH+BL), low part
+        stp     x3, x4, [x2,#0]
+        $mul_kc_low
+
+        // x15-x19, x7 <- (AH+AL) x (BH+BL), final step
+        adds    x15, x15, x21
+        adcs    x16, x16, x22
+        adcs    x17, x17, x23
+        adcs    x19, x19, x24
+        adc     x7, x7, xzr
+
+        // Load AL
+        ldp     x3, x4, [x0]
+        ldp     x5, x6, [x0,#16]
+        // Load BL
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+
+        // Temporarily store x8,x9 in x2
+        stp     x8,x9, [x2,#0]
+        // x21-x28 <- AL x BL
+        $mul_albl
+        // Restore x8,x9
+        ldp     x8,x9, [x2,#0]
+
+        // x8-x10,x20,x15-x17,x19 <- maskd (AH+AL) x (BH+BL) - ALxBL
+        subs    x8, x8, x21
+        sbcs    x9, x9, x22
+        sbcs    x10, x10, x23
+        sbcs    x20, x20, x24
+        sbcs    x15, x15, x25
+        sbcs    x16, x16, x26
+        sbcs    x17, x17, x27
+        sbcs    x19, x19, x28
+        sbc     x7, x7, xzr
+
+        // Store ALxBL, low
+        stp     x21, x22, [x2]
+        stp     x23, x24, [x2,#16]
+
+        // Load AH
+        ldp     x3, x4, [x0,#32]
+        ldp     x5, x6, [x0,#48]
+        // Load BH
+        ldp     x11, x12, [x1,#32]
+        ldp     x13, x14, [x1,#48]
+
+        adds    x8, x8, x25
+        adcs    x9, x9, x26
+        adcs    x10, x10, x27
+        adcs    x20, x20, x28
+        adc     x1, xzr, xzr
+
+        add     x0, x0, #32
+        // Temporarily store x8,x9 in x2
+        stp     x8,x9, [x2,#32]
+        // x21-x28 <- AH x BH
+        $mul_ahbh
+        // Restore x8,x9
+        ldp     x8,x9, [x2,#32]
+
+        neg     x1, x1
+
+        // x8-x10,x20,x15-x17,x19 <- (AH+AL) x (BH+BL) - ALxBL - AHxBH
+        subs    x8, x8, x21
+        sbcs    x9, x9, x22
+        sbcs    x10, x10, x23
+        sbcs    x20, x20, x24
+        sbcs    x15, x15, x25
+        sbcs    x16, x16, x26
+        sbcs    x17, x17, x27
+        sbcs    x19, x19, x28
+        sbc     x7, x7, xzr
+
+        // Store (AH+AL) x (BH+BL) - ALxBL - AHxBH, low
+        stp     x8, x9, [x2,#32]
+        stp     x10, x20, [x2,#48]
+
+        adds    x1, x1, #1
+        adcs    x15, x15, x21
+        adcs    x16, x16, x22
+        adcs    x17, x17, x23
+        adcs    x19, x19, x24
+        adcs    x25, x7, x25
+        adcs    x26, x26, xzr
+        adcs    x27, x27, xzr
+        adc     x28, x28, xzr
+
+        stp     x15, x16, [x2,#64]
+        stp     x17, x19, [x2,#80]
+        stp     x25, x26, [x2,#96]
+        stp     x27, x28, [x2,#112]
+
+        ldp     x19, x20, [x29,#16]
+        ldp     x21, x22, [x29,#32]
+        ldp     x23, x24, [x29,#48]
+        ldp     x25, x26, [x29,#64]
+        ldp     x27, x28, [x29,#80]
+        ldp     x29, x30, [sp],#96
+        ret
+___
+    return $body;
+}
+$code.=&mul();
+
+# Computes C0-C4 = (A0-A1) * (B0-B3)
+# Inputs remain intact
+sub mul128x256_comba {
+    my ($A0,$A1,$B0,$B1,$B2,$B3,$C0,$C1,$C2,$C3,$C4,$T0,$T1,$T2,$T3)=@_;
+    my $body=<<___;
+        mul     $T0, $A1, $B0
+        umulh   $T1, $A1, $B0
+        adds    $C0, $C0, $C2
+        adc     $C1, $C1, xzr
+
+        mul     $T2, $A0, $B2
+        umulh   $T3, $A0, $B2
+        adds    $C0, $C0, $T0
+        adcs    $C1, $C1, $T1
+        adc     $C2, xzr, xzr
+
+        mul     $T0, $A1, $B1
+        umulh   $T1, $A1, $B1
+        adds    $C1, $C1, $T2
+        adcs    $C2, $C2, $T3
+        adc     $C3, xzr, xzr
+
+        mul     $T2, $A0, $B3
+        umulh   $T3, $A0, $B3
+        adds    $C1, $C1, $T0
+        adcs    $C2, $C2, $T1
+        adc     $C3, $C3, xzr
+
+        mul     $T0, $A1, $B2
+        umulh   $T1, $A1, $B2
+        adds    $C2, $C2, $T2
+        adcs    $C3, $C3, $T3
+        adc     $C4, xzr, xzr
+
+        mul     $T2, $A1, $B3
+        umulh   $T3, $A1, $B3
+        adds    $C2, $C2, $T0
+        adcs    $C3, $C3, $T1
+        adc     $C4, $C4, xzr
+        adds    $C3, $C3, $T2
+        adc     $C4, $C4, $T3
+
+___
+    return $body;
+}
+
+#  Montgomery reduction
+#  Based on method described in Faz-Hernandez et al. https://eprint.iacr.org/2017/1015
+#  Operation: mc [x1] = ma [x0]
+#  NOTE: ma=mc is not allowed
+sub rdc {
+    my $mul01=&mul128x256_comba(
+        "x2","x3",                  # A0-A1
+        "x24","x25","x26","x27",    # B0-B3
+        "x5","x6","x7","x8","x9",   # C0-B4
+        "x1","x10","x11","x19");    # TMP
+    my $mul23=&mul128x256_comba(
+        "x2","x3",                  # A0-A1
+        "x24","x25","x26","x27",    # B0-B3
+        "x5","x6","x7","x8","x9",   # C0-C4
+        "x1","x10","x11","x19");    # TMP
+    my $mul45=&mul128x256_comba(
+        "x12","x13",                # A0-A1
+        "x24","x25","x26","x27",    # B0-B3
+        "x5","x6","x7","x8","x9",   # C0-C4
+        "x1","x10","x11","x19");    # TMP
+    my $mul67=&mul128x256_comba(
+        "x14","x15",                # A0-A1
+        "x24","x25","x26","x27",    # B0-B3
+        "x5","x6","x7","x8","x9",   # C0-C4
+        "x1","x10","x11","x19");    # TMP
+    my $body=<<___;
+    .global ${PREFIX}_fprdc
+    .align 4
+    ${PREFIX}_fprdc:
+        stp     x29, x30, [sp, #-112]!
+        add     x29, sp, #0
+        stp     x19, x20, [sp,#16]
+        stp     x21, x22, [sp,#32]
+        stp     x23, x24, [sp,#48]
+        stp     x25, x26, [sp,#64]
+        stp     x27, x28, [sp,#80]
+        str     x1, [sp,#96]
+
+        ldp     x2, x3, [x0,#0]       // a[0-1]
+
+        // Load the prime constant
+        adrp    x23, :pg_hi21:.Lp503p1_nz_s8
+        add     x23, x23, :lo12:.Lp503p1_nz_s8
+        ldp     x24, x25, [x23, #0]
+        ldp     x26, x27, [x23, #16]
+
+        // a[0-1] x .Lp503p1_nz_s8 --> result: x4:x9
+        mul     x4, x2, x24           // a[0] x .Lp503p1_nz_s8[0]
+        umulh   x7, x2, x24
+        mul     x5, x2, x25           // a[0] x .Lp503p1_nz_s8[1]
+        umulh   x6, x2, x25
+
+        $mul01
+
+        ldp      x2,  x3, [x0,#16]     // a[2]
+        ldp     x12, x13, [x0,#32]
+        ldp     x14, x15, [x0,#48]
+
+        orr     x10, xzr, x9, lsr #8
+        lsl     x9, x9, #56
+        orr     x9, x9, x8, lsr #8
+        lsl     x8, x8, #56
+        orr     x8, x8, x7, lsr #8
+        lsl     x7, x7, #56
+        orr     x7, x7, x6, lsr #8
+        lsl     x6, x6, #56
+        orr     x6, x6, x5, lsr #8
+        lsl     x5, x5, #56
+        orr     x5, x5, x4, lsr #8
+        lsl     x4, x4, #56
+
+        adds     x3, x4,  x3          // a[3]
+        adcs    x12, x5, x12          // a[4]
+        adcs    x13, x6, x13
+        adcs    x14, x7, x14
+        adcs    x15, x8, x15
+        ldp     x16, x17, [x0,#64]
+        ldp     x28, x30, [x0,#80]
+        mul     x4,  x2, x24          // a[2] x .Lp503p1_nz_s8[0]
+        umulh   x7,  x2, x24
+        adcs    x16, x9, x16
+        adcs    x17, x10, x17
+        adcs    x28, xzr, x28
+        adcs    x30, xzr, x30
+        ldp     x20, x21, [x0,#96]
+        ldp     x22, x23, [x0,#112]
+        mul     x5,  x2, x25          // a[2] x .Lp503p1_nz_s8[1]
+        umulh   x6,  x2, x25
+        adcs    x20, xzr, x20
+        adcs    x21, xzr, x21
+        adcs    x22, xzr, x22
+        adc     x23, xzr, x23
+
+        // a[2-3] x .Lp503p1_nz_s8 --> result: x4:x9
+        $mul23
+
+        orr     x10, xzr, x9, lsr #8
+        lsl     x9, x9, #56
+        orr     x9, x9, x8, lsr #8
+        lsl     x8, x8, #56
+        orr     x8, x8, x7, lsr #8
+        lsl     x7, x7, #56
+        orr     x7, x7, x6, lsr #8
+        lsl     x6, x6, #56
+        orr     x6, x6, x5, lsr #8
+        lsl     x5, x5, #56
+        orr     x5, x5, x4, lsr #8
+        lsl     x4, x4, #56
+
+        adds    x13, x4, x13          // a[5]
+        adcs    x14, x5, x14          // a[6]
+        adcs    x15, x6, x15
+        adcs    x16, x7, x16
+        mul     x4, x12, x24          // a[4] x .Lp503p1_nz_s8[0]
+        umulh   x7, x12, x24
+        adcs    x17, x8, x17
+        adcs    x28, x9, x28
+        adcs    x30, x10, x30
+        adcs    x20, xzr, x20
+        mul     x5, x12, x25          // a[4] x .Lp503p1_nz_s8[1]
+        umulh   x6, x12, x25
+        adcs    x21, xzr, x21
+        adcs    x22, xzr, x22
+        adc     x23, xzr, x23
+
+        // a[4-5] x .Lp503p1_nz_s8 --> result: x4:x9
+        $mul45
+
+        orr     x10, xzr, x9, lsr #8
+        lsl     x9, x9, #56
+        orr     x9, x9, x8, lsr #8
+        lsl     x8, x8, #56
+        orr     x8, x8, x7, lsr #8
+        lsl     x7, x7, #56
+        orr     x7, x7, x6, lsr #8
+        lsl     x6, x6, #56
+        orr     x6, x6, x5, lsr #8
+        lsl     x5, x5, #56
+        orr     x5, x5, x4, lsr #8
+        lsl     x4, x4, #56
+
+        adds    x15, x4, x15          // a[7]
+        adcs    x16, x5, x16          // a[8]
+        adcs    x17, x6, x17
+        adcs    x28, x7, x28
+        mul     x4, x14, x24          // a[6] x .Lp503p1_nz_s8[0]
+        umulh   x7, x14, x24
+        adcs    x30, x8, x30
+        adcs    x20, x9, x20
+        adcs    x21, x10, x21
+        mul     x5, x14, x25          // a[6] x .Lp503p1_nz_s8[1]
+        umulh   x6, x14, x25
+        adcs    x22, xzr, x22
+        adc     x23, xzr, x23
+
+        // a[6-7] x .Lp503p1_nz_s8 --> result: x4:x9
+        $mul67
+
+        orr     x10, xzr, x9, lsr #8
+        lsl     x9, x9, #56
+        orr     x9, x9, x8, lsr #8
+        lsl     x8, x8, #56
+        orr     x8, x8, x7, lsr #8
+        lsl     x7, x7, #56
+        orr     x7, x7, x6, lsr #8
+        lsl     x6, x6, #56
+        orr     x6, x6, x5, lsr #8
+        lsl     x5, x5, #56
+        orr     x5, x5, x4, lsr #8
+        lsl     x4, x4, #56
+
+        adds    x17, x4, x17
+        adcs    x28, x5, x28
+        ldr     x1, [sp,#96]
+        adcs    x30, x6, x30
+        adcs    x20, x7, x20
+        stp     x16, x17, [x1,#0]     // Final result
+        stp     x28, x30, [x1,#16]
+        adcs    x21, x8, x21
+        adcs    x22, x9, x22
+        adc     x23, x10, x23
+        stp     x20, x21, [x1,#32]
+        stp     x22, x23, [x1,#48]
+
+        ldp     x19, x20, [x29,#16]
+        ldp     x21, x22, [x29,#32]
+        ldp     x23, x24, [x29,#48]
+        ldp     x25, x26, [x29,#64]
+        ldp     x27, x28, [x29,#80]
+        ldp     x29, x30, [sp],#112
+        ret
+
+___
+}
+
+$code.=&rdc();
+
+
+#  Field addition
+#  Operation: c [x2] = a [x0] + b [x1]
+$code.=<<___;
+    .global ${PREFIX}_fpadd
+    .align 4
+    ${PREFIX}_fpadd:
+        stp     x29,x30, [sp,#-16]!
+        add     x29, sp, #0
+
+        ldp     x3, x4,   [x0,#0]
+        ldp     x5, x6,   [x0,#16]
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+
+        // Add a + b
+        adds    x3, x3, x11
+        adcs    x4, x4, x12
+        adcs    x5, x5, x13
+        adcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#32]
+        ldp     x9, x10,  [x0,#48]
+        ldp     x11, x12, [x1,#32]
+        ldp     x13, x14, [x1,#48]
+        adcs    x7, x7, x11
+        adcs    x8, x8, x12
+        adcs    x9, x9, x13
+        adc     x10, x10, x14
+
+        //  Subtract 2xp503
+        adrp    x17, :pg_hi21:.Lp503x2
+        add     x17, x17, :lo12:.Lp503x2
+        ldp     x11, x12, [x17, #0]
+        ldp     x13, x14, [x17, #16]
+        subs    x3, x3, x11
+        sbcs    x4, x4, x12
+        sbcs    x5, x5, x12
+        sbcs    x6, x6, x13
+        sbcs    x7, x7, x14
+
+        ldp     x15, x16, [x17, #32]
+        ldr     x17,      [x17, #48]
+        sbcs    x8, x8, x15
+        sbcs    x9, x9, x16
+        sbcs    x10, x10, x17
+        sbc     x0, xzr, xzr    // x0 can be reused now
+
+        // Add 2xp503 anded with the mask in x0
+        and     x11, x11, x0
+        and     x12, x12, x0
+        and     x13, x13, x0
+        and     x14, x14, x0
+        and     x15, x15, x0
+        and     x16, x16, x0
+        and     x17, x17, x0
+
+        adds    x3, x3, x11
+        adcs    x4, x4, x12
+        adcs    x5, x5, x12
+        adcs    x6, x6, x13
+        adcs    x7, x7, x14
+        adcs    x8, x8, x15
+        adcs    x9, x9, x16
+        adc     x10, x10, x17
+
+        stp     x3, x4,  [x2,#0]
+        stp     x5, x6,  [x2,#16]
+        stp     x7, x8,  [x2,#32]
+        stp     x9, x10, [x2,#48]
+
+        ldp     x29, x30, [sp],#16
+        ret
+
+___
+
+#  Field subtraction
+#  Operation: c [x2] = a [x0] - b [x1]
+$code.=<<___;
+    .global ${PREFIX}_fpsub
+    .align 4
+    ${PREFIX}_fpsub:
+        stp     x29, x30, [sp,#-16]!
+        add     x29, sp, #0
+
+        ldp     x3, x4,   [x0,#0]
+        ldp     x5, x6,   [x0,#16]
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+
+        // Subtract a - b
+        subs    x3, x3, x11
+        sbcs    x4, x4, x12
+        sbcs    x5, x5, x13
+        sbcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#32]
+        ldp     x11, x12, [x1,#32]
+        sbcs    x7, x7, x11
+        sbcs    x8, x8, x12
+        ldp     x9, x10,  [x0,#48]
+        ldp     x11, x12, [x1,#48]
+        sbcs    x9, x9, x11
+        sbcs    x10, x10, x12
+        sbc     x17, xzr, xzr
+
+        // Add 2xp503 anded with the mask in x17
+        adrp    x16, :pg_hi21:.Lp503x2
+        add     x16, x16, :lo12:.Lp503x2
+
+        // First half
+        ldp     x11, x12, [x16, #0]
+        ldp     x13, x14, [x16, #16]
+        and     x11, x11, x17
+        and     x12, x12, x17
+        and     x13, x13, x17
+        adds    x3, x3, x11
+        adcs    x4, x4, x12
+        adcs    x5, x5, x12
+        adcs    x6, x6, x13
+        stp     x3, x4,  [x2,#0]
+        stp     x5, x6,  [x2,#16]
+
+        // Second half
+        ldp     x11, x12, [x16, #32]
+        ldr     x13,      [x16, #48]
+        and     x14, x14, x17
+        and     x11, x11, x17
+        and     x12, x12, x17
+        and     x13, x13, x17
+        adcs    x7, x7, x14
+        adcs    x8, x8, x11
+        adcs    x9, x9, x12
+        adc     x10, x10, x13
+        stp     x7, x8,  [x2,#32]
+        stp     x9, x10, [x2,#48]
+
+        ldp     x29, x30, [sp],#16
+        ret
+___
+
+# 503-bit multiprecision addition
+# Operation: c [x2] = a [x0] + b [x1]
+$code.=<<___;
+    .global ${PREFIX}_mpadd_asm
+    .align 4
+    ${PREFIX}_mpadd_asm:
+        stp     x29, x30, [sp,#-16]!
+        add     x29, sp, #0
+
+        ldp     x3, x4,   [x0,#0]
+        ldp     x5, x6,   [x0,#16]
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+
+        adds    x3, x3, x11
+        adcs    x4, x4, x12
+        adcs    x5, x5, x13
+        adcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#32]
+        ldp     x9, x10,  [x0,#48]
+        ldp     x11, x12, [x1,#32]
+        ldp     x13, x14, [x1,#48]
+        adcs    x7, x7, x11
+        adcs    x8, x8, x12
+        adcs    x9, x9, x13
+        adc     x10, x10, x14
+
+        stp     x3, x4,   [x2,#0]
+        stp     x5, x6,   [x2,#16]
+        stp     x7, x8,   [x2,#32]
+        stp     x9, x10,  [x2,#48]
+
+        ldp     x29, x30, [sp],#16
+        ret
+___
+
+
+# 2x503-bit multiprecision addition
+# Operation: c [x2] = a [x0] + b [x1]
+$code.=<<___;
+    .global ${PREFIX}_mpadd503x2_asm
+    .align 4
+    ${PREFIX}_mpadd503x2_asm:
+        stp     x29, x30, [sp,#-16]!
+        add     x29, sp, #0
+
+        ldp     x3, x4,   [x0,#0]
+        ldp     x5, x6,   [x0,#16]
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+        adds    x3, x3, x11
+        adcs    x4, x4, x12
+        adcs    x5, x5, x13
+        adcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#32]
+        ldp     x9, x10,  [x0,#48]
+        ldp     x11, x12, [x1,#32]
+        ldp     x13, x14, [x1,#48]
+        adcs    x7, x7, x11
+        adcs    x8, x8, x12
+        adcs    x9, x9, x13
+        adcs    x10, x10, x14
+
+        stp     x3, x4,   [x2,#0]
+        stp     x5, x6,   [x2,#16]
+        stp     x7, x8,   [x2,#32]
+        stp     x9, x10,  [x2,#48]
+
+        ldp     x3, x4,   [x0,#64]
+        ldp     x5, x6,   [x0,#80]
+        ldp     x11, x12, [x1,#64]
+        ldp     x13, x14, [x1,#80]
+        adcs    x3, x3, x11
+        adcs    x4, x4, x12
+        adcs    x5, x5, x13
+        adcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#96]
+        ldp     x9, x10,  [x0,#112]
+        ldp     x11, x12, [x1,#96]
+        ldp     x13, x14, [x1,#112]
+        adcs    x7, x7, x11
+        adcs    x8, x8, x12
+        adcs    x9, x9, x13
+        adc     x10, x10, x14
+
+        stp     x3, x4,   [x2,#64]
+        stp     x5, x6,   [x2,#80]
+        stp     x7, x8,   [x2,#96]
+        stp     x9, x10,  [x2,#112]
+
+        ldp     x29, x30, [sp],#16
+        ret
+___
+
+
+
+# 2x503-bit multiprecision subtraction
+# Operation: c [x2] = a [x0] - b [x1].
+# Returns borrow mask
+$code.=<<___;
+    .global ${PREFIX}_mpsubx2_asm
+    .align 4
+    ${PREFIX}_mpsubx2_asm:
+        stp     x29, x30, [sp,#-16]!
+        add     x29, sp, #0
+
+        ldp     x3, x4,   [x0,#0]
+        ldp     x5, x6,   [x0,#16]
+        ldp     x11, x12, [x1,#0]
+        ldp     x13, x14, [x1,#16]
+        subs    x3, x3, x11
+        sbcs    x4, x4, x12
+        sbcs    x5, x5, x13
+        sbcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#32]
+        ldp     x9, x10,  [x0,#48]
+        ldp     x11, x12, [x1,#32]
+        ldp     x13, x14, [x1,#48]
+        sbcs    x7, x7, x11
+        sbcs    x8, x8, x12
+        sbcs    x9, x9, x13
+        sbcs    x10, x10, x14
+
+        stp     x3, x4,   [x2,#0]
+        stp     x5, x6,   [x2,#16]
+        stp     x7, x8,   [x2,#32]
+        stp     x9, x10,  [x2,#48]
+
+        ldp     x3, x4,   [x0,#64]
+        ldp     x5, x6,   [x0,#80]
+        ldp     x11, x12, [x1,#64]
+        ldp     x13, x14, [x1,#80]
+        sbcs    x3, x3, x11
+        sbcs    x4, x4, x12
+        sbcs    x5, x5, x13
+        sbcs    x6, x6, x14
+        ldp     x7, x8,   [x0,#96]
+        ldp     x9, x10,  [x0,#112]
+        ldp     x11, x12, [x1,#96]
+        ldp     x13, x14, [x1,#112]
+        sbcs    x7, x7, x11
+        sbcs    x8, x8, x12
+        sbcs    x9, x9, x13
+        sbcs    x10, x10, x14
+        sbc     x0, xzr, xzr
+
+        stp     x3, x4,   [x2,#64]
+        stp     x5, x6,   [x2,#80]
+        stp     x7, x8,   [x2,#96]
+        stp     x9, x10,  [x2,#112]
+
+        ldp     x29, x30, [sp],#16
+        ret
+___
+
+
+# Double 2x503-bit multiprecision subtraction
+# Operation: c [x2] = c [x2] - a [x0] - b [x1]
+$code.=<<___;
+    .global ${PREFIX}_mpdblsubx2_asm
+    .align 4
+    ${PREFIX}_mpdblsubx2_asm:
+        stp     x29, x30, [sp, #-64]!
+        add     x29, sp, #0
+
+        stp     x20, x21, [sp, #16]
+        stp     x22, x23, [sp, #32]
+        str     x24,      [sp, #48]
+
+        ldp     x3, x4,   [x2,#0]
+        ldp     x5, x6,   [x2,#16]
+        ldp     x7, x8,   [x2,#32]
+        ldp     x9, x10,  [x2,#48]
+        ldp     x11, x12, [x2,#64]
+        ldp     x13, x14, [x2,#80]
+        ldp     x15, x16, [x2,#96]
+        ldp     x17, x24, [x2,#112]
+
+        ldp     x20, x21, [x0,#0]
+        ldp     x22, x23, [x0,#16]
+        subs    x3, x3, x20
+        sbcs    x4, x4, x21
+        sbcs    x5, x5, x22
+        sbcs    x6, x6, x23
+        ldp     x20, x21, [x0,#32]
+        ldp     x22, x23, [x0,#48]
+        sbcs    x7, x7, x20
+        sbcs    x8, x8, x21
+        sbcs    x9, x9, x22
+        sbcs    x10, x10, x23
+        ldp     x20, x21, [x0,#64]
+        ldp     x22, x23, [x0,#80]
+        sbcs    x11, x11, x20
+        sbcs    x12, x12, x21
+        sbcs    x13, x13, x22
+        sbcs    x14, x14, x23
+        ldp     x20, x21, [x0,#96]
+        ldp     x22, x23, [x0,#112]
+        sbcs    x15, x15, x20
+        sbcs    x16, x16, x21
+        sbcs    x17, x17, x22
+        sbc     x24, x24, x23
+
+        ldp     x20, x21, [x1,#0]
+        ldp     x22, x23, [x1,#16]
+        subs    x3, x3, x20
+        sbcs    x4, x4, x21
+        sbcs    x5, x5, x22
+        sbcs    x6, x6, x23
+        ldp     x20, x21, [x1,#32]
+        ldp     x22, x23, [x1,#48]
+        sbcs    x7, x7, x20
+        sbcs    x8, x8, x21
+        sbcs    x9, x9, x22
+        sbcs    x10, x10, x23
+        ldp     x20, x21, [x1,#64]
+        ldp     x22, x23, [x1,#80]
+        sbcs    x11, x11, x20
+        sbcs    x12, x12, x21
+        sbcs    x13, x13, x22
+        sbcs    x14, x14, x23
+        ldp     x20, x21, [x1,#96]
+        ldp     x22, x23, [x1,#112]
+        sbcs    x15, x15, x20
+        sbcs    x16, x16, x21
+        sbcs    x17, x17, x22
+        sbc     x24, x24, x23
+
+        stp     x3, x4,   [x2,#0]
+        stp     x5, x6,   [x2,#16]
+        stp     x7, x8,   [x2,#32]
+        stp     x9, x10,  [x2,#48]
+        stp     x11, x12, [x2,#64]
+        stp     x13, x14, [x2,#80]
+        stp     x15, x16, [x2,#96]
+        stp     x17, x24, [x2,#112]
+
+        ldp     x20, x21, [x29,#16]
+        ldp     x22, x23, [x29,#32]
+        ldr     x24,      [x29,#48]
+
+        ldp     x29, x30, [sp],#64
+        ret
+___
+
+foreach (split("\n",$code)) {
+  s/\`([^\`]*)\`/eval($1)/ge;
+  print $_,"\n";
+}
+
+close STDOUT;
diff --git a/third_party/sike/asm/fp-x86_64.pl b/third_party/sike/asm/fp-x86_64.pl
new file mode 100755
index 0000000..c093c20
--- /dev/null
+++ b/third_party/sike/asm/fp-x86_64.pl
@@ -0,0 +1,1850 @@
+#! /usr/bin/env perl
+#
+# April 2019
+#
+# Abstract: field arithmetic in x64 assembly for SIDH/p503
+
+$flavour = shift;
+$output  = shift;
+if ($flavour =~ /\./) { $output = $flavour; undef $flavour; }
+
+$0 =~ m/(.*[\/\\])[^\/\\]+$/; $dir=$1;
+( $xlate="${dir}x86_64-xlate.pl" and -f $xlate ) or
+( $xlate="${dir}../../../crypto/perlasm/x86_64-xlate.pl" and -f $xlate) or
+die "can't locate x86_64-xlate.pl";
+
+open OUT,"| \"$^X\" \"$xlate\" $flavour \"$output\"";
+*STDOUT=*OUT;
+
+$PREFIX="sike";
+$bmi2_adx = 1;
+
+$code.=<<___;
+.text
+
+# p503 x 2
+.Lp503x2:
+.quad   0xFFFFFFFFFFFFFFFE
+.quad   0xFFFFFFFFFFFFFFFF
+.quad   0x57FFFFFFFFFFFFFF
+.quad   0x2610B7B44423CF41
+.quad   0x3737ED90F6FCFB5E
+.quad   0xC08B8D7BB4EF49A0
+.quad   0x0080CDEA83023C3C
+
+# p503 + 1
+.Lp503p1:
+.quad   0xAC00000000000000
+.quad   0x13085BDA2211E7A0
+.quad   0x1B9BF6C87B7E7DAF
+.quad   0x6045C6BDDA77A4D0
+.quad   0x004066F541811E1E
+
+.Lp503p1_nz:
+.quad    0xAC00000000000000
+.quad    0x13085BDA2211E7A0
+.quad    0x1B9BF6C87B7E7DAF
+.quad    0x6045C6BDDA77A4D0
+.quad    0x004066F541811E1E
+
+.extern OPENSSL_ia32cap_P
+.hidden OPENSSL_ia32cap_P
+
+___
+
+# Performs schoolbook multiplication of 128-bit with 320-bit
+# number. Uses MULX, ADOX, ADCX instruction.
+sub mul128x320_school {
+  my ($idxM0,$M0,$M1,$T0,$T1,$T2,$T3,$T4,$T5,$T6,$T7,$T8,$T9)=@_;
+  my ($MUL0,$MUL8)=map("$idxM0+$_(%$M0)", (0,8));
+  my $body.=<<___;
+    mov    $MUL0, %rdx
+    mulx   0+$M1, %$T0, %$T1       # T0 <- C0_final
+    mulx   8+$M1, %$T4, %$T2
+
+    xor    %rax, %rax
+    mulx   16+$M1, %$T5, %$T3
+    adox   %$T4, %$T1
+    adox   %$T5, %$T2
+    mulx   24+$M1, %$T7, %$T4
+    adox   %$T7, %$T3
+    mulx   32+$M1, %$T6, %$T5
+    adox   %$T6, %$T4
+    adox   %rax, %$T5
+
+    mov    $MUL8, %rdx
+    mulx   0+$M1, %$T6, %$T7
+    adcx   %$T6, %$T1               # T1 <- C1_final
+    adcx   %$T7, %$T2
+    mulx   8+$M1, %$T8, %$T6
+    adcx   %$T6, %$T3
+    mulx   16+$M1, %$T7, %$T9
+    adcx   %$T9, %$T4
+    mulx   24+$M1, %$T9, %$T6
+    adcx   %$T6, %$T5
+    mulx   32+$M1, %rdx, %$T6
+    adcx   %rax, %$T6
+
+    xor    %rax, %rax
+    adox   %$T8, %$T2
+    adox   %$T7, %$T3
+    adox   %$T9, %$T4
+    adox   %rdx, %$T5
+    adox   %rax, %$T6
+
+___
+  return $body;
+}
+
+# Compute z = x + y (mod p).
+# Operation: c [rdx] = a [rdi] + b [rsi]
+$code.=<<___;
+.globl  ${PREFIX}_fpadd
+.type   ${PREFIX}_fpadd,\@function,3
+${PREFIX}_fpadd:
+.cfi_startproc
+  push %r12
+.cfi_adjust_cfa_offset  8
+.cfi_offset r12, -16
+  push %r13
+.cfi_adjust_cfa_offset  8
+.cfi_offset r13, -24
+  push %r14
+.cfi_adjust_cfa_offset  8
+.cfi_offset r14, -32
+  push %r15
+.cfi_adjust_cfa_offset  8
+.cfi_offset r15, -40
+
+  xor  %rax, %rax
+
+  mov  0x0(%rdi),  %r8
+  mov  0x8(%rdi),  %r9
+  mov 0x10(%rdi), %r10
+  mov 0x18(%rdi), %r11
+  mov 0x20(%rdi), %r12
+  mov 0x28(%rdi), %r13
+  mov 0x30(%rdi), %r14
+  mov 0x38(%rdi), %r15
+
+  add  0x0(%rsi),  %r8
+  adc  0x8(%rsi),  %r9
+  adc 0x10(%rsi), %r10
+  adc 0x18(%rsi), %r11
+  adc 0x20(%rsi), %r12
+  adc 0x28(%rsi), %r13
+  adc 0x30(%rsi), %r14
+  adc 0x38(%rsi), %r15
+
+  mov .Lp503x2(%rip), %rcx;
+  sub %rcx, %r8
+  mov 8+.Lp503x2(%rip), %rcx;
+  sbb %rcx, %r9
+  sbb %rcx, %r10
+  mov 16+.Lp503x2(%rip), %rcx;
+  sbb %rcx, %r11
+  mov 24+.Lp503x2(%rip), %rcx;
+  sbb %rcx, %r12
+  mov 32+.Lp503x2(%rip), %rcx;
+  sbb %rcx, %r13
+  mov 40+.Lp503x2(%rip), %rcx;
+  sbb %rcx, %r14
+  mov 48+.Lp503x2(%rip), %rcx;
+  sbb %rcx, %r15
+  sbb \$0, %rax
+
+  mov .Lp503x2(%rip), %rdi
+  and %rax, %rdi
+  mov 8+.Lp503x2(%rip), %rsi
+  and %rax, %rsi
+  mov 16+.Lp503x2(%rip), %rcx
+  and %rax, %rcx
+
+  add %rdi, %r8
+  mov %r8, 0x0(%rdx)
+  adc %rsi, %r9
+  mov %r9, 0x8(%rdx)
+  adc %rsi, %r10
+  mov %r10, 0x10(%rdx)
+  adc %rcx, %r11
+  mov %r11, 0x18(%rdx)
+
+  setc   %cl
+
+  mov 24+.Lp503x2(%rip), %r8
+  and %rax, %r8
+  mov 32+.Lp503x2(%rip), %r9
+  and %rax, %r9
+  mov 40+.Lp503x2(%rip), %r10
+  and %rax, %r10
+  mov 48+.Lp503x2(%rip), %r11
+  and %rax, %r11
+
+  bt    \$0, %rcx
+
+  adc  %r8, %r12
+  mov %r12, 0x20(%rdx)
+  adc  %r9, %r13
+  mov %r13, 0x28(%rdx)
+  adc %r10, %r14
+  mov %r14, 0x30(%rdx)
+  adc %r11, %r15
+  mov %r15, 0x38(%rdx)
+
+  pop %r15
+.cfi_adjust_cfa_offset  -8
+  pop %r14
+.cfi_adjust_cfa_offset  -8
+  pop %r13
+.cfi_adjust_cfa_offset  -8
+  pop %r12
+.cfi_adjust_cfa_offset  -8
+  ret
+.cfi_endproc
+___
+
+
+
+# Loads data to XMM0 and XMM1 and
+# conditionaly swaps depending on XMM3
+sub cswap_block16() {
+  my $idx = shift;
+  $idx *= 16;
+  ("
+    movdqu   $idx(%rdi), %xmm0
+    movdqu   $idx(%rsi), %xmm1
+    movdqa   %xmm1, %xmm2
+    pxor     %xmm0, %xmm2
+    pand     %xmm3, %xmm2
+    pxor     %xmm2, %xmm0
+    pxor     %xmm2, %xmm1
+    movdqu   %xmm0, $idx(%rdi)
+    movdqu   %xmm1, $idx(%rsi)
+  ");
+}
+
+# Conditionally swaps bits in x and y in constant time.
+# mask indicates bits to be swapped (set bits are swapped)
+# Operation: [rdi] <-> [rsi] if rdx==1
+sub cswap {
+  # P[0].X with Q[0].X
+  foreach ( 0.. 3){$BLOCKS.=eval "&cswap_block16($_)";}
+  # P[0].Z with Q[0].Z
+  foreach ( 4.. 7){$BLOCKS.=eval "&cswap_block16($_)";}
+  # P[1].X with Q[1].X
+  foreach ( 8..11){$BLOCKS.=eval "&cswap_block16($_)";}
+  # P[1].Z with Q[1].Z
+  foreach (12..15){$BLOCKS.=eval "&cswap_block16($_)";}
+
+  my $body =<<___;
+.globl  ${PREFIX}_cswap_asm
+.type   ${PREFIX}_cswap_asm,\@function,3
+${PREFIX}_cswap_asm:
+  # Fill XMM3. After this step first half of XMM3 is
+  # just zeros and second half is whatever in RDX
+  mov   %rdx, %xmm3
+
+  # Copy lower double word everywhere else. So that
+  # XMM3=RDX|RDX. As RDX has either all bits set
+  # or non result will be that XMM3 has also either
+  # all bits set or non of them. 68 = 01000100b
+  pshufd  \$68, %xmm3, %xmm3
+  $BLOCKS
+  ret
+___
+  ($body)
+}
+$code.=&cswap();
+
+# Field subtraction
+# Operation: c [rdx] = a [rdi] - b [rsi]
+$code.=<<___;
+.globl  ${PREFIX}_fpsub
+.type   ${PREFIX}_fpsub,\@function,3
+${PREFIX}_fpsub:
+.cfi_startproc
+  push   %r12
+.cfi_adjust_cfa_offset  8
+.cfi_offset r12, -16
+  push   %r13
+.cfi_adjust_cfa_offset  8
+.cfi_offset r13, -24
+  push   %r14
+.cfi_adjust_cfa_offset  8
+.cfi_offset r14, -32
+  push   %r15
+.cfi_adjust_cfa_offset  8
+.cfi_offset r15, -40
+
+  xor %rax, %rax
+
+  mov  0x0(%rdi), %r8
+  mov  0x8(%rdi), %r9
+  mov 0x10(%rdi), %r10
+  mov 0x18(%rdi), %r11
+  mov 0x20(%rdi), %r12
+  mov 0x28(%rdi), %r13
+  mov 0x30(%rdi), %r14
+  mov 0x38(%rdi), %r15
+
+  sub  0x0(%rsi), %r8
+  sbb  0x8(%rsi), %r9
+  sbb 0x10(%rsi), %r10
+  sbb 0x18(%rsi), %r11
+  sbb 0x20(%rsi), %r12
+  sbb 0x28(%rsi), %r13
+  sbb 0x30(%rsi), %r14
+  sbb 0x38(%rsi), %r15
+  sbb \$0x0, %rax
+
+  mov .Lp503x2(%rip), %rdi
+  and %rax, %rdi
+  mov 0x8+.Lp503x2(%rip), %rsi
+  and %rax, %rsi
+  mov 0x10+.Lp503x2(%rip), %rcx
+  and %rax, %rcx
+
+  add %rdi,        %r8
+  adc %rsi,        %r9
+  adc %rsi,       %r10
+  adc %rcx,       %r11
+  mov %r8,   0x0(%rdx)
+  mov %r9,   0x8(%rdx)
+  mov %r10, 0x10(%rdx)
+  mov %r11, 0x18(%rdx)
+
+  setc %cl
+
+  mov 0x18+.Lp503x2(%rip),  %r8
+  and %rax,  %r8
+  mov 0x20+.Lp503x2(%rip),  %r9
+  and %rax,  %r9
+  mov 0x28+.Lp503x2(%rip), %r10
+  and %rax, %r10
+  mov 0x30+.Lp503x2(%rip), %r11
+  and %rax, %r11
+
+  bt \$0x0, %rcx
+
+  adc %r8, %r12
+  adc %r9, %r13
+  adc %r10, %r14
+  adc %r11, %r15
+  mov %r12, 0x20(%rdx)
+  mov %r13, 0x28(%rdx)
+  mov %r14, 0x30(%rdx)
+  mov %r15, 0x38(%rdx)
+
+  pop %r15
+.cfi_adjust_cfa_offset  -8
+  pop %r14
+.cfi_adjust_cfa_offset  -8
+  pop %r13
+.cfi_adjust_cfa_offset  -8
+  pop %r12
+.cfi_adjust_cfa_offset  -8
+  ret
+.cfi_endproc
+___
+
+#  503-bit multiprecision addition
+#  Operation: c [rdx] = a [rdi] + b [rsi]
+$code.=<<___;
+.globl  ${PREFIX}_mpadd_asm
+.type   ${PREFIX}_mpadd_asm,\@function,3
+${PREFIX}_mpadd_asm:
+.cfi_startproc
+  mov  0x0(%rdi), %r8
+  mov  0x8(%rdi), %r9
+  mov 0x10(%rdi), %r10
+  mov 0x18(%rdi), %r11
+  add  0x0(%rsi), %r8
+  adc  0x8(%rsi), %r9
+  adc 0x10(%rsi), %r10
+  adc 0x18(%rsi), %r11
+  mov %r8,   0x0(%rdx)
+  mov %r9,   0x8(%rdx)
+  mov %r10, 0x10(%rdx)
+  mov %r11, 0x18(%rdx)
+
+  mov 0x20(%rdi), %r8
+  mov 0x28(%rdi), %r9
+  mov 0x30(%rdi), %r10
+  mov 0x38(%rdi), %r11
+  adc 0x20(%rsi), %r8
+  adc 0x28(%rsi), %r9
+  adc 0x30(%rsi), %r10
+  adc 0x38(%rsi), %r11
+  mov %r8,  0x20(%rdx)
+  mov %r9,  0x28(%rdx)
+  mov %r10, 0x30(%rdx)
+  mov %r11, 0x38(%rdx)
+  ret
+.cfi_endproc
+___
+
+#  2x503-bit multiprecision subtraction
+#  Operation: c [rdx] = a [rdi] - b [rsi].
+#  Returns borrow mask
+$code.=<<___;
+.globl  ${PREFIX}_mpsubx2_asm
+.type   ${PREFIX}_mpsubx2_asm,\@function,3
+${PREFIX}_mpsubx2_asm:
+.cfi_startproc
+  xor %rax, %rax
+
+  mov  0x0(%rdi), %r8
+  mov  0x8(%rdi), %r9
+  mov 0x10(%rdi), %r10
+  mov 0x18(%rdi), %r11
+  mov 0x20(%rdi), %rcx
+  sub  0x0(%rsi), %r8
+  sbb  0x8(%rsi), %r9
+  sbb 0x10(%rsi), %r10
+  sbb 0x18(%rsi), %r11
+  sbb 0x20(%rsi), %rcx
+  mov %r8,   0x0(%rdx)
+  mov %r9,   0x8(%rdx)
+  mov %r10, 0x10(%rdx)
+  mov %r11, 0x18(%rdx)
+  mov %rcx, 0x20(%rdx)
+
+  mov 0x28(%rdi), %r8
+  mov 0x30(%rdi), %r9
+  mov 0x38(%rdi), %r10
+  mov 0x40(%rdi), %r11
+  mov 0x48(%rdi), %rcx
+  sbb 0x28(%rsi), %r8
+  sbb 0x30(%rsi), %r9
+  sbb 0x38(%rsi), %r10
+  sbb 0x40(%rsi), %r11
+  sbb 0x48(%rsi), %rcx
+  mov %r8,  0x28(%rdx)
+  mov %r9,  0x30(%rdx)
+  mov %r10, 0x38(%rdx)
+  mov %r11, 0x40(%rdx)
+  mov %rcx, 0x48(%rdx)
+
+  mov 0x50(%rdi), %r8
+  mov 0x58(%rdi), %r9
+  mov 0x60(%rdi), %r10
+  mov 0x68(%rdi), %r11
+  mov 0x70(%rdi), %rcx
+  sbb 0x50(%rsi), %r8
+  sbb 0x58(%rsi), %r9
+  sbb 0x60(%rsi), %r10
+  sbb 0x68(%rsi), %r11
+  sbb 0x70(%rsi), %rcx
+  mov %r8,  0x50(%rdx)
+  mov %r9,  0x58(%rdx)
+  mov %r10, 0x60(%rdx)
+  mov %r11, 0x68(%rdx)
+  mov %rcx, 0x70(%rdx)
+
+  mov 0x78(%rdi), %r8
+  sbb 0x78(%rsi), %r8
+  sbb \$0x0, %rax
+  mov %r8, 0x78(%rdx)
+  ret
+.cfi_endproc
+___
+
+#  Double 2x503-bit multiprecision subtraction
+#  Operation: c [rdx] = c [rdx] - a [rdi] - b [rsi]
+$code.=<<___;
+.globl  ${PREFIX}_mpdblsubx2_asm
+.type   ${PREFIX}_mpdblsubx2_asm,\@function,3
+${PREFIX}_mpdblsubx2_asm:
+.cfi_startproc
+  push   %r12
+.cfi_adjust_cfa_offset 8
+.cfi_offset r12, -16
+  push   %r13
+.cfi_adjust_cfa_offset 8
+.cfi_offset r13, -24
+  push   %r14
+.cfi_adjust_cfa_offset 8
+.cfi_offset r14, -32
+
+  xor %rax, %rax
+
+  mov  0x0(%rdx), %r8
+  mov  0x8(%rdx), %r9
+  mov 0x10(%rdx), %r10
+  mov 0x18(%rdx), %r11
+  mov 0x20(%rdx), %r12
+  mov 0x28(%rdx), %r13
+  mov 0x30(%rdx), %r14
+  mov 0x38(%rdx), %rcx
+  sub  0x0(%rdi), %r8
+  sbb  0x8(%rdi), %r9
+  sbb 0x10(%rdi), %r10
+  sbb 0x18(%rdi), %r11
+  sbb 0x20(%rdi), %r12
+  sbb 0x28(%rdi), %r13
+  sbb 0x30(%rdi), %r14
+  sbb 0x38(%rdi), %rcx
+  adc \$0x0, %rax
+
+  sub  0x0(%rsi), %r8
+  sbb  0x8(%rsi), %r9
+  sbb 0x10(%rsi), %r10
+  sbb 0x18(%rsi), %r11
+  sbb 0x20(%rsi), %r12
+  sbb 0x28(%rsi), %r13
+  sbb 0x30(%rsi), %r14
+  sbb 0x38(%rsi), %rcx
+  adc \$0x0, %rax
+
+  mov %r8,   0x0(%rdx)
+  mov %r9,   0x8(%rdx)
+  mov %r10, 0x10(%rdx)
+  mov %r11, 0x18(%rdx)
+  mov %r12, 0x20(%rdx)
+  mov %r13, 0x28(%rdx)
+  mov %r14, 0x30(%rdx)
+  mov %rcx, 0x38(%rdx)
+
+  mov 0x40(%rdx), %r8
+  mov 0x48(%rdx), %r9
+  mov 0x50(%rdx), %r10
+  mov 0x58(%rdx), %r11
+  mov 0x60(%rdx), %r12
+  mov 0x68(%rdx), %r13
+  mov 0x70(%rdx), %r14
+  mov 0x78(%rdx), %rcx
+
+  sub %rax, %r8
+  sbb 0x40(%rdi), %r8
+  sbb 0x48(%rdi), %r9
+  sbb 0x50(%rdi), %r10
+  sbb 0x58(%rdi), %r11
+  sbb 0x60(%rdi), %r12
+  sbb 0x68(%rdi), %r13
+  sbb 0x70(%rdi), %r14
+  sbb 0x78(%rdi), %rcx
+  sub 0x40(%rsi), %r8
+  sbb 0x48(%rsi), %r9
+  sbb 0x50(%rsi), %r10
+  sbb 0x58(%rsi), %r11
+  sbb 0x60(%rsi), %r12
+  sbb 0x68(%rsi), %r13
+  sbb 0x70(%rsi), %r14
+  sbb 0x78(%rsi), %rcx
+
+  mov %r8,  0x40(%rdx)
+  mov %r9,  0x48(%rdx)
+  mov %r10, 0x50(%rdx)
+  mov %r11, 0x58(%rdx)
+  mov %r12, 0x60(%rdx)
+  mov %r13, 0x68(%rdx)
+  mov %r14, 0x70(%rdx)
+  mov %rcx, 0x78(%rdx)
+
+  pop %r14
+.cfi_adjust_cfa_offset -8
+  pop %r13
+.cfi_adjust_cfa_offset -8
+  pop %r12
+.cfi_adjust_cfa_offset -8
+  ret
+.cfi_endproc
+
+___
+
+# Performs schoolbook multiplication of 2 256-bit numbers. Uses
+# MULX instruction. Result is stored in 256 bits pointed by $DST.
+sub mul256_school {
+  my ($idxM0,$M0,$idxM1,$M1,$idxDST,$DST,$T0,$T1,$T2,$T3,$T4,$T5,$T6,$T7,$T8,$T9)=@_;
+  my ($ML0,$ML8,$ML16,$ML24)=map("$idxM0+$_(%$M0)",(0,8,16,24));
+  my ($MR0,$MR8,$MR16,$MR24)=map("$idxM1+$_(%$M1)",(0,8,16,24));
+  my ($D0,$D1,$D2,$D3,$D4,$D5,$D6,$D7)=map("$idxDST+$_(%$DST)",(0,8,16,24,32,40,48,56));
+
+  $body=<<___;
+  mov    $ML0, %rdx
+  mulx   $MR0, %$T1, %$T0   # T0:T1 = A0*B0
+  mov    %$T1, $D0          # DST0_final
+  mulx   $MR8, %$T2, %$T1   # T1:T2 = A0*B1
+  xor    %rax, %rax
+  adox   %$T2, %$T0
+  mulx   $MR16,%$T3, %$T2   # T2:T3 = A0*B2
+  adox   %$T3, %$T1
+  mulx   $MR24,%$T4, %$T3   # T3:T4 = A0*B3
+  adox   %$T4, %$T2
+
+  mov    $ML8, %rdx
+  mulx   $MR0, %$T4, %$T5   # T5:T4 = A1*B0
+  adox   %rax, %$T3
+  xor    %rax, %rax
+  mulx   $MR8, %$T7, %$T6   # T6:T7 = A1*B1
+  adox   %$T0, %$T4
+  mov    %$T4, $D1          # DST1_final
+  adcx   %$T7, %$T5
+  mulx   $MR16,%$T8, %$T7   # T7:T8 = A1*B2
+  adcx   %$T8, %$T6
+  adox   %$T1, %$T5
+  mulx   $MR24,%$T9, %$T8   # T8:T9 = A1*B3
+  adcx   %$T9, %$T7
+  adcx   %rax, %$T8
+  adox   %$T2, %$T6
+
+  mov    $ML16,%rdx
+  mulx   $MR0, %$T0, %$T1   # T1:T0 = A2*B0
+  adox   %$T3, %$T7
+  adox   %rax, %$T8
+  xor    %rax, %rax
+  mulx   $MR8, %$T3, %$T2   # T2:T3 = A2*B1
+  adox   %$T5, %$T0
+  mov    %$T0, $D2          # DST2_final
+  adcx   %$T3, %$T1
+  mulx   $MR16,%$T4, %$T3   # T3:T4 = A2*B2
+  adcx   %$T4, %$T2
+  adox   %$T6, %$T1
+  mulx   $MR24,%$T9, %$T4   # T3:T4 = A2*B3
+  adcx   %$T9, %$T3
+
+  adcx   %rax, %$T4
+  adox   %$T7, %$T2
+  adox   %$T8, %$T3
+  adox   %rax, %$T4
+
+  mov    $ML24, %rdx
+  mulx   $MR0,  %$T0, %$T5   # T5:T0 = A3*B0
+  xor    %rax,  %rax
+  mulx   $MR8,  %$T7, %$T6   # T6:T7 = A3*B1
+  adcx   %$T7,  %$T5
+  adox   %$T0,  %$T1
+  mulx   $MR16, %$T8, %$T7   # T7:T8 = A3*B2
+  adcx   %$T8,  %$T6
+  adox   %$T5,  %$T2
+  mulx   $MR24, %$T9, %$T8   # T8:T9 = A3*B3
+  adcx   %$T9,  %$T7
+  adcx   %rax,  %$T8
+  adox   %$T6,  %$T3
+  adox   %$T7,  %$T4
+  adox   %rax,  %$T8
+  mov    %$T1,  $D3          # DST3_final
+  mov    %$T2,  $D4          # DST4_final
+  mov    %$T3,  $D5          # DST5_final
+  mov    %$T4,  $D6          # DST6_final
+  mov    %$T8,  $D7          # DST7_final
+
+___
+  return $body;
+}
+
+# 503-bit multiplication using Karatsuba (one level),
+# schoolbook (one level).
+sub mul_mulx {
+  # [rcx+64] <- (AH+AL) x (BH+BL)
+  my $mul256_low=&mul256_school(0,"rsp",32,"rsp",64,"rcx",map("r$_",(8..15)),"rbx","rbp");
+  # [rcx] <- AL x BL
+  my $mul256_albl=&mul256_school(0,"rdi",0,"rsi",0,"rcx",map("r$_",(8..15)),"rbx","rbp");
+  # [rsp] <- AH x BH
+  my $mul256_ahbh=&mul256_school(32,"rdi",32,"rsi",0,"rsp",map("r$_",(8..15)),"rbx","rbp");
+
+  $body=<<___;
+  .Lmul_mulx:
+  .cfi_startproc
+    # sike_mpmul has already pushed r12--15 by this point.
+  .cfi_adjust_cfa_offset 32
+  .cfi_offset r12, -16
+  .cfi_offset r13, -24
+  .cfi_offset r14, -32
+  .cfi_offset r15, -40
+
+    mov %rdx, %rcx
+
+    # r8-r11 <- AH + AL, rax <- mask
+    xor %rax, %rax
+    mov (%rdi), %r8
+    mov 0x8(%rdi), %r9
+    mov 0x10(%rdi), %r10
+    mov 0x18(%rdi), %r11
+    push %rbx
+
+  .cfi_adjust_cfa_offset 8
+  .cfi_offset rbx, -48
+    push %rbp
+  .cfi_offset rbp, -56
+  .cfi_adjust_cfa_offset 8
+    sub \$96, %rsp
+  .cfi_adjust_cfa_offset 96
+    add 0x20(%rdi), %r8
+    adc 0x28(%rdi), %r9
+    adc 0x30(%rdi), %r10
+    adc 0x38(%rdi), %r11
+    sbb \$0x0, %rax
+    mov %r8, (%rsp)
+    mov %r9, 0x8(%rsp)
+    mov %r10, 0x10(%rsp)
+    mov %r11, 0x18(%rsp)
+
+    # r12-r15 <- BH + BL, rbx <- mask
+    xor %rbx, %rbx
+    mov (%rsi), %r12
+    mov 0x8(%rsi), %r13
+    mov 0x10(%rsi), %r14
+    mov 0x18(%rsi), %r15
+    add 0x20(%rsi), %r12
+    adc 0x28(%rsi), %r13
+    adc 0x30(%rsi), %r14
+    adc 0x38(%rsi), %r15
+    sbb \$0x0, %rbx
+    mov %r12, 0x20(%rsp)
+    mov %r13, 0x28(%rsp)
+    mov %r14, 0x30(%rsp)
+    mov %r15, 0x38(%rsp)
+
+    # r12-r15 <- masked (BH + BL)
+    and %rax, %r12
+    and %rax, %r13
+    and %rax, %r14
+    and %rax, %r15
+
+    # r8-r11 <- masked (AH + AL)
+    and %rbx, %r8
+    and %rbx, %r9
+    and %rbx, %r10
+    and %rbx, %r11
+
+    # r8-r11 <- masked (AH + AL) + masked (AH + AL)
+    add %r12, %r8
+    adc %r13, %r9
+    adc %r14, %r10
+    adc %r15, %r11
+    mov %r8, 0x40(%rsp)
+    mov %r9, 0x48(%rsp)
+    mov %r10, 0x50(%rsp)
+    mov %r11, 0x58(%rsp)
+
+    # [rcx+64] <- (AH+AL) x (BH+BL)
+    $mul256_low
+    # [rcx] <- AL x BL (Result c0-c3)
+    $mul256_albl
+    # [rsp] <- AH x BH
+    $mul256_ahbh
+
+    # r8-r11 <- (AH+AL) x (BH+BL), final step
+    mov 0x40(%rsp), %r8
+    mov 0x48(%rsp), %r9
+    mov 0x50(%rsp), %r10
+    mov 0x58(%rsp), %r11
+    mov 0x60(%rcx), %rax
+    add %rax, %r8
+    mov 0x68(%rcx), %rax
+    adc %rax, %r9
+    mov 0x70(%rcx), %rax
+    adc %rax, %r10
+    mov 0x78(%rcx), %rax
+    adc %rax, %r11
+
+    # [rcx+64], x3-x5 <- (AH+AL) x (BH+BL) - ALxBL
+    mov 0x40(%rcx), %r12
+    mov 0x48(%rcx), %r13
+    mov 0x50(%rcx), %r14
+    mov 0x58(%rcx), %r15
+    sub (%rcx), %r12
+    sbb 0x8(%rcx), %r13
+    sbb 0x10(%rcx), %r14
+    sbb 0x18(%rcx), %r15
+    sbb 0x20(%rcx), %r8
+    sbb 0x28(%rcx), %r9
+    sbb 0x30(%rcx), %r10
+    sbb 0x38(%rcx), %r11
+
+    # r8-r15 <- (AH+AL) x (BH+BL) - ALxBL - AHxBH
+    sub (%rsp), %r12
+    sbb 0x8(%rsp), %r13
+    sbb 0x10(%rsp), %r14
+    sbb 0x18(%rsp), %r15
+    sbb 0x20(%rsp), %r8
+    sbb 0x28(%rsp), %r9
+    sbb 0x30(%rsp), %r10
+    sbb 0x38(%rsp), %r11
+
+    add 0x20(%rcx), %r12
+    mov %r12, 0x20(%rcx)    # Result C4-C7
+    adc 0x28(%rcx), %r13
+    mov %r13, 0x28(%rcx)
+    adc 0x30(%rcx), %r14
+    mov %r14, 0x30(%rcx)
+    adc 0x38(%rcx), %r15
+    mov %r15, 0x38(%rcx)
+    mov (%rsp), %rax
+    adc %rax, %r8           # Result C8-C15
+    mov %r8, 0x40(%rcx)
+    mov 0x8(%rsp), %rax
+    adc %rax, %r9
+    mov %r9, 0x48(%rcx)
+    mov 0x10(%rsp), %rax
+    adc %rax, %r10
+    mov %r10, 0x50(%rcx)
+    mov 0x18(%rsp), %rax
+    adc %rax, %r11
+    mov %r11, 0x58(%rcx)
+    mov 0x20(%rsp), %r12
+    adc \$0x0, %r12
+    mov %r12, 0x60(%rcx)
+    mov 0x28(%rsp), %r13
+    adc \$0x0, %r13
+    mov %r13, 0x68(%rcx)
+    mov 0x30(%rsp), %r14
+    adc \$0x0, %r14
+    mov %r14, 0x70(%rcx)
+    mov 0x38(%rsp), %r15
+    adc \$0x0, %r15
+    mov %r15, 0x78(%rcx)
+
+    add \$96, %rsp
+  .cfi_adjust_cfa_offset -96
+    pop %rbp
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value rbp
+    pop %rbx
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value rbx
+    pop %r15
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r15
+    pop %r14
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r14
+    pop %r13
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r13
+    pop %r12
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r12
+      ret
+  .cfi_endproc
+
+___
+  return $body;
+}
+
+# Jump to alternative implemenatation provided as an
+# argument in case CPU supports ADOX/ADCX and MULX instructions.
+sub alt_impl {
+  $jmp_func = shift;
+
+  $body=<<___;
+  lea OPENSSL_ia32cap_P(%rip), %rcx
+  mov 8(%rcx), %rcx
+  and \$0x80100, %ecx
+  cmp \$0x80100, %ecx
+  je  $jmp_func
+
+___
+  return $body
+}
+
+#  Integer multiplication based on Karatsuba method
+#  Operation: c [rdx] = a [rdi] * b [rsi]
+#  NOTE: a=c or b=c are not allowed
+sub mul {
+  my $jump_optim.=&alt_impl(".Lmul_mulx") if ($bmi2_adx);
+  my $body.=&mul_mulx() if ($bmi2_adx);
+
+  $body.=<<___;
+  .globl  ${PREFIX}_mpmul
+  .type   ${PREFIX}_mpmul,\@function,3
+  ${PREFIX}_mpmul:
+  .cfi_startproc
+    push %r12
+  .cfi_adjust_cfa_offset 8
+  .cfi_offset r12, -16
+    push %r13
+  .cfi_adjust_cfa_offset 8
+  .cfi_offset r13, -24
+    push %r14
+  .cfi_adjust_cfa_offset 8
+  .cfi_offset r14, -32
+    push %r15
+  .cfi_adjust_cfa_offset 8
+  .cfi_offset r15, -40
+
+    $jump_optim
+
+    mov %rdx, %rcx
+
+    # rcx[0-3] <- AH+AL
+    xor %rax, %rax
+    mov 0x20(%rdi), %r8
+    mov 0x28(%rdi), %r9
+    mov 0x30(%rdi), %r10
+    mov 0x38(%rdi), %r11
+    add  0x0(%rdi), %r8
+    adc  0x8(%rdi), %r9
+    adc 0x10(%rdi), %r10
+    adc 0x18(%rdi), %r11
+    mov %r8,   0x0(%rcx)
+    mov %r9,   0x8(%rcx)
+    mov %r10, 0x10(%rcx)
+    mov %r11, 0x18(%rcx)
+    sbb  \$0,  %rax
+    sub \$80,  %rsp           # Allocating space in stack
+  .cfi_adjust_cfa_offset 80
+
+    # r12-r15 <- BH+BL
+    xor %rdx, %rdx
+    mov 0x20(%rsi), %r12
+    mov 0x28(%rsi), %r13
+    mov 0x30(%rsi), %r14
+    mov 0x38(%rsi), %r15
+    add  0x0(%rsi), %r12
+    adc  0x8(%rsi), %r13
+    adc 0x10(%rsi), %r14
+    adc 0x18(%rsi), %r15
+    sbb \$0x0, %rdx
+    mov %rax, 0x40(%rsp)
+    mov %rdx, 0x48(%rsp)
+
+    # (rsp[0-3],r8,r9,r10,r11) <- (AH+AL)*(BH+BL)
+    mov (%rcx), %rax
+    mul %r12
+    mov %rax, (%rsp)            # c0
+    mov %rdx, %r8
+
+    xor %r9, %r9
+    mov (%rcx), %rax
+    mul %r13
+    add %rax, %r8
+    adc %rdx, %r9
+
+    xor %r10, %r10
+    mov 0x8(%rcx), %rax
+    mul %r12
+    add %rax, %r8
+    mov %r8, 0x8(%rsp)          # c1
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    xor %r8, %r8
+    mov (%rcx), %rax
+    mul %r14
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x10(%rcx), %rax
+    mul %r12
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x8(%rcx), %rax
+    mul %r13
+    add %rax, %r9
+    mov %r9, 0x10(%rsp)         # c2
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    xor %r9, %r9
+    mov (%rcx), %rax
+    mul %r15
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x18(%rcx), %rax
+    mul %r12
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x8(%rcx), %rax
+    mul %r14
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x10(%rcx), %rax
+    mul %r13
+    add %rax, %r10
+    mov %r10, 0x18(%rsp)        # c3
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    xor %r10, %r10
+    mov 0x8(%rcx), %rax
+    mul %r15
+    add %rax, %r8
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    mov 0x18(%rcx), %rax
+    mul %r13
+    add %rax, %r8
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    mov 0x10(%rcx), %rax
+    mul %r14
+    add %rax, %r8
+    mov %r8, 0x20(%rsp)          # c4
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    xor %r11, %r11
+    mov 0x10(%rcx), %rax
+    mul %r15
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r11
+
+    mov 0x18(%rcx), %rax
+    mul %r14
+    add %rax, %r9               # c5
+    adc %rdx, %r10
+    adc \$0x0, %r11
+
+    mov 0x18(%rcx), %rax
+    mul %r15
+    add %rax, %r10              # c6
+    adc %rdx, %r11              # c7
+
+    mov 0x40(%rsp), %rax
+    and %rax, %r12
+    and %rax, %r13
+    and %rax, %r14
+    and %rax, %r15
+    add %r8, %r12
+    adc %r9, %r13
+    adc %r10, %r14
+    adc %r11, %r15
+
+    mov 0x48(%rsp), %rax
+    mov (%rcx), %r8
+    mov 0x8(%rcx), %r9
+    mov 0x10(%rcx), %r10
+    mov 0x18(%rcx), %r11
+    and %rax, %r8
+    and %rax, %r9
+    and %rax, %r10
+    and %rax, %r11
+    add %r12, %r8
+    adc %r13, %r9
+    adc %r14, %r10
+    adc %r15, %r11
+    mov %r8, 0x20(%rsp)
+    mov %r9, 0x28(%rsp)
+    mov %r10, 0x30(%rsp)
+    mov %r11, 0x38(%rsp)
+
+    mov (%rdi), %r11
+    mov (%rsi), %rax
+    mul %r11
+    xor %r9, %r9
+    mov %rax, (%rcx)              # c0
+    mov %rdx, %r8
+
+    mov 0x10(%rdi), %r14
+    mov 0x8(%rsi), %rax
+    mul %r11
+    xor %r10, %r10
+    add %rax, %r8
+    adc %rdx, %r9
+
+    mov 0x8(%rdi), %r12
+    mov (%rsi), %rax
+    mul %r12
+    add %rax, %r8
+    mov %r8, 0x8(%rcx)            # c1
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    xor %r8, %r8
+    mov 0x10(%rsi), %rax
+    mul %r11
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov (%rsi), %r13
+    mov %r14, %rax
+    mul %r13
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x8(%rsi), %rax
+    mul %r12
+    add %rax, %r9
+    mov %r9, 0x10(%rcx)           # c2
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    xor %r9, %r9
+    mov 0x18(%rsi), %rax
+    mul %r11
+    mov 0x18(%rdi), %r15
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov %r15, %rax
+    mul %r13
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x10(%rsi), %rax
+    mul %r12
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x8(%rsi), %rax
+    mul %r14
+    add %rax, %r10
+    mov %r10, 0x18(%rcx)           # c3
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    xor %r10, %r10
+    mov 0x18(%rsi), %rax
+    mul %r12
+    add %rax, %r8
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    mov 0x8(%rsi), %rax
+    mul %r15
+    add %rax, %r8
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    mov 0x10(%rsi), %rax
+    mul %r14
+    add %rax, %r8
+    mov %r8, 0x20(%rcx)           # c4
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    xor %r8, %r8
+    mov 0x18(%rsi), %rax
+    mul %r14
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x10(%rsi), %rax
+    mul %r15
+    add %rax, %r9
+    mov %r9, 0x28(%rcx)           # c5
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x18(%rsi), %rax
+    mul %r15
+    add %rax, %r10
+    mov %r10, 0x30(%rcx)          # c6
+    adc %rdx, %r8
+    mov %r8, 0x38(%rcx)           # c7
+
+    # rcx[8-15] <- AH*BH
+    mov 0x20(%rdi), %r11
+    mov 0x20(%rsi), %rax
+    mul %r11
+    xor %r9, %r9
+    mov %rax, 0x40(%rcx)          # c0
+    mov %rdx, %r8
+
+    mov 0x30(%rdi), %r14
+    mov 0x28(%rsi), %rax
+    mul %r11
+    xor %r10, %r10
+    add %rax, %r8
+    adc %rdx, %r9
+
+    mov 0x28(%rdi), %r12
+    mov 0x20(%rsi), %rax
+    mul %r12
+    add %rax, %r8
+    mov %r8, 0x48(%rcx)           # c1
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    xor %r8, %r8
+    mov 0x30(%rsi), %rax
+    mul %r11
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x20(%rsi), %r13
+    mov %r14, %rax
+    mul %r13
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x28(%rsi), %rax
+    mul %r12
+    add %rax, %r9
+    mov %r9, 0x50(%rcx)             # c2
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    xor %r9, %r9
+    mov 0x38(%rsi), %rax
+    mul %r11
+    mov 0x38(%rdi), %r15
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov %r15, %rax
+    mul %r13
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x30(%rsi), %rax
+    mul %r12
+    add %rax, %r10
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    mov 0x28(%rsi), %rax
+    mul %r14
+    add %rax, %r10
+    mov %r10, 0x58(%rcx)            # c3
+    adc %rdx, %r8
+    adc \$0x0, %r9
+
+    xor %r10, %r10
+    mov 0x38(%rsi), %rax
+    mul %r12
+    add %rax, %r8
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    mov 0x28(%rsi), %rax
+    mul %r15
+    add %rax, %r8
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    mov 0x30(%rsi), %rax
+    mul %r14
+    add %rax, %r8
+    mov %r8, 0x60(%rcx)             # c4
+    adc %rdx, %r9
+    adc \$0x0, %r10
+
+    xor %r8, %r8
+    mov 0x38(%rsi), %rax
+    mul %r14
+    add %rax, %r9
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x30(%rsi), %rax
+    mul %r15
+    add %rax, %r9
+    mov %r9, 0x68(%rcx)             # c5
+    adc %rdx, %r10
+    adc \$0x0, %r8
+
+    mov 0x38(%rsi), %rax
+    mul %r15
+    add %rax, %r10
+    mov %r10, 0x70(%rcx)            # c6
+    adc %rdx, %r8
+    mov %r8, 0x78(%rcx)             # c7
+
+    # [r8-r15] <- (AH+AL)*(BH+BL) - AL*BL
+    mov  0x0(%rsp), %r8
+    sub  0x0(%rcx), %r8
+    mov  0x8(%rsp), %r9
+    sbb  0x8(%rcx), %r9
+    mov 0x10(%rsp), %r10
+    sbb 0x10(%rcx), %r10
+    mov 0x18(%rsp), %r11
+    sbb 0x18(%rcx), %r11
+    mov 0x20(%rsp), %r12
+    sbb 0x20(%rcx), %r12
+    mov 0x28(%rsp), %r13
+    sbb 0x28(%rcx), %r13
+    mov 0x30(%rsp), %r14
+    sbb 0x30(%rcx), %r14
+    mov 0x38(%rsp), %r15
+    sbb 0x38(%rcx), %r15
+
+    # [r8-r15] <- (AH+AL)*(BH+BL) - AL*BL - AH*BH
+    mov 0x40(%rcx), %rax
+    sub %rax, %r8
+    mov 0x48(%rcx), %rax
+    sbb %rax, %r9
+    mov 0x50(%rcx), %rax
+    sbb %rax, %r10
+    mov 0x58(%rcx), %rax
+    sbb %rax, %r11
+    mov 0x60(%rcx), %rax
+    sbb %rax, %r12
+    mov 0x68(%rcx), %rdx
+    sbb %rdx, %r13
+    mov 0x70(%rcx), %rdi
+    sbb %rdi, %r14
+    mov 0x78(%rcx), %rsi
+    sbb %rsi, %r15
+
+    # Final result
+    add 0x20(%rcx),  %r8
+    mov %r8,  0x20(%rcx)
+    adc 0x28(%rcx),  %r9
+    mov %r9,  0x28(%rcx)
+    adc 0x30(%rcx), %r10
+    mov %r10, 0x30(%rcx)
+    adc 0x38(%rcx), %r11
+    mov %r11, 0x38(%rcx)
+    adc 0x40(%rcx), %r12
+    mov %r12, 0x40(%rcx)
+    adc 0x48(%rcx), %r13
+    mov %r13, 0x48(%rcx)
+    adc 0x50(%rcx), %r14
+    mov %r14, 0x50(%rcx)
+    adc 0x58(%rcx), %r15
+    mov %r15, 0x58(%rcx)
+    adc \$0x0, %rax
+    mov %rax, 0x60(%rcx)
+    adc \$0x0, %rdx
+    mov %rdx, 0x68(%rcx)
+    adc \$0x0, %rdi
+    mov %rdi, 0x70(%rcx)
+    adc \$0x0, %rsi
+    mov %rsi, 0x78(%rcx)
+
+    add \$80, %rsp           # Restoring space in stack
+  .cfi_adjust_cfa_offset -80
+    pop %r15
+  .cfi_adjust_cfa_offset -8
+    pop %r14
+  .cfi_adjust_cfa_offset -8
+    pop %r13
+  .cfi_adjust_cfa_offset -8
+    pop %r12
+  .cfi_adjust_cfa_offset -8
+    ret
+  .cfi_endproc
+
+___
+  return $body;
+}
+
+$code.=&mul();
+
+#  Optimized Montgomery reduction for CPUs with ADOX/ADCX and MULX
+#  Based on method described in Faz-Hernandez et al. https://eprint.iacr.org/2017/1015
+#  Operation: c [rsi] = a [rdi]
+#  NOTE: a=c is not allowed
+sub rdc_mulx {
+  # a[0-1] x .Lp503p1_nz --> result: r8:r14
+  my $mul01=&mul128x320_school(0,"rdi",".Lp503p1_nz(%rip)",map("r$_",(8..14)),"rbx","rcx","r15");
+  # a[2-3] x .Lp503p1_nz --> result: r8:r14
+  my $mul23=&mul128x320_school(16,"rdi",".Lp503p1_nz(%rip)",map("r$_",(8..14)),"rbx","rcx","r15");
+  # a[4-5] x .Lp503p1_nz --> result: r8:r14
+  my $mul45=&mul128x320_school(32,"rdi",".Lp503p1_nz(%rip)",map("r$_",(8..14)),"rbx","rcx","r15");
+  # a[6-7] x .Lp503p1_nz --> result: r8:r14
+  my $mul67=&mul128x320_school(48,"rdi",".Lp503p1_nz(%rip)",map("r$_", (8..14)),"rbx","rcx","r15");
+
+  my $body=<<___;
+    .Lrdc_mulx_asm:
+    .cfi_startproc
+      # sike_fprdc has already pushed r12--15 and rbx by this point.
+    .cfi_adjust_cfa_offset 32
+    .cfi_offset r12, -16
+    .cfi_offset r13, -24
+    .cfi_offset r14, -32
+    .cfi_offset r15, -40
+    .cfi_offset rbx, -48
+    .cfi_adjust_cfa_offset 8
+
+    $mul01
+
+    xor %r15, %r15
+    add 0x18(%rdi), %r8
+    adc 0x20(%rdi), %r9
+    adc 0x28(%rdi), %r10
+    adc 0x30(%rdi), %r11
+    adc 0x38(%rdi), %r12
+    adc 0x40(%rdi), %r13
+    adc 0x48(%rdi), %r14
+    adc 0x50(%rdi), %r15
+    mov %r8, 0x18(%rdi)
+    mov %r9, 0x20(%rdi)
+    mov %r10, 0x28(%rdi)
+    mov %r11, 0x30(%rdi)
+    mov %r12, 0x38(%rdi)
+    mov %r13, 0x40(%rdi)
+    mov %r14, 0x48(%rdi)
+    mov %r15, 0x50(%rdi)
+    mov 0x58(%rdi), %r8
+    mov 0x60(%rdi), %r9
+    mov 0x68(%rdi), %r10
+    mov 0x70(%rdi), %r11
+    mov 0x78(%rdi), %r12
+    adc \$0x0, %r8
+    adc \$0x0, %r9
+    adc \$0x0, %r10
+    adc \$0x0, %r11
+    adc \$0x0, %r12
+    mov %r8, 0x58(%rdi)
+    mov %r9, 0x60(%rdi)
+    mov %r10, 0x68(%rdi)
+    mov %r11, 0x70(%rdi)
+    mov %r12, 0x78(%rdi)
+
+    $mul23
+
+    xor %r15, %r15
+    add 0x28(%rdi), %r8
+    adc 0x30(%rdi), %r9
+    adc 0x38(%rdi), %r10
+    adc 0x40(%rdi), %r11
+    adc 0x48(%rdi), %r12
+    adc 0x50(%rdi), %r13
+    adc 0x58(%rdi), %r14
+    adc 0x60(%rdi), %r15
+    mov %r8, 0x28(%rdi)
+    mov %r9, 0x30(%rdi)
+    mov %r10, 0x38(%rdi)
+    mov %r11, 0x40(%rdi)
+    mov %r12, 0x48(%rdi)
+    mov %r13, 0x50(%rdi)
+    mov %r14, 0x58(%rdi)
+    mov %r15, 0x60(%rdi)
+    mov 0x68(%rdi), %r8
+    mov 0x70(%rdi), %r9
+    mov 0x78(%rdi), %r10
+    adc \$0x0, %r8
+    adc \$0x0, %r9
+    adc \$0x0, %r10
+    mov %r8, 0x68(%rdi)
+    mov %r9, 0x70(%rdi)
+    mov %r10, 0x78(%rdi)
+
+    $mul45
+
+    xor %r15, %r15
+    xor %rbx, %rbx
+    add 0x38(%rdi), %r8
+    adc 0x40(%rdi), %r9
+    adc 0x48(%rdi), %r10
+    adc 0x50(%rdi), %r11
+    adc 0x58(%rdi), %r12
+    adc 0x60(%rdi), %r13
+    adc 0x68(%rdi), %r14
+    adc 0x70(%rdi), %r15
+    adc 0x78(%rdi), %rbx
+    mov %r8, 0x38(%rdi)
+    mov %r9, (%rsi)         # Final result c0
+    mov %r10, 0x48(%rdi)
+    mov %r11, 0x50(%rdi)
+    mov %r12, 0x58(%rdi)
+    mov %r13, 0x60(%rdi)
+    mov %r14, 0x68(%rdi)
+    mov %r15, 0x70(%rdi)
+    mov %rbx, 0x78(%rdi)
+
+    $mul67
+
+    add 0x48(%rdi), %r8
+    adc 0x50(%rdi), %r9
+    adc 0x58(%rdi), %r10
+    adc 0x60(%rdi), %r11
+    adc 0x68(%rdi), %r12
+    adc 0x70(%rdi), %r13
+    adc 0x78(%rdi), %r14
+    mov %r8, 0x8(%rsi)
+    mov %r9, 0x10(%rsi)
+    mov %r10, 0x18(%rsi)
+    mov %r11, 0x20(%rsi)
+    mov %r12, 0x28(%rsi)
+    mov %r13, 0x30(%rsi)
+    mov %r14, 0x38(%rsi)
+
+    pop %rbx
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value rbx
+    pop %r15
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r15
+    pop %r14
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r14
+    pop %r13
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r13
+    pop %r12
+  .cfi_adjust_cfa_offset -8
+  .cfi_same_value r12
+    ret
+  .cfi_endproc
+___
+  return $body;
+}
+
+#  Montgomery reduction
+#  Based on comba method
+#  Operation: c [rsi] = a [rdi]
+#  NOTE: a=c is not allowed
+sub rdc {
+  my $jump_optim=&alt_impl(".Lrdc_mulx_asm") if ($bmi2_adx);
+  my $body=&rdc_mulx() if ($bmi2_adx);
+
+  $body.=<<___;
+    .globl  ${PREFIX}_fprdc
+    .type   ${PREFIX}_fprdc,\@function,3
+    ${PREFIX}_fprdc:
+    .cfi_startproc
+      push %r12
+    .cfi_adjust_cfa_offset  8
+    .cfi_offset r12, -16
+      push %r13
+    .cfi_adjust_cfa_offset  8
+    .cfi_offset r13, -24
+      push %r14
+    .cfi_adjust_cfa_offset  8
+    .cfi_offset r14, -32
+      push %r15
+    .cfi_adjust_cfa_offset  8
+    .cfi_offset r15, -40
+      push %rbx
+    .cfi_adjust_cfa_offset  8
+    .cfi_offset rbx, -48
+
+    $jump_optim
+
+    # Reduction, generic x86 implementation
+    lea .Lp503p1(%rip), %rbx
+
+    mov (%rdi), %r11
+    mov (%rbx), %rax
+    mul %r11
+    xor %r8, %r8
+    add 0x18(%rdi), %rax
+    mov %rax, 0x18(%rsi)  # z3
+    adc %rdx, %r8
+
+    xor %r9, %r9
+    mov 0x8(%rbx), %rax
+    mul %r11
+    xor %r10, %r10
+    add %rax, %r8
+    adc %rdx, %r9
+
+    mov 0x8(%rdi), %r12
+    mov (%rbx), %rax
+    mul %r12
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+    add 0x20(%rdi), %r8
+    mov %r8, 0x20(%rsi)   # z4
+    adc  \$0, %r9
+    adc  \$0, %r10
+
+    xor %r8, %r8
+    mov 0x10(%rbx), %rax
+    mul %r11
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 8(%rbx), %rax
+    mul %r12
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x10(%rdi), %r13
+    mov (%rbx), %rax
+    mul %r13
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+    add 0x28(%rdi), %r9
+    mov %r9, 0x28(%rsi)   # z5
+    adc  \$0, %r10
+    adc  \$0, %r8
+
+    xor %r9, %r9
+    mov 0x18(%rbx), %rax
+    mul %r11
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x10(%rbx), %rax
+    mul %r12
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x8(%rbx), %rax
+    mul %r13
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x18(%rsi), %r14
+    mov (%rbx), %rax
+    mul %r14
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+    add 0x30(%rdi), %r10
+    mov %r10, 0x30(%rsi)    # z6
+    adc  \$0, %r8
+    adc  \$0, %r9
+
+    xor %r10, %r10
+    mov 0x20(%rbx), %rax
+    mul %r11
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x18(%rbx), %rax
+    mul %r12
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x10(%rbx), %rax
+    mul %r13
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x8(%rbx), %rax
+    mul %r14
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x20(%rsi), %r15
+    mov (%rbx), %rax
+    mul %r15
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+    add 0x38(%rdi), %r8     # Z7
+    mov %r8, 0x38(%rsi)
+    adc  \$0, %r9
+    adc  \$0, %r10
+
+    xor %r8, %r8
+    mov 0x20(%rbx), %rax
+    mul %r12
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x18(%rbx), %rax
+    mul %r13
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x10(%rbx), %rax
+    mul %r14
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x8(%rbx), %rax
+    mul %r15
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x28(%rsi), %rcx
+    mov (%rbx), %rax
+    mul %rcx
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+    add 0x40(%rdi), %r9
+    mov %r9, (%rsi)       # Z9
+    adc  \$0, %r10
+    adc  \$0, %r8
+
+    xor %r9, %r9
+    mov 0x20(%rbx), %rax
+    mul %r13
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x18(%rbx), %rax
+    mul %r14
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x10(%rbx), %rax
+    mul %r15
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 8(%rbx), %rax
+    mul %rcx
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x30(%rsi), %r13
+    mov (%rbx), %rax
+    mul %r13
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+    add 0x48(%rdi), %r10
+    mov %r10, 0x8(%rsi)     # Z1
+    adc  \$0, %r8
+    adc  \$0, %r9
+
+    xor %r10, %r10
+    mov 0x20(%rbx), %rax
+    mul %r14
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x18(%rbx), %rax
+    mul %r15
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x10(%rbx), %rax
+    mul %rcx
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 8(%rbx), %rax
+    mul %r13
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x38(%rsi), %r14
+    mov (%rbx), %rax
+    mul %r14
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+    add 0x50(%rdi), %r8
+    mov %r8, 0x10(%rsi)     # Z2
+    adc  \$0, %r9
+    adc  \$0, %r10
+
+    xor %r8, %r8
+    mov 0x20(%rbx), %rax
+    mul %r15
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x18(%rbx), %rax
+    mul %rcx
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 0x10(%rbx), %rax
+    mul %r13
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+
+    mov 8(%rbx), %rax
+    mul %r14
+    add %rax, %r9
+    adc %rdx, %r10
+    adc  \$0, %r8
+    add 0x58(%rdi), %r9
+    mov %r9, 0x18(%rsi)     # Z3
+    adc  \$0, %r10
+    adc  \$0, %r8
+
+    xor %r9, %r9
+    mov 0x20(%rbx), %rax
+    mul %rcx
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x18(%rbx), %rax
+    mul %r13
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+
+    mov 0x10(%rbx), %rax
+    mul %r14
+    add %rax, %r10
+    adc %rdx, %r8
+    adc  \$0, %r9
+    add 0x60(%rdi), %r10
+    mov %r10, 0x20(%rsi)    # Z4
+    adc  \$0, %r8
+    adc  \$0, %r9
+
+    xor %r10, %r10
+    mov 0x20(%rbx), %rax
+    mul %r13
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+
+    mov 0x18(%rbx), %rax
+    mul %r14
+    add %rax, %r8
+    adc %rdx, %r9
+    adc  \$0, %r10
+    add 0x68(%rdi), %r8     # Z5
+    mov %r8, 0x28(%rsi)     # Z5
+    adc  \$0, %r9
+    adc  \$0, %r10
+
+    mov 0x20(%rbx), %rax
+    mul %r14
+    add %rax, %r9
+    adc %rdx, %r10
+    add 0x70(%rdi), %r9     # Z6
+    mov %r9, 0x30(%rsi)     # Z6
+    adc  \$0, %r10
+    add 0x78(%rdi), %r10    # Z7
+    mov %r10, 0x38(%rsi)    # Z7
+
+    pop %rbx
+  .cfi_adjust_cfa_offset -8
+    pop %r15
+  .cfi_adjust_cfa_offset -8
+    pop %r14
+  .cfi_adjust_cfa_offset -8
+    pop %r13
+  .cfi_adjust_cfa_offset -8
+    pop %r12
+  .cfi_adjust_cfa_offset -8
+    ret
+  .cfi_endproc
+___
+  return $body;
+}
+
+$code.=&rdc();
+
+foreach (split("\n",$code)) {
+  s/\`([^\`]*)\`/eval($1)/ge;
+  print $_,"\n";
+}
+
+close STDOUT;
diff --git a/third_party/sike/asm/fp_generic.c b/third_party/sike/asm/fp_generic.c
new file mode 100644
index 0000000..cb786ff
--- /dev/null
+++ b/third_party/sike/asm/fp_generic.c
@@ -0,0 +1,181 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: portable modular arithmetic for P503
+*********************************************************************************************/
+
+#include <openssl/base.h>
+
+#if defined(OPENSSL_NO_ASM) || \
+    (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
+
+#include "../utils.h"
+#include "../fpx.h"
+
+// Global constants
+extern const struct params_t p503;
+
+static void digit_x_digit(const crypto_word_t a, const crypto_word_t b, crypto_word_t* c)
+{ // Digit multiplication, digit * digit -> 2-digit result
+    crypto_word_t al, ah, bl, bh, temp;
+    crypto_word_t albl, albh, ahbl, ahbh, res1, res2, res3, carry;
+    crypto_word_t mask_low = (crypto_word_t)(-1) >> (sizeof(crypto_word_t)*4);
+    crypto_word_t mask_high = (crypto_word_t)(-1) << (sizeof(crypto_word_t)*4);
+
+    al = a & mask_low;                              // Low part
+    ah = a >> (sizeof(crypto_word_t) * 4);          // High part
+    bl = b & mask_low;
+    bh = b >> (sizeof(crypto_word_t) * 4);
+
+    albl = al*bl;
+    albh = al*bh;
+    ahbl = ah*bl;
+    ahbh = ah*bh;
+    c[0] = albl & mask_low;                         // C00
+
+    res1 = albl >> (sizeof(crypto_word_t) * 4);
+    res2 = ahbl & mask_low;
+    res3 = albh & mask_low;
+    temp = res1 + res2 + res3;
+    carry = temp >> (sizeof(crypto_word_t) * 4);
+    c[0] ^= temp << (sizeof(crypto_word_t) * 4);    // C01
+
+    res1 = ahbl >> (sizeof(crypto_word_t) * 4);
+    res2 = albh >> (sizeof(crypto_word_t) * 4);
+    res3 = ahbh & mask_low;
+    temp = res1 + res2 + res3 + carry;
+    c[1] = temp & mask_low;                         // C10
+    carry = temp & mask_high;
+    c[1] ^= (ahbh & mask_high) + carry;             // C11
+}
+
+void sike_fpadd(const felm_t a, const felm_t b, felm_t c)
+{ // Modular addition, c = a+b mod p503.
+  // Inputs: a, b in [0, 2*p503-1]
+  // Output: c in [0, 2*p503-1]
+    unsigned int i, carry = 0;
+    crypto_word_t mask;
+
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        ADDC(carry, a[i], b[i], carry, c[i]);
+    }
+
+    carry = 0;
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        SUBC(carry, c[i], ((crypto_word_t*)p503.prime_x2)[i], carry, c[i]);
+    }
+    mask = 0 - (crypto_word_t)carry;
+
+    carry = 0;
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        ADDC(carry, c[i], ((crypto_word_t*)p503.prime_x2)[i] & mask, carry, c[i]);
+    }
+}
+
+void sike_fpsub(const felm_t a, const felm_t b, felm_t c)
+{ // Modular subtraction, c = a-b mod p503.
+  // Inputs: a, b in [0, 2*p503-1]
+  // Output: c in [0, 2*p503-1]
+    unsigned int i, borrow = 0;
+    crypto_word_t mask;
+
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        SUBC(borrow, a[i], b[i], borrow, c[i]);
+    }
+    mask = 0 - (crypto_word_t)borrow;
+
+    borrow = 0;
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        ADDC(borrow, c[i], ((crypto_word_t*)p503.prime_x2)[i] & mask, borrow, c[i]);
+    }
+}
+
+void sike_mpmul(const felm_t a, const felm_t b, dfelm_t c)
+{ // Multiprecision comba multiply, c = a*b, where lng(a) = lng(b) = NWORDS_FIELD.
+    unsigned int i, j;
+    crypto_word_t t = 0, u = 0, v = 0, UV[2];
+    unsigned int carry = 0;
+
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        for (j = 0; j <= i; j++) {
+            MUL(a[j], b[i-j], UV+1, UV[0]);
+            ADDC(0, UV[0], v, carry, v);
+            ADDC(carry, UV[1], u, carry, u);
+            t += carry;
+        }
+        c[i] = v;
+        v = u;
+        u = t;
+        t = 0;
+    }
+
+    for (i = NWORDS_FIELD; i < 2*NWORDS_FIELD-1; i++) {
+        for (j = i-NWORDS_FIELD+1; j < NWORDS_FIELD; j++) {
+            MUL(a[j], b[i-j], UV+1, UV[0]);
+            ADDC(0, UV[0], v, carry, v);
+            ADDC(carry, UV[1], u, carry, u);
+            t += carry;
+        }
+        c[i] = v;
+        v = u;
+        u = t;
+        t = 0;
+    }
+    c[2*NWORDS_FIELD-1] = v;
+}
+
+void sike_fprdc(const felm_t ma, felm_t mc)
+{ // Efficient Montgomery reduction using comba and exploiting the special form of the prime p503.
+  // mc = ma*R^-1 mod p503x2, where R = 2^512.
+  // If ma < 2^512*p503, the output mc is in the range [0, 2*p503-1].
+  // ma is assumed to be in Montgomery representation.
+    unsigned int i, j, carry, count = p503_ZERO_WORDS;
+    crypto_word_t UV[2], t = 0, u = 0, v = 0;
+
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        mc[i] = 0;
+    }
+
+    for (i = 0; i < NWORDS_FIELD; i++) {
+        for (j = 0; j < i; j++) {
+            if (j < (i-p503_ZERO_WORDS+1)) {
+                MUL(mc[j], ((crypto_word_t*)p503.prime_p1)[i-j], UV+1, UV[0]);
+                ADDC(0, UV[0], v, carry, v);
+                ADDC(carry, UV[1], u, carry, u);
+                t += carry;
+            }
+        }
+        ADDC(0, v, ma[i], carry, v);
+        ADDC(carry, u, 0, carry, u);
+        t += carry;
+        mc[i] = v;
+        v = u;
+        u = t;
+        t = 0;
+    }
+
+    for (i = NWORDS_FIELD; i < 2*NWORDS_FIELD-1; i++) {
+        if (count > 0) {
+            count -= 1;
+        }
+        for (j = i-NWORDS_FIELD+1; j < NWORDS_FIELD; j++) {
+            if (j < (NWORDS_FIELD-count)) {
+                MUL(mc[j], ((crypto_word_t*)p503.prime_p1)[i-j], UV+1, UV[0]);
+                ADDC(0, UV[0], v, carry, v);
+                ADDC(carry, UV[1], u, carry, u);
+                t += carry;
+            }
+        }
+        ADDC(0, v, ma[i], carry, v);
+        ADDC(carry, u, 0, carry, u);
+        t += carry;
+        mc[i-NWORDS_FIELD] = v;
+        v = u;
+        u = t;
+        t = 0;
+    }
+    ADDC(0, v, ma[2*NWORDS_FIELD-1], carry, v);
+    mc[NWORDS_FIELD-1] = v;
+}
+
+#endif  // NO_ASM || (!X86_64 && !AARCH64)
diff --git a/third_party/sike/fpx.c b/third_party/sike/fpx.c
new file mode 100644
index 0000000..abb1ca2
--- /dev/null
+++ b/third_party/sike/fpx.c
@@ -0,0 +1,305 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: core functions over GF(p) and GF(p^2)
+*********************************************************************************************/
+#include <openssl/base.h>
+
+#include "utils.h"
+#include "fpx.h"
+
+extern const struct params_t p503;
+
+// Multiprecision squaring, c = a^2 mod p.
+static void fpsqr_mont(const felm_t ma, felm_t mc)
+{
+    dfelm_t temp = {0};
+    sike_mpmul(ma, ma, temp);
+    sike_fprdc(temp, mc);
+}
+
+// Chain to compute a^(p-3)/4 using Montgomery arithmetic.
+static void fpinv_chain_mont(felm_t a)
+{
+    unsigned int i, j;
+    felm_t t[15], tt;
+
+    // Precomputed table
+    fpsqr_mont(a, tt);
+    sike_fpmul_mont(a, tt, t[0]);
+    for (i = 0; i <= 13; i++) sike_fpmul_mont(t[i], tt, t[i+1]);
+
+    sike_fpcopy(a, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(a, tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[8], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[9], tt, tt);
+    for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[0], tt, tt);
+    for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(a, tt, tt);
+    for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[2], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[8], tt, tt);
+    for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(a, tt, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[10], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[0], tt, tt);
+    for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[10], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[10], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[5], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[2], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[3], tt, tt);
+    for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[5], tt, tt);
+    for (i = 0; i < 12; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[12], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[8], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[12], tt, tt);
+    for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[11], tt, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[5], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[14], tt, tt);
+    for (i = 0; i < 7; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[14], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[5], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[8], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(a, tt, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[4], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[6], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[5], tt, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[7], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(a, tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[0], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[11], tt, tt);
+    for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[13], tt, tt);
+    for (i = 0; i < 8; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[1], tt, tt);
+    for (i = 0; i < 6; i++) fpsqr_mont(tt, tt);
+    sike_fpmul_mont(t[10], tt, tt);
+    for (j = 0; j < 49; j++) {
+        for (i = 0; i < 5; i++) fpsqr_mont(tt, tt);
+        sike_fpmul_mont(t[14], tt, tt);
+    }
+    sike_fpcopy(tt, a);
+}
+
+// Field inversion using Montgomery arithmetic, a = a^(-1)*R mod p.
+static void fpinv_mont(felm_t a)
+{
+    felm_t tt = {0};
+    sike_fpcopy(a, tt);
+    fpinv_chain_mont(tt);
+    fpsqr_mont(tt, tt);
+    fpsqr_mont(tt, tt);
+    sike_fpmul_mont(a, tt, a);
+}
+
+// Multiprecision addition, c = a+b, where lng(a) = lng(b) = nwords. Returns the carry bit.
+#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
+inline static unsigned int mp_add(const felm_t a, const felm_t b, felm_t c, const unsigned int nwords) {
+    uint8_t carry = 0;
+    for (size_t i = 0; i < nwords; i++) {
+        ADDC(carry, a[i], b[i], carry, c[i]);
+    }
+    return carry;
+}
+
+// Multiprecision subtraction, c = a-b, where lng(a) = lng(b) = nwords. Returns the borrow bit.
+inline static unsigned int mp_sub(const felm_t a, const felm_t b, felm_t c, const unsigned int nwords) {
+    uint32_t borrow = 0;
+    for (size_t i = 0; i < nwords; i++) {
+        SUBC(borrow, a[i], b[i], borrow, c[i]);
+    }
+    return borrow;
+}
+#endif
+
+// Multiprecision addition, c = a+b.
+inline static void mp_addfast(const felm_t a, const felm_t b, felm_t c)
+{
+#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
+    mp_add(a, b, c, NWORDS_FIELD);
+#else
+    sike_mpadd_asm(a, b, c);
+#endif
+}
+
+// Multiprecision subtraction, c = a-b, where lng(a) = lng(b) = 2*NWORDS_FIELD.
+// If c < 0 then returns mask = 0xFF..F, else mask = 0x00..0
+inline static crypto_word_t mp_subfast(const dfelm_t a, const dfelm_t b, dfelm_t c) {
+#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
+    return (0 - (crypto_word_t)mp_sub(a, b, c, 2*NWORDS_FIELD));
+#else
+    return sike_mpsubx2_asm(a, b, c);
+#endif
+}
+
+// Multiprecision subtraction, c = c-a-b, where lng(a) = lng(b) = 2*NWORDS_FIELD.
+// Inputs should be s.t. c > a and c > b
+inline static void mp_dblsubfast(const dfelm_t a, const dfelm_t b, dfelm_t c) {
+#if defined(OPENSSL_NO_ASM) || (!defined(OPENSSL_X86_64) && !defined(OPENSSL_AARCH64))
+    mp_sub(c, a, c, 2*NWORDS_FIELD);
+    mp_sub(c, b, c, 2*NWORDS_FIELD);
+#else
+    sike_mpdblsubx2_asm(a, b, c);
+#endif
+}
+
+// Copy a field element, c = a.
+void sike_fpcopy(const felm_t a, felm_t c) {
+    for (size_t i = 0; i < NWORDS_FIELD; i++) {
+        c[i] = a[i];
+    }
+}
+
+// Field multiplication using Montgomery arithmetic, c = a*b*R^-1 mod p503, where R=2^768
+void sike_fpmul_mont(const felm_t ma, const felm_t mb, felm_t mc)
+{
+    dfelm_t temp = {0};
+    sike_mpmul(ma, mb, temp);
+    sike_fprdc(temp, mc);
+}
+
+// Conversion from Montgomery representation to standard representation,
+// c = ma*R^(-1) mod p = a mod p, where ma in [0, p-1].
+void sike_from_mont(const felm_t ma, felm_t c)
+{
+    felm_t one = {0};
+    one[0] = 1;
+
+    sike_fpmul_mont(ma, one, c);
+    sike_fpcorrection(c);
+}
+
+// GF(p^2) squaring using Montgomery arithmetic, c = a^2 in GF(p^2).
+// Inputs: a = a0+a1*i, where a0, a1 are in [0, 2*p-1]
+// Output: c = c0+c1*i, where c0, c1 are in [0, 2*p-1]
+void sike_fp2sqr_mont(const f2elm_t a, f2elm_t c) {
+    felm_t t1, t2, t3;
+
+    mp_addfast(a->c0, a->c1, t1);                      // t1 = a0+a1
+    sike_fpsub(a->c0, a->c1, t2);                      // t2 = a0-a1
+    mp_addfast(a->c0, a->c0, t3);                      // t3 = 2a0
+    sike_fpmul_mont(t1, t2, c->c0);                    // c0 = (a0+a1)(a0-a1)
+    sike_fpmul_mont(t3, a->c1, c->c1);                 // c1 = 2a0*a1
+}
+
+// Modular negation, a = -a mod p503.
+// Input/output: a in [0, 2*p503-1]
+void sike_fpneg(felm_t a) {
+  uint32_t borrow = 0;
+  for (size_t i = 0; i < NWORDS_FIELD; i++) {
+    SUBC(borrow, ((crypto_word_t*)p503.prime_x2)[i], a[i], borrow, a[i]);
+  }
+}
+
+// Modular division by two, c = a/2 mod p503.
+// Input : a in [0, 2*p503-1]
+// Output: c in [0, 2*p503-1]
+void sike_fpdiv2(const felm_t a, felm_t c) {
+  uint32_t carry = 0;
+  crypto_word_t mask;
+
+  mask = 0 - (crypto_word_t)(a[0] & 1);    // If a is odd compute a+p503
+  for (size_t i = 0; i < NWORDS_FIELD; i++) {
+    ADDC(carry, a[i], ((crypto_word_t*)p503.prime)[i] & mask, carry, c[i]);
+  }
+
+  // Multiprecision right shift by one.
+  for (size_t i = 0; i < NWORDS_FIELD-1; i++) {
+    c[i] = (c[i] >> 1) ^ (c[i+1] << (RADIX - 1));
+  }
+  c[NWORDS_FIELD-1] >>= 1;
+}
+
+// Modular correction to reduce field element a in [0, 2*p503-1] to [0, p503-1].
+void sike_fpcorrection(felm_t a) {
+  uint32_t borrow = 0;
+  crypto_word_t mask;
+
+  for (size_t i = 0; i < NWORDS_FIELD; i++) {
+    SUBC(borrow, a[i], ((crypto_word_t*)p503.prime)[i], borrow, a[i]);
+  }
+  mask = 0 - (crypto_word_t)borrow;
+
+  borrow = 0;
+  for (size_t i = 0; i < NWORDS_FIELD; i++) {
+    ADDC(borrow, a[i], ((crypto_word_t*)p503.prime)[i] & mask, borrow, a[i]);
+  }
+}
+
+// GF(p^2) multiplication using Montgomery arithmetic, c = a*b in GF(p^2).
+// Inputs: a = a0+a1*i and b = b0+b1*i, where a0, a1, b0, b1 are in [0, 2*p-1]
+// Output: c = c0+c1*i, where c0, c1 are in [0, 2*p-1]
+void sike_fp2mul_mont(const f2elm_t a, const f2elm_t b, f2elm_t c) {
+    felm_t t1, t2;
+    dfelm_t tt1, tt2, tt3;
+    crypto_word_t mask;
+
+    mp_addfast(a->c0, a->c1, t1);                      // t1 = a0+a1
+    mp_addfast(b->c0, b->c1, t2);                      // t2 = b0+b1
+    sike_mpmul(a->c0, b->c0, tt1);                     // tt1 = a0*b0
+    sike_mpmul(a->c1, b->c1, tt2);                     // tt2 = a1*b1
+    sike_mpmul(t1, t2, tt3);                           // tt3 = (a0+a1)*(b0+b1)
+    mp_dblsubfast(tt1, tt2, tt3);                      // tt3 = (a0+a1)*(b0+b1) - a0*b0 - a1*b1
+    mask = mp_subfast(tt1, tt2, tt1);                  // tt1 = a0*b0 - a1*b1. If tt1 < 0 then mask = 0xFF..F, else if tt1 >= 0 then mask = 0x00..0
+
+    for (size_t i = 0; i < NWORDS_FIELD; i++) {
+        t1[i] = ((crypto_word_t*)p503.prime)[i] & mask;
+    }
+
+    sike_fprdc(tt3, c->c1);                             // c[1] = (a0+a1)*(b0+b1) - a0*b0 - a1*b1
+    mp_addfast(&((crypto_word_t*)&tt1)[NWORDS_FIELD], t1, &((crypto_word_t*)&tt1)[NWORDS_FIELD]);
+    sike_fprdc(tt1, c->c0);                             // c[0] = a0*b0 - a1*b1
+}
+
+// GF(p^2) inversion using Montgomery arithmetic, a = (a0-i*a1)/(a0^2+a1^2).
+void sike_fp2inv_mont(f2elm_t a) {
+    f2elm_t t1;
+
+    fpsqr_mont(a->c0, t1->c0);                         // t10 = a0^2
+    fpsqr_mont(a->c1, t1->c1);                         // t11 = a1^2
+    sike_fpadd(t1->c0, t1->c1, t1->c0);                // t10 = a0^2+a1^2
+    fpinv_mont(t1->c0);                                // t10 = (a0^2+a1^2)^-1
+    sike_fpneg(a->c1);                                 // a = a0-i*a1
+    sike_fpmul_mont(a->c0, t1->c0, a->c0);
+    sike_fpmul_mont(a->c1, t1->c0, a->c1);             // a = (a0-i*a1)*(a0^2+a1^2)^-1
+}
diff --git a/third_party/sike/fpx.h b/third_party/sike/fpx.h
new file mode 100644
index 0000000..ed67768
--- /dev/null
+++ b/third_party/sike/fpx.h
@@ -0,0 +1,112 @@
+#ifndef FPX_H_
+#define FPX_H_
+
+#include "utils.h"
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+// Modular addition, c = a+b mod p503.
+void sike_fpadd(const felm_t a, const felm_t b, felm_t c);
+// Modular subtraction, c = a-b mod p503.
+void sike_fpsub(const felm_t a, const felm_t b, felm_t c);
+// Modular division by two, c = a/2 mod p503.
+void sike_fpdiv2(const felm_t a, felm_t c);
+// Modular correction to reduce field element a in [0, 2*p503-1] to [0, p503-1].
+void sike_fpcorrection(felm_t a);
+// Multiprecision multiply, c = a*b, where lng(a) = lng(b) = nwords.
+void sike_mpmul(const felm_t a, const felm_t b, dfelm_t c);
+// 503-bit Montgomery reduction, c = a mod p
+void sike_fprdc(const dfelm_t a, felm_t c);
+// Double 2x503-bit multiprecision subtraction, c = c-a-b
+void sike_mpdblsubx2_asm(const felm_t a, const felm_t b, felm_t c);
+// Multiprecision subtraction, c = a-b
+crypto_word_t sike_mpsubx2_asm(const dfelm_t a, const dfelm_t b, dfelm_t c);
+// 503-bit multiprecision addition, c = a+b
+void sike_mpadd_asm(const felm_t a, const felm_t b, felm_t c);
+// Modular negation, a = -a mod p503.
+void sike_fpneg(felm_t a);
+// Copy of a field element, c = a
+void sike_fpcopy(const felm_t a, felm_t c);
+// Copy a field element, c = a.
+void sike_fpzero(felm_t a);
+// If option = 0xFF...FF x=y; y=x, otherwise swap doesn't happen. Constant time.
+void sike_cswap_asm(point_proj_t x, point_proj_t y, const crypto_word_t option);
+// Conversion from Montgomery representation to standard representation,
+// c = ma*R^(-1) mod p = a mod p, where ma in [0, p-1].
+void sike_from_mont(const felm_t ma, felm_t c);
+// Field multiplication using Montgomery arithmetic, c = a*b*R^-1 mod p503, where R=2^768
+void sike_fpmul_mont(const felm_t ma, const felm_t mb, felm_t mc);
+// GF(p503^2) multiplication using Montgomery arithmetic, c = a*b in GF(p503^2)
+void sike_fp2mul_mont(const f2elm_t a, const f2elm_t b, f2elm_t c);
+// GF(p503^2) inversion using Montgomery arithmetic, a = (a0-i*a1)/(a0^2+a1^2)
+void sike_fp2inv_mont(f2elm_t a);
+// GF(p^2) squaring using Montgomery arithmetic, c = a^2 in GF(p^2).
+void sike_fp2sqr_mont(const f2elm_t a, f2elm_t c);
+// Modular correction, a = a in GF(p^2).
+void sike_fp2correction(f2elm_t a);
+
+#if defined(__cplusplus)
+}  // extern C
+#endif
+
+// GF(p^2) addition, c = a+b in GF(p^2).
+#define sike_fp2add(a, b, c)             \
+do {                                     \
+    sike_fpadd(a->c0, b->c0, c->c0);     \
+    sike_fpadd(a->c1, b->c1, c->c1);     \
+} while(0)
+
+// GF(p^2) subtraction, c = a-b in GF(p^2).
+#define sike_fp2sub(a,b,c)               \
+do {                                     \
+    sike_fpsub(a->c0, b->c0, c->c0);     \
+    sike_fpsub(a->c1, b->c1, c->c1);     \
+} while(0)
+
+// Copy a GF(p^2) element, c = a.
+#define sike_fp2copy(a, c)               \
+do {                                     \
+    sike_fpcopy(a->c0, c->c0);           \
+    sike_fpcopy(a->c1, c->c1);           \
+} while(0)
+
+// GF(p^2) negation, a = -a in GF(p^2).
+#define sike_fp2neg(a)                   \
+do {                                     \
+    sike_fpneg(a->c0);                   \
+    sike_fpneg(a->c1);                   \
+} while(0)
+
+// GF(p^2) division by two, c = a/2  in GF(p^2).
+#define sike_fp2div2(a, c)               \
+do {                                     \
+    sike_fpdiv2(a->c0, c->c0);           \
+    sike_fpdiv2(a->c1, c->c1);           \
+} while(0)
+
+// Modular correction, a = a in GF(p^2).
+#define sike_fp2correction(a)            \
+do {                                     \
+    sike_fpcorrection(a->c0);            \
+    sike_fpcorrection(a->c1);            \
+} while(0)
+
+// Conversion of a GF(p^2) element to Montgomery representation,
+// mc_i = a_i*R^2*R^(-1) = a_i*R in GF(p^2).
+#define sike_to_fp2mont(a, mc)           \
+do {                                     \
+    sike_fpmul_mont(a->c0, (crypto_word_t*)&p503.mont_R2, mc->c0);   \
+    sike_fpmul_mont(a->c1, (crypto_word_t*)&p503.mont_R2, mc->c1);   \
+} while(0)
+
+// Conversion of a GF(p^2) element from Montgomery representation to standard representation,
+// c_i = ma_i*R^(-1) = a_i in GF(p^2).
+#define sike_from_fp2mont(ma, c)         \
+do {                                     \
+    sike_from_mont(ma->c0, c->c0);       \
+    sike_from_mont(ma->c1, c->c1);       \
+} while(0)
+
+#endif // FPX_H_
diff --git a/third_party/sike/isogeny.c b/third_party/sike/isogeny.c
new file mode 100644
index 0000000..45f9c40
--- /dev/null
+++ b/third_party/sike/isogeny.c
@@ -0,0 +1,260 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: elliptic curve and isogeny functions
+*********************************************************************************************/
+#include "utils.h"
+#include "isogeny.h"
+#include "fpx.h"
+
+static void xDBL(const point_proj_t P, point_proj_t Q, const f2elm_t A24plus, const f2elm_t C24)
+{ // Doubling of a Montgomery point in projective coordinates (X:Z).
+  // Input: projective Montgomery x-coordinates P = (X1:Z1), where x1=X1/Z1 and Montgomery curve constants A+2C and 4C.
+  // Output: projective Montgomery x-coordinates Q = 2*P = (X2:Z2).
+    f2elm_t t0, t1;
+
+    sike_fp2sub(P->X, P->Z, t0);                         // t0 = X1-Z1
+    sike_fp2add(P->X, P->Z, t1);                         // t1 = X1+Z1
+    sike_fp2sqr_mont(t0, t0);                            // t0 = (X1-Z1)^2
+    sike_fp2sqr_mont(t1, t1);                            // t1 = (X1+Z1)^2
+    sike_fp2mul_mont(C24, t0, Q->Z);                     // Z2 = C24*(X1-Z1)^2
+    sike_fp2mul_mont(t1, Q->Z, Q->X);                    // X2 = C24*(X1-Z1)^2*(X1+Z1)^2
+    sike_fp2sub(t1, t0, t1);                             // t1 = (X1+Z1)^2-(X1-Z1)^2
+    sike_fp2mul_mont(A24plus, t1, t0);                   // t0 = A24plus*[(X1+Z1)^2-(X1-Z1)^2]
+    sike_fp2add(Q->Z, t0, Q->Z);                         // Z2 = A24plus*[(X1+Z1)^2-(X1-Z1)^2] + C24*(X1-Z1)^2
+    sike_fp2mul_mont(Q->Z, t1, Q->Z);                    // Z2 = [A24plus*[(X1+Z1)^2-(X1-Z1)^2] + C24*(X1-Z1)^2]*[(X1+Z1)^2-(X1-Z1)^2]
+}
+
+void xDBLe(const point_proj_t P, point_proj_t Q, const f2elm_t A24plus, const f2elm_t C24, size_t e)
+{ // Computes [2^e](X:Z) on Montgomery curve with projective constant via e repeated doublings.
+  // Input: projective Montgomery x-coordinates P = (XP:ZP), such that xP=XP/ZP and Montgomery curve constants A+2C and 4C.
+  // Output: projective Montgomery x-coordinates Q <- (2^e)*P.
+
+    memmove(Q, P, sizeof(*P));
+    for (size_t i = 0; i < e; i++) {
+        xDBL(Q, Q, A24plus, C24);
+    }
+}
+
+void get_4_isog(const point_proj_t P, f2elm_t A24plus, f2elm_t C24, f2elm_t* coeff)
+{ // Computes the corresponding 4-isogeny of a projective Montgomery point (X4:Z4) of order 4.
+  // Input:  projective point of order four P = (X4:Z4).
+  // Output: the 4-isogenous Montgomery curve with projective coefficients A+2C/4C and the 3 coefficients
+  //         that are used to evaluate the isogeny at a point in eval_4_isog().
+
+    sike_fp2sub(P->X, P->Z, coeff[1]);                   // coeff[1] = X4-Z4
+    sike_fp2add(P->X, P->Z, coeff[2]);                   // coeff[2] = X4+Z4
+    sike_fp2sqr_mont(P->Z, coeff[0]);                    // coeff[0] = Z4^2
+    sike_fp2add(coeff[0], coeff[0], coeff[0]);           // coeff[0] = 2*Z4^2
+    sike_fp2sqr_mont(coeff[0], C24);                     // C24 = 4*Z4^4
+    sike_fp2add(coeff[0], coeff[0], coeff[0]);           // coeff[0] = 4*Z4^2
+    sike_fp2sqr_mont(P->X, A24plus);                     // A24plus = X4^2
+    sike_fp2add(A24plus, A24plus, A24plus);              // A24plus = 2*X4^2
+    sike_fp2sqr_mont(A24plus, A24plus);                  // A24plus = 4*X4^4
+}
+
+void eval_4_isog(point_proj_t P, f2elm_t* coeff)
+{ // Evaluates the isogeny at the point (X:Z) in the domain of the isogeny, given a 4-isogeny phi defined
+  // by the 3 coefficients in coeff (computed in the function get_4_isog()).
+  // Inputs: the coefficients defining the isogeny, and the projective point P = (X:Z).
+  // Output: the projective point P = phi(P) = (X:Z) in the codomain.
+    f2elm_t t0, t1;
+
+    sike_fp2add(P->X, P->Z, t0);                         // t0 = X+Z
+    sike_fp2sub(P->X, P->Z, t1);                         // t1 = X-Z
+    sike_fp2mul_mont(t0, coeff[1], P->X);                // X = (X+Z)*coeff[1]
+    sike_fp2mul_mont(t1, coeff[2], P->Z);                // Z = (X-Z)*coeff[2]
+    sike_fp2mul_mont(t0, t1, t0);                        // t0 = (X+Z)*(X-Z)
+    sike_fp2mul_mont(t0, coeff[0], t0);                  // t0 = coeff[0]*(X+Z)*(X-Z)
+    sike_fp2add(P->X, P->Z, t1);                         // t1 = (X-Z)*coeff[2] + (X+Z)*coeff[1]
+    sike_fp2sub(P->X, P->Z, P->Z);                       // Z = (X-Z)*coeff[2] - (X+Z)*coeff[1]
+    sike_fp2sqr_mont(t1, t1);                            // t1 = [(X-Z)*coeff[2] + (X+Z)*coeff[1]]^2
+    sike_fp2sqr_mont(P->Z, P->Z);                        // Z = [(X-Z)*coeff[2] - (X+Z)*coeff[1]]^2
+    sike_fp2add(t1, t0, P->X);                           // X = coeff[0]*(X+Z)*(X-Z) + [(X-Z)*coeff[2] + (X+Z)*coeff[1]]^2
+    sike_fp2sub(P->Z, t0, t0);                           // t0 = [(X-Z)*coeff[2] - (X+Z)*coeff[1]]^2 - coeff[0]*(X+Z)*(X-Z)
+    sike_fp2mul_mont(P->X, t1, P->X);                    // Xfinal
+    sike_fp2mul_mont(P->Z, t0, P->Z);                    // Zfinal
+}
+
+
+void xTPL(const point_proj_t P, point_proj_t Q, const f2elm_t A24minus, const f2elm_t A24plus)
+{ // Tripling of a Montgomery point in projective coordinates (X:Z).
+  // Input: projective Montgomery x-coordinates P = (X:Z), where x=X/Z and Montgomery curve constants A24plus = A+2C and A24minus = A-2C.
+  // Output: projective Montgomery x-coordinates Q = 3*P = (X3:Z3).
+    f2elm_t t0, t1, t2, t3, t4, t5, t6;
+
+    sike_fp2sub(P->X, P->Z, t0);                         // t0 = X-Z
+    sike_fp2sqr_mont(t0, t2);                            // t2 = (X-Z)^2
+    sike_fp2add(P->X, P->Z, t1);                         // t1 = X+Z
+    sike_fp2sqr_mont(t1, t3);                            // t3 = (X+Z)^2
+    sike_fp2add(t0, t1, t4);                             // t4 = 2*X
+    sike_fp2sub(t1, t0, t0);                             // t0 = 2*Z
+    sike_fp2sqr_mont(t4, t1);                            // t1 = 4*X^2
+    sike_fp2sub(t1, t3, t1);                             // t1 = 4*X^2 - (X+Z)^2
+    sike_fp2sub(t1, t2, t1);                             // t1 = 4*X^2 - (X+Z)^2 - (X-Z)^2
+    sike_fp2mul_mont(t3, A24plus, t5);                   // t5 = A24plus*(X+Z)^2
+    sike_fp2mul_mont(t3, t5, t3);                        // t3 = A24plus*(X+Z)^3
+    sike_fp2mul_mont(A24minus, t2, t6);                  // t6 = A24minus*(X-Z)^2
+    sike_fp2mul_mont(t2, t6, t2);                        // t2 = A24minus*(X-Z)^3
+    sike_fp2sub(t2, t3, t3);                             // t3 = A24minus*(X-Z)^3 - coeff*(X+Z)^3
+    sike_fp2sub(t5, t6, t2);                             // t2 = A24plus*(X+Z)^2 - A24minus*(X-Z)^2
+    sike_fp2mul_mont(t1, t2, t1);                        // t1 = [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2]
+    sike_fp2add(t3, t1, t2);                             // t2 = [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2] + A24minus*(X-Z)^3 - coeff*(X+Z)^3
+    sike_fp2sqr_mont(t2, t2);                            // t2 = t2^2
+    sike_fp2mul_mont(t4, t2, Q->X);                      // X3 = 2*X*t2
+    sike_fp2sub(t3, t1, t1);                             // t1 = A24minus*(X-Z)^3 - A24plus*(X+Z)^3 - [4*X^2 - (X+Z)^2 - (X-Z)^2]*[A24plus*(X+Z)^2 - A24minus*(X-Z)^2]
+    sike_fp2sqr_mont(t1, t1);                            // t1 = t1^2
+    sike_fp2mul_mont(t0, t1, Q->Z);                      // Z3 = 2*Z*t1
+}
+
+void xTPLe(const point_proj_t P, point_proj_t Q, const f2elm_t A24minus, const f2elm_t A24plus, size_t e)
+{ // Computes [3^e](X:Z) on Montgomery curve with projective constant via e repeated triplings.
+  // Input: projective Montgomery x-coordinates P = (XP:ZP), such that xP=XP/ZP and Montgomery curve constants A24plus = A+2C and A24minus = A-2C.
+  // Output: projective Montgomery x-coordinates Q <- (3^e)*P.
+    memmove(Q, P, sizeof(*P));
+    for (size_t i = 0; i < e; i++) {
+        xTPL(Q, Q, A24minus, A24plus);
+    }
+}
+
+void get_3_isog(const point_proj_t P, f2elm_t A24minus, f2elm_t A24plus, f2elm_t* coeff)
+{ // Computes the corresponding 3-isogeny of a projective Montgomery point (X3:Z3) of order 3.
+  // Input:  projective point of order three P = (X3:Z3).
+  // Output: the 3-isogenous Montgomery curve with projective coefficient A/C.
+    f2elm_t t0, t1, t2, t3, t4;
+
+    sike_fp2sub(P->X, P->Z, coeff[0]);                   // coeff0 = X-Z
+    sike_fp2sqr_mont(coeff[0], t0);                      // t0 = (X-Z)^2
+    sike_fp2add(P->X, P->Z, coeff[1]);                   // coeff1 = X+Z
+    sike_fp2sqr_mont(coeff[1], t1);                      // t1 = (X+Z)^2
+    sike_fp2add(t0, t1, t2);                             // t2 = (X+Z)^2 + (X-Z)^2
+    sike_fp2add(coeff[0], coeff[1], t3);                 // t3 = 2*X
+    sike_fp2sqr_mont(t3, t3);                            // t3 = 4*X^2
+    sike_fp2sub(t3, t2, t3);                             // t3 = 4*X^2 - (X+Z)^2 - (X-Z)^2
+    sike_fp2add(t1, t3, t2);                             // t2 = 4*X^2 - (X-Z)^2
+    sike_fp2add(t3, t0, t3);                             // t3 = 4*X^2 - (X+Z)^2
+    sike_fp2add(t0, t3, t4);                             // t4 = 4*X^2 - (X+Z)^2 + (X-Z)^2
+    sike_fp2add(t4, t4, t4);                             // t4 = 2(4*X^2 - (X+Z)^2 + (X-Z)^2)
+    sike_fp2add(t1, t4, t4);                             // t4 = 8*X^2 - (X+Z)^2 + 2*(X-Z)^2
+    sike_fp2mul_mont(t2, t4, A24minus);                  // A24minus = [4*X^2 - (X-Z)^2]*[8*X^2 - (X+Z)^2 + 2*(X-Z)^2]
+    sike_fp2add(t1, t2, t4);                             // t4 = 4*X^2 + (X+Z)^2 - (X-Z)^2
+    sike_fp2add(t4, t4, t4);                             // t4 = 2(4*X^2 + (X+Z)^2 - (X-Z)^2)
+    sike_fp2add(t0, t4, t4);                             // t4 = 8*X^2 + 2*(X+Z)^2 - (X-Z)^2
+    sike_fp2mul_mont(t3, t4, t4);                        // t4 = [4*X^2 - (X+Z)^2]*[8*X^2 + 2*(X+Z)^2 - (X-Z)^2]
+    sike_fp2sub(t4, A24minus, t0);                       // t0 = [4*X^2 - (X+Z)^2]*[8*X^2 + 2*(X+Z)^2 - (X-Z)^2] - [4*X^2 - (X-Z)^2]*[8*X^2 - (X+Z)^2 + 2*(X-Z)^2]
+    sike_fp2add(A24minus, t0, A24plus);                  // A24plus = 8*X^2 - (X+Z)^2 + 2*(X-Z)^2
+}
+
+
+void eval_3_isog(point_proj_t Q, f2elm_t* coeff)
+{ // Computes the 3-isogeny R=phi(X:Z), given projective point (X3:Z3) of order 3 on a Montgomery curve and
+  // a point P with 2 coefficients in coeff (computed in the function get_3_isog()).
+  // Inputs: projective points P = (X3:Z3) and Q = (X:Z).
+  // Output: the projective point Q <- phi(Q) = (X3:Z3).
+    f2elm_t t0, t1, t2;
+
+    sike_fp2add(Q->X, Q->Z, t0);                       // t0 = X+Z
+    sike_fp2sub(Q->X, Q->Z, t1);                       // t1 = X-Z
+    sike_fp2mul_mont(t0, coeff[0], t0);                // t0 = coeff0*(X+Z)
+    sike_fp2mul_mont(t1, coeff[1], t1);                // t1 = coeff1*(X-Z)
+    sike_fp2add(t0, t1, t2);                           // t2 = coeff0*(X+Z) + coeff1*(X-Z)
+    sike_fp2sub(t1, t0, t0);                           // t0 = coeff1*(X-Z) - coeff0*(X+Z)
+    sike_fp2sqr_mont(t2, t2);                          // t2 = [coeff0*(X+Z) + coeff1*(X-Z)]^2
+    sike_fp2sqr_mont(t0, t0);                          // t0 = [coeff1*(X-Z) - coeff0*(X+Z)]^2
+    sike_fp2mul_mont(Q->X, t2, Q->X);                  // X3final = X*[coeff0*(X+Z) + coeff1*(X-Z)]^2
+    sike_fp2mul_mont(Q->Z, t0, Q->Z);                  // Z3final = Z*[coeff1*(X-Z) - coeff0*(X+Z)]^2
+}
+
+
+void inv_3_way(f2elm_t z1, f2elm_t z2, f2elm_t z3)
+{ // 3-way simultaneous inversion
+  // Input:  z1,z2,z3
+  // Output: 1/z1,1/z2,1/z3 (override inputs).
+    f2elm_t t0, t1, t2, t3;
+
+    sike_fp2mul_mont(z1, z2, t0);                      // t0 = z1*z2
+    sike_fp2mul_mont(z3, t0, t1);                      // t1 = z1*z2*z3
+    sike_fp2inv_mont(t1);                              // t1 = 1/(z1*z2*z3)
+    sike_fp2mul_mont(z3, t1, t2);                      // t2 = 1/(z1*z2)
+    sike_fp2mul_mont(t2, z2, t3);                      // t3 = 1/z1
+    sike_fp2mul_mont(t2, z1, z2);                      // z2 = 1/z2
+    sike_fp2mul_mont(t0, t1, z3);                      // z3 = 1/z3
+    sike_fp2copy(t3, z1);                              // z1 = 1/z1
+}
+
+
+void get_A(const f2elm_t xP, const f2elm_t xQ, const f2elm_t xR, f2elm_t A)
+{ // Given the x-coordinates of P, Q, and R, returns the value A corresponding to the Montgomery curve E_A: y^2=x^3+A*x^2+x such that R=Q-P on E_A.
+  // Input:  the x-coordinates xP, xQ, and xR of the points P, Q and R.
+  // Output: the coefficient A corresponding to the curve E_A: y^2=x^3+A*x^2+x.
+    f2elm_t t0, t1, one = F2ELM_INIT;
+
+    extern const struct params_t p503;
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, one->c0);
+    sike_fp2add(xP, xQ, t1);                           // t1 = xP+xQ
+    sike_fp2mul_mont(xP, xQ, t0);                      // t0 = xP*xQ
+    sike_fp2mul_mont(xR, t1, A);                       // A = xR*t1
+    sike_fp2add(t0, A, A);                             // A = A+t0
+    sike_fp2mul_mont(t0, xR, t0);                      // t0 = t0*xR
+    sike_fp2sub(A, one, A);                            // A = A-1
+    sike_fp2add(t0, t0, t0);                           // t0 = t0+t0
+    sike_fp2add(t1, xR, t1);                           // t1 = t1+xR
+    sike_fp2add(t0, t0, t0);                           // t0 = t0+t0
+    sike_fp2sqr_mont(A, A);                            // A = A^2
+    sike_fp2inv_mont(t0);                              // t0 = 1/t0
+    sike_fp2mul_mont(A, t0, A);                        // A = A*t0
+    sike_fp2sub(A, t1, A);                             // Afinal = A-t1
+}
+
+
+void j_inv(const f2elm_t A, const f2elm_t C, f2elm_t jinv)
+{ // Computes the j-invariant of a Montgomery curve with projective constant.
+  // Input: A,C in GF(p^2).
+  // Output: j=256*(A^2-3*C^2)^3/(C^4*(A^2-4*C^2)), which is the j-invariant of the Montgomery curve B*y^2=x^3+(A/C)*x^2+x or (equivalently) j-invariant of B'*y^2=C*x^3+A*x^2+C*x.
+    f2elm_t t0, t1;
+
+    sike_fp2sqr_mont(A, jinv);                           // jinv = A^2
+    sike_fp2sqr_mont(C, t1);                             // t1 = C^2
+    sike_fp2add(t1, t1, t0);                             // t0 = t1+t1
+    sike_fp2sub(jinv, t0, t0);                           // t0 = jinv-t0
+    sike_fp2sub(t0, t1, t0);                             // t0 = t0-t1
+    sike_fp2sub(t0, t1, jinv);                           // jinv = t0-t1
+    sike_fp2sqr_mont(t1, t1);                            // t1 = t1^2
+    sike_fp2mul_mont(jinv, t1, jinv);                    // jinv = jinv*t1
+    sike_fp2add(t0, t0, t0);                             // t0 = t0+t0
+    sike_fp2add(t0, t0, t0);                             // t0 = t0+t0
+    sike_fp2sqr_mont(t0, t1);                            // t1 = t0^2
+    sike_fp2mul_mont(t0, t1, t0);                        // t0 = t0*t1
+    sike_fp2add(t0, t0, t0);                             // t0 = t0+t0
+    sike_fp2add(t0, t0, t0);                             // t0 = t0+t0
+    sike_fp2inv_mont(jinv);                              // jinv = 1/jinv
+    sike_fp2mul_mont(jinv, t0, jinv);                    // jinv = t0*jinv
+}
+
+
+void xDBLADD(point_proj_t P, point_proj_t Q, const f2elm_t xPQ, const f2elm_t A24)
+{ // Simultaneous doubling and differential addition.
+  // Input: projective Montgomery points P=(XP:ZP) and Q=(XQ:ZQ) such that xP=XP/ZP and xQ=XQ/ZQ, affine difference xPQ=x(P-Q) and Montgomery curve constant A24=(A+2)/4.
+  // Output: projective Montgomery points P <- 2*P = (X2P:Z2P) such that x(2P)=X2P/Z2P, and Q <- P+Q = (XQP:ZQP) such that = x(Q+P)=XQP/ZQP.
+    f2elm_t t0, t1, t2;
+
+    sike_fp2add(P->X, P->Z, t0);                         // t0 = XP+ZP
+    sike_fp2sub(P->X, P->Z, t1);                         // t1 = XP-ZP
+    sike_fp2sqr_mont(t0, P->X);                          // XP = (XP+ZP)^2
+    sike_fp2sub(Q->X, Q->Z, t2);                         // t2 = XQ-ZQ
+    sike_fp2correction(t2);
+    sike_fp2add(Q->X, Q->Z, Q->X);                       // XQ = XQ+ZQ
+    sike_fp2mul_mont(t0, t2, t0);                        // t0 = (XP+ZP)*(XQ-ZQ)
+    sike_fp2sqr_mont(t1, P->Z);                          // ZP = (XP-ZP)^2
+    sike_fp2mul_mont(t1, Q->X, t1);                      // t1 = (XP-ZP)*(XQ+ZQ)
+    sike_fp2sub(P->X, P->Z, t2);                         // t2 = (XP+ZP)^2-(XP-ZP)^2
+    sike_fp2mul_mont(P->X, P->Z, P->X);                  // XP = (XP+ZP)^2*(XP-ZP)^2
+    sike_fp2mul_mont(t2, A24, Q->X);                     // XQ = A24*[(XP+ZP)^2-(XP-ZP)^2]
+    sike_fp2sub(t0, t1, Q->Z);                           // ZQ = (XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)
+    sike_fp2add(Q->X, P->Z, P->Z);                       // ZP = A24*[(XP+ZP)^2-(XP-ZP)^2]+(XP-ZP)^2
+    sike_fp2add(t0, t1, Q->X);                           // XQ = (XP+ZP)*(XQ-ZQ)+(XP-ZP)*(XQ+ZQ)
+    sike_fp2mul_mont(P->Z, t2, P->Z);                    // ZP = [A24*[(XP+ZP)^2-(XP-ZP)^2]+(XP-ZP)^2]*[(XP+ZP)^2-(XP-ZP)^2]
+    sike_fp2sqr_mont(Q->Z, Q->Z);                        // ZQ = [(XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)]^2
+    sike_fp2sqr_mont(Q->X, Q->X);                        // XQ = [(XP+ZP)*(XQ-ZQ)+(XP-ZP)*(XQ+ZQ)]^2
+    sike_fp2mul_mont(Q->Z, xPQ, Q->Z);                   // ZQ = xPQ*[(XP+ZP)*(XQ-ZQ)-(XP-ZP)*(XQ+ZQ)]^2
+}
diff --git a/third_party/sike/isogeny.h b/third_party/sike/isogeny.h
new file mode 100644
index 0000000..460c8c6
--- /dev/null
+++ b/third_party/sike/isogeny.h
@@ -0,0 +1,49 @@
+#ifndef ISOGENY_H_
+#define ISOGENY_H_
+
+// Computes [2^e](X:Z) on Montgomery curve with projective
+// constant via e repeated doublings.
+void xDBLe(
+    const point_proj_t P, point_proj_t Q, const f2elm_t A24plus,
+    const f2elm_t C24, size_t e);
+// Simultaneous doubling and differential addition.
+void xDBLADD(
+    point_proj_t P, point_proj_t Q, const f2elm_t xPQ,
+    const f2elm_t A24);
+// Tripling of a Montgomery point in projective coordinates (X:Z).
+void xTPL(
+    const point_proj_t P, point_proj_t Q, const f2elm_t A24minus,
+    const f2elm_t A24plus);
+// Computes [3^e](X:Z) on Montgomery curve with projective constant
+// via e repeated triplings.
+void xTPLe(
+    const point_proj_t P, point_proj_t Q, const f2elm_t A24minus,
+    const f2elm_t A24plus, size_t e);
+// Given the x-coordinates of P, Q, and R, returns the value A
+// corresponding to the Montgomery curve E_A: y^2=x^3+A*x^2+x such that R=Q-P on E_A.
+void get_A(
+    const f2elm_t xP, const f2elm_t xQ, const f2elm_t xR, f2elm_t A);
+// Computes the j-invariant of a Montgomery curve with projective constant.
+void j_inv(
+    const f2elm_t A, const f2elm_t C, f2elm_t jinv);
+// Computes the corresponding 4-isogeny of a projective Montgomery
+// point (X4:Z4) of order 4.
+void get_4_isog(
+    const point_proj_t P, f2elm_t A24plus, f2elm_t C24, f2elm_t* coeff);
+// Computes the corresponding 3-isogeny of a projective Montgomery
+// point (X3:Z3) of order 3.
+void get_3_isog(
+    const point_proj_t P, f2elm_t A24minus, f2elm_t A24plus,
+    f2elm_t* coeff);
+// Computes the 3-isogeny R=phi(X:Z), given projective point (X3:Z3)
+// of order 3 on a Montgomery curve and a point P with coefficients given in coeff.
+void eval_3_isog(
+    point_proj_t Q, f2elm_t* coeff);
+// Evaluates the isogeny at the point (X:Z) in the domain of the isogeny.
+void eval_4_isog(
+    point_proj_t P, f2elm_t* coeff);
+// 3-way simultaneous inversion
+void inv_3_way(
+    f2elm_t z1, f2elm_t z2, f2elm_t z3);
+
+#endif // ISOGENY_H_
diff --git a/third_party/sike/sike.c b/third_party/sike/sike.c
new file mode 100644
index 0000000..01add94
--- /dev/null
+++ b/third_party/sike/sike.c
@@ -0,0 +1,571 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: supersingular isogeny key encapsulation (SIKE) protocol
+*********************************************************************************************/
+
+#include <assert.h>
+#include <stdint.h>
+#include <string.h>
+#include <openssl/bn.h>
+#include <openssl/base.h>
+#include <openssl/rand.h>
+#include <openssl/mem.h>
+#include <openssl/hmac.h>
+#include <openssl/sha.h>
+
+#include "utils.h"
+#include "isogeny.h"
+#include "fpx.h"
+
+extern const struct params_t p503;
+
+// Domain separation parameters for HMAC
+static const uint8_t G[2] = {0,0};
+static const uint8_t H[2] = {1,0};
+static const uint8_t F[2] = {2,0};
+
+// SIDHp503_JINV_BYTESZ is a number of bytes used for encoding j-invariant.
+#define SIDHp503_JINV_BYTESZ    126U
+// SIDHp503_PRV_A_BITSZ is a number of bits of SIDH private key (2-isogeny)
+#define SIDHp503_PRV_A_BITSZ    250U
+// SIDHp503_PRV_A_BITSZ is a number of bits of SIDH private key (3-isogeny)
+#define SIDHp503_PRV_B_BITSZ    253U
+// MAX_INT_POINTS_ALICE is a number of points used in 2-isogeny tree computation
+#define MAX_INT_POINTS_ALICE    7U
+// MAX_INT_POINTS_ALICE is a number of points used in 3-isogeny tree computation
+#define MAX_INT_POINTS_BOB      8U
+
+// Produces HMAC-SHA256 of data |S| mac'ed with the key |key|. Result is stored in |out|
+// which must have size of at least |outsz| bytes and must be not bigger than
+// SHA256_DIGEST_LENGTH. The output of a HMAC may be truncated.
+// The |key| buffer is reused by the hmac_sum and hence, it's size must be equal
+// to SHA256_CBLOCK. The HMAC key provided in |key| buffer must be smaller or equal
+// to SHA256_DIGHEST_LENTH. |key| can overlap |out|.
+static void hmac_sum(
+    uint8_t *out, size_t outsz, const uint8_t S[2], uint8_t key[SHA256_CBLOCK]) {
+    for(size_t i=0; i<SHA256_DIGEST_LENGTH; i++) {
+        key[i] = key[i] ^ 0x36;
+    }
+    // set rest of the buffer to ipad = 0x36
+    memset(&key[SHA256_DIGEST_LENGTH], 0x36, SHA256_CBLOCK - SHA256_DIGEST_LENGTH);
+
+    SHA256_CTX ctx;
+    SHA256_Init(&ctx);
+    SHA256_Update(&ctx, key, SHA256_CBLOCK);
+    SHA256_Update(&ctx, S, 2);
+    uint8_t digest[SHA256_DIGEST_LENGTH];
+    SHA256_Final(digest, &ctx);
+
+    // XOR key with an opad = 0x5C
+    for(size_t i=0; i<SHA256_CBLOCK; i++) {
+        key[i] = key[i] ^ 0x36 ^ 0x5C;
+    }
+
+    SHA256_Init(&ctx);
+    SHA256_Update(&ctx, key, SHA256_CBLOCK);
+    SHA256_Update(&ctx, digest, SHA256_DIGEST_LENGTH);
+    SHA256_Final(digest, &ctx);
+    assert(outsz <= sizeof(digest));
+    memcpy(out, digest, outsz);
+}
+
+// Swap points.
+// If option = 0 then P <- P and Q <- Q, else if option = 0xFF...FF then P <- Q and Q <- P
+#if !defined(OPENSSL_X86_64) || defined(OPENSSL_NO_ASM)
+static void sike_cswap(point_proj_t P, point_proj_t Q, const crypto_word_t option)
+{
+    crypto_word_t temp;
+    for (size_t i = 0; i < NWORDS_FIELD; i++) {
+        temp = option & (P->X->c0[i] ^ Q->X->c0[i]);
+        P->X->c0[i] = temp ^ P->X->c0[i];
+        Q->X->c0[i] = temp ^ Q->X->c0[i];
+        temp = option & (P->Z->c0[i] ^ Q->Z->c0[i]);
+        P->Z->c0[i] = temp ^ P->Z->c0[i];
+        Q->Z->c0[i] = temp ^ Q->Z->c0[i];
+        temp = option & (P->X->c1[i] ^ Q->X->c1[i]);
+        P->X->c1[i] = temp ^ P->X->c1[i];
+        Q->X->c1[i] = temp ^ Q->X->c1[i];
+        temp = option & (P->Z->c1[i] ^ Q->Z->c1[i]);
+        P->Z->c1[i] = temp ^ P->Z->c1[i];
+        Q->Z->c1[i] = temp ^ Q->Z->c1[i];
+    }
+}
+#endif
+
+// Swap points.
+// If option = 0 then P <- P and Q <- Q, else if option = 0xFF...FF then P <- Q and Q <- P
+static inline void sike_fp2cswap(point_proj_t P, point_proj_t Q, const crypto_word_t option)
+{
+#if defined(OPENSSL_X86_64) && !defined(OPENSSL_NO_ASM)
+    sike_cswap_asm(P, Q, option);
+#else
+    sike_cswap(P, Q, option);
+#endif
+}
+
+static void LADDER3PT(
+    const f2elm_t xP, const f2elm_t xQ, const f2elm_t xPQ, const crypto_word_t* m,
+    int is_A, point_proj_t R, const f2elm_t A) {
+    point_proj_t R0 = POINT_PROJ_INIT, R2 = POINT_PROJ_INIT;
+    f2elm_t A24 = F2ELM_INIT;
+    crypto_word_t mask;
+    int bit, swap, prevbit = 0;
+
+    const size_t nbits = is_A?SIDHp503_PRV_A_BITSZ:SIDHp503_PRV_B_BITSZ;
+
+    // Initializing constant
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, A24[0].c0);
+    sike_fp2add(A24, A24, A24);
+    sike_fp2add(A, A24, A24);
+    sike_fp2div2(A24, A24);
+    sike_fp2div2(A24, A24); // A24 = (A+2)/4
+
+    // Initializing points
+    sike_fp2copy(xQ, R0->X);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, R0->Z[0].c0);
+    sike_fp2copy(xPQ, R2->X);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, R2->Z[0].c0);
+    sike_fp2copy(xP, R->X);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, R->Z[0].c0);
+    memset(R->Z->c1, 0, sizeof(R->Z->c1));
+
+    // Main loop
+    for (size_t i = 0; i < nbits; i++) {
+        bit = (m[i >> LOG2RADIX] >> (i & (RADIX-1))) & 1;
+        swap = bit ^ prevbit;
+        prevbit = bit;
+        mask = 0 - (crypto_word_t)swap;
+
+        sike_fp2cswap(R, R2, mask);
+        xDBLADD(R0, R2, R->X, A24);
+        sike_fp2mul_mont(R2->X, R->Z, R2->X);
+    }
+}
+
+// Initialization of basis points
+static inline void sike_init_basis(const crypto_word_t *gen, f2elm_t XP, f2elm_t XQ, f2elm_t XR) {
+    sike_fpcopy(gen,                  XP->c0);
+    sike_fpcopy(gen +   NWORDS_FIELD, XP->c1);
+    sike_fpcopy(gen + 2*NWORDS_FIELD, XQ->c0);
+    memset(XQ->c1, 0, sizeof(XQ->c1));
+    sike_fpcopy(gen + 3*NWORDS_FIELD, XR->c0);
+    sike_fpcopy(gen + 4*NWORDS_FIELD, XR->c1);
+}
+
+// Conversion of GF(p^2) element from Montgomery to standard representation.
+static inline void sike_fp2_encode(const f2elm_t x, uint8_t *enc) {
+    f2elm_t t;
+    sike_from_fp2mont(x, t);
+
+    // convert to bytes in little endian form
+    for (size_t i=0; i<FIELD_BYTESZ; i++) {
+        enc[i+           0] = (t[0].c0[i/LSZ] >> (8*(i%LSZ))) & 0xFF;
+        enc[i+FIELD_BYTESZ] = (t[0].c1[i/LSZ] >> (8*(i%LSZ))) & 0xFF;
+    }
+}
+
+// Parse byte sequence back into GF(p^2) element, and conversion to Montgomery representation.
+// Elements over GF(p503) are encoded in 63 octets in little endian format
+// (i.e., the least significant octet is located in the lowest memory address).
+static inline void fp2_decode(const uint8_t *enc, f2elm_t t) {
+    memset(t[0].c0, 0, sizeof(t[0].c0));
+    memset(t[0].c1, 0, sizeof(t[0].c1));
+    // convert bytes in little endian form to f2elm_t
+    for (size_t i = 0; i < FIELD_BYTESZ; i++) {
+        t[0].c0[i/LSZ] |= ((crypto_word_t)enc[i+           0]) << (8*(i%LSZ));
+        t[0].c1[i/LSZ] |= ((crypto_word_t)enc[i+FIELD_BYTESZ]) << (8*(i%LSZ));
+    }
+    sike_to_fp2mont(t, t);
+}
+
+// Alice's ephemeral public key generation
+// Input:  a private key prA in the range [0, 2^250 - 1], stored in 32 bytes.
+// Output: the public key pkA consisting of 3 GF(p503^2) elements encoded in 378 bytes.
+static void gen_iso_A(const uint8_t* skA, uint8_t* pkA)
+{
+    point_proj_t R, pts[MAX_INT_POINTS_ALICE];
+    point_proj_t phiP = POINT_PROJ_INIT;
+    point_proj_t phiQ = POINT_PROJ_INIT;
+    point_proj_t phiR = POINT_PROJ_INIT;
+    f2elm_t XPA, XQA, XRA, coeff[3];
+    f2elm_t A24plus = F2ELM_INIT;
+    f2elm_t C24 = F2ELM_INIT;
+    f2elm_t A = F2ELM_INIT;
+    unsigned int m, index = 0, pts_index[MAX_INT_POINTS_ALICE], npts = 0, ii = 0;
+
+    // Initialize basis points
+    sike_init_basis((crypto_word_t*)p503.A_gen, XPA, XQA, XRA);
+    sike_init_basis((crypto_word_t*)p503.B_gen, phiP->X, phiQ->X, phiR->X);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, (phiP->Z)->c0);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, (phiQ->Z)->c0);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, (phiR->Z)->c0);
+
+    // Initialize constants
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, A24plus->c0);
+    sike_fp2add(A24plus, A24plus, C24);
+
+    // Retrieve kernel point
+    LADDER3PT(XPA, XQA, XRA, (crypto_word_t*)skA, 1, R, A);
+
+    // Traverse tree
+    index = 0;
+    for (size_t row = 1; row < A_max; row++) {
+        while (index < A_max-row) {
+            sike_fp2copy(R->X, pts[npts]->X);
+            sike_fp2copy(R->Z, pts[npts]->Z);
+            pts_index[npts++] = index;
+            m = p503.A_strat[ii++];
+            xDBLe(R, R, A24plus, C24, (2*m));
+            index += m;
+        }
+        get_4_isog(R, A24plus, C24, coeff);
+
+        for (size_t i = 0; i < npts; i++) {
+            eval_4_isog(pts[i], coeff);
+        }
+        eval_4_isog(phiP, coeff);
+        eval_4_isog(phiQ, coeff);
+        eval_4_isog(phiR, coeff);
+
+        sike_fp2copy(pts[npts-1]->X, R->X);
+        sike_fp2copy(pts[npts-1]->Z, R->Z);
+        index = pts_index[npts-1];
+        npts -= 1;
+    }
+
+    get_4_isog(R, A24plus, C24, coeff);
+    eval_4_isog(phiP, coeff);
+    eval_4_isog(phiQ, coeff);
+    eval_4_isog(phiR, coeff);
+
+    inv_3_way(phiP->Z, phiQ->Z, phiR->Z);
+    sike_fp2mul_mont(phiP->X, phiP->Z, phiP->X);
+    sike_fp2mul_mont(phiQ->X, phiQ->Z, phiQ->X);
+    sike_fp2mul_mont(phiR->X, phiR->Z, phiR->X);
+
+    // Format public key
+    sike_fp2_encode(phiP->X, pkA);
+    sike_fp2_encode(phiQ->X, pkA + SIDHp503_JINV_BYTESZ);
+    sike_fp2_encode(phiR->X, pkA + 2*SIDHp503_JINV_BYTESZ);
+}
+
+// Bob's ephemeral key-pair generation
+// It produces a private key skB and computes the public key pkB.
+// The private key is an integer in the range [0, 2^Floor(Log(2,3^159)) - 1], stored in 32 bytes.
+// The public key consists of 3 GF(p503^2) elements encoded in 378 bytes.
+static void gen_iso_B(const uint8_t* skB, uint8_t* pkB)
+{
+    point_proj_t R, pts[MAX_INT_POINTS_BOB];
+    point_proj_t phiP = POINT_PROJ_INIT;
+    point_proj_t phiQ = POINT_PROJ_INIT;
+    point_proj_t phiR = POINT_PROJ_INIT;
+    f2elm_t XPB, XQB, XRB, coeff[3];
+    f2elm_t A24plus = F2ELM_INIT;
+    f2elm_t A24minus = F2ELM_INIT;
+    f2elm_t A = F2ELM_INIT;
+    unsigned int m, index = 0, pts_index[MAX_INT_POINTS_BOB], npts = 0, ii = 0;
+
+    // Initialize basis points
+    sike_init_basis((crypto_word_t*)p503.B_gen, XPB, XQB, XRB);
+    sike_init_basis((crypto_word_t*)p503.A_gen, phiP->X, phiQ->X, phiR->X);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, (phiP->Z)->c0);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, (phiQ->Z)->c0);
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, (phiR->Z)->c0);
+
+    // Initialize constants
+    sike_fpcopy((crypto_word_t*)&p503.mont_one, A24plus->c0);
+    sike_fp2add(A24plus, A24plus, A24plus);
+    sike_fp2copy(A24plus, A24minus);
+    sike_fp2neg(A24minus);
+
+    // Retrieve kernel point
+    LADDER3PT(XPB, XQB, XRB, (crypto_word_t*)skB, 0, R, A);
+
+    // Traverse tree
+    index = 0;
+    for (size_t row = 1; row < B_max; row++) {
+        while (index < B_max-row) {
+            sike_fp2copy(R->X, pts[npts]->X);
+            sike_fp2copy(R->Z, pts[npts]->Z);
+            pts_index[npts++] = index;
+            m = p503.B_strat[ii++];
+            xTPLe(R, R, A24minus, A24plus, m);
+            index += m;
+        }
+        get_3_isog(R, A24minus, A24plus, coeff);
+
+        for (size_t i = 0; i < npts; i++) {
+            eval_3_isog(pts[i], coeff);
+        }
+        eval_3_isog(phiP, coeff);
+        eval_3_isog(phiQ, coeff);
+        eval_3_isog(phiR, coeff);
+
+        sike_fp2copy(pts[npts-1]->X, R->X);
+        sike_fp2copy(pts[npts-1]->Z, R->Z);
+        index = pts_index[npts-1];
+        npts -= 1;
+    }
+
+    get_3_isog(R, A24minus, A24plus, coeff);
+    eval_3_isog(phiP, coeff);
+    eval_3_isog(phiQ, coeff);
+    eval_3_isog(phiR, coeff);
+
+    inv_3_way(phiP->Z, phiQ->Z, phiR->Z);
+    sike_fp2mul_mont(phiP->X, phiP->Z, phiP->X);
+    sike_fp2mul_mont(phiQ->X, phiQ->Z, phiQ->X);
+    sike_fp2mul_mont(phiR->X, phiR->Z, phiR->X);
+
+    // Format public key
+    sike_fp2_encode(phiP->X, pkB);
+    sike_fp2_encode(phiQ->X, pkB + SIDHp503_JINV_BYTESZ);
+    sike_fp2_encode(phiR->X, pkB + 2*SIDHp503_JINV_BYTESZ);
+}
+
+// Alice's ephemeral shared secret computation
+// It produces a shared secret key ssA using her secret key skA and Bob's public key pkB
+// Inputs: Alice's skA is an integer in the range [0, 2^250 - 1], stored in 32 bytes.
+//         Bob's pkB consists of 3 GF(p503^2) elements encoded in 378 bytes.
+// Output: a shared secret ssA that consists of one element in GF(p503^2) encoded in 126 bytes.
+static void ex_iso_A(const uint8_t* skA, const uint8_t* pkB, uint8_t* ssA)
+{
+    point_proj_t R, pts[MAX_INT_POINTS_ALICE];
+    f2elm_t coeff[3], PKB[3], jinv;
+    f2elm_t A24plus = F2ELM_INIT;
+    f2elm_t C24 = F2ELM_INIT;
+    f2elm_t A = F2ELM_INIT;
+    unsigned int m, index = 0, pts_index[MAX_INT_POINTS_ALICE], npts = 0, ii = 0;
+
+    // Initialize images of Bob's basis
+    fp2_decode(pkB, PKB[0]);
+    fp2_decode(pkB + SIDHp503_JINV_BYTESZ, PKB[1]);
+    fp2_decode(pkB + 2*SIDHp503_JINV_BYTESZ, PKB[2]);
+
+    // Initialize constants
+    get_A(PKB[0], PKB[1], PKB[2], A); // TODO: Can return projective A?
+    sike_fpadd((crypto_word_t*)&p503.mont_one, (crypto_word_t*)&p503.mont_one, C24->c0);
+    sike_fp2add(A, C24, A24plus);
+    sike_fpadd(C24->c0, C24->c0, C24->c0);
+
+    // Retrieve kernel point
+    LADDER3PT(PKB[0], PKB[1], PKB[2], (crypto_word_t*)skA, 1, R, A);
+
+    // Traverse tree
+    index = 0;
+    for (size_t row = 1; row < A_max; row++) {
+        while (index < A_max-row) {
+            sike_fp2copy(R->X, pts[npts]->X);
+            sike_fp2copy(R->Z, pts[npts]->Z);
+            pts_index[npts++] = index;
+            m = p503.A_strat[ii++];
+            xDBLe(R, R, A24plus, C24, (2*m));
+            index += m;
+        }
+        get_4_isog(R, A24plus, C24, coeff);
+
+        for (size_t i = 0; i < npts; i++) {
+            eval_4_isog(pts[i], coeff);
+        }
+
+        sike_fp2copy(pts[npts-1]->X, R->X);
+        sike_fp2copy(pts[npts-1]->Z, R->Z);
+        index = pts_index[npts-1];
+        npts -= 1;
+    }
+
+    get_4_isog(R, A24plus, C24, coeff);
+    sike_fp2div2(C24, C24);
+    sike_fp2sub(A24plus, C24, A24plus);
+    sike_fp2div2(C24, C24);
+    j_inv(A24plus, C24, jinv);
+    sike_fp2_encode(jinv, ssA);
+}
+
+// Bob's ephemeral shared secret computation
+// It produces a shared secret key ssB using his secret key skB and Alice's public key pkA
+// Inputs: Bob's skB is an integer in the range [0, 2^Floor(Log(2,3^159)) - 1], stored in 32 bytes.
+//         Alice's pkA consists of 3 GF(p503^2) elements encoded in 378 bytes.
+// Output: a shared secret ssB that consists of one element in GF(p503^2) encoded in 126 bytes.
+static void ex_iso_B(const uint8_t* skB, const uint8_t* pkA, uint8_t* ssB)
+{
+    point_proj_t R, pts[MAX_INT_POINTS_BOB];
+    f2elm_t coeff[3], PKB[3], jinv;
+    f2elm_t A24plus = F2ELM_INIT;
+    f2elm_t A24minus = F2ELM_INIT;
+    f2elm_t A = F2ELM_INIT;
+    unsigned int m, index = 0, pts_index[MAX_INT_POINTS_BOB], npts = 0, ii = 0;
+
+    // Initialize images of Alice's basis
+    fp2_decode(pkA, PKB[0]);
+    fp2_decode(pkA + SIDHp503_JINV_BYTESZ, PKB[1]);
+    fp2_decode(pkA + 2*SIDHp503_JINV_BYTESZ, PKB[2]);
+
+    // Initialize constants
+    get_A(PKB[0], PKB[1], PKB[2], A);
+    sike_fpadd((crypto_word_t*)&p503.mont_one, (crypto_word_t*)&p503.mont_one, A24minus->c0);
+    sike_fp2add(A, A24minus, A24plus);
+    sike_fp2sub(A, A24minus, A24minus);
+
+    // Retrieve kernel point
+    LADDER3PT(PKB[0], PKB[1], PKB[2], (crypto_word_t*)skB, 0, R, A);
+
+    // Traverse tree
+    index = 0;
+    for (size_t row = 1; row < B_max; row++) {
+        while (index < B_max-row) {
+            sike_fp2copy(R->X, pts[npts]->X);
+            sike_fp2copy(R->Z, pts[npts]->Z);
+            pts_index[npts++] = index;
+            m = p503.B_strat[ii++];
+            xTPLe(R, R, A24minus, A24plus, m);
+            index += m;
+        }
+        get_3_isog(R, A24minus, A24plus, coeff);
+
+        for (size_t i = 0; i < npts; i++) {
+            eval_3_isog(pts[i], coeff);
+        }
+
+        sike_fp2copy(pts[npts-1]->X, R->X);
+        sike_fp2copy(pts[npts-1]->Z, R->Z);
+        index = pts_index[npts-1];
+        npts -= 1;
+    }
+
+    get_3_isog(R, A24minus, A24plus, coeff);
+    sike_fp2add(A24plus, A24minus, A);
+    sike_fp2add(A, A, A);
+    sike_fp2sub(A24plus, A24minus, A24plus);
+    j_inv(A, A24plus, jinv);
+    sike_fp2_encode(jinv, ssB);
+}
+
+int SIKE_keypair(uint8_t out_priv[SIKEp503_PRV_BYTESZ],
+                 uint8_t out_pub[SIKEp503_PUB_BYTESZ]) {
+  int ret = 0;
+
+  // Calculate private key for Alice. Needs to be in range [0, 2^0xFA - 1] and <
+  // 253 bits
+  BIGNUM *bn_sidh_prv = BN_new();
+  if (!bn_sidh_prv ||
+      !BN_rand(bn_sidh_prv, SIDHp503_PRV_B_BITSZ, BN_RAND_TOP_ONE,
+               BN_RAND_BOTTOM_ANY) ||
+      !BN_bn2le_padded(out_priv, BITS_TO_BYTES(SIDHp503_PRV_B_BITSZ),
+                       bn_sidh_prv)) {
+    goto end;
+  }
+
+  gen_iso_B(out_priv, out_pub);
+  ret = 1;
+
+end:
+  BN_free(bn_sidh_prv);
+  return ret;
+}
+
+void SIKE_encaps(uint8_t out_shared_key[SIKEp503_SS_BYTESZ],
+                 uint8_t out_ciphertext[SIKEp503_CT_BYTESZ],
+                 const uint8_t pub_key[SIKEp503_PUB_BYTESZ]) {
+  // Secret buffer is reused by the function to store some ephemeral
+  // secret data. It's size must be maximum of SHA256_CBLOCK,
+  // SIKEp503_MSG_BYTESZ and SIDHp503_PRV_A_BITSZ in bytes.
+  uint8_t secret[SHA256_CBLOCK];
+  uint8_t j[SIDHp503_JINV_BYTESZ];
+  uint8_t temp[SIKEp503_MSG_BYTESZ + SIKEp503_CT_BYTESZ];
+  SHA256_CTX ctx;
+
+  // Generate secret key for A
+  // secret key A = HMAC({0,1}^n || pub_key), G) mod SIDHp503_PRV_A_BITSZ
+  RAND_bytes(temp, SIKEp503_MSG_BYTESZ);
+
+  SHA256_Init(&ctx);
+  SHA256_Update(&ctx, temp, SIKEp503_MSG_BYTESZ);
+  SHA256_Update(&ctx, pub_key, SIKEp503_PUB_BYTESZ);
+  SHA256_Final(secret, &ctx);
+  hmac_sum(secret, BITS_TO_BYTES(SIDHp503_PRV_A_BITSZ), G, secret);
+  secret[BITS_TO_BYTES(SIDHp503_PRV_A_BITSZ) - 1] &=
+      (1 << (SIDHp503_PRV_A_BITSZ % 8)) - 1;
+
+  // Generate public key for A - first part of the ciphertext
+  gen_iso_A(secret, out_ciphertext);
+
+  // Generate c1:
+  //  h = HMAC(j-invariant(secret key A, public key B), F)
+  // c1 = h ^ m
+  ex_iso_A(secret, pub_key, j);
+  SHA256_Init(&ctx);
+  SHA256_Update(&ctx, j, sizeof(j));
+  SHA256_Final(secret, &ctx);
+  hmac_sum(secret, SIKEp503_MSG_BYTESZ, F, secret);
+
+  // c1 = h ^ m
+  uint8_t *c1 = &out_ciphertext[SIKEp503_PUB_BYTESZ];
+  for (size_t i = 0; i < SIKEp503_MSG_BYTESZ; i++) {
+    c1[i] = temp[i] ^ secret[i];
+  }
+
+  SHA256_Init(&ctx);
+  SHA256_Update(&ctx, temp, SIKEp503_MSG_BYTESZ);
+  SHA256_Update(&ctx, out_ciphertext, SIKEp503_CT_BYTESZ);
+  SHA256_Final(secret, &ctx);
+  // Generate shared secret out_shared_key = HMAC(m||out_ciphertext, F)
+  hmac_sum(out_shared_key, SIKEp503_SS_BYTESZ, H, secret);
+}
+
+void SIKE_decaps(uint8_t out_shared_key[SIKEp503_SS_BYTESZ],
+                 const uint8_t ciphertext[SIKEp503_CT_BYTESZ],
+                 const uint8_t pub_key[SIKEp503_PUB_BYTESZ],
+                 const uint8_t priv_key[SIKEp503_PRV_BYTESZ]) {
+  // Secret buffer is reused by the function to store some ephemeral
+  // secret data. It's size must be maximum of SHA256_CBLOCK,
+  // SIKEp503_MSG_BYTESZ and SIDHp503_PRV_A_BITSZ in bytes.
+  uint8_t secret[SHA256_CBLOCK];
+  uint8_t j[SIDHp503_JINV_BYTESZ];
+  uint8_t c0[SIKEp503_PUB_BYTESZ];
+  uint8_t temp[SIKEp503_MSG_BYTESZ];
+  uint8_t shared_nok[SIKEp503_MSG_BYTESZ];
+  SHA256_CTX ctx;
+
+  RAND_bytes(shared_nok, SIKEp503_MSG_BYTESZ);
+
+  // Recover m
+  // Let ciphertext = c0 || c1 - both have fixed sizes
+  // m = F(j-invariant(c0, priv_key)) ^ c1
+  ex_iso_B(priv_key, ciphertext, j);
+
+  SHA256_Init(&ctx);
+  SHA256_Update(&ctx, j, sizeof(j));
+  SHA256_Final(secret, &ctx);
+  hmac_sum(secret, SIKEp503_MSG_BYTESZ, F, secret);
+
+  const uint8_t *c1 = &ciphertext[sizeof(c0)];
+  for (size_t i = 0; i < SIKEp503_MSG_BYTESZ; i++) {
+    temp[i] = c1[i] ^ secret[i];
+  }
+
+  SHA256_Init(&ctx);
+  SHA256_Update(&ctx, temp, SIKEp503_MSG_BYTESZ);
+  SHA256_Update(&ctx, pub_key, SIKEp503_PUB_BYTESZ);
+  SHA256_Final(secret, &ctx);
+  hmac_sum(secret, BITS_TO_BYTES(SIDHp503_PRV_A_BITSZ), G, secret);
+
+  // Recover secret key A = G(m||pub_key) mod
+  secret[BITS_TO_BYTES(SIDHp503_PRV_A_BITSZ) - 1] &=
+      (1 << (SIDHp503_PRV_A_BITSZ % 8)) - 1;
+
+  // Recover c0 = public key A
+  gen_iso_A(secret, c0);
+  crypto_word_t ok = constant_time_is_zero_w(
+      CRYPTO_memcmp(c0, ciphertext, SIKEp503_PUB_BYTESZ));
+  for (size_t i = 0; i < SIKEp503_MSG_BYTESZ; i++) {
+    temp[i] = constant_time_select_8(ok, temp[i], shared_nok[i]);
+  }
+
+  SHA256_Init(&ctx);
+  SHA256_Update(&ctx, temp, SIKEp503_MSG_BYTESZ);
+  SHA256_Update(&ctx, ciphertext, SIKEp503_CT_BYTESZ);
+  SHA256_Final(secret, &ctx);
+  hmac_sum(out_shared_key, SIKEp503_SS_BYTESZ, H, secret);
+}
diff --git a/third_party/sike/sike.h b/third_party/sike/sike.h
new file mode 100644
index 0000000..09093cd
--- /dev/null
+++ b/third_party/sike/sike.h
@@ -0,0 +1,64 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: API header file for SIKE
+*********************************************************************************************/
+
+#ifndef SIKE_H_
+#define SIKE_H_
+
+#include <stdint.h>
+#include <openssl/base.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* SIKEp503
+ *
+ * SIKE is a isogeny based post-quantum key encapsulation mechanism. Description of the
+ * algorithm is provided in [SIKE]. This implementation uses 503-bit field size. The code
+ * is based on "Additional_Implementations" from PQC NIST submission package which can
+ * be found here:
+ * https://csrc.nist.gov/CSRC/media/Projects/Post-Quantum-Cryptography/documents/round-1/submissions/SIKE.zip
+ *
+ * [SIKE] https://sike.org/files/SIDH-spec.pdf
+ */
+
+// SIKEp503_PUB_BYTESZ is the number of bytes in a public key.
+#define SIKEp503_PUB_BYTESZ 378
+// SIKEp503_PRV_BYTESZ is the number of bytes in a private key.
+#define SIKEp503_PRV_BYTESZ 32
+// SIKEp503_SS_BYTESZ is the number of bytes in a shared key.
+#define SIKEp503_SS_BYTESZ  16
+// SIKEp503_MSG_BYTESZ is the number of bytes in a random bit string concatenated
+// with the public key (see 1.4 of SIKE).
+#define SIKEp503_MSG_BYTESZ 24
+// SIKEp503_SS_BYTESZ is the number of bytes in a ciphertext.
+#define SIKEp503_CT_BYTESZ  (SIKEp503_PUB_BYTESZ + SIKEp503_MSG_BYTESZ)
+
+// SIKE_keypair outputs a public and secret key. Internally it uses BN_rand() as
+// an entropy source. In case of success function returns 1, otherwise 0.
+OPENSSL_EXPORT int SIKE_keypair(
+    uint8_t out_priv[SIKEp503_PRV_BYTESZ],
+    uint8_t out_pub[SIKEp503_PUB_BYTESZ]);
+
+// SIKE_encaps generates and encrypts a random session key, writing those values to
+// |out_shared_key| and |out_ciphertext|, respectively.
+OPENSSL_EXPORT void SIKE_encaps(
+    uint8_t out_shared_key[SIKEp503_SS_BYTESZ],
+    uint8_t out_ciphertext[SIKEp503_CT_BYTESZ],
+    const uint8_t pub_key[SIKEp503_PUB_BYTESZ]);
+
+// SIKE_decaps outputs a random session key, writing it to |out_shared_key|.
+OPENSSL_EXPORT void SIKE_decaps(
+    uint8_t out_shared_key[SIKEp503_SS_BYTESZ],
+    const uint8_t ciphertext[SIKEp503_CT_BYTESZ],
+    const uint8_t pub_key[SIKEp503_PUB_BYTESZ],
+    const uint8_t priv_key[SIKEp503_PRV_BYTESZ]);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/third_party/sike/sike_test.cc b/third_party/sike/sike_test.cc
new file mode 100644
index 0000000..a1426ef
--- /dev/null
+++ b/third_party/sike/sike_test.cc
@@ -0,0 +1,208 @@
+/* Copyright (c) 2018, Google Inc.
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+ * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
+ * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
+ * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */
+
+#include <stdint.h>
+
+#include <gtest/gtest.h>
+
+#include "../../crypto/test/abi_test.h"
+#include "sike.h"
+#include "fpx.h"
+
+TEST(SIKE, RoundTrip) {
+    uint8_t sk[SIKEp503_PRV_BYTESZ] = {0};
+    uint8_t pk[SIKEp503_PUB_BYTESZ] = {0};
+    uint8_t ct[SIKEp503_CT_BYTESZ] = {0};
+    uint8_t ss_enc[SIKEp503_SS_BYTESZ] = {0};
+    uint8_t ss_dec[SIKEp503_SS_BYTESZ] = {0};
+
+    EXPECT_EQ(SIKE_keypair(sk, pk), 1);
+    SIKE_encaps(ss_enc, ct, pk);
+    SIKE_decaps(ss_dec, ct, pk, sk);
+
+    EXPECT_EQ(memcmp(ss_enc, ss_dec, SIKEp503_SS_BYTESZ), 0);
+}
+
+TEST(SIKE, Decapsulation) {
+    const uint8_t sk[SIKEp503_PRV_BYTESZ] = {
+        0xDB, 0xAF, 0x2C, 0x89, 0xCA, 0x5A, 0xD4, 0x9D, 0x4F, 0x13,
+        0x40, 0xDF, 0x2D, 0xB1, 0x5F, 0x4C, 0x91, 0xA7, 0x1F, 0x0B,
+        0x29, 0x15, 0x01, 0x59, 0xBC, 0x5F, 0x0B, 0x4A, 0x03, 0x27,
+        0x6F, 0x18};
+
+    const uint8_t pk[SIKEp503_PUB_BYTESZ] = {
+        0x07, 0xAA, 0x51, 0x45, 0x3E, 0x1F, 0x53, 0x2A, 0x0A, 0x05,
+        0x46, 0xF6, 0x54, 0x7F, 0x5D, 0x56, 0xD6, 0x76, 0xD3, 0xEA,
+        0x4B, 0x6B, 0x01, 0x9B, 0x11, 0x72, 0x6F, 0x75, 0xEA, 0x34,
+        0x3C, 0x28, 0x2C, 0x36, 0xFD, 0x77, 0xDA, 0xBE, 0xB6, 0x20,
+        0x18, 0xC1, 0x93, 0x98, 0x18, 0x86, 0x30, 0x2F, 0x2E, 0xD2,
+        0x00, 0x61, 0xFF, 0xAE, 0x78, 0xAE, 0xFB, 0x6F, 0x32, 0xAC,
+        0x06, 0xBF, 0x35, 0xF6, 0xF7, 0x5B, 0x98, 0x26, 0x95, 0xC2,
+        0xD8, 0xD6, 0x1C, 0x0E, 0x47, 0xDA, 0x76, 0xCE, 0xB5, 0xF1,
+        0x19, 0xCC, 0x01, 0xE1, 0x17, 0xA9, 0x62, 0xF7, 0x82, 0x6C,
+        0x25, 0x51, 0x25, 0xAE, 0xFE, 0xE3, 0xE2, 0xE1, 0x35, 0xAE,
+        0x2E, 0x8F, 0x38, 0xE0, 0x7C, 0x74, 0x3C, 0x1D, 0x39, 0x91,
+        0x1B, 0xC7, 0x9F, 0x8E, 0x33, 0x4E, 0x84, 0x19, 0xB8, 0xD9,
+        0xC2, 0x71, 0x35, 0x02, 0x47, 0x3E, 0x79, 0xEF, 0x47, 0xE1,
+        0xD8, 0x21, 0x96, 0x1F, 0x11, 0x59, 0x39, 0x34, 0x76, 0xEF,
+        0x3E, 0xB7, 0x4E, 0xFB, 0x7C, 0x55, 0xA1, 0x85, 0xAA, 0xAB,
+        0xAD, 0xF0, 0x09, 0xCB, 0xD1, 0xE3, 0x7C, 0x4F, 0x5D, 0x2D,
+        0xE1, 0x13, 0xF0, 0x71, 0xD9, 0xE5, 0xF6, 0xAF, 0x7F, 0xC1,
+        0x27, 0x95, 0x8D, 0x52, 0xD5, 0x96, 0x42, 0x38, 0x41, 0xF7,
+        0x24, 0x3F, 0x3A, 0xB5, 0x7E, 0x11, 0xE4, 0xF9, 0x33, 0xEE,
+        0x4D, 0xBE, 0x74, 0x48, 0xF9, 0x98, 0x04, 0x01, 0x16, 0xEB,
+        0xA9, 0x0D, 0x61, 0xC6, 0xFD, 0x4C, 0xCF, 0x98, 0x84, 0x4A,
+        0x94, 0xAC, 0x69, 0x2C, 0x02, 0x8B, 0xE3, 0xD1, 0x41, 0x0D,
+        0xF2, 0x2D, 0x46, 0x1F, 0x57, 0x1C, 0x77, 0x86, 0x18, 0xE3,
+        0x63, 0xDE, 0xF3, 0xE3, 0x02, 0x30, 0x54, 0x73, 0xAE, 0xC2,
+        0x32, 0xA2, 0xCE, 0xEB, 0xCF, 0x81, 0x46, 0x54, 0x5C, 0xF4,
+        0x5D, 0x2A, 0x03, 0x5D, 0x9C, 0xAE, 0xE0, 0x60, 0x03, 0x80,
+        0x11, 0x30, 0xA5, 0xAA, 0xD1, 0x75, 0x67, 0xE0, 0x1C, 0x2B,
+        0x6B, 0x5D, 0x83, 0xDE, 0x92, 0x9B, 0x0E, 0xD7, 0x11, 0x0F,
+        0x00, 0xC4, 0x59, 0xE4, 0x81, 0x04, 0x3B, 0xEE, 0x5C, 0x04,
+        0xD1, 0x0E, 0xD0, 0x67, 0xF5, 0xCC, 0xAA, 0x72, 0x73, 0xEA,
+        0xC4, 0x76, 0x99, 0x3B, 0x4C, 0x90, 0x2F, 0xCB, 0xD8, 0x0A,
+        0x5B, 0xEC, 0x0E, 0x0E, 0x1F, 0x59, 0xEA, 0x14, 0x8D, 0x34,
+        0x53, 0x65, 0x4C, 0x1A, 0x59, 0xA8, 0x95, 0x66, 0x60, 0xBB,
+        0xC4, 0xCC, 0x32, 0xA9, 0x8D, 0x2A, 0xAA, 0x14, 0x6F, 0x0F,
+        0x81, 0x4D, 0x32, 0x02, 0xFD, 0x33, 0x58, 0x42, 0xCF, 0xF3,
+        0x67, 0xD0, 0x9F, 0x0B, 0xB1, 0xCC, 0x18, 0xA5, 0xC4, 0x19,
+        0xB6, 0x00, 0xED, 0xFA, 0x32, 0x1A, 0x5F, 0x67, 0xC8, 0xC3,
+        0xEB, 0x0D, 0xB5, 0x9A, 0x36, 0x47, 0x82, 0x00};
+
+    const uint8_t ct_exp[SIKEp503_CT_BYTESZ] = {
+        0xE6, 0xB7, 0xE5, 0x7B, 0xA9, 0x19, 0xD1, 0x2C, 0xB8, 0x5C,
+        0x7B, 0x66, 0x74, 0xB0, 0x71, 0xA1, 0xFF, 0x71, 0x7F, 0x4B,
+        0xB5, 0xA6, 0xAF, 0x48, 0x32, 0x52, 0xD5, 0x82, 0xEE, 0x8A,
+        0xBB, 0x08, 0x1E, 0xF6, 0xAC, 0x91, 0xA2, 0xCB, 0x6B, 0x6A,
+        0x09, 0x2B, 0xD9, 0xC6, 0x27, 0xD6, 0x3A, 0x6B, 0x8D, 0xFC,
+        0xB8, 0x90, 0x8F, 0x72, 0xB3, 0xFA, 0x7D, 0x34, 0x7A, 0xC4,
+        0x7E, 0xE3, 0x30, 0xC5, 0xA0, 0xFE, 0x3D, 0x43, 0x14, 0x4E,
+        0x3A, 0x14, 0x76, 0x3E, 0xFB, 0xDF, 0xE3, 0xA8, 0xE3, 0x5E,
+        0x38, 0xF2, 0xE0, 0x39, 0x67, 0x60, 0xFD, 0xFB, 0xB4, 0x19,
+        0xCD, 0xE1, 0x93, 0xA2, 0x06, 0xCC, 0x65, 0xCD, 0x6E, 0xC8,
+        0xB4, 0x5E, 0x41, 0x4B, 0x6C, 0xA5, 0xF4, 0xE4, 0x9D, 0x52,
+        0x8C, 0x25, 0x60, 0xDD, 0x3D, 0xA9, 0x7F, 0xF2, 0x88, 0xC1,
+        0x0C, 0xEE, 0x97, 0xE0, 0xE7, 0x3B, 0xB7, 0xD3, 0x6F, 0x28,
+        0x79, 0x2F, 0x50, 0xB2, 0x4F, 0x74, 0x3A, 0x0C, 0x88, 0x27,
+        0x98, 0x3A, 0x27, 0xD3, 0x26, 0x83, 0x59, 0x49, 0x81, 0x5B,
+        0x0D, 0xA7, 0x0C, 0x4F, 0xEF, 0xFB, 0x1E, 0xAF, 0xE9, 0xD2,
+        0x1C, 0x10, 0x25, 0xEC, 0x9E, 0xFA, 0x57, 0x36, 0xAA, 0x3F,
+        0xC1, 0xA3, 0x2C, 0xE9, 0xB5, 0xC9, 0xED, 0x72, 0x51, 0x4C,
+        0x02, 0xB4, 0x7B, 0xB3, 0xED, 0x9F, 0x45, 0x03, 0x34, 0xAC,
+        0x9A, 0x9E, 0x62, 0x5F, 0x82, 0x7A, 0x77, 0x34, 0xF9, 0x21,
+        0x94, 0xD2, 0x38, 0x3D, 0x05, 0xF0, 0x8A, 0x60, 0x1C, 0xB7,
+        0x1D, 0xF5, 0xB7, 0x53, 0x77, 0xD3, 0x9D, 0x3D, 0x70, 0x6A,
+        0xCB, 0x18, 0x20, 0x6B, 0x29, 0x17, 0x3A, 0x6D, 0xA1, 0xB2,
+        0x64, 0xDB, 0x6C, 0xE6, 0x1A, 0x95, 0xA7, 0xF4, 0x1A, 0x78,
+        0x1D, 0xA2, 0x40, 0x15, 0x41, 0x59, 0xDD, 0xEE, 0x23, 0x57,
+        0xCE, 0x36, 0x0D, 0x55, 0xBD, 0xB8, 0xFD, 0x0F, 0x35, 0xBD,
+        0x5B, 0x92, 0xD6, 0x1C, 0x84, 0x8C, 0x32, 0x64, 0xA6, 0x5C,
+        0x45, 0x18, 0x07, 0x6B, 0xF9, 0xA9, 0x43, 0x9A, 0x83, 0xCD,
+        0xB5, 0xB3, 0xD9, 0x17, 0x99, 0x2C, 0x2A, 0x8B, 0xE0, 0x8E,
+        0xAF, 0xA6, 0x4C, 0x95, 0xBB, 0x70, 0x60, 0x1A, 0x3A, 0x97,
+        0xAA, 0x2F, 0x3D, 0x22, 0x83, 0xB7, 0x4F, 0x59, 0xED, 0x3F,
+        0x4E, 0xF4, 0x19, 0xC6, 0x25, 0x0B, 0x0A, 0x5E, 0x21, 0xB9,
+        0x91, 0xB8, 0x19, 0x84, 0x48, 0x78, 0xCE, 0x27, 0xBF, 0x41,
+        0x89, 0xF6, 0x30, 0xFD, 0x6B, 0xD9, 0xB8, 0x1D, 0x72, 0x8A,
+        0x56, 0xCC, 0x2F, 0x82, 0xE4, 0x46, 0x4D, 0x75, 0xD8, 0x92,
+        0xE6, 0x9C, 0xCC, 0xD2, 0xCD, 0x35, 0xE4, 0xFC, 0x2A, 0x85,
+        0x6B, 0xA9, 0xB2, 0x27, 0xC9, 0xA1, 0xFF, 0xB3, 0x96, 0x3E,
+        0x59, 0xF6, 0x4C, 0x66, 0x56, 0x2E, 0xF5, 0x1B, 0x97, 0x32,
+        0xB0, 0x71, 0x5A, 0x9C, 0x50, 0x4B, 0x6F, 0xC4, 0xCA, 0x94,
+        0x75, 0x37, 0x46, 0x10, 0x12, 0x2F, 0x4F, 0xA3, 0x82, 0xCD,
+        0xBD, 0x7C};
+
+    const uint8_t ss_exp[SIKEp503_SS_BYTESZ] = {
+        0x74, 0x3D, 0x25, 0x36, 0x00, 0x24, 0x63, 0x1A, 0x39, 0x1A,
+        0xB4, 0xAD, 0x01, 0x17, 0x78, 0xE9};
+
+    uint8_t ss_dec[SIKEp503_SS_BYTESZ] = {0};
+    SIKE_decaps(ss_dec, ct_exp, pk, sk);
+    EXPECT_EQ(memcmp(ss_dec, ss_exp, sizeof(ss_exp)), 0);
+}
+
+// SIKE_encaps and SIKE_keypair doesn't return zeros.
+TEST(SIKE, NonZero) {
+    uint8_t sk[SIKEp503_PRV_BYTESZ] = {0};
+    uint8_t pk[SIKEp503_PUB_BYTESZ] = {0};
+    uint8_t ct[SIKEp503_CT_BYTESZ] = {0};
+    uint8_t ss[SIKEp503_SS_BYTESZ] = {0};
+
+    // Check secret and public key returned by SIKE_keypair
+    EXPECT_EQ(SIKE_keypair(sk, pk), 1);
+    uint8_t tmp = 0;
+    for (size_t i=0; i<sizeof(sk); i++) tmp|=sk[i];
+    EXPECT_NE(tmp, 0);
+
+    tmp = 0;
+    for (size_t i=0; i<sizeof(pk); i++) tmp|=pk[i];
+    EXPECT_NE(tmp, 0);
+
+    // Check shared secret and ciphertext returned by SIKE_encaps
+    SIKE_encaps(ss, ct, pk);
+    tmp = 0;
+    for (size_t i=0; i<sizeof(ct); i++) tmp|=ct[i];
+    EXPECT_NE(tmp, 0);
+
+    tmp = 0;
+    for (size_t i=0; i<sizeof(ss); i++) tmp|=ss[i];
+    EXPECT_NE(tmp, 0);
+}
+
+TEST(SIKE, Negative) {
+    uint8_t sk[SIKEp503_PRV_BYTESZ] = {0};
+    uint8_t pk[SIKEp503_PUB_BYTESZ] = {0};
+    uint8_t ct[SIKEp503_CT_BYTESZ] = {0};
+    uint8_t ss_enc[SIKEp503_SS_BYTESZ] = {0};
+    uint8_t ss_dec[SIKEp503_SS_BYTESZ] = {0};
+
+    EXPECT_EQ(SIKE_keypair(sk, pk), 1);
+    SIKE_encaps(ss_enc, ct, pk);
+
+    // Change cipertext
+    uint8_t ct_tmp[SIKEp503_CT_BYTESZ] = {0};
+    memcpy(ct_tmp, ct, sizeof(ct));
+    ct_tmp[0] = ~ct_tmp[0];
+    SIKE_decaps(ss_dec, ct_tmp, pk, sk);
+    EXPECT_NE(memcmp(ss_enc, ss_dec, SIKEp503_SS_BYTESZ), 0);
+
+    // Change secret key
+    uint8_t sk_tmp[SIKEp503_PRV_BYTESZ] = {0};
+    memcpy(sk_tmp, sk, sizeof(sk));
+    sk_tmp[0] = ~sk_tmp[0];
+    SIKE_decaps(ss_dec, ct, pk, sk_tmp);
+    EXPECT_NE(memcmp(ss_enc, ss_dec, SIKEp503_SS_BYTESZ), 0);
+
+    // Change public key
+    uint8_t pk_tmp[SIKEp503_PUB_BYTESZ] = {0};
+    memcpy(pk_tmp, pk, sizeof(pk));
+    pk_tmp[0] = ~pk_tmp[0];
+    SIKE_decaps(ss_dec, ct, pk_tmp, sk);
+    EXPECT_NE(memcmp(ss_enc, ss_dec, SIKEp503_SS_BYTESZ), 0);
+}
+
+#if defined(SUPPORTS_ABI_TEST) && (defined(OPENSSL_X86_64) || defined(OPENSSL_AARCH64))
+TEST(SIKE, ABI) {
+  felm_t a, b, c;
+  dfelm_t d, e, f;
+  CHECK_ABI(sike_fpadd, a, b, c);
+  CHECK_ABI(sike_fpsub, a, b, c);
+  CHECK_ABI(sike_mpmul, a, b, d);
+  CHECK_ABI(sike_fprdc, d, a);
+  CHECK_ABI(sike_mpadd_asm, a, b, c);
+  CHECK_ABI(sike_mpsubx2_asm, d, e, f);
+  CHECK_ABI(sike_mpdblsubx2_asm, d, e, f);
+}
+#endif  // SUPPORTS_ABI_TEST && (X86_64 || AARCH64)
diff --git a/third_party/sike/utils.h b/third_party/sike/utils.h
new file mode 100644
index 0000000..ee9e685
--- /dev/null
+++ b/third_party/sike/utils.h
@@ -0,0 +1,144 @@
+/********************************************************************************************
+* SIDH: an efficient supersingular isogeny cryptography library
+*
+* Abstract: internal header file for P503
+*********************************************************************************************/
+
+#ifndef UTILS_H_
+#define UTILS_H_
+
+#include <stdbool.h>
+#include "openssl/base.h"
+#include "../crypto/internal.h"
+#include "sike.h"
+
+// Conversion macro from number of bits to number of bytes
+#define BITS_TO_BYTES(nbits)      (((nbits)+7)/8)
+
+// Bit size of the field
+#define BITS_FIELD             503
+// Byte size of the field
+#define FIELD_BYTESZ            BITS_TO_BYTES(BITS_FIELD)
+// Number of 64-bit words of a 503-bit field element
+#define NWORDS64_FIELD          ((BITS_FIELD+63)/64)
+// Number of 64-bit words of a 256-bit element
+#define NBITS_ORDER             256
+#define NWORDS64_ORDER          ((NBITS_ORDER+63)/64)
+// Number of elements in Alice's strategy
+#define A_max                   125
+// Number of elements in Bob's strategy
+#define B_max                   159
+// Word size size
+#define RADIX                   sizeof(crypto_word_t)*8
+// Byte size of a limb
+#define LSZ                     sizeof(crypto_word_t)
+
+#if defined(OPENSSL_64_BIT)
+    // Number of words of a 503-bit field element
+    #define NWORDS_FIELD    8
+    // Number of "0" digits in the least significant part of p503 + 1
+    #define p503_ZERO_WORDS 3
+    // log_2(RADIX)
+    #define LOG2RADIX       6
+#else
+    // Number of words of a 503-bit field element
+    #define NWORDS_FIELD    16
+    // Number of "0" digits in the least significant part of p503 + 1
+    #define p503_ZERO_WORDS 7
+    // log_2(RADIX)
+    #define LOG2RADIX       5
+#endif
+
+// Extended datatype support
+#if !defined(BORINGSSL_HAS_UINT128)
+    typedef uint64_t uint128_t[2];
+#endif
+
+// The following functions return 1 (TRUE) if condition is true, 0 (FALSE) otherwise
+// Digit multiplication
+#define MUL(multiplier, multiplicand, hi, lo) digit_x_digit((multiplier), (multiplicand), &(lo));
+
+// If mask |x|==0xff.ff set |x| to 1, otherwise 0
+#define M2B(x) ((x)>>(RADIX-1))
+
+// Digit addition with carry
+#define ADDC(carryIn, addend1, addend2, carryOut, sumOut)                   \
+do {                                                                        \
+  crypto_word_t tempReg = (addend1) + (crypto_word_t)(carryIn);             \
+  (sumOut) = (addend2) + tempReg;                                           \
+  (carryOut) = M2B(constant_time_lt_w(tempReg, (crypto_word_t)(carryIn)) |  \
+                   constant_time_lt_w((sumOut), tempReg));                  \
+} while(0)
+
+// Digit subtraction with borrow
+#define SUBC(borrowIn, minuend, subtrahend, borrowOut, differenceOut)           \
+do {                                                                            \
+    crypto_word_t tempReg = (minuend) - (subtrahend);                           \
+    crypto_word_t borrowReg = M2B(constant_time_lt_w((minuend), (subtrahend))); \
+    borrowReg |= ((borrowIn) & constant_time_is_zero_w(tempReg));               \
+    (differenceOut) = tempReg - (crypto_word_t)(borrowIn);                      \
+    (borrowOut) = borrowReg;                                                    \
+} while(0)
+
+/* Old GCC 4.9 (jessie) doesn't implement {0} initialization properly,
+   which violates C11 as described in 6.7.9, 21 (similarily C99, 6.7.8).
+   Defines below are used to work around the bug, and provide a way
+   to initialize f2elem_t and point_proj_t structs.
+   Bug has been fixed in GCC6 (debian stretch).
+*/
+#define F2ELM_INIT {{ {0}, {0} }}
+#define POINT_PROJ_INIT {{ F2ELM_INIT, F2ELM_INIT }}
+
+// Datatype for representing 503-bit field elements (512-bit max.)
+// Elements over GF(p503) are encoded in 63 octets in little endian format
+// (i.e., the least significant octet is located in the lowest memory address).
+typedef crypto_word_t felm_t[NWORDS_FIELD];
+
+// An element in F_{p^2}, is composed of two coefficients from F_p, * i.e.
+// Fp2 element = c0 + c1*i in F_{p^2}
+// Datatype for representing double-precision 2x503-bit field elements (512-bit max.)
+// Elements (a+b*i) over GF(p503^2), where a and b are defined over GF(p503), are
+// encoded as {a, b}, with a in the lowest memory portion.
+typedef struct {
+    felm_t c0;
+    felm_t c1;
+} fp2;
+
+// Our F_{p^2} element type is a pointer to the struct.
+typedef fp2 f2elm_t[1];
+
+// Datatype for representing double-precision 2x503-bit
+// field elements in contiguous memory.
+typedef crypto_word_t dfelm_t[2*NWORDS_FIELD];
+
+// Constants used during SIKEp503 computation.
+struct params_t {
+    // Stores P503 prime
+    const uint64_t prime[NWORDS64_FIELD];
+    // Stores P503 + 1
+    const uint64_t prime_p1[NWORDS64_FIELD];
+    // Stores P503 * 2
+    const uint64_t prime_x2[NWORDS64_FIELD];
+    // Alice's generator values {XPA0 + XPA1*i, XQA0, XRA0 + XRA1*i}
+    // in GF(p503^2), expressed in Montgomery representation
+    const uint64_t A_gen[5*NWORDS64_FIELD];
+    // Bob's generator values {XPB0 + XPB1*i, XQB0, XRB0 + XRB1*i}
+    // in GF(p503^2), expressed in Montgomery representation
+    const uint64_t B_gen[5*NWORDS64_FIELD];
+    // Montgomery constant mont_R2 = (2^512)^2 mod p503
+    const uint64_t mont_R2[NWORDS64_FIELD];
+    // Value 'one' in Montgomery representation
+    const uint64_t mont_one[NWORDS64_FIELD];
+    // Fixed parameters for isogeny tree computation
+    const unsigned int A_strat[A_max-1];
+    const unsigned int B_strat[B_max-1];
+};
+
+// Point representation in projective XZ Montgomery coordinates.
+typedef struct {
+    f2elm_t X;
+    f2elm_t Z;
+} point_proj;
+typedef point_proj point_proj_t[1];
+
+#endif // UTILS_H_
diff --git a/tool/speed.cc b/tool/speed.cc
index b26470e..47edc75 100644
--- a/tool/speed.cc
+++ b/tool/speed.cc
@@ -51,6 +51,8 @@
 #include "../crypto/internal.h"
 #include "internal.h"
 
+#include "../third_party/sike/sike.h"
+
 
 // TimeResults represents the results of benchmarking a function.
 struct TimeResults {
@@ -294,6 +296,64 @@
   return true;
 }
 
+static bool SpeedSIKEP503(const std::string &selected) {
+  if (!selected.empty() && selected.find("SIKE") == std::string::npos) {
+    return true;
+  }
+  // speed generation
+  uint8_t public_SIKE[SIKEp503_PUB_BYTESZ];
+  uint8_t private_SIKE[SIKEp503_PRV_BYTESZ];
+  uint8_t ct[SIKEp503_CT_BYTESZ];
+  bool res;
+
+  {
+    TimeResults results;
+    res = TimeFunction(&results,
+                [&private_SIKE, &public_SIKE]() -> bool {
+      return (SIKE_keypair(private_SIKE, public_SIKE) == 1);
+    });
+    results.Print("SIKE/P503 generate");
+  }
+
+  if (!res) {
+    fprintf(stderr, "Failed to time SIKE_keypair.\n");
+    return false;
+  }
+
+  {
+    TimeResults results;
+    TimeFunction(&results,
+                [&ct, &public_SIKE]() -> bool {
+      uint8_t ss[SIKEp503_SS_BYTESZ];
+      SIKE_encaps(ss, ct, public_SIKE);
+      return true;
+    });
+    results.Print("SIKE/P503 encap");
+  }
+
+  if (!res) {
+    fprintf(stderr, "Failed to time SIKE_encaps.\n");
+    return false;
+  }
+
+  {
+    TimeResults results;
+    TimeFunction(&results,
+                [&ct, &public_SIKE, &private_SIKE]() -> bool {
+      uint8_t ss[SIKEp503_SS_BYTESZ];
+      SIKE_decaps(ss, ct, public_SIKE, private_SIKE);
+      return true;
+    });
+    results.Print("SIKE/P503 decap");
+  }
+
+  if (!res) {
+    fprintf(stderr, "Failed to time SIKE_decaps.\n");
+    return false;
+  }
+  return true;
+}
+
 static uint8_t *align(uint8_t *in, unsigned alignment) {
   return reinterpret_cast<uint8_t *>(
       (reinterpret_cast<uintptr_t>(in) + alignment) &
@@ -938,6 +998,7 @@
       !SpeedECDH(selected) ||
       !SpeedECDSA(selected) ||
       !Speed25519(selected) ||
+      !SpeedSIKEP503(selected) ||
       !SpeedSPAKE2(selected) ||
       !SpeedScrypt(selected) ||
       !SpeedRSAKeyGen(selected) ||