--- /dev/null
+/// -*- mode: asm; asm-comment-char: ?/ -*-
+///
+/// GCM acceleration for x86 processors
+///
+/// (c) 2018 Straylight/Edgeware
+///
+
+///----- Licensing notice ---------------------------------------------------
+///
+/// This file is part of Catacomb.
+///
+/// Catacomb is free software: you can redistribute it and/or modify it
+/// under the terms of the GNU Library General Public License as published
+/// by the Free Software Foundation; either version 2 of the License, or
+/// (at your option) any later version.
+///
+/// Catacomb is distributed in the hope that it will be useful, but
+/// WITHOUT ANY WARRANTY; without even the implied warranty of
+/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+/// Library General Public License for more details.
+///
+/// You should have received a copy of the GNU Library General Public
+/// License along with Catacomb. If not, write to the Free Software
+/// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+/// USA.
+
+///--------------------------------------------------------------------------
+/// Preliminaries.
+
+#include "config.h"
+#include "asm-common.h"
+
+ .arch .pclmul
+
+ .text
+
+///--------------------------------------------------------------------------
+/// Common register allocation.
+
+#if CPUFAM_X86
+# define A eax
+# define K edx
+#elif CPUFAM_AMD64 && ABI_SYSV
+# define A rdi
+# define K rsi
+#elif CPUFAM_AMD64 && ABI_WIN
+# define A rcx
+# define K rdx
+#endif
+
+///--------------------------------------------------------------------------
+/// Multiplication macros.
+
+ // The good news is that we have a fancy instruction to do the
+ // multiplications. The bad news is that it's not particularly well-
+ // suited to the job.
+ //
+ // For one thing, it only does a 64-bit multiplication, so in general
+ // we'll need to synthesize the full-width multiply by hand. For
+ // another thing, it doesn't help with the reduction, so we have to
+ // do that by hand too. And, finally, GCM has crazy bit ordering,
+ // and the instruction does nothing useful for that at all.
+ //
+ // Focusing on that last problem first: the bits aren't in monotonic
+ // significance order unless we permute them. If we reverse the byte
+ // order, then we'll have the bits in monotonic order, but backwards,
+ // so the degree-0 coefficient will be in the most-significant bit.
+ //
+ // This is less of a difficulty than it seems at first, because
+ // algebra. Suppose we are given u = SUM_{0<=i<n} u_i t^i and v =
+ // SUM_{0<=j<n} v_j t^j; then
+ //
+ // u v = SUM_{0<=i,j<n} u_i v_j t^{i+j}
+ //
+ // Suppose instead that we're given ũ = SUM_{0<=i<n} u_{n-i-1} t^i
+ // and ṽ = SUM_{0<=j<n} v_{n-j-1} t^j, so the bits are backwards.
+ // Then
+ //
+ // ũ ṽ = SUM_{0<=i,j<n} u_{n-i-1} v_{n-j-1} t^{i+j}
+ // = SUM_{0<=i,j<n} u_i v_j t^{2n-2-(i+j)}
+ //
+ // which is almost the bit-reversal of u v, only it's shifted right
+ // by one place. Oh, well: we'll have to shift it back later.
+ //
+ // That was important to think about, but there's not a great deal to
+ // do about it yet other than to convert what we've got from the
+ // blockcipher's byte-ordering convention to our big-endian
+ // convention. Since this depends on the blockcipher convention,
+ // we'll leave the caller to cope with this: the macros here will
+ // assume that the operands are in `register' format, which is the
+ // byte-reversal of the external representation, padded at the
+ // most-significant end except for 96-bit blocks, which are
+ // zero-padded at the least-significant end (see `mul96' for the
+ // details). In the commentary, pieces of polynomial are numbered
+ // according to the degree of the coefficients, so the unit
+ // coefficient of some polynomial a is in a_0.
+ //
+ // The commentary for `mul128' is the most detailed. The other
+ // macros assume that you've already read and understood that.
+
+.macro mul128
+ // Enter with u and v in xmm0 and xmm1 respectively; leave with z =
+ // u v in xmm0. Clobbers xmm1--xmm4.
+
+ // First for the double-precision multiplication. It's tempting to
+ // use Karatsuba's identity here, but I suspect that loses more in
+ // the shifting, bit-twiddling, and dependency chains that it gains
+ // in saving a multiplication which otherwise pipelines well.
+ // xmm0 = // (u_1; u_0)
+ // xmm1 = // (v_1; v_0)
+ movdqa xmm2, xmm1 // (v_1; v_0) again
+ movdqa xmm3, xmm0 // (u_1; u_0) again
+ movdqa xmm4, xmm0 // (u_1; u_0) yet again
+ pclmulhqlqdq xmm2, xmm0 // u_1 v_0
+ pclmullqlqdq xmm0, xmm1 // u_1 v_1
+ pclmulhqlqdq xmm3, xmm1 // u_0 v_1
+ pclmulhqhqdq xmm4, xmm1 // u_0 v_0
+
+ // Arrange the pieces to form a double-precision polynomial.
+ pxor xmm2, xmm3 // (m_1; m_0) = u_1 v_0 + u_0 v_1
+ movdqa xmm1, xmm2 // (m_1; m_0) again
+ pslldq xmm2, 8 // (0; m_1)
+ psrldq xmm1, 8 // (m_0; 0)
+ pxor xmm0, xmm2 // x_1 = u_1 v_1 + m_1
+ pxor xmm1, xmm4 // x_0 = u_0 v_0 + t^64 m_0
+
+ // Two problems remain. The first is that this product is shifted
+ // left (from GCM's backwards perspective) by one place, which is
+ // annoying. Let's take care of that now. Once this is done, we'll
+ // be properly in GCM's backwards bit-ordering, so xmm1 will hold the
+ // low half of the product and xmm0 the high half. (The following
+ // diagrams show bit 0 consistently on the right.)
+ //
+ // xmm1
+ // ,-------------.-------------.-------------.-------------.
+ // | 0 x_0-x_30 | x_31-x_62 | x_63-x_94 | x_95-x_126 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // xmm0
+ // ,-------------.-------------.-------------.-------------.
+ // | x_127-x_158 | x_159-x_190 | x_191-x_222 | x_223-x_254 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // We start by shifting each 32-bit lane right (from GCM's point of
+ // view -- physically, left) by one place, which gives us this:
+ //
+ // low (xmm3)
+ // ,-------------.-------------.-------------.-------------.
+ // | x_0-x_30 0 | x_32-x_62 0 | x_64-x_94 0 | x_96-x_126 0|
+ // `-------------^-------------^-------------^-------------'
+ //
+ // high (xmm2)
+ // ,-------------.-------------.-------------.-------------.
+ // |x_128-x_158 0|x_160-x_190 0|x_192-x_222 0|x_224-x_254 0|
+ // `-------------^-------------^-------------^-------------'
+ //
+ // but we've lost a bunch of bits. We separately shift each lane
+ // left by 31 places to give us the bits we lost.
+ //
+ // low (xmm1)
+ // ,-------------.-------------.-------------.-------------.
+ // | 0...0 | 0...0 x_31 | 0...0 x_63 | 0...0 x_95 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // high (xmm0)
+ // ,-------------.-------------.-------------.-------------.
+ // | 0...0 x_127 | 0...0 x_159 | 0...0 x_191 | 0...0 x_223 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // Which is close, but we don't get a cigar yet. To get the missing
+ // bits into position, we shift each of these right by a lane, but,
+ // alas, the x_127 falls off, so, separately, we shift the high
+ // register left by three lanes, so that everything is lined up
+ // properly when we OR them all together:
+ //
+ // low (xmm1)
+ // ,-------------.-------------.-------------.-------------.
+ // ? 0...0 x_31 | 0...0 x_63 | 0...0 x_95 | 0...0 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // wrap (xmm4)
+ // ,-------------.-------------.-------------.-------------.
+ // | 0...0 | 0...0 | 0...0 | 0...0 x_127 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // high (xmm0)
+ // ,-------------.-------------.-------------.-------------.
+ // | 0...0 x_159 | 0...0 x_191 | 0...0 x_223 | 0...0 |
+ // `-------------^-------------^-------------^-------------'
+ //
+ // The `low' and `wrap' registers (xmm1, xmm3, xmm4) then collect the
+ // low 128 coefficients, while the `high' registers (xmm0, xmm2)
+ // collect the high 127 registers, leaving a zero bit at the most
+ // significant end as we expect.
+
+ // xmm0 = // (x_7, x_6; x_5, x_4)
+ // xmm1 = // (x_3, x_2; x_1, x_0)
+ movdqa xmm3, xmm1 // (x_3, x_2; x_1, x_0) again
+ movdqa xmm2, xmm0 // (x_7, x_6; x_5, x_4) again
+ psrld xmm1, 31 // shifted left; just the carries
+ psrld xmm0, 31
+ pslld xmm3, 1 // shifted right, but dropped carries
+ pslld xmm2, 1
+ movdqa xmm4, xmm0 // another copy for the carry around
+ pslldq xmm1, 4 // move carries over
+ pslldq xmm0, 4
+ psrldq xmm4, 12 // the big carry wraps around
+ por xmm1, xmm3
+ por xmm0, xmm2 // (y_7, y_6; y_5, y_4)
+ por xmm1, xmm4 // (y_3, y_2; y_1, y_0)
+
+ // And the other problem is that the result needs to be reduced
+ // modulo p(t) = t^128 + t^7 + t^2 + t + 1. Let R = t^128 = t^7 +
+ // t^2 + t + 1 in our field. So far, we've calculated z_0 and z_1
+ // such that z_0 + z_1 R = u v using the identity R = t^128: now we
+ // must collapse the two halves of z together using the other
+ // identity R = t^7 + t^2 + t + 1.
+ //
+ // We do this by working on each 32-bit word of the high half of z
+ // separately, so consider y_i, for some 4 <= i < 8. Certainly, y_i
+ // t^{32i} = y_i R t^{32(i-4)} = (t^7 + t^2 + t + 1) y_i t^{32(i-4)},
+ // but we can't use that directly without breaking up the 32-bit word
+ // structure. Instead, we start by considering just y_i t^7
+ // t^{32(i-4)}, which again looks tricky. Now, split y_i = a_i +
+ // t^25 b_i, with deg a_i < 25; then
+ //
+ // y_i t^7 t^{32(i-4)} = a_i t^7 t^{32(i-4)} + b_i t^{32(i-3)}
+ //
+ // We can similarly decompose y_i t^2 and y_i t into a pair of 32-bit
+ // contributions to the t^{32(i-4)} and t^{32(i-3)} words, but the
+ // splits are different. This is lovely, with one small snag: when
+ // we do this to y_7, we end up with a contribution back into the
+ // t^128 coefficient word. But notice that only the low seven bits
+ // of this word are affected, so there's no knock-on contribution
+ // into the t^32 word. Therefore, if we handle the high bits of each
+ // word together, and then the low bits, everything will be fine.
+
+ // First, shift the high bits down.
+ movdqa xmm2, xmm0 // (y_7, y_6; y_5, y_4) again
+ movdqa xmm3, xmm0 // (y_7, y_6; y_5, y_4) yet again
+ movdqa xmm4, xmm0 // (y_7, y_6; y_5, y_4) again again
+ pslld xmm2, 31 // the b_i for t
+ pslld xmm3, 30 // the b_i for t^2
+ pslld xmm4, 25 // the b_i for t^7
+ pxor xmm2, xmm3 // add them all together
+ pxor xmm2, xmm4
+ movdqa xmm3, xmm2 // and a copy for later
+ psrldq xmm2, 4 // contribution into low half
+ pslldq xmm3, 12 // and high half
+ pxor xmm1, xmm2
+ pxor xmm0, xmm3
+
+ // And then shift the low bits up.
+ movdqa xmm2, xmm0
+ movdqa xmm3, xmm0
+ pxor xmm1, xmm0 // mix in the unit contribution
+ psrld xmm0, 1
+ psrld xmm2, 2
+ psrld xmm3, 7
+ pxor xmm1, xmm2 // low half, unit, and t^2 contribs
+ pxor xmm0, xmm3 // t and t^7 contribs
+ pxor xmm0, xmm1 // mix them together and we're done
+.endm
+
+.macro mul64
+ // Enter with u and v in the low halves of xmm0 and xmm1
+ // respectively; leave with z = u v in xmm0. Clobbers xmm1--xmm4.
+
+ // The multiplication is thankfully easy.
+ pclmullqlqdq xmm0, xmm1 // u v
+
+ // Shift the product up by one place. After this, we're in GCM
+ // bizarro-world.
+ movdqa xmm1, xmm0 // u v again
+ psrld xmm0, 31 // shifted left; just the carries
+ pslld xmm1, 1 // shifted right, but dropped carries
+ pslldq xmm0, 4 // move carries over
+ por xmm1, xmm0 // (y_3, y_2; y_1, y_0)
+
+ // Now we must reduce. This is essentially the same as the 128-bit
+ // case above, but mostly simpler because everything is smaller. The
+ // polynomial this time is p(t) = t^64 + t^4 + t^3 + t + 1.
+
+ // First, we must detach the top (`low'!) half of the result.
+ movdqa xmm0, xmm1 // (y_3, y_2; y_1, y_0) again
+ psrldq xmm1, 8 // (y_1, y_0; 0, 0)
+
+ // Next, shift the high bits down.
+ movdqa xmm2, xmm0 // (y_3, y_2; ?, ?) again
+ movdqa xmm3, xmm0 // (y_3, y_2; ?, ?) yet again
+ movdqa xmm4, xmm0 // (y_3, y_2; ?, ?) again again
+ pslld xmm2, 31 // b_i for t
+ pslld xmm3, 29 // b_i for t^3
+ pslld xmm4, 28 // b_i for t^4
+ pxor xmm2, xmm3 // add them all together
+ pxor xmm2, xmm4
+ movdqa xmm3, xmm2 // and a copy for later
+ movq xmm2, xmm2 // zap high half
+ pslldq xmm3, 4 // contribution into high half
+ psrldq xmm2, 4 // and low half
+ pxor xmm0, xmm3
+ pxor xmm1, xmm2
+
+ // And then shift the low bits up.
+ movdqa xmm2, xmm0
+ movdqa xmm3, xmm0
+ pxor xmm1, xmm0 // mix in the unit contribution
+ psrld xmm0, 1
+ psrld xmm2, 3
+ psrld xmm3, 4
+ pxor xmm1, xmm2 // low half, unit, and t^3 contribs
+ pxor xmm0, xmm3 // t and t^4 contribs
+ pxor xmm0, xmm1 // mix them together and we're done
+.endm
+
+.macro mul96
+ // Enter with u and v in the /high/ three words of xmm0 and xmm1
+ // respectively (and zero in the low word); leave with z = u v in the
+ // high three words of xmm0, and /junk/ in the low word. Clobbers
+ // xmm1--xmm4.
+
+ // This is an inconvenient size. There's nothing for it but to do
+ // four multiplications, as if for the 128-bit case. It's possible
+ // that there's cruft in the top 32 bits of the input registers, so
+ // shift both of them up by four bytes before we start. This will
+ // mean that the high 64 bits of the result (from GCM's viewpoint)
+ // will be zero.
+ // xmm0 = // (0, u_2; u_1, u_0)
+ // xmm1 = // (0, v_2; v_1, v_0)
+ movdqa xmm2, xmm1 // (0, v_2; v_1, v_0) again
+ movdqa xmm3, xmm0 // (0, u_2; u_1, u_0) again
+ movdqa xmm4, xmm0 // (0, u_2; u_1, u_0) yet again
+ pclmulhqlqdq xmm2, xmm0 // u_2 (v_1 t^32 + v_0) = e_0
+ pclmullqlqdq xmm0, xmm1 // u_2 v_2 = d = (0; d)
+ pclmulhqlqdq xmm3, xmm1 // v_2 (u_1 t^32 + u_0) = e_1
+ pclmulhqhqdq xmm4, xmm1 // u_0 v_0 + (u_1 v_0 + u_0 v_1) t^32
+ // + u_1 v_1 t^64 = f
+
+ // Extract the high and low halves of the 192-bit result. We don't
+ // need be too picky about the unused high words of the result
+ // registers. The answer we want is d t^128 + e t^64 + f, where e =
+ // e_0 + e_1.
+ //
+ // The place values for the two halves are (t^160, t^128; t^96, ?)
+ // and (?, t^64; t^32, 1).
+ psrldq xmm0, 8 // (d; 0) = d t^128
+ pxor xmm2, xmm3 // e = (e_0 + e_1)
+ movdqa xmm1, xmm4 // f again
+ pxor xmm0, xmm2 // d t^128 + e t^64
+ psrldq xmm2, 12 // e[31..0] t^64
+ psrldq xmm1, 4 // f[95..0]
+ pslldq xmm4, 8 // f[127..96]
+ pxor xmm1, xmm2 // low 96 bits of result
+ pxor xmm0, xmm4 // high 96 bits of result
+
+ // Next, shift everything one bit to the left to compensate for GCM's
+ // strange ordering. This will be easier if we shift up the high
+ // half by a word before we start. After this we're in GCM bizarro-
+ // world.
+ movdqa xmm3, xmm1 // low half again
+ pslldq xmm0, 4 // shift high half
+ psrld xmm1, 31 // shift low half down: just carries
+ movdqa xmm2, xmm0 // copy high half
+ pslld xmm3, 1 // shift low half down: drop carries
+ psrld xmm0, 31 // shift high half up: just carries
+ pslld xmm2, 1 // shift high half down: drop carries
+ movdqa xmm4, xmm0 // copy high carries for carry-around
+ pslldq xmm0, 4 // shift carries down
+ pslldq xmm1, 4
+ psrldq xmm4, 12 // the big carry wraps around
+ por xmm1, xmm3
+ por xmm0, xmm2
+ por xmm1, xmm4
+
+ // Finally, the reduction. This is essentially the same as the
+ // 128-bit case, except that the polynomial is p(t) = t^96 + t^10 +
+ // t^9 + t^6 + 1. The degrees are larger but not enough to cause
+ // trouble for the general approach.
+
+ // First, shift the high bits down.
+ movdqa xmm2, xmm0 // copies of the high part
+ movdqa xmm3, xmm0
+ movdqa xmm4, xmm0
+ pslld xmm2, 26 // b_i for t^6
+ pslld xmm3, 23 // b_i for t^9
+ pslld xmm4, 22 // b_i for t^10
+ pxor xmm2, xmm3 // add them all together
+ pslldq xmm1, 4 // shift low part up to match
+ pxor xmm2, xmm4
+ movdqa xmm3, xmm2 // and a copy for later
+ pslldq xmm2, 8 // contribution to high half
+ psrldq xmm3, 4 // contribution to low half
+ pxor xmm1, xmm3
+ pxor xmm0, xmm2
+
+ // And then shift the low bits up.
+ movdqa xmm2, xmm0 // copies of the high part
+ movdqa xmm3, xmm0
+ pxor xmm1, xmm0 // mix in the unit contribution
+ psrld xmm0, 6
+ psrld xmm2, 9
+ psrld xmm3, 10
+ pxor xmm1, xmm2 // low half, unit, and t^9 contribs
+ pxor xmm0, xmm3 // t^6 and t^10 contribs
+ pxor xmm0, xmm1 // mix them together and we're done
+.endm
+
+.macro mul192
+ // Enter with u and v in xmm0/xmm1 and xmm2/xmm3 respectively; leave
+ // with z = u v in xmm0/xmm1 -- the top halves of the high registers
+ // are unimportant. Clobbers xmm2--xmm7.
+
+ // Start multiplying and accumulating pieces of product.
+ // xmm0 = // (u_2; u_1)
+ // xmm1 = // (u_0; ?)
+ // xmm2 = // (v_2; v_1)
+ // xmm3 = // (v_0; ?)
+ movdqa xmm4, xmm0 // (u_2; u_1) again
+ movdqa xmm5, xmm0 // (u_2; u_1) yet again
+ movdqa xmm6, xmm0 // (u_2; u_1) again again
+ movdqa xmm7, xmm1 // (u_0; ?) again
+ punpcklqdq xmm1, xmm3 // (u_0; v_0)
+ pclmulhqhqdq xmm4, xmm2 // u_1 v_1
+ pclmullqlqdq xmm3, xmm0 // u_2 v_0
+ pclmullqhqdq xmm5, xmm2 // u_2 v_1
+ pclmulhqlqdq xmm6, xmm2 // u_1 v_2
+ pxor xmm4, xmm3 // u_2 v_0 + u_1 v_1
+ pclmullqlqdq xmm7, xmm2 // u_0 v_2
+ pxor xmm5, xmm6 // b = u_2 v_1 + u_1 v_2
+ movdqa xmm6, xmm0 // (u_2; u_1) like a bad penny
+ pxor xmm4, xmm7 // c = u_0 v_2 + u_1 v_1 + u_2 v_0
+ pclmullqlqdq xmm0, xmm2 // a = u_2 v_2
+ pclmulhqhqdq xmm6, xmm1 // u_1 v_0
+ pclmulhqlqdq xmm2, xmm1 // u_0 v_1
+ pclmullqhqdq xmm1, xmm1 // e = u_0 v_0
+ pxor xmm2, xmm6 // d = u_1 v_0 + u_0 v_1
+
+ // Next, the piecing together of the product.
+ // xmm0 = // (a_1; a_0) = a = u_2 v_2
+ // xmm5 = // (b_1; b_0) = b = u_1 v_2 + u_2 v_1
+ // xmm4 = // (c_1; c_0) = c = u_0 v_2 +
+ // u_1 v_1 + u_2 v_0
+ // xmm2 = // (d_1; d_0) = d = u_0 v_1 + u_1 v_0
+ // xmm1 = // (e_1; e_0) = e = u_0 v_0
+ // xmm3, xmm6, xmm7 spare
+ movdqa xmm3, xmm2 // (d_1; d_0) again
+ movdqa xmm6, xmm5 // (b_1; b_0) again
+ pslldq xmm2, 8 // (0; d_1)
+ psrldq xmm5, 8 // (b_0; 0)
+ psrldq xmm3, 8 // (d_0; 0)
+ pslldq xmm6, 8 // (0; b_1)
+ pxor xmm5, xmm2 // (b_0; d_1)
+ pxor xmm0, xmm6 // x_2 = (a_1; a_0 + b_1)
+ pxor xmm3, xmm1 // x_0 = (e_1 + d_0; e_0)
+ pxor xmm4, xmm5 // x_1 = (b_0 + c_1; c_0 + d_1)
+
+ // Now, shift it right (from GCM's point of view) by one bit, and try
+ // to leave the result in less random registers. After this, we'll
+ // be in GCM bizarro-world.
+ // xmm1, xmm2, xmm5, xmm6, xmm7 spare
+ movdqa xmm5, xmm0 // copy x_2
+ movdqa xmm1, xmm4 // copy x_1
+ movdqa xmm2, xmm3 // copy x_0
+ psrld xmm0, 31 // x_2 carries
+ psrld xmm4, 31 // x_1 carries
+ psrld xmm3, 31 // x_0 carries
+ pslld xmm5, 1 // x_2 shifted
+ pslld xmm1, 1 // x_1 shifted
+ pslld xmm2, 1 // x_0 shifted
+ movdqa xmm6, xmm0 // x_2 carry copy
+ movdqa xmm7, xmm4 // x_1 carry copy
+ pslldq xmm0, 4 // x_2 carry shifted
+ pslldq xmm4, 4 // x_1 carry shifted
+ pslldq xmm3, 4 // x_0 carry shifted
+ psrldq xmm6, 12 // x_2 carry out
+ psrldq xmm7, 12 // x_1 carry out
+ por xmm0, xmm5 // (y_5; y_4)
+ por xmm1, xmm4
+ por xmm2, xmm3
+ por xmm1, xmm6 // (y_3; y_2)
+ por xmm2, xmm7 // (y_1; y_0)
+
+ // Next, the reduction. Our polynomial this time is p(x) = t^192 +
+ // t^7 + t^2 + t + 1. Yes, the magic numbers are the same as the
+ // 128-bit case. I don't know why.
+
+ // First, shift the high bits down.
+ // xmm0 = // (y_5; y_4)
+ // xmm1 = // (y_3; y_2)
+ // xmm2 = // (y_1; y_0)
+ // xmm3--xmm7 spare
+ movdqa xmm3, xmm0 // (y_5; y_4) copy
+ movdqa xmm4, xmm0 // (y_5; y_4) copy
+ movdqa xmm5, xmm0 // (y_5; y_4) copy
+ pslld xmm3, 31 // (y_5; y_4) b_i for t
+ pslld xmm4, 30 // (y_5; y_4) b_i for t^2
+ pslld xmm5, 25 // (y_5; y_4) b_i for t^7
+ movq xmm6, xmm1 // (y_3; 0) copy
+ pxor xmm3, xmm4
+ movq xmm7, xmm1 // (y_3; 0) copy
+ pxor xmm3, xmm5
+ movq xmm5, xmm1 // (y_3; 0) copy
+ movdqa xmm4, xmm3 // (y_5; y_4) b_i combined
+ pslld xmm6, 31 // (y_3; 0) b_i for t
+ pslld xmm7, 30 // (y_3; 0) b_i for t^2
+ pslld xmm5, 25 // (y_3; 0) b_i for t^7
+ psrldq xmm3, 12 // (y_5; y_4) low contrib
+ pslldq xmm4, 4 // (y_5; y_4) high contrib
+ pxor xmm6, xmm7
+ pxor xmm2, xmm3
+ pxor xmm6, xmm5
+ pxor xmm1, xmm4
+ pslldq xmm6, 4
+ pxor xmm2, xmm6
+
+ // And finally shift the low bits up. Unfortunately, we also have to
+ // split the low bits out.
+ // xmm0 = // (y'_5; y'_4)
+ // xmm1 = // (y'_3; y'_2)
+ // xmm2 = // (y'_1; y'_0)
+ movdqa xmm5, xmm1 // copies of (y'_3; y'_2)
+ movdqa xmm6, xmm1
+ movdqa xmm7, xmm1
+ psrldq xmm1, 8 // bring down (y'_2; ?)
+ movdqa xmm3, xmm0 // copies of (y'_5; y'_4)
+ movdqa xmm4, xmm0
+ punpcklqdq xmm1, xmm2 // (y'_2; y'_1)
+ psrldq xmm2, 8 // (y'_0; ?)
+ pxor xmm2, xmm5 // low half and unit contrib
+ pxor xmm1, xmm0
+ psrld xmm5, 1
+ psrld xmm0, 1
+ psrld xmm6, 2
+ psrld xmm3, 2
+ psrld xmm7, 7
+ psrld xmm4, 7
+ pxor xmm2, xmm6 // low half, unit, t^2 contribs
+ pxor xmm1, xmm3
+ pxor xmm5, xmm7 // t and t^7 contribs
+ pxor xmm0, xmm4
+ pxor xmm5, xmm2 // mix everything together
+ pxor xmm0, xmm1
+ movq xmm1, xmm5 // shunt (z_0; ?) into proper place
+.endm
+
+.macro mul256
+ // Enter with u and v in xmm0/xmm1 and xmm2/xmm3 respectively; leave
+ // with z = u v in xmm0/xmm1. Clobbers xmm2--xmm7. On 32-bit x86,
+ // requires 16 bytes aligned space at SP; on amd64, also clobbers
+ // xmm8.
+
+ // Now it's starting to look worthwhile to do Karatsuba. Suppose
+ // u = u_0 + u_1 B and v = v_0 + v_1 B. Then
+ //
+ // u v = (u_0 v_0) + (u_0 v_1 + u_1 v_0) B + (u_1 v_1) B^2
+ //
+ // Name these coefficients of B^i be a, b, and c, respectively, and
+ // let r = u_0 + u_1 and s = v_0 + v_1. Then observe that
+ //
+ // q = r s = (u_0 + u_1) (v_0 + v_1)
+ // = (u_0 v_0) + (u1 v_1) + (u_0 v_1 + u_1 v_0)
+ // = a + d + c
+ //
+ // The first two terms we've already calculated; the last is the
+ // remaining one we want. We'll set B = t^128. We know how to do
+ // 128-bit multiplications already, and Karatsuba is too annoying
+ // there, so there'll be 12 multiplications altogether, rather than
+ // the 16 we'd have if we did this the naïve way.
+ //
+ // On x86, there aren't quite enough registers, so spill one for a
+ // bit. On AMD64, we can keep on going, so it's all good.
+
+ // xmm0 = // u_1 = (u_11; u_10)
+ // xmm1 = // u_0 = (u_01; u_00)
+ // xmm2 = // v_1 = (v_11; v_10)
+ // xmm3 = // v_0 = (v_01; v_00)
+ movdqa xmm4, xmm0 // u_1 again
+#if CPUFAM_X86
+ movdqa [esp + 0], xmm3
+#elif CPUFAM_AMD64
+ movdqa xmm8, xmm3
+# define V0 xmm8
+#endif
+ pxor xmm4, xmm1 // u_* = (u_01 + u_11; u_00 + u_10)
+ pxor xmm3, xmm2 // v_* = (v_01 + v_11; v_00 + v_10)
+
+ // Start by building the cross product, q = u_* v_*.
+ movdqa xmm7, xmm4 // more copies of u_*
+ movdqa xmm5, xmm4
+ movdqa xmm6, xmm4
+ pclmullqhqdq xmm4, xmm3 // u_*1 v_*0
+ pclmulhqlqdq xmm7, xmm3 // u_*0 v_*1
+ pclmullqlqdq xmm5, xmm3 // u_*1 v_*1
+ pclmulhqhqdq xmm6, xmm3 // u_*0 v_*0
+ pxor xmm4, xmm7 // u_*1 v_*0 + u_*0 v_*1
+ movdqa xmm7, xmm4
+ pslldq xmm4, 8
+ psrldq xmm7, 8
+ pxor xmm5, xmm4 // q_1
+ pxor xmm6, xmm7 // q_0
+
+ // Next, work on the high half, a = u_1 v_1.
+ movdqa xmm3, xmm0 // more copies of u_1
+ movdqa xmm4, xmm0
+ movdqa xmm7, xmm0
+ pclmullqhqdq xmm0, xmm2 // u_11 v_10
+ pclmulhqlqdq xmm3, xmm2 // u_10 v_11
+ pclmullqlqdq xmm4, xmm2 // u_11 v_11
+ pclmulhqhqdq xmm7, xmm2 // u_10 v_10
+#if CPUFAM_X86
+ movdqa xmm2, [esp + 0]
+# define V0 xmm2
+#endif
+ pxor xmm0, xmm3 // u_10 v_11 + u_11 v_10
+ movdqa xmm3, xmm0
+ pslldq xmm0, 8
+ psrldq xmm3, 8
+ pxor xmm4, xmm0 // x_1 = a_1
+ pxor xmm7, xmm3 // a_0
+
+ // Mix that into the product now forming in xmm4--xmm7.
+ pxor xmm5, xmm4 // a_1 + q_1
+ pxor xmm6, xmm7 // a_0 + q_0
+ pxor xmm5, xmm7 // a_0 + (a_1 + q_1)
+
+ // Finally, the low half, c = u_0 v_0.
+ movdqa xmm0, xmm1 // more copies of u_0
+ movdqa xmm3, xmm1
+ movdqa xmm7, xmm1
+ pclmullqhqdq xmm1, V0 // u_01 v_00
+ pclmulhqlqdq xmm0, V0 // u_00 v_01
+ pclmullqlqdq xmm3, V0 // u_01 v_01
+ pclmulhqhqdq xmm7, V0 // u_00 v_00
+ pxor xmm0, xmm1 // u_10 v_11 + u_11 v_10
+ movdqa xmm1, xmm0
+ pslldq xmm0, 8
+ psrldq xmm1, 8
+ pxor xmm3, xmm0 // c_1
+ pxor xmm7, xmm1 // x_0 = c_0
+
+ // And mix that in to complete the product.
+ pxor xmm6, xmm3 // (a_0 + q_0) + c_1
+ pxor xmm5, xmm3 // x_2 = a_0 + (a_1 + c_1 + q_1) = a_0 + b_1
+ pxor xmm6, xmm7 // x_1 = (a_0 + c_0 + q_0) + c_1 = b_0 + c_1
+
+#undef V0
+
+ // Now we need to shift that whole lot one bit to the left. This
+ // will also give us an opportunity to put the product back in
+ // xmm0--xmm3. This is a slightly merry dance because it's nearly
+ // pipelined but we don't have enough registers.
+ //
+ // After this, we'll be in GCM bizarro-world.
+ movdqa xmm0, xmm4 // x_3 again
+ psrld xmm4, 31 // x_3 carries
+ pslld xmm0, 1 // x_3 shifted left
+ movdqa xmm3, xmm4 // x_3 copy carries
+ movdqa xmm1, xmm5 // x_2 again
+ pslldq xmm4, 4 // x_3 carries shifted up
+ psrld xmm5, 31 // x_2 carries
+ psrldq xmm3, 12 // x_3 big carry out
+ pslld xmm1, 1 // x_2 shifted left
+ por xmm0, xmm4 // x_3 mixed together
+ movdqa xmm4, xmm5 // x_2 copy carries
+ movdqa xmm2, xmm6 // x_1 again
+ pslldq xmm5, 4 // x_2 carries shifted up
+ psrld xmm6, 31 // x_1 carries
+ psrldq xmm4, 12 // x_2 big carry out
+ pslld xmm2, 1 // x_1 shifted
+ por xmm1, xmm5 // x_2 mixed together
+ movdqa xmm5, xmm6 // x_1 copy carries
+ por xmm1, xmm3 // x_2 with carry from x_3
+ movdqa xmm3, xmm7 // x_0 again
+ pslldq xmm6, 4 // x_1 carries shifted up
+ psrld xmm7, 31 // x_2 carries
+ psrldq xmm5, 12 // x_1 big carry out
+ pslld xmm3, 1 // x_0 shifted
+ por xmm2, xmm6 // x_1 mixed together
+ pslldq xmm7, 4 // x_0 carries shifted up
+ por xmm2, xmm4 // x_1 with carry from x_2
+ por xmm3, xmm7 // x_0 mixed together
+ por xmm3, xmm5 // x_0 with carry from x_1
+
+ // Now we must reduce. This is essentially the same as the 128-bit
+ // case above, but more complicated because everything is bigger.
+ // The polynomial this time is p(t) = t^256 + t^10 + t^5 + t^2 + 1.
+
+ // First, shift the high bits down.
+ movdqa xmm4, xmm0 // y_3 again
+ movdqa xmm5, xmm0 // y_3 yet again
+ movdqa xmm6, xmm0 // y_3 again again
+ pslld xmm4, 30 // y_3: b_i for t^2
+ pslld xmm5, 27 // y_3: b_i for t^5
+ pslld xmm6, 22 // y_3: b_i for t^10
+ movdqa xmm7, xmm1 // y_2 again
+ pxor xmm4, xmm5
+ movdqa xmm5, xmm1 // y_2 again
+ pxor xmm4, xmm6
+ movdqa xmm6, xmm1 // y_2 again
+ pslld xmm7, 30 // y_2: b_i for t^2
+ pslld xmm5, 27 // y_2: b_i for t^5
+ pslld xmm6, 22 // y_2: b_i for t^10
+ pxor xmm7, xmm5
+ movdqa xmm5, xmm4
+ pxor xmm7, xmm6
+ psrldq xmm4, 4
+ movdqa xmm6, xmm7
+ pslldq xmm5, 12
+ psrldq xmm7, 4
+ pxor xmm2, xmm4
+ pslldq xmm6, 12
+ pxor xmm3, xmm7
+ pxor xmm1, xmm5
+ pxor xmm2, xmm6
+
+ // And then shift the low bits up.
+ movdqa xmm4, xmm0 // y_3 again
+ movdqa xmm5, xmm1 // y_2 again
+ movdqa xmm6, xmm0 // y_3 yet again
+ movdqa xmm7, xmm1 // y_2 yet again
+ pxor xmm2, xmm0 // y_1 and unit contrib from y_3
+ pxor xmm3, xmm1 // y_0 and unit contrib from y_2
+ psrld xmm0, 2
+ psrld xmm1, 2
+ psrld xmm4, 5
+ psrld xmm5, 5
+ psrld xmm6, 10
+ psrld xmm7, 10
+ pxor xmm0, xmm2 // y_1, with y_3 units and t^2
+ pxor xmm1, xmm3 // y_0, with y_2 units and t^2
+ pxor xmm4, xmm6 // y_3 t^5 and t^10 contribs
+ pxor xmm5, xmm7 // y_2 t^5 and t^10 contribs
+ pxor xmm0, xmm4 // high half of reduced result
+ pxor xmm1, xmm5 // low half; all done
+.endm
+
+///--------------------------------------------------------------------------
+/// Main code.
+
+// There are a number of representations of field elements in this code and
+// it can be confusing.
+//
+// * The `external format' consists of a sequence of contiguous bytes in
+// memory called a `block'. The GCM spec explains how to interpret this
+// block as an element of a finite field. As discussed extensively, this
+// representation is very annoying for a number of reasons. On the other
+// hand, this code never actually deals with it directly.
+//
+// * The `register format' consists of one or more XMM registers, depending
+// on the block size. The bytes in these registers are in reverse order
+// -- so the least-significant byte of the lowest-numbered register holds
+// the /last/ byte in the block. If the block size is not a multiple of
+// 16 bytes, then there must be padding. 96-bit blocks are weird: the
+// padding is inserted at the /least/ significant end, so the register
+// holds (0, x_0; x_1, x_2); otherwise, the padding goes at the most
+// significant end.
+//
+// * The `words' format consists of a sequence of bytes, as in the
+// `external format', but, according to the blockcipher in use, the bytes
+// within each 32-bit word may be reversed (`big-endian') or not
+// (`little-endian'). Accordingly, there are separate entry points for
+// each variant, identified with `b' or `l'.
+
+#define SSEFUNC(f) \
+ FUNC(f##_avx); vzeroupper; endprologue; ENDFUNC; \
+ FUNC(f)
+
+SSEFUNC(gcm_mulk_128b_x86ish_pclmul)
+ // On entry, A points to a 128-bit field element in big-endian words
+ // format; K points to a field-element in register format. On exit,
+ // A is updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+#endif
+ endprologue
+ movdqu xmm0, [A]
+ movdqu xmm1, [K]
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ mul128
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ movdqu [A], xmm0
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_128l_x86ish_pclmul)
+ // On entry, A points to a 128-bit field element in little-endian
+ // words format; K points to a field-element in register format. On
+ // exit, A is updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+ ldgot ecx
+#endif
+ endprologue
+ movdqa xmm7, [INTADDR(swaptab_128l, ecx)]
+ movdqu xmm0, [A]
+ movdqu xmm1, [K]
+ pshufb xmm0, xmm7
+ mul128
+ pshufb xmm0, xmm7
+ movdqu [A], xmm0
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_64b_x86ish_pclmul)
+ // On entry, A points to a 64-bit field element in big-endian words
+ // format; K points to a field-element in register format. On exit,
+ // A is updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+#endif
+ endprologue
+ movq xmm0, [A]
+ movq xmm1, [K]
+ pshufd xmm0, xmm0, SHUF(1, 0, 3, 3)
+ mul64
+ pshufd xmm0, xmm0, SHUF(1, 0, 3, 3)
+ movq [A], xmm0
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_64l_x86ish_pclmul)
+ // On entry, A points to a 64-bit field element in little-endian
+ // words format; K points to a field-element in register format. On
+ // exit, A is updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+ ldgot ecx
+#endif
+ endprologue
+ movdqa xmm7, [INTADDR(swaptab_64l, ecx)]
+ movq xmm0, [A]
+ movq xmm1, [K]
+ pshufb xmm0, xmm7
+ mul64
+ pshufb xmm0, xmm7
+ movq [A], xmm0
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_96b_x86ish_pclmul)
+ // On entry, A points to a 96-bit field element in big-endian words
+ // format; K points to a field-element in register format (i.e., 16
+ // bytes, with the first four bytes zero). On exit, A is updated
+ // with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+#endif
+ endprologue
+ movq xmm0, [A + 0]
+ movd xmm2, [A + 8]
+ movdqu xmm1, [K]
+ punpcklqdq xmm0, xmm2
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ mul96
+ pshufd xmm1, xmm0, SHUF(3, 2, 1, 0)
+ psrldq xmm0, 4
+ movq [A + 0], xmm1
+ movd [A + 8], xmm0
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_96l_x86ish_pclmul)
+ // On entry, A points to a 96-bit field element in little-endian
+ // words format; K points to a field-element in register format
+ // (i.e., 16 bytes, with the first four bytes zero). On exit, A is
+ // updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+ ldgot ecx
+#endif
+ endprologue
+ movdqa xmm7, [INTADDR(swaptab_128l, ecx)]
+ movq xmm0, [A + 0]
+ movd xmm2, [A + 8]
+ movdqu xmm1, [K]
+ punpcklqdq xmm0, xmm2
+ pshufb xmm0, xmm7
+ mul96
+ pshufb xmm0, xmm7
+ movq [A + 0], xmm0
+ psrldq xmm0, 8
+ movd [A + 8], xmm0
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_192b_x86ish_pclmul)
+ // On entry, A points to a 192-bit field element in big-endian words
+ // format; K points to a field-element in register format. On exit,
+ // A is updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+#endif
+#if CPUFAM_AMD64 && ABI_WIN
+ stalloc 2*16 + 8
+ savexmm xmm6, 0
+ savexmm xmm7, 16
+#endif
+ endprologue
+ movdqu xmm0, [A + 8]
+ movq xmm1, [A + 0]
+ movdqu xmm2, [K + 0]
+ movq xmm3, [K + 16]
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ pshufd xmm1, xmm1, SHUF(1, 0, 3, 3)
+ mul192
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ pshufd xmm1, xmm1, SHUF(1, 0, 3, 3)
+ movdqu [A + 8], xmm0
+ movq [A + 0], xmm1
+#if CPUFAM_AMD64 && ABI_WIN
+ rstrxmm xmm6, 0
+ rstrxmm xmm7, 16
+ stfree 2*16 + 8
+#endif
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_192l_x86ish_pclmul)
+ // On entry, A points to a 192-bit field element in little-endian
+ // words format; K points to a field-element in register format. On
+ // exit, A is updated with the product A K.
+
+#if CPUFAM_X86
+ mov A, [esp + 4]
+ mov K, [esp + 8]
+ ldgot ecx
+#endif
+#if CPUFAM_AMD64 && ABI_WIN
+ stalloc 2*16 + 8
+ savexmm xmm6, 0
+ savexmm xmm7, 16
+#endif
+ endprologue
+ movdqu xmm0, [A + 8]
+ movq xmm1, [A + 0]
+ movdqu xmm2, [K + 0]
+ movq xmm3, [K + 16]
+ pshufb xmm0, [INTADDR(swaptab_128l, ecx)]
+ pshufb xmm1, [INTADDR(swaptab_64l, ecx)]
+ mul192
+ pshufb xmm0, [INTADDR(swaptab_128l, ecx)]
+ pshufb xmm1, [INTADDR(swaptab_64l, ecx)]
+ movdqu [A + 8], xmm0
+ movq [A + 0], xmm1
+#if CPUFAM_AMD64 && ABI_WIN
+ rstrxmm xmm6, 0
+ rstrxmm xmm7, 16
+ stfree 2*16 + 8
+#endif
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_256b_x86ish_pclmul)
+ // On entry, A points to a 256-bit field element in big-endian words
+ // format; K points to a field-element in register format. On exit,
+ // A is updated with the product A K.
+
+#if CPUFAM_X86
+ pushreg ebp
+ setfp
+ mov A, [esp + 8]
+ mov K, [esp + 12]
+ and esp, ~15
+ sub esp, 16
+#endif
+#if CPUFAM_AMD64 && ABI_WIN
+ stalloc 3*16 + 8
+ savexmm xmm6, 0
+ savexmm xmm7, 16
+ savexmm xmm8, 32
+#endif
+ endprologue
+ movdqu xmm0, [A + 16]
+ movdqu xmm1, [A + 0]
+ movdqu xmm2, [K + 0]
+ movdqu xmm3, [K + 16]
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ pshufd xmm1, xmm1, SHUF(3, 2, 1, 0)
+ mul256
+ pshufd xmm0, xmm0, SHUF(3, 2, 1, 0)
+ pshufd xmm1, xmm1, SHUF(3, 2, 1, 0)
+ movdqu [A + 16], xmm0
+ movdqu [A + 0], xmm1
+#if CPUFAM_X86
+ dropfp
+ popreg ebp
+#endif
+#if CPUFAM_AMD64 && ABI_WIN
+ rstrxmm xmm6, 0
+ rstrxmm xmm7, 16
+ rstrxmm xmm8, 32
+ stfree 3*16 + 8
+#endif
+ ret
+ENDFUNC
+
+SSEFUNC(gcm_mulk_256l_x86ish_pclmul)
+ // On entry, A points to a 256-bit field element in little-endian
+ // words format; K points to a field-element in register format. On
+ // exit, A is updated with the product A K.
+
+#if CPUFAM_X86
+ pushreg ebp
+ setfp
+ mov A, [esp + 8]
+ mov K, [esp + 12]
+ and esp, ~15
+ ldgot ecx
+ sub esp, 16
+#endif
+#if CPUFAM_AMD64 && ABI_WIN
+ stalloc 3*16 + 8
+ savexmm xmm6, 0
+ savexmm xmm7, 16
+ savexmm xmm8, 32
+#endif
+ endprologue
+ movdqa xmm7, [INTADDR(swaptab_128l, ecx)]
+ movdqu xmm0, [A + 16]
+ movdqu xmm1, [A + 0]
+ movdqu xmm2, [K + 0]
+ movdqu xmm3, [K + 16]
+ pshufb xmm0, xmm7
+ pshufb xmm1, xmm7
+ mul256
+ movdqa xmm7, [INTADDR(swaptab_128l, ecx)]
+ pshufb xmm0, xmm7
+ pshufb xmm1, xmm7
+ movdqu [A + 16], xmm0
+ movdqu [A + 0], xmm1
+#if CPUFAM_X86
+ dropfp
+ popreg ebp
+#endif
+#if CPUFAM_AMD64 && ABI_WIN
+ rstrxmm xmm6, 0
+ rstrxmm xmm7, 16
+ rstrxmm xmm8, 32
+ stfree 3*16 + 8
+#endif
+ ret
+ENDFUNC
+
+ RODATA
+
+ .balign 16
+swaptab_128l:
+ // Table for byte-swapping little-endian words-format blocks larger
+ // than 64 bits.
+ .byte 15, 14, 13, 12, 11, 10, 9, 8
+ .byte 7, 6, 5, 4, 3, 2, 1, 0
+
+ .balign 16
+swaptab_64l:
+ // Table for byte-swapping 64-bit little-endian words-format blocks.
+ .byte 7, 6, 5, 4, 3, 2, 1, 0
+ .byte 255, 255, 255, 255, 255, 255, 255, 255
+
+///----- That's all, folks --------------------------------------------------
~mdw / catacomb / blobdiff