curve25519-donna-x86-64.s

# 2008, Google Inc.
# All rights reserved.
#
# Code released into the public domain

################################################################################
# curve25519-donna.s - an x86-64 bit implementation of curve25519. See the
# comments at the top of curve25519-donna.c
#
# Adam Langley <agl@imperialviolet.org>
#
# Derived from public domain C code by Daniel J. Bernstein <djb@cr.yp.to>
#
# More information about curve25519 can be found here
#   http://cr.yp.to/ecdh.html
################################################################################
.text

.extern fmonty

.globl fmul
.globl fsquare
.globl fexpand
.globl fcontract
.globl freduce_coefficients
.globl fscalar
.globl fdifference_backwards
.globl cmult

################################################################################
# fmul - multiply two 256-bit numbers
# Registers: RDI (output): uint64_t[5] product
#            RSI (input): uint64_t[5] input 1
#            RDX (input): uint64_t[5] input 2
################################################################################
fmul:

# Input pointers: rdi (output), rsi (in1), rdx (in2)
# Spill: rdi, rbx, r12..15

push %rbx
push %r12
push %r13
push %r14
push %r15
push %rdi

# Load 5 64-bit values from *rsi to rsi, r8..11

mov %rsi,%rcx
mov (%rcx),%rsi
mov 8(%rcx),%r8
mov 16(%rcx),%r9
mov 24(%rcx),%r10
mov 32(%rcx),%r11

# Load 5 64-bit values from *rdx to rdi, r12..15

mov (%rdx),%rdi
mov 8(%rdx),%r12
mov 16(%rdx),%r13
mov 24(%rdx),%r14
mov 32(%rdx),%r15

# We are going to perform a polynomial multiplication of two, five element
# polynomials. I and J and the polynomials and I2 would be the coefficient of
# x^2 etc.

#define I0 %rsi
#define I1 %r8
#define I2 %r9
#define I3 %r10
#define I4 %r11

#define J0 %rdi
#define J1 %r12
#define J2 %r13
#define J3 %r14
#define J4 %r15

#define MUL(a, b) \
mov a,%rax ; \
mul b

# We accumululate results in RCX:RBX

#define FIRSTMUL(a,b) \
MUL(a,b) ; \
mov %rax,%rbx ; \
mov %rdx,%rcx

#define MAC(a,b) \
MUL (a,b) ; \
add %rax,%rbx ; \
adc %rdx,%rcx

# p[0] = i[0] * j[0]
# p[0] stored in xmm0, xmm1

MUL(I0,J0)
movq %rax,%xmm0
movq %rdx,%xmm1

# p[1] = i[0] * j[1] + i[1] * j[0]

FIRSTMUL(I0,J1)
MAC(I1, J0)

movq %rbx,%xmm2
movq %rcx,%xmm3

# p[2] = i[1] * j[1] + i[0] * j[2] + i[2] * j[0]

FIRSTMUL(I1,J1)
MAC(I0,J2)
MAC(I2,J0)

movq %rbx,%xmm4
movq %rcx,%xmm5

# p[3]

FIRSTMUL(I0,J3)
MAC(I3,J0)
MAC(I1,J2)
MAC(I2,J1)

movq %rbx,%xmm6
movq %rcx,%xmm7

# p[4]

FIRSTMUL(I0,J4)
MAC(I4,J0)
MAC(I3,J1)
MAC(I1,J3)
MAC(I2,J2)

movq %rbx,%xmm8
movq %rcx,%xmm9

# p[5]

FIRSTMUL(I4,J1)
MAC(I1,J4)
MAC(I2,J3)
MAC(I3,J2)

movq %rbx,%xmm10
movq %rcx,%xmm11

# p[6]

FIRSTMUL(I4,J2)
MAC(I2,J4)
MAC(I3,J3)

movq %rbx,%xmm12
movq %rcx,%xmm13

# p[7]

