diff --git a/docs/avx2-notes.md b/docs/avx2-notes.md
index fb8085d..d902650 100644
--- a/docs/avx2-notes.md
+++ b/docs/avx2-notes.md
@@ -1,42 +1,48 @@
-An implementation of group operations on the twisted Edwards form of
-Curve25519, using AVX2 to implement the 4-way parallel formulas of
-Hisil, Wong, Carter, and Dawson (HWCD).
-Their 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08], which
-introduced the extended coordinates used in other parts of `-dalek`,
-also describes 4-way parallel formulas for point addition and
-doubling:
+A vectorized implementation of group operations on the twisted Edwards
+form of Curve25519, using a modification of the 4-way parallel
+formulas of Hisil, Wong, Carter, and Dawson.
 
-* a unified addition algorithm taking an effective \\(2\mathbf M +
-1\mathbf D\\);
+# Overview
 
-* a doubling algorithm taking an effective \\(1\mathbf M + 1\mathbf
-S\\);
+The 2008 paper [_Twisted Edwards Curves Revisited_][hwcd08] by Hisil,
+Wong, Carter, and Dawson (HWCD) introduced the “extended coordinates”
+and mixed-model representations which are used by most Edwards curve
+implementations.
 
-* a dedicated (i.e., for distinct points) addition algorithm taking
-an effective \\(2 \mathbf M \\).
+However, they also describe 4-way parallel formulas for point addition
+and doubling: a unified addition algorithm taking an effective
+\\(2\mathbf M + 1\mathbf D\\), a doubling algorithm taking an
+effective \\(1\mathbf M + 1\mathbf S\\), and a dedicated (i.e., for
+distinct points) addition algorithm taking an effective \\(2 \mathbf M
+\\).  They compare these formulas with a 2-way parallel variant of the
+Montgomery ladder.
 
-Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of
-multiplication and squaring of generic field elements and \\(\mathbf
-D\\) represents the cost of multiplication by a curve constant.
+Unlike their serial formulas, which are used widely, their parallel
+formulas do not seem to have been implemented in software before.  The
+2-way parallel Montgomery ladder was used in 2015 by Tung Chou's
+`sandy2x` implementation.  Curiously, however, although the [`sandy2x`
+paper][sandy2x] also implements Edwards arithmetic, and cites HWCD08,
+it doesn't mention their parallel Edwards formulas.
+A 2015 paper by Hernández and López describes an AVX2 implementation
+of X25519. Neither the paper nor the code are publicly available, but
+it apparently gives only a [slight speedup][avx2trac], suggesting that
+it uses a 4-way parallel Montgomery ladder rather than parallel
+Edwards formulas.
 
-These formulas do not seem to have been implemented using SIMD before.
-A 2015 paper by Hernández and López mentions using AVX2 for the X25519
-Montgomery ladder, but neither the paper nor the code are publicly
-available, and it apparently gives only a [slight speedup][avx2trac].
-The 2008 HWCD paper also describes and analyzes a 2-wide variant of the
-Montgomery ladder (for comparison with parallel Edwards formulas); this
-strategy was used in 2015 by Tung Chou's `sandy2x` implementation, which
-used a 2-wide field implementation in 128-bit vector registers.
-Curiously, however, although the [`sandy2x` paper][sandy2x] also
-implements Edwards arithmetic, and cites the HWCD paper, it doesn't
-mention the parallel formulas from HWCD, suggesting that they have been
-overlooked for software implementations.
+The reason may be that HWCD08 describe their formulas as operating on
+four independent processors, which would make a software
+implementation impractical: all of the operations are too low-latency
+to effectively synchronize.  But a closer inspection reveals that the
+(more expensive) multiplication and squaring steps are uniform, while
+the instruction divergence occurs in the (much cheaper) addition and
+subtraction steps.  This means that a SIMD implementation can perform
+the expensive steps uniformly, and handle divergence in the
+inexpensive steps using masking.
 
-The notes below describe a tweak to the \\( 2\mathbf M + 1\mathbf D \\)
-unified addition formulas to give \\( 2\mathbf M \\) readdition with
-\\(1\mathbf D\\) precomputation, and a tweak to the doubling formulas to
-avoid an extra reduction.  These tweaked formulas are the ones used by
-the `avx2` backend of `curve25519-dalek`.
+These notes describe modifications to the original parallel formulas
+to allow a SIMD implementation, and this module contains an
+implementation of the modified formulas using 256-bit AVX2 vector
+operations.
 
 # Parallel formulas in HWCD'08
 
@@ -60,218 +66,153 @@ and the unified addition algorithm is presented as follows:
 |                  | \\( R\_1 \gets R\_6 - R\_5 \\) | \\( R\_2 \gets R\_8 - R\_7 \\) | \\( R\_3 \gets R\_8 + R\_7 \\) | \\( R\_4 \gets R\_6 + R\_5 \\) |
 | \\(1\mathbf M\\) | \\( X\_3 \gets R\_1 R\_2 \\)   | \\( Y\_3 \gets R\_3 R\_4 \\)   | \\( T\_3 \gets R\_1 R\_4 \\)   | \\( Z\_3 \gets R\_2 R\_3 \\)   |
 