FIRSTMUL(I3,J4)
MAC(I4,J3)

movq %rbx,%xmm14
movq %rcx,%xmm15

# p[8], keeping it in RDX:RAX

MUL(I4,J4)

donna_reduce:
# We done with the original inputs now, so we start reusing them
# At this point we have a degree 8 resulting polynomial and we need to reduce
# mod 2**255-19. Since 2**255 is in our polynomial, we can multiply the
# coefficients of the higher powers and add them to the lower powers. The limb
# size (51-bits) is chosen to avoid overflows.

mov $19,%r15

# p[8] *= 19, store in R13:R12

mov %rdx,%r13
mul %r15
imul %r15,%r13
add %rdx,%r13
mov %rax,%r12

# p[3] += p[8] * 19

movq %xmm7,%rcx
movq %xmm6,%rbx
add %rbx,%r12
adc %rcx,%r13

#define MUL19ADD(srchi,srclo,desthi,destlo,src2hi,src2lo) \
movq %srclo,%rax ; \
mul %r15 ; \
movq %srchi,%desthi ; \
imul %r15,%desthi ; \
add %rdx,%desthi ; \
mov %rax,%destlo ; \
movq %src2hi,%rcx ; \
movq %src2lo,%rbx ; \
add %rbx,%destlo ; \
adc %rcx,%desthi

# p[2] += p[7] * 19, store in R11:R10

MUL19ADD(xmm15,xmm14,r11,r10,xmm5,xmm4)

# p[1] += p[6] * 19, store in R9:R8

MUL19ADD(xmm13,xmm12,r9,r8,xmm3,xmm2)

# p[0] += p[5] * 19, store in RDI:RSI

MUL19ADD(xmm11,xmm10,rdi,rsi,xmm1,xmm0)

# p[4], store in R15:R14

movq %xmm9,%r15
movq %xmm8,%r14

# Coefficient reduction

# Bottom 51-bits set
mov $0x7ffffffffffff,%rbx
mov $19,%rcx

coeffreduction:

# The carry chain takes the excess bits from a 128-bit result (excess are
# anything over 51-bits and above) and adds them to the next value. If the top
# value spills over, we reduce mod 2**255-19 again by multipling by 19 and
# adding onto the bottom.

#define CARRYCHAIN(srchi,srclo,desthi,destlo) \
mov %srclo,%rax ; \
shr $51,%srclo ; \
shl $13,%srchi ; \
or %srclo,%srchi ; \
add %srchi,%destlo ; \
adc $0,%desthi ; \
xor %srchi,%srchi ; \
mov %rax,%srclo ; \
and %rbx,%srclo

#define CARRYCHAIN19(srchi,srclo,desthi,destlo) \
mov %srclo,%rax ; \
shr $51,%srclo ; \
shl $13,%srchi ; \
or %srclo,%srchi ; \
imul $19,%srchi ; \
add %srchi,%destlo ; \
adc $0,%desthi ; \
xor %srchi,%srchi ; \
mov %rax,%srclo ; \
and %rbx,%srclo

CARRYCHAIN(rdi,rsi,r9,r8)
CARRYCHAIN(r9,r8,r11,r10)
CARRYCHAIN(r11,r10,r13,r12)
CARRYCHAIN(r13,r12,r15,r14)
CARRYCHAIN19(r15,r14,rdi,rsi)
CARRYCHAIN(rdi,rsi,r9,r8)

# write out results, which are in rsi, r8, r10, r12, rax
# output pointer is on top of the stack

pop %rdi

mov %rsi,(%rdi)
mov %r8,8(%rdi)
mov %r10,16(%rdi)
mov %r12,24(%rdi)
mov %r14,32(%rdi)

pop %r15
pop %r14
pop %r13
pop %r12
pop %rbx

ret