-Here \\( k = 2d \\) is a curve constant.
+Here \\(\mathbf M\\) and \\(\mathbf S\\) represent the cost of
+multiplication and squaring of generic field elements, \\(\mathbf D\\)
+represents the cost of multiplication by a curve constant (in this
+case \\( k = 2d \\)).
 
-For a software implementation, each processor's operations are too
-low-latency to parallelize across threads.  However, the main cost
-is in the multiplication and squaring steps, which are uniform, while
-the divergent steps involve inexpensive additions and subtractions.
+Notice that the \\(1\mathbf M\\) and \\(1\mathbf S\\) steps are
+uniform.  The non-uniform steps are all inexpensive additions or
+subtractions, with the exception of the multiplication by the curve
+constant \\(k = 2d\\):
+$$
+R\_7 \gets 2 d R\_7.
+$$
 
-This means we can use SIMD to implement the expensive portions in
-parallel, and handle the instruction divergence on the inexpensive parts
-using masking.
-
-The remaining obstacle to parallelism is the multiplication by the curve
-constant \\(k = 2d\\).  In the Curve25519 case, this is
-
-$$ k \equiv 2 \frac{-121665}{121666} \\ \equiv 16295367250680780974490674513165176452449235426866156013048779062215315747161 \pmod p. $$
-
-HWCD suggest parallelising this step by breaking \\(k\\) into four
+HWCD suggest parallelising this step by breaking \\(k = 2d\\) into four
 parts as \\(k = k_0 + 2\^n k_1 + 2\^{2n} k_2 + 2\^{3n} k_3 \\) and
-computing \\(k_i R_7 \\) in parallel.  However, this would be
-somewhat awkward in our case, since we would normally represent
-\\(k\\) as \\( 10 \\) 32-bit limbs, and \\(10 \\) is not divisible
-by \\(4\\), so we would need a specialized routine to perform a
-vectorized multiplication by 64-bit constants.
+computing \\(k_i R_7 \\) in parallel.  This is quite awkward, but if
+the curve constant is a ratio \\( d = d\_1/d\_2 \\), then projective
+coordinates allow us to instead compute 
+$$
+(R\_5, R\_6, R\_7, R\_8) \gets (d\_2 R\_5, d\_2 R\_6, 2d\_1 R\_7, d\_2 R\_8).
+$$
+This can be performed as a uniform multiplication by a vector of
+constants, and if \\(d\_1, d\_2\\) are small, it is relatively
+inexpensive.  (This trick was suggested by Mike Hamburg).
+In the Curve25519 case, we have 
+$$
+d = \frac{d\_1}{d\_2} = \frac{-121665}{121666};
+$$
+Since \\(2 \cdot 121666 < 2\^{18}\\), all the constants above fit (up
+to sign) in 32 bits, so this can be done in parallel as four
+multiplications by small constants \\( (121666, 121666, 2\cdot 121665,
+2\cdot 121666) \\), followed by a negation to compute \\( - 2\cdot 121665\\).
 
-Instead, since we are working projectively, we can multiply
-\\(R_7\\) by \\( -2\cdot 121665 \\) and multiply the other three
-variables by \\(121666\\).  This trick was suggested by Mike
-Hamburg.  Ignoring the sign for the moment, since
-\\(2 \cdot 121666 < 2\^{18}\\), all these constants fit in 32 bits,
-so (up to sign) this can be done in parallel as four multiplications
-by small constants \\( (121666, 121666, 2\cdot 121665, 2\cdot 121666) \\).
+# Modified parallel formulas
 
-How do we handle the sign?
-Since we're primarily interested in Ristretto performance, not
-Curve25519 performance, we could alternately work on the
-\\(4\\)-isogenous "IsoEd25519" curve, which has \\(d = 121665\\).
-However, this would only save the negation step, since multiplying
-one field element by a 32-bit constant is not much easier than
-multiplying four field elements by 32-bit constants, and it would
-prevent accelerating Curve25519, so we don't make this choice.
-Instead, we just negate one lane, and move the \\(1 \mathbf D\\)
-into precomputation (see below).
-
-# Tweaked formulas
-
-After tweaking the formulas as described above, we obtain the
-following.  To avoid confusion with the original HWCD formulas,
-temporary variables are named \\(S\\) instead of \\(R\\) and are in
-static single-assignment form.
+Using the modifications sketched above, we can write SIMD-friendly
+versions of the parallel formulas as follows.  To avoid confusion with
+the original formulas, temporary variables are named \\(S\\) instead
+of \\(R\\) and are in static single-assignment form.
 
 ## Addition
 
-This implementation only implements readdition, but the tweaked addition
-formulas are described first.  To add points \\(P_1 = (X_1 : Y_1 : Z_1 :
-T_1) \\) and \\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\), we compute
-
+To add points
+\\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\)
+and
+\\(P_2 = (X_2 : Y_2 : Z_2 : T_2 ) \\),
+we compute
 $$
 \begin{aligned}