################################################################################
# fsquare - square a 256-bit number
# Registers: RDI (output): uint64_t[5] product
#            RSI (input): uint64_t[5] input
#
# This is very similar to fmul, above, however when squaring a number we can
# save some multiplications and replace them with doublings.
################################################################################

fsquare:

push %rbx
push %r12
push %r13
push %r14
push %r15
push %rdi

# Load 5 64-bit values from *rsi to rsi, r8..11

mov %rsi,%rcx
mov (%rcx),%rsi
mov 8(%rcx),%r8
mov 16(%rcx),%r9
mov 24(%rcx),%r10
mov 32(%rcx),%r11

# p[0] = i[0] * j[0]
# p[0] stored in xmm0, xmm1

MUL(I0,I0)
movq %rax,%xmm0
movq %rdx,%xmm1

# p[1] = i[0] * j[1] + i[1] * j[0]

MUL(I0,I1)
sal $1,%rax
rcl $1,%rdx

movq %rax,%xmm2
movq %rdx,%xmm3

# p[2] = i[1] * j[1] + i[0] * j[2] + i[2] * j[0]

#define FIRSTMULDOUBLE(x,y) \
FIRSTMUL(x,y) ; \
sal $1,%rbx ; \
rcl $1,%rcx

#define DOUBLEMAC(x,y) \
MUL(x,y) ; \
sal $1,%rax ; \
rcl $1,%rdx ; \
add %rax,%rbx ; \
adc %rdx,%rcx

FIRSTMUL(I1,I1)
DOUBLEMAC(I0,I2)

movq %rbx,%xmm4
movq %rcx,%xmm5

# p[3]

FIRSTMULDOUBLE(I0,I3)
DOUBLEMAC(I1,I2)

movq %rbx,%xmm6
movq %rcx,%xmm7

# p[4]

FIRSTMULDOUBLE(I0,I4)
DOUBLEMAC(I3,I1)
MAC(I2,I2)

movq %rbx,%xmm8
movq %rcx,%xmm9

# p[5]

FIRSTMULDOUBLE(I4,I1)
DOUBLEMAC(I2,I3)

movq %rbx,%xmm10
movq %rcx,%xmm11

# p[6]

FIRSTMULDOUBLE(I4,I2)
MAC(I3,I3)

movq %rbx,%xmm12
movq %rcx,%xmm13

# p[7]

MUL(I3,I4)
sal $1,%rax
rcl $1,%rdx

movq %rax,%xmm14
movq %rdx,%xmm15

# p[8], keeping it in RDX:RAX

MUL(I4,I4)

jmp donna_reduce

################################################################################
# fdifference_backwards - set output to in - output (note order)
#
# /*
# static const felem twotothe51 = (1ul << 51);
#
# static void fdifference_backwards(felem *ooutput, const felem *oin) {
#   const int64_t *in = (const int64_t *) oin;
#   int64_t *out = (int64_t *) ooutput;
#
#   out[0] = in[0] - out[0];
#   out[1] = in[1] - out[1];
#   out[2] = in[2] - out[2];
#   out[3] = in[3] - out[3];
#   out[4] = in[4] - out[4];
#
#   do {
#     if (out[0] < 0) {
#       out[0] += twotothe51;
#       out[1]--;
#     }
#     if (out[1] < 0) {
#       out[1] += twotothe51;
#       out[2]--;
#     }
#     if (out[2] < 0) {
#       out[2] += twotothe51;
#       out[3]--;
#     }
#     if (out[3] < 0) {
#       out[3] += twotothe51;
#       out[4]--;
#     }
#     if (out[4] < 0) {
#       out[4] += twotothe51;
#       out[0] -= 19;
#     }
#   } while (out[0] < 0);
# }
# */
################################################################################

fdifference_backwards:

mov (%rsi),%rax
mov 8(%rsi),%r8
mov 16(%rsi),%r9
mov 24(%rsi),%r10
mov 32(%rsi),%r11

sub (%rdi),%rax
sub 8(%rdi),%r8
sub 16(%rdi),%r9
sub 24(%rdi),%r10
sub 32(%rdi),%r11

# 2**51
mov $0x8000000000000,%rdx

fdifference_backwards_loop:

# In the C code, above, we have lots of branches. We replace these branches
# with a trick. An arithmetic shift right of 63-bits turns a positive number to
# 0, but a negative number turns to all ones. This gives us a bit-mask that we
# can AND against to add 2**51, conditionally.

#define NEGCHAIN(src,dest) \
mov %src,%rcx ; \
sar $63,%rcx ; \
and %rdx,%rcx ; \
add %rcx,%src ; \
shr $51,%rcx ; \
sub %rcx,%dest

#define NEGCHAIN19(src,dest) \
mov %src,%rcx ; \
sar $63,%rcx ; \
mov %rcx,%rsi ; \
and %rdx,%rcx ; \
add %rcx,%src ; \
and $19,%rsi ; \
sub %rsi,%dest

NEGCHAIN(rax,r8)
NEGCHAIN(r8,r9)
NEGCHAIN(r9,r10)
NEGCHAIN(r10,r11)
NEGCHAIN19(r11,rax)
NEGCHAIN(rax,r8)
NEGCHAIN(r8,r9)
NEGCHAIN(r9,r10)
NEGCHAIN(r10,r11)

mov %rax,(%rdi)
mov %r8,8(%rdi)
mov %r9,16(%rdi)
mov %r10,24(%rdi)
mov %r11,32(%rdi)

ret


################################################################################
# fscalar - multiply by 121665
#
# Registers: RDI: (out) pointer to uint64_t[5]
#            RSI: (in) pointer to uint64_t[5]
#
# Since we only have 13-bits of space at the top of our limbs, this is a full,
# cascading multiplication.
################################################################################
fscalar:

mov $121665,%rcx

mov (%rsi),%rax
mul %rcx
shl $13,%rdx
mov %rdx,%r8
mov %rax,%r9

mov 8(%rsi),%rax
mul %rcx
add %r8,%rax
shl $13,%rdx
mov %rdx,%r8
mov %rax,8(%rdi)

mov 16(%rsi),%rax
mul %rcx
add %r8,%rax
shl $13,%rdx
mov %rdx,%r8
mov %rax,16(%rdi)

mov 24(%rsi),%rax
mul %rcx
add %r8,%rax
shl $13,%rdx
mov %rdx,%r8
mov %rax,24(%rdi)

mov 32(%rsi),%rax
mul %rcx
add %r8,%rax
mov %rax,32(%rdi)
shl $13,%rdx
mov $19,%rcx
mov %rdx,%rax
mul %rcx
add %rax,%r9
mov %r9,0(%rdi)

ret

################################################################################
# freduce_coefficients
#
# Registers: RDI: (in/out) pointer to uint64_t[5]
################################################################################
freduce_coefficients:

push %r12

mov $0x7ffffffffffff,%rcx
mov $19,%rdx

mov (%rdi),%r8
mov 8(%rdi),%r9
mov 16(%rdi),%r10
mov 24(%rdi),%r11
mov 32(%rdi),%r12

carrychain_:

mov %r8,%rax
shr $51,%rax
add %rax,%r9
and %rcx,%r8

mov %r9,%rax
shr $51,%rax
add %rax,%r10
and %rcx,%r9

mov %r10,%rax
shr $51,%rax
add %rax,%r11
and %rcx,%r10

mov %r11,%rax
shr $51,%rax
add %rax,%r12
and %rcx,%r11

mov %r12,%rax
shr $51,%rax
imul $19,%rax
add %rax,%r8
and %rcx,%r12

mov %r8,(%rdi)
mov %r9,8(%rdi)
mov %r10,16(%rdi)
mov %r11,24(%rdi)
mov %r12,32(%rdi)

pop %r12
ret