-S\_0 &\gets Y\_1 - X\_1 \\\\
-S\_1 &\gets Y\_1 + X\_1 \\\\
-S\_2 &\gets Y\_2 - X\_2 \\\\
-S\_3 &\gets Y\_2 + X\_2
+(S\_0 &&,&& S\_1 &&,&& S\_2 &&,&& S\_3 )
+&\gets
+(Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Y\_2 - X\_2&&,&& Y\_2 + X\_2)
+\\\\
+(S\_4 &&,&& S\_5 &&,&& S\_6 &&,&& S\_7 )
+&\gets
+(S\_0 \cdot S\_2&&,&& S\_1 \cdot S\_3&&,&& Z\_1 \cdot Z\_2&&,&& T\_1 \cdot T\_2)
+\\\\
+(S\_8    &&,&& S\_9    &&,&& S\_{10} &&,&& S\_{11} )
+&\gets
+(d\_2 \cdot S\_4 &&,&& d\_2 \cdot S\_5 &&,&& 2 d\_2 \cdot S\_6 &&,&& 2 d\_1 \cdot S\_7 )
+\\\\
+(S\_{12} &&,&& S\_{13} &&,&& S\_{14} &&,&& S\_{15})
+&\gets
+(S\_9 - S\_8&&,&& S\_9 + S\_8&&,&& S\_{10} - S\_{11}&&,&& S\_{10} + S\_{11})
+\\\\
+(X\_3&&,&& Y\_3&&,&& Z\_3&&,&& T\_3)
+&\gets
+(S\_{12} \cdot S\_{14}&&,&& S\_{15} \cdot S\_{13}&&,&& S\_{15} \cdot S\_{14}&&,&& S\_{12} \cdot S\_{13})
 \end{aligned}
 $$
-
-$$
-\begin{aligned}
-S\_4 &\gets S\_0 S\_2 \\\\
-S\_5 &\gets S\_1 S\_3 \\\\
-S\_6 &\gets Z\_1 Z\_2 \\\\
-S\_7 &\gets T\_1 T\_2
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-S\_8    &\gets S\_4 \cdot 121666 \\\\
-S\_9    &\gets S\_5 \cdot 121666 \\\\
-S\_{10} &\gets S\_6 \cdot 2 \cdot 121666 \\\\
-S\_{11} &\gets S\_7 \cdot -2 \cdot 121665
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-S\_{12} &\gets S\_9 - S\_8 \\\\
-S\_{13} &\gets S\_9 + S\_8 \\\\
-S\_{14} &\gets S\_{10} - S\_{11} \\\\
-S\_{15} &\gets S\_{10} + S\_{11}
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-X\_3 &\gets S\_{12} S\_{14} \\\\
-Y\_3 &\gets S\_{15} S\_{13} \\\\
-Z\_3 &\gets S\_{15} S\_{14} \\\\
-T\_3 &\gets S\_{12} S\_{13}
-\end{aligned}
-$$
-
 to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\).
+This costs \\( 2\mathbf M + 1 \mathbf D\\).
 
 ## Readdition
 
-If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we can precompute
-
+If the point \\( P_2 = (X\_2 : Y\_2 : Z\_2 : T\_2) \\) is fixed, we
+can cache the multiplication of the curve constants by computing
 $$
 \begin{aligned}
-S\_2 &\gets Y\_2 - X\_2 \\\\
-S\_3 &\gets Y\_2 + X\_2
+(S\_2' &&,&& S\_3' &&,&& Z\_2' &&,&& T\_2' )
+&\gets
+(d\_2 \cdot (Y\_2 - X\_2)&&,&& d\_2 \cdot (Y\_1 + X\_1)&&,&& 2d\_2 \cdot Z\_2 &&,&& 2d\_1 \cdot T\_2).
 \end{aligned}
 $$
-
+This costs \\( 1\mathbf D\\); with \\( (S\_2', S\_3', Z\_2', T\_2')\\)
+in hand, the addition formulas above become
 $$
 \begin{aligned}
-S\_2'    &\gets S\_2 \cdot 121666 \\\\
-S\_3'    &\gets S\_3 \cdot 121666 \\\\
-Z\_2'    &\gets Z\_2 \cdot 2 \cdot 121666 \\\\
-T\_2'    &\gets T\_2 \cdot -2 \cdot 121665 \\\\
+(S\_0 &&,&& S\_1 &&,&& Z\_1 &&,&& T\_1 )
+&\gets
+(Y\_1 - X\_1&&,&& Y\_1 + X\_1&&,&& Z\_1 &&,&& T\_1)
+\\\\
+(S\_8    &&,&& S\_9    &&,&& S\_{10} &&,&& S\_{11} )
+&\gets
+(S\_0 \cdot S\_2' &&,&& S\_1 \cdot S\_3'&&,&& Z\_1 \cdot Z\_2' &&,&& T\_1 \cdot T\_2')
+\\\\
+(S\_{12} &&,&& S\_{13} &&,&& S\_{14} &&,&& S\_{15})
+&\gets
+(S\_9 - S\_8&&,&& S\_9 + S\_8&&,&& S\_{10} - S\_{11}&&,&& S\_{10} + S\_{11})
+\\\\
+(X\_3&&,&& Y\_3&&,&& Z\_3&&,&& T\_3)
+&\gets
+(S\_{12} \cdot S\_{14}&&,&& S\_{15} \cdot S\_{13}&&,&& S\_{15} \cdot S\_{14}&&,&& S\_{12} \cdot S\_{13})
 \end{aligned}
 $$
-
-to obtain the `CachedPoint` \\( (S\_2', S\_3', Z\_2', T\_2') \\).
-This precomputation is essentially the same as that suggested in
-§3.1 of HWCD, with the difference that the multiplication by the curve
-constant \\( -121665 / 121666 \\) is spread over all four
-coordinates, to allow a vectorized computation of four
-multiplications of small constants instead of a serial computation
-of multiplication by a large constant.
-
-To perform readdition of \\(P_1 = (X_1 : Y_1 : Z_1 : T_1) \\) and
-\\(P_2 = (S\_2', S\_3', Z\_2', T\_2') \\), we compute
-
-$$
-\begin{aligned}
-S\_0 &\gets Y\_1 - X\_1 \\\\
-S\_1 &\gets Y\_1 + X\_1
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-S\_8    &\gets S\_0 S\_2' \\\\
-S\_9    &\gets S\_1 S\_3' \\\\
-S\_{10} &\gets Z\_1 Z\_2' \\\\
-S\_{11} &\gets T\_1 T\_2'
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-S\_{12} &\gets S\_9 - S\_8 \\\\
-S\_{13} &\gets S\_9 + S\_8 \\\\
-S\_{14} &\gets S\_{10} - S\_{11} \\\\
-S\_{15} &\gets S\_{10} + S\_{11}
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-X\_3 &\gets S\_{12} S\_{14} \\\\
-Y\_3 &\gets S\_{15} S\_{13} \\\\
-Z\_3 &\gets S\_{15} S\_{14} \\\\
-T\_3 &\gets S\_{12} S\_{13}
-\end{aligned}
-$$
-
-to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = P\_1 + P\_2 \\).
-
-Compared to the addition formulas above, this saves \\( 1\mathbf D \\).
+which costs only \\( 2\mathbf M \\).  This precomputation is
+essentially similar to the precomputation that HWCD suggest for their
+serial formulas.  Because the cost of precomputation and then
+readdition is the same as addition, it's sufficient to only
+implement caching and readdition.
 
 ## Doubling
 
+The non-uniform portions of the (re)addition formulas have a fairly
+regular structure.  Unfortunately, this is not the case for the
+doubling formulas, which are much less nice.
+
 To double a point \\( P = (X\_1 : Y\_1 : Z\_1 : T\_1) \\), we compute
-
-$$ S\_0 \gets X\_1 + Y\_1 $$
-
 $$
 \begin{aligned}
-S\_1 &\gets X\_1\^2 \\\\
-S\_2 &\gets Y\_1\^2 \\\\
-S\_3 &\gets Z\_1\^2 \\\\
-S\_4 &\gets S\_0\^2
+(X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& S\_0)
+&\gets
+(X\_1 &&,&& Y\_1 &&,&& Z\_1 &&,&& X\_1 + Y\_1)
+\\\\
+(S\_1 &&,&& S\_2 &&,&& S\_3 &&,&& S\_4 )
+&\gets
+(X\_1\^2 &&,&& Y\_1\^2&&,&& Z\_1\^2 &&,&& S\_0\^2)
+\\\\
+(S\_5 &&,&& S\_6 &&,&& S\_8 &&,&& S\_9 )
+&\gets
+(S\_1 + S\_2 &&,&& S\_1 - S\_2 &&,&& S\_1 + 2S\_3 - S\_2 &&,&& S\_1 + S\_2 - S\_4)
+\\\\
+(X\_3 &&,&& Y\_3 &&,&& Z\_3 &&,&& T\_3 )
+&\gets
+(S\_8 \cdot S\_9 &&,&& S\_5 \cdot S\_6 &&,&& S\_8 \cdot S\_6 &&,&& S\_5 \cdot S\_9)
 \end{aligned}
 $$
-
-$$
-\begin{aligned}
-S\_5 &\gets S\_1 + S\_2 \\\\
-S\_6 &\gets S\_1 - S\_2 \\\\
-S\_7 &\gets 2S\_3 \\\\
-S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\
-S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4
-\end{aligned}
-$$
-
-$$
-\begin{aligned}
-X\_3 &\gets S\_8 S\_9 \\\\
-Y\_3 &\gets S\_5 S\_6 \\\\
-Z\_3 &\gets S\_8 S\_6 \\\\
-T\_3 &\gets S\_5 S\_9
-\end{aligned}
-$$
-
 to obtain \\( P\_3 = (X\_3 : Y\_3 : Z\_3 : T\_3) = [2]P\_1 \\).
 
-Unlike the (re)addition formulas, the divergent parts of these formulas
-are less nice.  However, with some careful bounds-juggling, it is
-possible to implement them without inserting extra carry chains, as
-described below.
+The intermediate step between the squaring and multiplication requires
+a long chain of additions, but with some care and finesse,
+described below, it is possible (in our case) to arrange this
+computation without requiring an intermediate reduction.
+
+However, it does mean that the doubling formulas have proportionately
+more vectorization overhead than the (re)addition formulas.  The
+effects of this are discussed in the comparison section below.
 
 # Field element representation
 
@@ -329,39 +270,36 @@ much difficulty.  Going the other direction, to extend this to AVX512,
 we could either run two point operations in parallel in lower and upper
 halves of the registers, or use 2-way parallelism within a field operation.
 
-# Handling the Doubling Formulas
+# Avoiding Overflow in Doubling
 