################################################################################
# fexpand - convert a packed (32 byte) representation to 5 uint64_t's
# Registers: RDI: (output) pointer to uint64_t[5]
#            RSI: (input) pointer to uint8_t[32]
################################################################################
fexpand:

mov $0x7ffffffffffff,%rdx

mov (%rsi),%rax
and %rdx,%rax
mov %rax,(%rdi)

mov 6(%rsi),%rax
shr $3,%rax
and %rdx,%rax
mov %rax,8(%rdi)

mov 12(%rsi),%rax
shr $6,%rax
and %rdx,%rax
mov %rax,16(%rdi)

mov 19(%rsi),%rax
shr $1,%rax
and %rdx,%rax
mov %rax,24(%rdi)

mov 25(%rsi),%rax
shr $4,%rax
and %rdx,%rax
mov %rax,32(%rdi)

ret

################################################################################
# fcontract - convert 5 uint64_t's to a packed (32 byte) representation
#
# Registers: RDI: (output) pointer to uint8_t[32]
#            RSI: (input) pointer to uint64_t[5]
################################################################################
fcontract:

mov (%rsi),%rax
mov 8(%rsi),%rdx
mov 16(%rsi),%r8
mov 24(%rsi),%r9
mov 32(%rsi),%r10

mov %rdx,%rcx
shl $51,%rcx
or %rcx,%rax
mov %rax,(%rdi)

shr $13,%rdx
mov %r8,%rcx
shl $38,%rcx
or %rcx,%rdx
mov %rdx,8(%rdi)

shr $26,%r8
mov %r9,%rcx
shl $25,%rcx
or %rcx,%r8
mov %r8,16(%rdi)

shr $39,%r9
shl $12,%r10
or %r10,%r9
mov %r9,24(%rdi)

ret

################################################################################
# cmult - calculates nQ wher Q is the x-coordinate of a point on the curve
#
# Registers: RDI: (output) final x
#            RSI: (output) final z
#            RDX: (input) n (big-endian)
#            RCX: (input) q (big-endian)
#
 /* Calculates nQ where Q is the x-coordinate of a point on the curve
  *
  *   resultx/resultz: the x coordinate of the resulting curve point (short form)
  *   n: a little endian, 32-byte number
  *   q: a point of the curve (short form)
  */
/*
 static void
 cmult(felem *resultx, felem *resultz, const u8 *n, const felem *q) {
   felem a[5] = {0}, b[5] = {1}, c[5] = {1}, d[5] = {0};
   felem *nqpqx = a, *nqpqz = b, *nqx = c, *nqz = d, *t;
   felem e[5] = {0}, f[5] = {1}, g[5] = {0}, h[5] = {1};
   felem *nqpqx2 = e, *nqpqz2 = f, *nqx2 = g, *nqz2 = h;

   unsigned i, j;

   memcpy(nqpqx, q, sizeof(felem) * 5);

   for (i = 0; i < 32; ++i) {
     u8 byte = n[31 - i];
     for (j = 0; j < 8; ++j) {
       if (byte & 0x80) {
         fmonty(nqpqx2, nqpqz2,
                nqx2, nqz2,
                nqpqx, nqpqz,
                nqx, nqz,
                q);
       } else {
         fmonty(nqx2, nqz2,
                nqpqx2, nqpqz2,
                nqx, nqz,
                nqpqx, nqpqz,
                q);
       }

       t = nqx;
       nqx = nqx2;
       nqx2 = t;
       t = nqz;
       nqz = nqz2;
       nqz2 = t;
       t = nqpqx;
       nqpqx = nqpqx2;
       nqpqx2 = t;
       t = nqpqz;
       nqpqz = nqpqz2;
       nqpqz2 = t;

       byte <<= 1;
     }
   }

   memcpy(resultx, nqx, sizeof(felem) * 5);
   memcpy(resultz, nqz, sizeof(felem) * 5);
 }
*/
################################################################################
cmult:

push %rbp
push %r13
push %r14

mov %rsp,%rbp
mov $63,%r8
not %r8
and %r8,%rsp

mov %rdx,%r13
mov %rcx,%r14

sub $512,%rsp

# value nQ+Q (x)
movq (%rcx),%rax
movq %rax,(%rsp)
movq 8(%rcx),%r8
movq %r8,8(%rsp)
movq 16(%rcx),%r9
movq %r9,16(%rsp)
movq 24(%rcx),%r10
movq %r10,24(%rsp)
movq 32(%rcx),%r11
movq %r11,32(%rsp)

# value nQ+Q (z)
movq $1,64(%rsp)
movq $0,72(%rsp)
movq $0,80(%rsp)
movq $0,88(%rsp)
movq $0,96(%rsp)

# value nQ (x)
movq $1,128(%rsp)
movq $0,136(%rsp)
movq $0,144(%rsp)
movq $0,152(%rsp)
movq $0,160(%rsp)

# value nQ (z)
movq $0,192(%rsp)
movq $0,200(%rsp)
movq $0,208(%rsp)
movq $0,216(%rsp)
movq $0,224(%rsp)

push %rbx
push %r12
push %r15
push %rdi
push %rsi

# The stack looks like
# (nQ)'
# (nQ+Q)'
# nQ
# nQ+Q
# saved registers  (40-bytes)   <-- %rsp
# We switch between the two banks with an offset in %r12, starting by writing
# into the prime bank and reading from the non-prime bank.
# Based on the current MSB of the operand, we flip the two values over based
# on an offset in %r8 for the first first member and %r9 for the second

mov $256,%r12
mov $32,%rbx

cmult_loop_outer:

# On entry to the loop, the word offset is kept in %rbx. We dec 8 bytes and
# then store the outer loop counter in the top 32-bits of %rbx. The inner loop
# counter is kept in %ebx

sub $8,%rbx
movq (%r13,%rbx),%r15
shl $32,%rbx

or $64,%rbx

cmult_loop_inner:

# Register allocation:
#  r8: the element switch offset
#  r9: complement r8
#  r11: complement r12
# Preserved by fmonty:
#  rbx: loop counters
#  r12: bank switch offset
#  r13: (input) n
#  r14: (input) q
#  r15: the current qword, getting left shifted

mov $128,%r8
xor %r9,%r9
bt $63,%r15
cmovc %r8,%r9
mov %r9,%r8
xor $128,%r8

shl $1,%r15

mov %r12,%r11
xor $256,%r11

lea 40(%rsp,%r12),%rdi
mov %rdi,%rsi
lea 40(%rsp,%r11),%rdx
mov %rdx,%rcx
add %r8,%rdi
add %r9,%rsi
add %r8,%rdx
add %r9,%rcx
mov %r14,%r8
call fmonty

xor $256,%r12

dec %rbx
cmp $0,%ebx
jnz cmult_loop_inner

shr $32,%rbx
cmp $0,%rbx
jnz cmult_loop_outer

pop %rsi
pop %rdi
pop %r15
pop %r12
pop %rbx

lea 128(%rsp),%r8

movq (%r8),%rax
movq %rax,(%rdi)
movq 8(%r8),%rax
movq %rax,8(%rdi)
movq 16(%r8),%rax
movq %rax,16(%rdi)
movq 24(%r8),%rax
movq %rax,24(%rdi)
movq 32(%r8),%rax
movq %rax,32(%rdi)

movq 64(%r8),%rax
movq %rax,(%rsi)
movq 72(%r8),%rax
movq %rax,8(%rsi)
movq 80(%r8),%rax
movq %rax,16(%rsi)
movq 88(%r8),%rax
movq %rax,24(%rsi)
movq 96(%r8),%rax
movq %rax,32(%rsi)

mov %rbp,%rsp
pop %r14
pop %r13
pop %rbp

ret