-The non-parallel portion of the doubling formulas is
+To analyze the size of the field element coefficients during the
+computations, we can parameterize the bounds on the limbs of each
+field element by \\( b \in \mathbb R \\) representing the excess bits
+above that limb's radix, so that each limb is bounded by either
+\\(2\^{25+b} \\) or \\( 2\^{26+b} \\), as appropriate.
 
-$$
-\begin{aligned}
-S\_5 &\gets S\_1 + S\_2 \\\\
-S\_6 &\gets S\_1 - S\_2 \\\\
-S\_7 &\gets 2S\_3 \\\\
-S\_8 &\gets S\_7 + S\_6 = S\_1 + 2S\_3 - S\_2 \\\\
-S\_9 &\gets S\_5 - S\_4 = S\_1 + S\_2 - S\_4
-\end{aligned}
-$$
-
-Performing too many intermediate additions and subtractions grows
-the bounds beyond what is allowed as input to multiplication,
-forcing an extra carry pass.  However, it is just possible to avoid
-this by rearranging signs.
-
-Assume that the bounds on the limbs of each field element are
-parameterized by \\( b \in \mathbb R \\) representing the excess
-bits, so that each limb is bounded by either 
-\\( 2\^{25+b} \\) or \\( 2\^{26+b} \\).
-
-The multiplication routine requires that its inputs are bounded by
+The multiplication routine requires that its inputs are bounded with
 \\( b < 1.75 \\), in order to fit a multiplication by \\( 19 \\)
 into 32 bits.  Since \\( \lg 19 < 4.25 \\), \\( 19x < 2\^{32} \\)
 when \\( x < 2\^{27.75} = 2\^{26 + 1.75} \\).  However, this is only
 required for one of the inputs; the other can grow up to \\( b < 2.5
 \\).
 
-Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as
+In addition, the multiplication and squaring routines do not
+canonically reduce their outputs, but can leave some small uncarried
+excesses, so that their reduced outputs are bounded with
+\\( b < 0.007 \\).
 
+The non-parallel portion of the doubling formulas is
+$$
+\begin{aligned}
+(S\_5 &&,&& S\_6 &&,&& S\_8 &&,&& S\_9 )
+&\gets
+(S\_1 + S\_2 &&,&& S\_1 - S\_2 &&,&& S\_1 + 2S\_3 - S\_2 &&,&& S\_1 + S\_2 - S\_4)
+\end{aligned}
+$$
+
+Computing \\( (S\_5, S\_6, S\_8, S\_9 ) \\) as
 $$
 \begin{matrix}
  & S\_1 & S\_1 & S\_1 & S\_1 \\\\
@@ -374,24 +312,22 @@ $$
 =& S\_5 & S\_6 & S\_8 & S\_9
 \end{matrix}
 $$
-
-results in bit-excesses \\( (1.00, 1.59, 2.33, 2.00)\\) for
+results in bit-excesses \\( < (1.01, 1.60, 2.33, 2.01)\\) for
 \\( (S\_5, S\_6, S\_8, S\_9 ) \\).  The products we want to compute
 are then
-
 $$
 \begin{aligned}
-X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.00) \\\\
-Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\
-Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\
-T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 2.00)
+X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 2.01) \\\\
+Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\
+Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\
+T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 2.01)
 \end{aligned}
 $$
-
-which are too large.  However, if we flip the sign of \\( S\_4 =
-S\_0\^2 \\) during squaring, so that we output \\(S\_4' = -S\_4
-\pmod p\\), then we can compute
-
+which are too large: it's not possible to arrange the multiplicands so
+that one vector has \\(b < 2.5\\) and the other has \\( b < 1.75 \\).
+However, if we flip the sign of \\( S\_4 = S\_0\^2 \\) during
+squaring, so that we output \\(S\_4' = -S\_4 \pmod p\\), then we can
+compute
 $$
 \begin{matrix}
  & S\_1 & S\_1 & S\_1 & S\_1 \\\\
@@ -404,61 +340,120 @@ $$
 =& S\_5 & S\_6 & S\_8 & S\_9
 \end{matrix}
 $$
-
-resulting in bit-excesses \\( (1.00, 1.59, 2.33, 1.59)\\) for
+resulting in bit-excesses \\( < (1.01, 1.60, 2.33, 1.60)\\) for
 \\( (S\_5, S\_6, S\_8, S\_9 ) \\).  The products we want to compute
 are then
-
 $$
 \begin{aligned}
-X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.59) \\\\
-Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.00, 1.59) \\\\
-Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.59) \\\\
-T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.00, 1.59)
+X\_3 &\gets S\_8 S\_9 \leftrightarrow (2.33, 1.60) \\\\
+Y\_3 &\gets S\_5 S\_6 \leftrightarrow (1.01, 1.60) \\\\
+Z\_3 &\gets S\_8 S\_6 \leftrightarrow (2.33, 1.60) \\\\
+T\_3 &\gets S\_5 S\_9 \leftrightarrow (1.01, 1.60)
 \end{aligned}
 $$
-
 whose right-hand sides are all bounded with \\( b < 1.75 \\) and
-whose left-hand sides are all bounded with \\( b < 2.5 \\).
+whose left-hand sides are all bounded with \\( b < 2.5 \\),
+so that we can avoid any intermediate reductions.
 
 # Comparison to non-vectorized formulas
 
-HWCD also suggest using a mixed representation, passing between \\(
-\mathbb P\^3 \\) "extended" coordinates and \\( \mathbb P\^2 \\)
-"projective" coordinates, where doubling is slightly cheaper (saving
-about \\(\mathbf 1M\\).  This approach is used for the
-non-vectorized `u32` and `u64` backends, and more
-details on the different coordinate systems can be found in the
-`curve_models` module documentation.
+In theory, the parallel Edwards formulas seem to allow a \\(4\\)-way
+speedup from parallelism.  However, an actual vectorized
+implementation has several slowdowns that cut into this speedup.
 
-This optimization is not compatible with the parallel formulas, which are
-therefore slightly less efficient when counting the total number of
-field multiplications and squarings.  In particular, vectorized doublings
-are less efficient than serial doublings.
-
-In addition, the parallel formulas can only use a \\( 32 \times 32
+First, the parallel formulas can only use a \\( 32 \times 32
 \rightarrow 64 \\)-bit integer multiplier, so the speedup from
 vectorization must overcome the disadvantage of losing the \\( 64
-\times 64 \rightarrow 128\\)-bit (serial) integer multiplier.
+\times 64 \rightarrow 128\\)-bit (serial) integer multiplier.  The
+effect of this slowdown is microarchitecture-dependent, since it
+requires accounting for the total number of multiplications and
+additions and their relative costs.  In the future, it will probably
+be possible to avoid this slowdown by using the `IFMA52` instructions,
+whose parallelism is perfectly suited to these formulas.
 
-When compiling with AVX512VL, LLVM is able to use the extra
-`ymm16..ymm31` registers to reduce register pressure, and avoid
-spills during field multiplication.  This gives a small but
-noticeable speedup.
+Second, the parallel doubling formulas incur both a theoretical and
+practical slowdown.  The parallel formulas described above work on the
+\\( \mathbb P\^3 \\) “extended” coordinates.  The \\( \mathbb P\^2 \\)
+model introduced earlier by [Bernstein, Birkner, Joye, Lange, and
+Peters][bbjlp08] allows slightly faster doublings, so HWCD suggest
+mixing coordinate systems while performing scalar multiplication
+(attributing the idea to [a 1998 paper][cmo98] by Cohen, Miyagi, and
+Ono).  The \\( T \\) coordinate is not required for doublings, so when
+doublings are followed by doublings, its computation can be skipped.
+More details on this approach and the different coordinate systems can
+be found in the [`curve_models` module documentation][curve_models].
 
-Another concern with AVX2 is that currently-available Intel processors
-(particularly Skylake and Skylake-X microarchitectures) perform thermal
-throttling when using wide vector instructions.  For a mixed workload,
+Unfortunately, this optimization is not compatible with the parallel
+formulas, which cannot save time by skipping a single variable, so the
+parallel doubling formulas do slightly more work when counting the
+total number of field multiplications and squarings.
+
+In addition, the parallel doubling formulas have a less regular
+pattern of additions and subtractions than the parallel addition
+formulas, so the vectorization overhead is proportionately greater.
+Both the parallel addition and parallel doubling formulas also require
+some shuffling to rearrange data within the vectors, which places more
+pressure on the shuffle unit than is desirable.
+
+This means that the speedup from using a vectorized implementation of
+parallel Edwards formulas is likely to be greatest in applications
+that do fewer doublings and more additions (like a large multiscalar
+multiplication) rather than applications that do fewer additions and
+more doublings (like a double-base scalar multiplication).
+
+Third, current Intel CPUs perform thermal throttling when using wide
+vector instructions.  A detailed description can be found in §15.26 of
+[the Intel Optimization Manual][intel], but using wide vector
+instructions prevents the core from operating at higher frequencies.
+The core can return to the higher-frequency state after 2
+milliseconds, but this timer is reset every time high-power
+instructions are used.
+
+Any speedup from vectorization therefore has to be weighed against a
+slowdown for the next few million instructions.  For a mixed workload,
 where point operations are interspersed with other tasks, this can
-reduce overall performance.  This probably means that this
-implementation is not suitable for basic applications, like signatures,
-but could still be worthwhile for complex applications, like
-zero-knowledge proofs, which do enough work to make it worthwhile.
+reduce overall performance.  This implementation is therefore probably
+not suitable for basic applications, like signatures, but is
+worthwhile for complex applications, like zero-knowledge proofs, which
+do sustained work.
+
+For this reason, the AVX2 backend is not enabled by default, but can
+be selected using the `avx2_backend` feature.
+
+# Future work
+
+There are several directions for future improvement:
+
+* Using the vectorized field arithmetic code to parallelize across
+  point operations rather than within a single point operation.  This
+  is less flexible, but would give a speedup both from allowing use of
+  the faster mixed-model arithmetic and from reducing shuffle
+  pressure.  One approach in this direction would be to implement
+  batched scalar-point operations using vectors of points (AoSoA
+  layout).  This less generally useful but would give a speedup for
+  Bulletproofs.
+
+* Extending the implementation to use the full width of AVX512, either
+  handling the extra parallelism internally to a single point
+  operation (by using a 2-way parallel implementation of field
+  arithmetic instead of a wordsliced one), or externally,
+  parallelizing across point operations.  Internal parallelism would
+  be preferable but might require too much shuffle pressure.
+
+* Generalizing the implementation to non-AVX2 instructions,
+  particularly NEON.  The current point arithmetic code is written in
+  terms of field element vectors, which are in turn implemented using
+  platform SIMD vectors.  It should be possible to write an alternate
+  implementation of the `FieldElement32x4` using NEON without changing
+  the point arithmetic.  NEON has 128-bit vectors rather than 256-bit
+  vectors, but this may still be worthwhile compared to a serial
+  implementation.
 
-On AMD's Zen microarchitecture, thermal throttling is not a concern,
-since AVX2 is implemented at half rate, so there is no penalty for mixed
-workloads (but also no speedup).
 
 [sandy2x]: https://eprint.iacr.org/2015/943.pdf
 [avx2trac]: https://trac.torproject.org/projects/tor/ticket/8897#comment:28
 [hwcd08]: https://www.iacr.org/archive/asiacrypt2008/53500329/53500329.pdf
+[curve_models]: https://doc-internal.dalek.rs/curve25519_dalek/curve_models/index.html
+[bbjlp08]: https://eprint.iacr.org/2008/013
+[cmo98]: https://link.springer.com/content/pdf/10.1007%2F3-540-49649-1_6.pdf
+[intel]: https://software.intel.com/sites/default/files/managed/9e/bc/64-ia-32-architectures-optimization-manual.pdf
\ No newline at end of file
diff --git a/src/backend/avx2/edwards.rs b/src/backend/avx2/edwards.rs
index a053414..43f02f7 100644
--- a/src/backend/avx2/edwards.rs
+++ b/src/backend/avx2/edwards.rs
@@ -8,7 +8,29 @@
 // - Isis Agora Lovecruft <isis@patternsinthevoid.net>
 // - Henry de Valence <hdevalence@hdevalence.ca>
 
-//! Extended Twisted Edwards for Curve25519, using AVX2.
+//! Parallel Edwards Arithmetic for Curve25519.
+//!
+//! This module currently has two point types:
+//!
+//! * `ExtendedPoint`: a point stored in vector-friendly format, with
+//! vectorized doubling and addition;
+//!
+//! * `CachedPoint`: used for readdition.
+//!
+//! Details on the formulas can be found in the documentation for the
+//! parent `avx2` module.
+//!
+//! This API is designed to be safe: vectorized points can only be
+//! created from serial points (which do validation on decompression),
+//! and operations on valid points return valid points, so invalid
+//! point states should be unrepresentable.
+//!
+//! This design goal is met, with one exception: the `Neg`
+//! implementation for the `CachedPoint` performs a lazy negation, so
+//! that subtraction can be efficiently implemented as a negation and
+//! an addition.  Repeatedly negating a `CachedPoint` will cause its
+//! coefficients to grow and eventually overflow.  Repeatedly negating
+//! a point should not be necessary anyways.
 
 #![allow(non_snake_case)]
 
@@ -26,7 +48,13 @@ use traits::Identity;
 use backend::avx2::field::{FieldElement32x4, Lanes, Shuffle};
 use backend::avx2::constants;
 
-/// A point on Curve25519, represented in an AVX2-friendly format.
+/// A point on Curve25519, using parallel Edwards formulas for curve
+/// operations.
+///
+/// # Invariant
+///
+/// The coefficients of an `ExtendedPoint` are bounded with
+/// \\( b < 0.007 \\).
 #[derive(Copy, Clone, Debug)]
 pub struct ExtendedPoint(pub(super) FieldElement32x4);
 
@@ -67,6 +95,7 @@ impl Identity for ExtendedPoint {
 }
 
 impl ExtendedPoint {
+    /// Compute the double of this point.
     pub fn double(&self) -> ExtendedPoint {
         // Want to compute (X1 Y1 Z1 X1+Y1).
         // Not sure how to do this less expensively than computing
@@ -138,6 +167,15 @@ impl ExtendedPoint {
 }
 
 /// A cached point with some precomputed variables used for readdition.
+///
+/// # Warning
+///
+/// It is not safe to negate this point more than once.
+///
+/// # Invariant
+///
+/// As long as the `CachedPoint` is not repeatedly negated, its
+/// coefficients will be bounded with \\( b < 1.0 \\).
 #[derive(Copy, Clone, Debug)]
 pub struct CachedPoint(pub(super) FieldElement32x4);
 
@@ -196,6 +234,12 @@ impl<'a, 'b> Add<&'b CachedPoint> for &'a ExtendedPoint {
 
     /// Add an `ExtendedPoint` and a `CachedPoint`.
     fn add(self, other: &'b CachedPoint) -> ExtendedPoint {
+        // The coefficients of an `ExtendedPoint` are reduced after
+        // every operation.  If the `CachedPoint` was negated, its
+        // coefficients grow by one bit.  So on input, `self` is
+        // bounded with `b < 0.007` and `other` is bounded with
+        // `b < 1.0`.
+        
         let mut tmp = self.0;
 
         tmp = tmp.blend(tmp.diff_sum(), Lanes::AB);
diff --git a/src/backend/avx2/field.rs b/src/backend/avx2/field.rs
index b5bf4b3..175ff14 100644
--- a/src/backend/avx2/field.rs
+++ b/src/backend/avx2/field.rs
@@ -8,7 +8,20 @@
 // - Isis Agora Lovecruft <isis@patternsinthevoid.net>
 // - Henry de Valence <hdevalence@hdevalence.ca>
 
-//! 4-way vectorized 32bit field arithmetic using AVX2.
+//! An implementation of 4-way vectorized 32bit field arithmetic using
+//! AVX2.
+//!
+//! The `FieldElement32x4` struct provides a vector of four field
+//! elements, implemented using AVX2 operations.  Its API is designed
+//! to abstract away the platform-dependent details, so that point
+//! arithmetic can be implemented only in terms of a vector of field
+//! elements.
+//!
+//! At this level, the API is optimized for speed and not safety.  The
+//! `FieldElement32x4` does not always perform reductions.  The pre-
+//! and post-conditions on the bounds of the coefficients are
+//! documented for each method, but it is the caller's responsibility
+//! to ensure that there are no overflows.
 
 #![allow(non_snake_case)]
 
@@ -32,7 +45,7 @@ use core::simd::{i32x8, u32x8, u64x4, IntoBits};
 use backend::avx2::constants::{P_TIMES_16_HI, P_TIMES_16_LO, P_TIMES_2_HI, P_TIMES_2_LO};
 use backend::u64::field::FieldElement64;
 
-/// Unpack
+/// Unpack 32-bit lanes into 64-bit lanes:
 /// ```
 /// (a0, b0, a1, b1, c0, d0, c1, d1)
 /// ```
@@ -55,7 +68,7 @@ fn unpack_pair(src: u32x8) -> (u32x8, u32x8) {
     (a, b)
 }
 
-/// Repack
+/// Repack 64-bit lanes into 32-bit lanes:
 /// ```
 /// (a0, 0, b0, 0, c0, 0, d0, 0)
 /// (a1, 0, b1, 0, c1, 0, d1, 0)
@@ -102,6 +115,10 @@ pub enum Lanes {
 }
 
 /// The `Shuffle` enum represents a shuffle of a `FieldElement32x4`.
+///
+/// The enum variants are named by what they do to a vector \\(
+/// (A,B,C,D) \\); for instance, `Shuffle::BADC` turns \\( (A, B, C,
+/// D) \\) into \\( (B, A, D, C) \\).
 #[derive(Copy, Clone, Debug)]
 pub enum Shuffle {
     AAAA,
@@ -116,7 +133,13 @@ pub enum Shuffle {
     ABDC,
 }
 
-/// A vector of four field elements, in an AVX2-friendly format.
+/// A vector of four field elements.
+///
+/// Each operation on a `FieldElement32x4` has documented effects on
+/// the bounds of the coefficients.  This API is designed for speed
+/// and not safety; it is the caller's responsibility to ensure that
+/// the post-conditions of one operation are compatible with the
+/// pre-conditions of the next.
 #[derive(Clone, Copy, Debug)]
 pub struct FieldElement32x4(pub(crate) [u32x8; 5]);
 
@@ -134,6 +157,8 @@ impl ConditionallyAssignable for FieldElement32x4 {
 }
 
 impl FieldElement32x4 {
+    /// Split this vector into an array of four (serial) field
+    /// elements.
     pub fn split(&self) -> [FieldElement64; 4] {
         let mut out = [FieldElement64::zero(); 4];
         for i in 0..5 {
@@ -155,6 +180,11 @@ impl FieldElement32x4 {
         out
     }
 
+    /// Rearrange the elements of this vector according to `control`.
+    ///
+    /// The `control` parameter should be a compile-time constant, so
+    /// that when this function is inlined, LLVM is able to lower the
+    /// shuffle using an immediate.
     #[inline]
     pub fn shuffle(&self, control: Shuffle) -> FieldElement32x4 {
         #[inline(always)]
@@ -190,8 +220,13 @@ impl FieldElement32x4 {
         ])
     }
 
+    /// Blend `self` with `other`, taking lanes specified in `control` from `other`.
+    ///
+    /// The `control` parameter should be a compile-time constant, so
+    /// that this function can be inlined and LLVM can lower it to a
+    /// blend instruction using an immediate.
     #[inline]
-    pub fn blend(&self, b: FieldElement32x4, control: Lanes) -> FieldElement32x4 {
+    pub fn blend(&self, other: FieldElement32x4, control: Lanes) -> FieldElement32x4 {
         #[inline(always)]
         fn blend_lanes(x: u32x8, y: u32x8, control: Lanes) -> u32x8 {
             unsafe {
@@ -254,11 +289,11 @@ impl FieldElement32x4 {
         }
 
         FieldElement32x4([
-            blend_lanes(self.0[0], b.0[0], control),
-            blend_lanes(self.0[1], b.0[1], control),
-            blend_lanes(self.0[2], b.0[2], control),
-            blend_lanes(self.0[3], b.0[3], control),
-            blend_lanes(self.0[4], b.0[4], control),
+            blend_lanes(self.0[0], other.0[0], control),
+            blend_lanes(self.0[1], other.0[1], control),
+            blend_lanes(self.0[2], other.0[2], control),
+            blend_lanes(self.0[3], other.0[3], control),
+            blend_lanes(self.0[4], other.0[4], control),
         ])
     }