diff --git a/Cargo.toml b/Cargo.toml
index c164472..484a117 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "curve25519-dalek"
-version = "0.6.0"
+version = "0.7.0"
 authors = ["Isis Lovecruft <isis@patternsinthevoid.net>",
            "Henry de Valence <hdevalence@hdevalence.ca>"]
 readme = "README.md"
@@ -36,12 +36,10 @@ version = "^0.6"
 version = "0.4"
 
 [features]
-nightly = ["basepoint_table_creation", "radix_51"]
+nightly = ["radix_51"]
 default = ["std"]
 std = ["rand"]
 yolocrypto = []
-# Needs nightly for placement new
-basepoint_table_creation = []
 bench = []
 # Radix-51 arithmetic using u128
 radix_51 = []
diff --git a/README.md b/README.md
index ef10f59..16bad4e 100644
--- a/README.md
+++ b/README.md
@@ -44,7 +44,7 @@ Extensive documentation is available [here](https://docs.rs/curve25519-dalek).
 To install, add the following to the dependencies section of your project's
 `Cargo.toml`:
 
-    curve25519-dalek = "^0.6"
+    curve25519-dalek = "^0.7"
 
 Then, in your library or executable source, add:
 
diff --git a/src/constants.rs b/src/constants.rs
index 9ff4d4d..2d3c790 100644
--- a/src/constants.rs
+++ b/src/constants.rs
@@ -24,6 +24,8 @@ use curve::ExtendedPoint;
 use curve::AffineNielsPoint;
 use curve::CompressedEdwardsY;
 use curve::EdwardsBasepointTable;
+#[cfg(feature = "yolocrypto")]
+use decaf::{DecafPoint, DecafBasepointTable};
 use scalar::Scalar;
 
 #[cfg(feature="radix_51")]
@@ -70,7 +72,7 @@ pub const SQRT_M1: FieldElement = FieldElement([
 pub const SQRT_M1: FieldElement = FieldElement([1718705420411056, 234908883556509, 2233514472574048, 2117202627021982, 765476049583133]);
 
 /// Precomputed value of the other square root of -1 (mod p),
-/// i.e., MSQRT_M1 = -SQRT_M1.
+/// i.e., `MSQRT_M1 = -SQRT_M1`.
 #[cfg(not(feature="radix_51"))]
 pub const MSQRT_M1: FieldElement = FieldElement([
     32595792,    7943725,  -9377950,  -3500415, -12389472,
@@ -92,7 +94,7 @@ pub const A: FieldElement       = FieldElement([
 #[cfg(feature="radix_51")]
 pub const A: FieldElement = FieldElement([486662, 0, 0, 0, 0]);
 
-/// SQRT_MINUS_A is sqrt(-486662)
+/// `SQRT_MINUS_A` is sqrt(-486662)
 // XXX I think that this was used in Adam's code for his elligator
 // implementation, but that should maybe be using sqrt(-486664)
 // instead...?  - hdevalence
@@ -103,7 +105,7 @@ pub const SQRT_MINUS_A: FieldElement = FieldElement([ // sqrtMinusA
 #[cfg(feature="radix_51")]
 pub const SQRT_MINUS_A: FieldElement = FieldElement([557817479725543, 1643290402203250, 16226468853936, 1304118542701054, 1985241807451647]);
 
-/// SQRT_MINUS_APLUS2 is sqrt(-486664)
+/// `SQRT_MINUS_APLUS2` is sqrt(-486664)
 #[cfg(not(feature="radix_51"))]
 pub const SQRT_MINUS_APLUS2: FieldElement = FieldElement([
     -12222970, -8312128, -11511410, 9067497, -15300785,
@@ -111,7 +113,7 @@ pub const SQRT_MINUS_APLUS2: FieldElement = FieldElement([
 #[cfg(feature="radix_51")]
 pub const SQRT_MINUS_APLUS2: FieldElement = FieldElement([1693982333959686, 608509411481997, 2235573344831311, 947681270984193, 266558006233600]);
 
-/// SQRT_MINUS_HALF is sqrt(-1/2)
+/// `SQRT_MINUS_HALF` is sqrt(-1/2)
 #[cfg(not(feature="radix_51"))]
 pub const SQRT_MINUS_HALF: FieldElement = FieldElement([ // sqrtMinusHalf
     -17256545,   3971863,  28865457,  -1750208,  27359696,
@@ -119,7 +121,7 @@ pub const SQRT_MINUS_HALF: FieldElement = FieldElement([ // sqrtMinusHalf
 #[cfg(feature="radix_51")]
 pub const SQRT_MINUS_HALF: FieldElement = FieldElement([266547196637087, 2134345371906993, 1135042577398223, 67298593331632, 743161882051057]);
 
-/// HALF_Q_MINUS_1_BYTES is (2^255-20)/2 expressed in little endian form.
+/// `HALF_Q_MINUS_1_BYTES` is (2^255-20)/2 expressed in little endian form.
 pub const HALF_Q_MINUS_1_BYTES: [u8; 32] = [ // halfQMinus1Bytes
     0xf6, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
     0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
@@ -136,6 +138,10 @@ pub const BASE_CMPRSSD: CompressedEdwardsY =
                         0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66,
                         0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66, 0x66]);
 
+/// The Ed25519 basepoint, as a `DecafPoint`.
+#[cfg(feature = "yolocrypto")]
+pub const DECAF_ED25519_BASEPOINT: DecafPoint = DecafPoint(ED25519_BASEPOINT);
+
 /// Basepoint has y = 4/5.
 #[cfg(not(feature="radix_51"))]
 pub const ED25519_BASEPOINT: ExtendedPoint = ExtendedPoint{
@@ -159,12 +165,20 @@ pub const l: Scalar = Scalar([ 0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58,
                                0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
                                0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]);
 
-/// `lminus1` is the order of base point minus one, i.e. 2^252 +
+/// `l_minus_1` is the order of base point minus one, i.e. 2^252 +
 /// 27742317777372353535851937790883648493 - 1, in little-endian form
-pub const lminus1: Scalar = Scalar([ 0xec, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58,
-                                     0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,
-                                     0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-                                     0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]);
+pub const l_minus_1: Scalar = Scalar([ 0xec, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58,
+                                       0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,
+                                       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+                                       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]);
+
+/// `lminus1` is the order of base point minus two, i.e. 2^252 +
+/// 27742317777372353535851937790883648493 - 2, in little-endian form
+pub const l_minus_2: Scalar = Scalar([ 0xeb, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58,
+                                       0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,
+                                       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+                                       0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10 ]);
+
 /// The 8-torsion subgroup Ɛ[8].
 ///
 /// In the case of Curve25519, it is cyclic; the `i`th element of the
@@ -376,6 +390,11 @@ pub const bi: [AffineNielsPoint; 8] = [
     }
 ];
 
+#[cfg(feature = "yolocrypto")]
+/// The Ed25519 basepoint
+pub const DECAF_ED25519_BASEPOINT_TABLE: DecafBasepointTable
+    = DecafBasepointTable(ED25519_BASEPOINT_TABLE);
+
 /// Table containing precomputed multiples of the basepoint `B = (x,4/5)`.
 ///
 /// The table is defined so `constants::base[i][j-1] = j*(16^2i)*B`,
diff --git a/src/curve.rs b/src/curve.rs
index aefecac..46358ad 100644
--- a/src/curve.rs
+++ b/src/curve.rs
@@ -59,15 +59,15 @@
 //! implementation for [Ed25519](https://ed25519.cr.yp.to/ed25519-20110926.pdf),
 //! we use several different models for curve points:
 //!
-//! * CompletedPoint: points in 𝗣^1 x 𝗣^1;
-//! * ExtendedPoint: points in 𝗣^3;
-//! * ProjectivePoint: points in 𝗣^2.
+//! * `CompletedPoint`: points in 𝗣^1 x 𝗣^1;
+//! * `ExtendedPoint`: points in 𝗣^3;
+//! * `ProjectivePoint`: points in 𝗣^2.
 //!
 //! Finally, to accelerate additions, we use two cached point formats,
 //! one for the affine model and one for the 𝗣^3 model:
 //!
-//! * AffineNielsPoint: `(y+x, y-x, 2dxy)`
-//! * ProjectiveNielsPoint: `(Y+X, Y-X, Z, 2dXY)`
+//! * `AffineNielsPoint`: `(y+x, y-x, 2dxy)`
+//! * `ProjectiveNielsPoint`: `(Y+X, Y-X, Z, 2dXY)`
 //!
 //! [1]: https://moderncrypto.org/mail-archive/curves/2016/000807.html
 
@@ -80,8 +80,12 @@
 use core::fmt::Debug;
 use core::iter::Iterator;
 use core::ops::{Add, Sub, Neg};
+use core::ops::{Mul, MulAssign};
+use core::ops::Index;
 
 use constants;
+#[cfg(feature = "yolocrypto")]
+use decaf::DecafPoint;
 use field::FieldElement;
 use scalar::Scalar;
 use subtle::arrays_equal_ct;
@@ -90,11 +94,6 @@ use subtle::CTAssignable;
 use subtle::CTEq;
 use subtle::CTNegatable;
 
-#[cfg(all(not(feature = "std"), feature = "basepoint_table_creation"))]
-use collections::boxed::Box;
-#[cfg(all(feature = "std", feature = "basepoint_table_creation"))]
-use std::boxed::Box;
-
 // ------------------------------------------------------------------------
 // Compressed points
 // ------------------------------------------------------------------------
@@ -103,7 +102,7 @@ use std::boxed::Box;
 /// determined by the `y`-coordinate and the sign of `x`, marshalled
 /// into a 32-byte array.
 ///
-/// The first 255 bits of a CompressedEdwardsY represent the
+/// The first 255 bits of a `CompressedEdwardsY` represent the
 /// y-coordinate. The high bit of the 32nd byte gives the sign of `x`.
 #[derive(Copy, Clone, Eq, PartialEq)]
 pub struct CompressedEdwardsY(pub [u8; 32]);
@@ -116,7 +115,7 @@ impl Debug for CompressedEdwardsY {
 
 impl CompressedEdwardsY {
     /// View this `CompressedEdwardsY` as an array of bytes.
-    pub fn as_bytes<'a>(&'a self) -> &'a [u8; 32] {
+    pub fn as_bytes(&self) -> &[u8; 32] {
         &self.0
     }
 
@@ -232,7 +231,7 @@ impl CompressedMontgomeryU {
     /// Montgomery `v` corresponding to this `u`.
     pub fn to_montgomery_v(u: &FieldElement) -> (u8, FieldElement) {
         let one:       FieldElement = FieldElement::one();
-        let v_squared: FieldElement = u * &(&(&u.square() + &(&(&constants::A * u) + &one)));
+        let v_squared: FieldElement = u * &(&u.square() + &(&(&constants::A * u) + &one));
 
         let (okay, v_inv) = v_squared.invsqrt();
         let v = &v_inv * &v_squared;
@@ -303,7 +302,7 @@ pub struct ProjectivePoint {
     Z: FieldElement,
 }
 
-/// A CompletedPoint is a point ((X:Z), (Y:T)) in 𝗣¹(𝔽ₚ)×𝗣¹(𝔽ₚ).
+/// A `CompletedPoint` is a point ((X:Z), (Y:T)) in 𝗣¹(𝔽ₚ)×𝗣¹(𝔽ₚ).
 /// A point (x,y) in the affine model corresponds to ((x:1),(y:1)).
 #[derive(Copy, Clone)]
 pub struct CompletedPoint {
@@ -448,6 +447,15 @@ impl CTAssignable for AffineNielsPoint {
     }
 }
 
+impl CTAssignable for ExtendedPoint {
+    fn conditional_assign(&mut self, other: &ExtendedPoint, choice: u8) {
+        self.X.conditional_assign(&other.X, choice);
+        self.Y.conditional_assign(&other.Y, choice);
+        self.Z.conditional_assign(&other.Z, choice);
+        self.T.conditional_assign(&other.T, choice);
+    }
+}
+
 // ------------------------------------------------------------------------
 // Constant-time Equality
 // ------------------------------------------------------------------------
@@ -470,13 +478,7 @@ pub trait IsIdentity {
 /// constructor.
 impl<T> IsIdentity for T where T: CTEq + Identity {
     fn is_identity(&self) -> bool {
-        let identity: T = T::identity();
-
-        if self.ct_eq(&identity) == 1u8 {
-            return true;
-        } else {
-            return false;
-        }
+        self.ct_eq(&T::identity()) == 1u8
     }
 }
 
@@ -795,18 +797,20 @@ impl<'a> Neg for &'a AffineNielsPoint {
 // Scalar multiplication
 // ------------------------------------------------------------------------
 
-/// Trait for scalar multiplication of an arbitrary point.
-pub trait ScalarMult<S> {
-    /// Compute `scalar * self`.
-    fn scalar_mult(&self, scalar: &S) -> Self;
+impl<'b> MulAssign<&'b Scalar> for ExtendedPoint {
+    fn mul_assign(&mut self, scalar: &'b Scalar) {
+        let result = (self as &ExtendedPoint) * scalar;
+        *self = result;
+    }
 }
 
-impl ScalarMult<Scalar> for ExtendedPoint {
+impl<'a, 'b> Mul<&'b Scalar> for &'a ExtendedPoint {
+    type Output = ExtendedPoint;
     /// Scalar multiplication: compute `scalar * self`.
     ///
     /// Uses a window of size 4.  Note: for scalar multiplication of
     /// the basepoint, `basepoint_mult` is approximately 4x faster.
-    fn scalar_mult(&self, scalar: &Scalar) -> ExtendedPoint {
+    fn mul(self, scalar: &'b Scalar) -> ExtendedPoint {
         let A = self.to_projective_niels();
         let mut As: [ProjectiveNielsPoint; 8] = [A; 8];
         for i in 0..7 {
@@ -824,31 +828,35 @@ impl ScalarMult<Scalar> for ExtendedPoint {
     }
 }
 
+impl<'a, 'b> Mul<&'b ExtendedPoint> for &'a Scalar {
+    type Output = ExtendedPoint;
+
+    /// Scalar multiplication: compute `self * point`.
+    ///
+    /// Uses a window of size 4.  Note: for scalar multiplication of
+    /// the basepoint, `basepoint_mult` is approximately 4x faster.
+    fn mul(self, point: &'b ExtendedPoint) -> ExtendedPoint {
+        point * &self
+    }
+}
+
+#[cfg(feature = "yolocrypto")]
+impl<'a, 'b> Mul<&'b DecafPoint> for &'a Scalar {
+    type Output = DecafPoint;
+
+    /// Scalar multiplication: compute `self * scalar`.
+    fn mul(self, point: &'b DecafPoint) -> DecafPoint {
+        DecafPoint(self * &point.0)
+    }
+}
+
+
 /// Precomputation
 #[derive(Clone)]
 pub struct EdwardsBasepointTable(pub [[AffineNielsPoint; 8]; 32]);
 
-impl EdwardsBasepointTable {
-    /// Create a table of precomputed multiples of `basepoint`.
-    #[cfg(feature="basepoint_table_creation")]
-    pub fn create(basepoint: &ExtendedPoint) -> Box<EdwardsBasepointTable> {
-        // Create the table storage
-        // XXX can we be assured that this is not allocated on the stack?
-        // XXX can we skip the initialization without too much unsafety?
-        let mut table = box EdwardsBasepointTable([[AffineNielsPoint::identity(); 8]; 32]);
-        let mut P = basepoint.clone();
-        for i in 0..32 {
-            // P = (16^2)^i * B
-            let mut jP = P.to_affine_niels();
-            for j in 1..9 {
-                // table[i][j-1] is supposed to be j*(16^2)^i*B
-                table.0[i][j-1] = jP;
-                jP = (&P + &jP).to_extended().to_affine_niels();
-            }
-            P = P.mult_by_pow_2(8);
-        }
-        return table
-    }
+impl<'a, 'b> Mul<&'b Scalar> for &'a EdwardsBasepointTable {
+    type Output = ExtendedPoint;
 
     /// Construct an `ExtendedPoint` from a `Scalar`, `scalar`, by
     /// computing the multiple `aB` of the basepoint `B`.
@@ -875,7 +883,7 @@ impl EdwardsBasepointTable {
     /// We then use the `select_precomputed_point` function, which
     /// takes `-8 ≤ x < 8` and `[16^2i * B, ..., 8 * 16^2i * B]`,
     /// and returns `x * 16^2i * B` in constant time.
-    pub fn basepoint_mult(&self, scalar: &Scalar) -> ExtendedPoint {
+    fn mul(self, scalar: &'b Scalar) -> ExtendedPoint {
         let e = scalar.to_radix_16();
         let mut h = ExtendedPoint::identity();
         let mut t: CompletedPoint;
@@ -896,21 +904,65 @@ impl EdwardsBasepointTable {
     }
 }
 
-/// Trait for scalar multiplication of a distinguished basepoint.
-pub trait BasepointMult<S> {
-    /// Return the basepoint `B`.
-    fn basepoint() -> Self;
-    /// Compute `scalar * B`.
-    fn basepoint_mult(scalar: &S) -> Self;
+impl<'a, 'b> Mul<&'a EdwardsBasepointTable> for &'b Scalar {
+    type Output = ExtendedPoint;
+
+    /// Construct an `ExtendedPoint` by via this `Scalar` times
+    /// a the basepoint, `B` included in a precomputed `basepoint_table`.
+    ///
+    /// Precondition: this scalar must be reduced.
+    ///
+    /// The computation proceeds as follows, as described on page 13
+    /// of the Ed25519 paper.  Write this scalar `a` in radix 16 with
+    /// coefficients in [-8,8), i.e.,
+    ///
+    ///    a = a_0 + a_1*16^1 + ... + a_63*16^63,
+    ///
+    /// with -8 ≤ a_i < 8.  Then
+    ///
+    ///    a*B = a_0*B + a_1*16^1*B + ... + a_63*16^63*B.
+    ///
+    /// Grouping even and odd coefficients gives
+    ///
+    ///    a*B =       a_0*16^0*B + a_2*16^2*B + ... + a_62*16^62*B
+    ///              + a_1*16^1*B + a_3*16^3*B + ... + a_63*16^63*B
+    ///        =      (a_0*16^0*B + a_2*16^2*B + ... + a_62*16^62*B)
+    ///          + 16*(a_1*16^0*B + a_3*16^2*B + ... + a_63*16^62*B).
+    ///
+    /// We then use the `select_precomputed_point` function, which
+    /// takes `-8 ≤ x < 8` and `[16^2i * B, ..., 8 * 16^2i * B]`,
+    /// and returns `x * 16^2i * B` in constant time.
+    fn mul(self, basepoint_table: &'a EdwardsBasepointTable) -> ExtendedPoint {
+        basepoint_table * &self
+    }
 }
 
-impl BasepointMult<Scalar> for ExtendedPoint {
-    fn basepoint() -> ExtendedPoint {
-        constants::ED25519_BASEPOINT
+impl EdwardsBasepointTable {
+    /// Create a table of precomputed multiples of `basepoint`.
+    pub fn create(basepoint: &ExtendedPoint) -> EdwardsBasepointTable {
+        // Create the table storage
+        // XXX can we skip the initialization without too much unsafety?
+        // stick 30K on the stack and call it a day.
+        let mut table = EdwardsBasepointTable([[AffineNielsPoint::identity(); 8]; 32]);
+        let mut P = basepoint.clone();
+        for i in 0..32 {
+            // P = (16^2)^i * B
+            let mut jP = P.to_affine_niels();
+            for j in 1..9 {
+                // table[i][j-1] is supposed to be j*(16^2)^i*B
+                table.0[i][j-1] = jP;
+                jP = (&P + &jP).to_extended().to_affine_niels();
+            }
+            P = P.mult_by_pow_2(8);
+        }
+        table
     }
 
-    fn basepoint_mult(scalar: &Scalar) -> ExtendedPoint {
-        constants::ED25519_BASEPOINT_TABLE.basepoint_mult(scalar)
+    /// Get the basepoint for this table as an `ExtendedPoint`.
+    pub fn basepoint(&self) -> ExtendedPoint {
+        // self.0[0][0] has 1*(16^2)^0*B, but as an `AffineNielsPoint`
+        // Add identity to convert to extended.
+        (&ExtendedPoint::identity() + &self.0[0][0]).to_extended()
     }
 }
 
@@ -933,8 +985,7 @@ impl ExtendedPoint {
             r = s.double(); s = r.to_projective();
         }
         // Unroll last iteration so we can go directly to_extended()
-        r = s.double();
-        return r.to_extended();
+        s.double().to_extended()
     }
 
     /// Determine if this point is of small order.
@@ -947,73 +998,10 @@ impl ExtendedPoint {
     ///
     /// True if it is of small order; false otherwise.
     pub fn is_small_order(&self) -> bool {
-        let p8: ExtendedPoint = self.mult_by_pow_2(3);
-
-        if p8.is_identity() {
-            return true;
-        } else {
-            return false;
-        }
+        self.mult_by_cofactor().is_identity()
     }
 }
 
-/// Given a point `A` and scalars `a` and `b`, compute the point
-/// `aA+bB`, where `B` is the Ed25519 basepoint (i.e., `B = (x,4/5)`
-/// with x positive).
-///
-/// # Warning
-///
-/// This function is *not* constant time, hence its name.
-// XXX should return ExtendedPoint?
-pub fn double_scalar_mult_vartime(a: &Scalar, A: &ExtendedPoint, b: &Scalar) -> ProjectivePoint {
-    let a_naf = a.non_adjacent_form();
-    let b_naf = b.non_adjacent_form();
-
-    // Build a lookup table of odd multiples of A
-    let mut Ai = [ProjectiveNielsPoint::identity(); 8];
-    let A2 = A.double();
-    Ai[0]  = A.to_projective_niels();
-    for i in 0..7 {
-        Ai[i+1] = (&A2 + &Ai[i]).to_extended().to_projective_niels();
-    }
-    // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A]
-
-    // Find starting index
-    let mut i: usize = 255;
-    for j in (0..255).rev() {
-        i = j;
-        if a_naf[i] != 0 || b_naf[i] != 0 {
-            break;
-        }
-    }
-
-    let mut r = ProjectivePoint::identity();
-    loop {
-        let mut t = r.double();
-
-        if a_naf[i] > 0 {
-            t = &t.to_extended() + &Ai[( a_naf[i]/2) as usize];
-        } else if a_naf[i] < 0 {
-            t = &t.to_extended() - &Ai[(-a_naf[i]/2) as usize];
-        }
-
-        if b_naf[i] > 0 {
-            t = &t.to_extended() + &constants::bi[( b_naf[i]/2) as usize];
-        } else if b_naf[i] < 0 {
-            t = &t.to_extended() - &constants::bi[(-b_naf[i]/2) as usize];
-        }
-
-        r = t.to_projective();
-
-        if i == 0 {
-            break;
-        }
-        i -= 1;
-    }
-
-    r
-}
-
 /// Given precomputed points `[P, 2P, 3P, ..., 8P]`, as well as `-8 ≤
 /// x ≤ 8`, compute `x * B` in constant time, i.e., without branching
 /// on x or using it as an array index.
@@ -1053,14 +1041,14 @@ impl ExtendedPoint {
     /// Returns `Some<[u8;32]>` if `self` is in the image of the
     /// Elligator2 map.  For a random point on the curve, this happens
     /// with probability 1/2.  Otherwise, returns `None`.
-    pub fn to_uniform_representative(&self) -> Option<[u8;32]> {
+    pub fn to_uniform_representative(&self) -> Option<[u8; 32]> {
         unimplemented!();
     }
 
     /// Use Elligator2 to convert a uniformly random string to a curve
     /// point.
     #[allow(unused_variables)] // REMOVE WHEN IMPLEMENTED
-    pub fn from_uniform_representative(bytes: &[u8;32]) -> ExtendedPoint {
+    pub fn from_uniform_representative(bytes: &[u8; 32]) -> ExtendedPoint {
         unimplemented!();
     }
 }
@@ -1104,12 +1092,133 @@ impl Debug for ProjectiveNielsPoint {
     }
 }
 
+// ------------------------------------------------------------------------
+// Variable-time functions
+// ------------------------------------------------------------------------
+
+pub mod vartime {
+    //! Variable-time operations on curve points, useful for non-secret data.
+    use super::*;
+
+    /// Holds odd multiples 1A, 3A, ..., 15A of a point A.
+    struct OddMultiples([ProjectiveNielsPoint; 8]);
+
+    impl OddMultiples {
+        fn create(A: &ExtendedPoint) -> OddMultiples {
+            let mut Ai = [ProjectiveNielsPoint::identity(); 8];
+            let A2 = A.double();
+            Ai[0]  = A.to_projective_niels();
+            for i in 0..7 {
+                Ai[i+1] = (&A2 + &Ai[i]).to_extended().to_projective_niels();
+            }
+            // Now Ai = [A, 3A, 5A, 7A, 9A, 11A, 13A, 15A]
+            OddMultiples(Ai)
+        }
+    }
+
+    impl Index<usize> for OddMultiples {
+        type Output = ProjectiveNielsPoint;
+
+        fn index<'a>(&'a self, _index: usize) -> &'a ProjectiveNielsPoint {
+            &(self.0[_index])
+        }
+    }
+
+    /// Given a vector of public scalars and a vector of (possibly secret)
+    /// points, compute
+    ///
+    ///    c_1 P_1 + ... + c_n P_n.
+    ///
+    /// # Input
+    ///
+    /// A vector of `Scalar`s and a vector of `ExtendedPoints`.  It is an
+    /// error to call this function with two vectors of different lengths.
+    pub fn k_fold_scalar_mult<'a,'b,I,J>(scalars: I, points: J) -> ExtendedPoint
+        where I: IntoIterator<Item=&'a Scalar>, J: IntoIterator<Item=&'b ExtendedPoint>
+    {
+        //assert_eq!(scalars.len(), points.len());
+
+        let nafs: Vec<_> = scalars.into_iter()
+            .map(|c| c.non_adjacent_form()).collect();
+        let odd_multiples: Vec<_> = points.into_iter()
+            .map(|P| OddMultiples::create(P)).collect();
+
+        let mut r = ProjectivePoint::identity();
+
+        for i in (0..255).rev() {
+            let mut t = r.double();
+
+            for (naf, odd_multiple) in nafs.iter().zip(odd_multiples.iter()) {
+                if naf[i] > 0 {
+                    t = &t.to_extended() + &odd_multiple[( naf[i]/2) as usize];
+                } else if naf[i] < 0 {
+                    t = &t.to_extended() - &odd_multiple[(-naf[i]/2) as usize];
+                }
+            }
+
+            r = t.to_projective();
+        }
+
+        r.to_extended()
+    }
+
+    /// Given a point `A` and scalars `a` and `b`, compute the point
+    /// `aA+bB`, where `B` is the Ed25519 basepoint (i.e., `B = (x,4/5)`
+    /// with x positive).
+    pub fn double_scalar_mult_basepoint(a: &Scalar,
+                                        A: &ExtendedPoint,
+                                        b: &Scalar) -> ProjectivePoint {
+        let a_naf = a.non_adjacent_form();
+        let b_naf = b.non_adjacent_form();
+
+        // Find starting index
+        let mut i: usize = 255;
+        for j in (0..255).rev() {
+            i = j;
+            if a_naf[i] != 0 || b_naf[i] != 0 {
+                break;
+            }
+        }
+
+        let odd_multiples_of_A = OddMultiples::create(A);
+
+        let mut r = ProjectivePoint::identity();
+        loop {
+            let mut t = r.double();
+
+            if a_naf[i] > 0 {
+                t = &t.to_extended() + &odd_multiples_of_A[( a_naf[i]/2) as usize];
+            } else if a_naf[i] < 0 {
+                t = &t.to_extended() - &odd_multiples_of_A[(-a_naf[i]/2) as usize];
+            }
+
+            if b_naf[i] > 0 {
+                t = &t.to_extended() + &constants::bi[( b_naf[i]/2) as usize];
+            } else if b_naf[i] < 0 {
+                t = &t.to_extended() - &constants::bi[(-b_naf[i]/2) as usize];
+            }
+
+            r = t.to_projective();
+
+            if i == 0 {
+                break;
+            }
+            i -= 1;
+        }
+
+        r
+    }
+
+}
+
 // ------------------------------------------------------------------------
 // Tests
 // ------------------------------------------------------------------------
 
 #[cfg(test)]
 mod test {
+    #[cfg(feature = "yolocrypto")]
+    use decaf::DecafPoint;
     use field::FieldElement;
     use scalar::Scalar;
     use subtle::CTAssignable;
@@ -1236,11 +1345,18 @@ mod test {
     /// Test that computing 1*basepoint gives the correct basepoint.
     #[test]
     fn basepoint_mult_one_vs_basepoint() {
-        let bp = ExtendedPoint::basepoint_mult(&Scalar::one());
+        let bp = &constants::ED25519_BASEPOINT_TABLE * &Scalar::one();
         let compressed = bp.compress_edwards();
         assert_eq!(compressed, constants::BASE_CMPRSSD);
     }
 
+    /// Test that `EdwardsBasepointTable::basepoint()` gives the correct basepoint.
+    #[test]
+    fn basepoint_table_basepoint_function_correct() {
+        let bp = constants::ED25519_BASEPOINT_TABLE.basepoint();
+        assert_eq!(bp.compress_edwards(), constants::BASE_CMPRSSD);
+    }
+
     /// Test `impl Add<ExtendedPoint> for ExtendedPoint`
     /// using basepoint + basepoint versus the 2*basepoint constant.
     #[test]
@@ -1288,7 +1404,7 @@ mod test {
     #[test]
     fn to_affine_niels_clears_denominators() {
         // construct a point as aB so it has denominators (ie. Z != 1)
-        let aB = ExtendedPoint::basepoint_mult(&A_SCALAR);
+        let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
         let aB_affine_niels = aB.to_affine_niels();
         let also_aB = (&ExtendedPoint::identity() + &aB_affine_niels).to_extended();
         assert_eq!(     aB.compress_edwards(),
@@ -1298,14 +1414,15 @@ mod test {
     /// Test basepoint_mult versus a known scalar multiple from ed25519.py
     #[test]
     fn basepoint_mult_vs_ed25519py() {
-        let aB = ExtendedPoint::basepoint_mult(&A_SCALAR);
+        let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
         assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT);
     }
 
     /// Test that multiplication by the basepoint order kills the basepoint
     #[test]
     fn basepoint_mult_by_basepoint_order() {
-        let should_be_id = ExtendedPoint::basepoint_mult(&constants::l);
+        let B = &constants::ED25519_BASEPOINT_TABLE;
+        let should_be_id = B * &constants::l;
         assert!(should_be_id.is_identity());
     }
 
@@ -1314,27 +1431,18 @@ mod test {
     #[cfg(feature="basepoint_table_creation")]
     fn test_precomputed_basepoint_mult() {
         let table = EdwardsBasepointTable::create(&constants::ED25519_BASEPOINT);
-        let aB_1 = ExtendedPoint::basepoint_mult(&A_SCALAR);
-        let aB_2 = table.basepoint_mult(&A_SCALAR);
-        assert_eq!(aB_1.compress_edwards(),
-                   aB_2.compress_edwards());
+        let aB_1 = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
+        let aB_2 = &table * &A_SCALAR;
+        assert_eq!(aB_1.compress_edwards(), aB_2.compress_edwards());
     }
 
     /// Test scalar_mult versus a known scalar multiple from ed25519.py
     #[test]
     fn scalar_mult_vs_ed25519py() {
-        let aB = constants::ED25519_BASEPOINT.scalar_mult(&A_SCALAR);
+        let aB = &constants::ED25519_BASEPOINT * &A_SCALAR;
         assert_eq!(aB.compress_edwards(), A_TIMES_BASEPOINT);
     }
 
-    /// Test double_scalar_mult_vartime vs ed25519.py
-    #[test]
-    fn double_scalar_mult_vartime_vs_ed25519py() {
-        let A = A_TIMES_BASEPOINT.decompress().unwrap();
-        let result = double_scalar_mult_vartime(&A_SCALAR, &A, &B_SCALAR);
-        assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
-    }
-
     /// Test basepoint.double() versus the 2*basepoint constant.
     #[test]
     fn basepoint_double_vs_basepoint2() {
@@ -1346,7 +1454,7 @@ mod test {
     #[test]
     fn basepoint_mult_two_vs_basepoint2() {
         let mut two_bytes = [0u8; 32]; two_bytes[0] = 2;
-        let bp2 = ExtendedPoint::basepoint_mult(&Scalar(two_bytes));
+        let bp2 = &constants::ED25519_BASEPOINT_TABLE * &Scalar(two_bytes);
         assert_eq!(bp2.compress_edwards(), BASE2_CMPRSSD);
     }
 
@@ -1412,11 +1520,57 @@ mod test {
     /// the type system and prove correctness).
     #[test]
     fn monte_carlo_overflow_underflow_debug_assert_test() {
-        let mut P = ExtendedPoint::basepoint();
+        let mut P = constants::ED25519_BASEPOINT;
         // N.B. each scalar_mult does 1407 field mults, 1024 field squarings,
         // so this does ~ 1M of each operation.
         for _ in 0..1_000 {
-            P = P.scalar_mult(&A_SCALAR);
+            P *= &A_SCALAR;
+        }
+    }
+
+    #[test]
+    fn scalarmult_extended_point_works_both_ways() {
+        let G: ExtendedPoint = constants::ED25519_BASEPOINT;
+        let s: Scalar = A_SCALAR;
+
+        let P1 = &G * &s;
+        let P2 = &s * &G;
+
+        assert!(P1.compress_edwards().to_bytes() == P2.compress_edwards().to_bytes());
+    }
+
+    #[test]
+    #[cfg(feature = "yolocrypto")]
+    fn scalarmult_decafpoint_works_both_ways() {
+        let P: DecafPoint = DecafPoint(constants::ED25519_BASEPOINT);
+        let s: Scalar = A_SCALAR;
+
+        let P1 = &P * &s;
+        let P2 = &s * &P;
+
+        assert!(P1.compress().as_bytes() == P2.compress().as_bytes());
+    }
+
+    mod vartime {
+        use super::super::*;
+        use super::{A_SCALAR, B_SCALAR, A_TIMES_BASEPOINT, DOUBLE_SCALAR_MULT_RESULT};
+        
+        /// Test double_scalar_mult_vartime vs ed25519.py
+        #[test]
+        fn double_scalar_mult_basepoint_vs_ed25519py() {
+            let A = A_TIMES_BASEPOINT.decompress().unwrap();
+            let result = vartime::double_scalar_mult_basepoint(&A_SCALAR, &A, &B_SCALAR);
+            assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
+        }
+
+        #[test]
+        fn k_fold_scalar_mult_vs_ed25519py() {
+            let A = A_TIMES_BASEPOINT.decompress().unwrap();
+            let result = vartime::k_fold_scalar_mult(
+                &[A_SCALAR, B_SCALAR],
+                &[A, constants::ED25519_BASEPOINT]
+            );
+            assert_eq!(result.compress_edwards(), DOUBLE_SCALAR_MULT_RESULT);
         }
     }
 }
@@ -1427,6 +1581,7 @@ mod test {
 
 #[cfg(all(test, feature = "bench"))]
 mod bench {
+    use rand::OsRng;
     use test::Bencher;
     use constants;
     use super::*;
@@ -1434,13 +1589,14 @@ mod bench {
 
     #[bench]
     fn basepoint_mult(b: &mut Bencher) {
-        b.iter(|| ExtendedPoint::basepoint_mult(&A_SCALAR));
+        let B = &constants::ED25519_BASEPOINT_TABLE;
+        b.iter(|| B * &A_SCALAR);
     }
 
     #[bench]
     fn scalar_mult(b: &mut Bencher) {
-        let bp = constants::ED25519_BASEPOINT;
-        b.iter(|| bp.scalar_mult(&A_SCALAR));
+        let B = &constants::ED25519_BASEPOINT;
+        b.iter(|| B * &A_SCALAR);
     }
 
     #[bench]
@@ -1448,12 +1604,6 @@ mod bench {
         b.iter(|| select_precomputed_point(0, &constants::ED25519_BASEPOINT_TABLE.0[0]));
     }
 
-    #[bench]
-    fn bench_double_scalar_mult_vartime(b: &mut Bencher) {
-        let A = A_TIMES_BASEPOINT.decompress().unwrap();
-        b.iter(|| double_scalar_mult_vartime(&A_SCALAR, &A, &B_SCALAR));
-    }
-
     #[bench]
     fn add_extended_and_projective_niels_output_completed(b: &mut Bencher) {
         let p1 = constants::ED25519_BASEPOINT;
@@ -1510,7 +1660,37 @@ mod bench {
     #[cfg(feature="basepoint_table_creation")]
     #[bench]
     fn create_basepoint_table(b: &mut Bencher) {
-        let aB = ExtendedPoint::basepoint_mult(&A_SCALAR);
+        let aB = &constants::ED25519_BASEPOINT_TABLE * &A_SCALAR;
         b.iter(|| EdwardsBasepointTable::create(&aB));
     }
+
+    mod vartime {
+        use super::super::*;
+        use super::super::test::{A_SCALAR, B_SCALAR, A_TIMES_BASEPOINT};
+        use super::{Bencher, OsRng};
+
+        #[bench]
+        fn bench_double_scalar_mult_basepoint(b: &mut Bencher) {
+            let A = A_TIMES_BASEPOINT.decompress().unwrap();
+            b.iter(|| vartime::double_scalar_mult_basepoint(&A_SCALAR, &A, &B_SCALAR));
+        }
+
+        #[bench]
+        fn ten_fold_scalar_mult(b: &mut Bencher) {
+            let mut csprng: OsRng = OsRng::new().unwrap();
+            // Create 10 random scalars
+            let scalars: Vec<_> = (0..10).map(|_| Scalar::random(&mut csprng)).collect();
+            // Create 10 points (by doing scalar mults)
+            let B = &constants::ED25519_BASEPOINT_TABLE;
+            let points: Vec<_> = scalars.iter().map(|s| B * &s).collect();
+
+            // XXX Currently Rust's benchmarking implementation doesn't
+            // allow you to specify a sequence of random inputs, but only
+            // many trials of the same input.
+            //
+            // Since this is a variable-time function, this means the
+            // benchmark is only useful as a ballpark measurement.
+            b.iter(|| vartime::k_fold_scalar_mult(&scalars, &points));
+        }
+    }
 }
diff --git a/src/decaf.rs b/src/decaf.rs
index 19d4977..9d36c10 100644
--- a/src/decaf.rs
+++ b/src/decaf.rs
@@ -30,16 +30,11 @@ use subtle::CTAssignable;
 use subtle::CTNegatable;
 
 use core::ops::{Add, Sub, Neg};
+use core::ops::{Mul, MulAssign};
 
-#[cfg(all(not(feature = "std"), feature = "basepoint_table_creation"))]
-use collections::boxed::Box;
-#[cfg(all(feature = "std", feature = "basepoint_table_creation"))]
-use std::boxed::Box;
-
+use curve;
 use curve::ExtendedPoint;
 use curve::EdwardsBasepointTable;
-use curve::BasepointMult;
-use curve::ScalarMult;
 use curve::Identity;
 use scalar::Scalar;
 
@@ -249,40 +244,82 @@ impl<'a> Neg for &'a DecafPoint {
     }
 }
 
-impl ScalarMult<Scalar> for DecafPoint {
-    fn scalar_mult(&self, scalar: &Scalar) -> DecafPoint {
-        DecafPoint(self.0.scalar_mult(scalar))
+impl<'b> MulAssign<&'b Scalar> for DecafPoint {
+    fn mul_assign(&mut self, scalar: &'b Scalar) {
+        let result = (self as &DecafPoint) * scalar;
+        *self = result;
     }
 }
 
-impl BasepointMult<Scalar> for DecafPoint {
-    // XXX is this actually in the image of the isogeny,
-    // or do we need a different basepoint?
-    fn basepoint() -> DecafPoint {
-        DecafPoint(ExtendedPoint::basepoint())
-    }
-
-    fn basepoint_mult(scalar: &Scalar) -> DecafPoint {
-        DecafPoint(ExtendedPoint::basepoint_mult(scalar))
+impl<'a, 'b> Mul<&'b Scalar> for &'a DecafPoint {
+    type Output = DecafPoint;
+    /// Scalar multiplication: compute `scalar * self`.
+    fn mul(self, scalar: &'b Scalar) -> DecafPoint {
+        DecafPoint(&self.0 * scalar)
     }
 }
 
-
 /// Precomputation
 #[derive(Clone)]
-pub struct DecafBasepointTable(EdwardsBasepointTable);
+pub struct DecafBasepointTable(pub EdwardsBasepointTable);
+
+impl<'a, 'b> Mul<&'b Scalar> for &'a DecafBasepointTable {
+    type Output = DecafPoint;
+
+    fn mul(self, scalar: &'b Scalar) -> DecafPoint {
+        DecafPoint(&self.0 * scalar)
+    }
+}
+
+impl<'a, 'b> Mul<&'a DecafBasepointTable> for &'b Scalar {
+    type Output = DecafPoint;
+
+    fn mul(self, basepoint_table: &'a DecafBasepointTable) -> DecafPoint {
+        DecafPoint(self * &basepoint_table.0)
+    }
+}
 
 impl DecafBasepointTable {
     /// Create a precomputed table of multiples of the given `basepoint`.
-    #[cfg(feature = "basepoint_table_creation")]
-    pub fn create(basepoint: &DecafPoint) -> Box<DecafBasepointTable> {
-        let edwards_table = EdwardsBasepointTable::create(&basepoint.0);
-        box DecafBasepointTable(*edwards_table)
+    pub fn create(basepoint: &DecafPoint) -> DecafBasepointTable {
+        DecafBasepointTable(EdwardsBasepointTable::create(&basepoint.0))
     }
 
-    /// Use the precomputed table to quickly compute `scalar * basepoint`
-    pub fn basepoint_mult(&self, scalar: &Scalar) -> DecafPoint {
-        DecafPoint(self.0.basepoint_mult(scalar))
+    /// Get the basepoint for this table as a `DecafPoint`.
+    pub fn basepoint(&self) -> DecafPoint {
+        DecafPoint(self.0.basepoint())
+    }
+}
+
+// ------------------------------------------------------------------------
+// Constant-time conditional assignment
+// ------------------------------------------------------------------------
+
+impl CTAssignable for DecafPoint {
+    /// Conditionally assign `other` to `self`, if `choice == 1u8`.
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// # use curve25519_dalek::curve::Identity;
+    /// # use curve25519_dalek::decaf::DecafPoint;
+    /// # use curve25519_dalek::subtle::CTAssignable;
+    /// # use curve25519_dalek::constants;
+    /// let A = DecafPoint::identity();
+    /// let B = constants::DECAF_ED25519_BASEPOINT;
+    ///
+    /// let mut P = A;
+    ///
+    /// P.conditional_assign(&B, 0u8);
+    /// assert!(P == A);
+    /// P.conditional_assign(&B, 1u8);
+    /// assert!(P == B);
+    /// ```
+    fn conditional_assign(&mut self, other: &DecafPoint, choice: u8) {
+        self.0.X.conditional_assign(&other.0.X, choice);
+        self.0.Y.conditional_assign(&other.0.Y, choice);
+        self.0.Z.conditional_assign(&other.0.Z, choice);
+        self.0.T.conditional_assign(&other.0.T, choice);
     }
 }
 
@@ -304,6 +341,31 @@ impl Debug for DecafPoint {
     }
 }
 
+// ------------------------------------------------------------------------
+// Variable-time functions
+// ------------------------------------------------------------------------
+
+pub mod vartime {
+    //! Variable-time operations on decaf points, useful for non-secret data.
+    use super::*;
+
+    /// Given a vector of public scalars and a vector of (possibly secret)
+    /// points, compute
+    ///
+    ///    c_1 P_1 + ... + c_n P_n.
+    ///
+    /// # Input
+    ///
+    /// A vector of `Scalar`s and a vector of `ExtendedPoints`.  It is an
+    /// error to call this function with two vectors of different lengths.
+    pub fn k_fold_scalar_mult<'a,'b,I,J>(scalars: I, points: J) -> DecafPoint
+        where I: IntoIterator<Item=&'a Scalar>, J: IntoIterator<Item=&'b DecafPoint>
+    {
+        let extended_points = points.into_iter().map(|P| &P.0);
+        DecafPoint(curve::vartime::k_fold_scalar_mult(scalars, extended_points))
+    }
+}
+
 // ------------------------------------------------------------------------
 // Tests
 // ------------------------------------------------------------------------
@@ -316,7 +378,6 @@ mod test {
     use constants;
     use curve::CompressedEdwardsY;
     use curve::ExtendedPoint;
-    use curve::BasepointMult;
     use curve::Identity;
     use super::*;
 
@@ -342,17 +403,17 @@ mod test {
 
     #[test]
     fn decaf_basepoint_roundtrip() {
-        let bp_compressed_decaf = DecafPoint::basepoint().compress();
+        let bp_compressed_decaf = constants::DECAF_ED25519_BASEPOINT.compress();
         let bp_recaf = bp_compressed_decaf.decompress().unwrap().0;
         // Check that bp_recaf differs from bp by a point of order 4
-        let diff = &ExtendedPoint::basepoint() - &bp_recaf;
-        let diff4 = diff.mult_by_pow_2(4);
+        let diff = &constants::ED25519_BASEPOINT - &bp_recaf;
+        let diff4 = diff.mult_by_pow_2(4); // XXX this is wrong
         assert_eq!(diff4.compress_edwards(), CompressedEdwardsY::identity());
     }
 
     #[test]
     fn decaf_four_torsion_basepoint() {
-        let bp = DecafPoint::basepoint();
+        let bp = constants::DECAF_ED25519_BASEPOINT;
         let bp_coset = bp.coset4();
         for i in 0..4 {
             assert_eq!(bp, DecafPoint(bp_coset[i]));
@@ -362,8 +423,8 @@ mod test {
     #[test]
     fn decaf_four_torsion_random() {
         let mut rng = OsRng::new().unwrap();
-        let s = Scalar::random(&mut rng);
-        let P = DecafPoint::basepoint_mult(&s);
+        let B = &constants::DECAF_ED25519_BASEPOINT_TABLE;
+        let P = B * &Scalar::random(&mut rng);
         let P_coset = P.coset4();
         for i in 0..4 {
             assert_eq!(P, DecafPoint(P_coset[i]));
@@ -373,27 +434,14 @@ mod test {
     #[test]
     fn decaf_random_roundtrip() {
         let mut rng = OsRng::new().unwrap();
+        let B = &constants::DECAF_ED25519_BASEPOINT_TABLE;
         for _ in 0..100 {
-            let s = Scalar::random(&mut rng);
-            let P = DecafPoint::basepoint_mult(&s);
+            let P = B * &Scalar::random(&mut rng);
             let compressed_P = P.compress();
             let Q = compressed_P.decompress().unwrap();
             assert_eq!(P, Q);
         }
     }
-
-    /// Test basepoint_mult versus a newly-generated DecafBasepointTable
-    #[test]
-    #[cfg(feature = "basepoint_table_creation")]
-    fn basepoint_mult_vs_decafbasepointtable() {
-        let table = DecafBasepointTable::create(&DecafPoint::basepoint());
-        let mut rng = OsRng::new().unwrap();
-        let s = Scalar::random(&mut rng);
-        let basepoint_mult_s = DecafPoint::basepoint_mult(&s);
-        let table_basepoint_mult_s = table.basepoint_mult(&s);
-
-        assert_eq!(basepoint_mult_s, table_basepoint_mult_s);
-    }
 }
 
 #[cfg(all(test, feature = "bench"))]
@@ -406,8 +454,8 @@ mod bench {
     #[bench]
     fn decompression(b: &mut Bencher) {
         let mut rng = OsRng::new().unwrap();
-        let s = Scalar::random(&mut rng);
-        let P = DecafPoint::basepoint_mult(&s);
+        let B = &constants::DECAF_ED25519_BASEPOINT_TABLE;
+        let P = B * &Scalar::random(&mut rng);
         let P_compressed = P.compress();
         b.iter(|| P_compressed.decompress().unwrap());
     }
@@ -415,8 +463,8 @@ mod bench {
     #[bench]
     fn compression(b: &mut Bencher) {
         let mut rng = OsRng::new().unwrap();
-        let s = Scalar::random(&mut rng);
-        let P = DecafPoint::basepoint_mult(&s);
+        let B = &constants::DECAF_ED25519_BASEPOINT_TABLE;
+        let P = B * &Scalar::random(&mut rng);
         b.iter(|| P.compress());
     }
 }
diff --git a/src/field.rs b/src/field.rs
index 78dcc9b..2d43545 100644
--- a/src/field.rs
+++ b/src/field.rs
@@ -15,7 +15,6 @@
 //! Based on Adam Langley's curve25519-donna and (Golang) ed25519
 //! implementations.
 
-use core::clone::Clone;
 use core::fmt::Debug;
 use core::ops::{Add, AddAssign};
 use core::ops::{Sub, SubAssign};
@@ -41,25 +40,27 @@ use constants;
 #[cfg(feature="radix_51")]
 pub type Limb = u64;
 
-/// FieldElement represents an element of the field GF(2^255 - 19).  An element
-/// t, entries t[0]...t[9], represents the integer t[0]+2^26 t[1]+2^51 t[2]+2^77
-/// t[3]+2^102 t[4]+...+2^230 t[9].  Bounds on each t[i] vary depending on
-/// context.
+/// A `FieldElement` represents an element of the field GF(2^255 - 19).
+///
+/// With the `radix_51` feature, a `FieldElement` is represented in
+/// radix 2^51 as five `u64`s; the coefficients are allowed to grow up
+/// to 2^54 between reductions mod `p`.
 #[cfg(feature="radix_51")]
 #[derive(Copy, Clone)]
 pub struct FieldElement(pub [u64; 5]);
 
-/// FieldElements are represented as an array of ten "Limbs", which are radix
-/// 25.5, that is, each Limb of a FieldElement alternates between being
-/// represented as a factor of 2^25 or 2^26 more than the last corresponding
-/// integer.
+/// Without the `radix51` feature enabled, `FieldElements` are represented
+/// in radix 2^25.5 as ten `i32`s.
 #[cfg(not(feature="radix_51"))]
 pub type Limb = i32;
 
-/// FieldElement represents an element of the field GF(2^255 - 19).  An element
-/// t, entries t[0]...t[9], represents the integer t[0]+2^26 t[1]+2^51 t[2]+2^77
-/// t[3]+2^102 t[4]+...+2^230 t[9].  Bounds on each t[i] vary depending on
-/// context.
+/// A `FieldElement` represents an element of the field GF(2^255 - 19).
+///
+/// With the `radix_51` feature, a `FieldElement` is represented in
+/// radix 2^25.5 as ten `i32`s, so that an element t, entries
+/// t[0],...,t[9], represents the integer t[0]+2^26 t[1]+2^51
+/// t[2]+2^77 t[3]+2^102 t[4]+...+2^230 t[9].  Bounds on each t[i]
+/// vary depending on context.
 #[cfg(not(feature="radix_51"))]
 #[derive(Copy, Clone)]
 pub struct FieldElement(pub [i32; 10]);
@@ -84,7 +85,7 @@ impl PartialEq for FieldElement {
         for i in 0..32 {
             are_equal &= self_bytes[i] == other_bytes[i];
         }
-        return are_equal;
+        are_equal
     }
 }
 
@@ -108,16 +109,14 @@ impl Debug for FieldElement {
 impl Index<usize> for FieldElement {
     type Output = Limb;
 
-    fn index<'a>(&'a self, _index: usize) -> &'a Limb {
-        let ret: &'a Limb = &(self.0[_index]);
-        ret
+    fn index(&self, _index: usize) -> &Limb {
+        &(self.0[_index])
     }
 }
 
 impl IndexMut<usize> for FieldElement {
-    fn index_mut<'a>(&'a mut self, _index: usize) -> &'a mut Limb {
-        let ret: &'a mut Limb = &mut(self.0[_index]);
-        ret
+    fn index_mut(&mut self, _index: usize) -> &mut Limb {
+        &mut(self.0[_index])
     }
 }
 
@@ -132,7 +131,7 @@ impl<'b> AddAssign<&'b FieldElement> for FieldElement {
 impl<'a, 'b> Add<&'b FieldElement> for &'a FieldElement {
     type Output = FieldElement;
     fn add(self, _rhs: &'b FieldElement) -> FieldElement {
-        let mut output = self.clone();
+        let mut output = *self;
         output += _rhs;
         output
     }
@@ -158,7 +157,7 @@ impl<'a, 'b> Sub<&'b FieldElement> for &'a FieldElement {
     type Output = FieldElement;
     #[cfg(not(feature="radix_51"))]
     fn sub(self, _rhs: &'b FieldElement) -> FieldElement {
-        let mut output = self.clone();
+        let mut output = *self;
         output -= _rhs;
         output
     }
@@ -319,14 +318,14 @@ impl<'a, 'b> Mul<&'b FieldElement> for &'a FieldElement {
         let h8 = f0*g8 + f1_2*g7 + f2*g6 + f3_2*g5 + f4*g4 + f5_2*g3 + f6*g2 + f7_2*g1 + f8*g0 + f9_2*g9_19;
         let h9 = f0*g9 + f1*g8 + f2*g7 + f3*g6 + f4*g5 + f5*g4 + f6*g3 + f7*g2 + f8*g1 + f9*g0;
 
-        FieldElement::reduce(&[h0, h1, h2, h3, h4, h5, h6, h7, h8, h9])
+        FieldElement::reduce([h0, h1, h2, h3, h4, h5, h6, h7, h8, h9])
     }
 }
 
 impl<'a> Neg for &'a FieldElement {
     type Output = FieldElement;
     fn neg(self) -> FieldElement {
-        let mut output = self.clone();
+        let mut output = *self;
         output.negate();
         output
     }
@@ -454,10 +453,10 @@ impl FieldElement {
 
         FieldElement(limbs)
     }
+
     #[cfg(not(feature="radix_51"))]
-    fn reduce(input: &[i64;10]) -> FieldElement { //FeCombine
+    fn reduce(mut h: [i64; 10]) -> FieldElement { //FeCombine
         let mut c = [0i64;10];
-        let mut h = input.clone();
 
         /*
           |h[0]| <= (1.1*1.1*2^52*(1+19+19+19+19)+1.1*1.1*2^50*(38+38+38+38+38))
@@ -533,7 +532,7 @@ impl FieldElement {
         /* |h[0]| <= 2^25; from now on fits into int32 unchanged */
         /* |h[1]| <= 1.01*2^24 */
 
-        let mut output = FieldElement([0i32;10]);
+        let mut output = FieldElement([0i32; 10]);
         output[0] = h[0] as i32;
         output[1] = h[1] as i32;
         output[2] = h[2] as i32;
@@ -581,11 +580,11 @@ impl FieldElement {
         h[8] =  load3(&data[26..]) << 4;
         h[9] = (load3(&data[29..]) & 8388607) << 2;
 
-        FieldElement::reduce(&h)
+        FieldElement::reduce(h)
     }
     /// Parse a `FieldElement` from 32 bytes.
     #[cfg(feature="radix_51")]
-    pub fn from_bytes(bytes: &[u8;32]) -> FieldElement {
+    pub fn from_bytes(bytes: &[u8; 32]) -> FieldElement {
         let low_51_bit_mask = (1u64 << 51) - 1;
         FieldElement(
         // load bits [  0, 64), no shift
@@ -621,7 +620,7 @@ impl FieldElement {
     /// assert!(data == bytes);
     /// ```
     #[cfg(not(feature="radix_51"))]
-    pub fn to_bytes(&self) -> [u8;32] { //FeToBytes
+    pub fn to_bytes(&self) -> [u8; 32] { //FeToBytes
         // Comment preserved from ed25519.go (presumably originally from ref10):
         //
         // # Preconditions
@@ -656,7 +655,7 @@ impl FieldElement {
         // so floor(2^-255 * (h + 19 * 2^-25 * h9 + 2^-1)) = q.
         //
         let mut carry = [0i32; 10];
-        let mut h = self.clone();
+        let mut h: [i32; 10] = self.0;
 
         let mut q:i32 = (19*h[9] + (1 << 24)) >> 25;
         q = (h[0] + q) >> 26;
@@ -752,7 +751,7 @@ impl FieldElement {
     }
     /// Serialize this `FieldElement` to bytes.
     #[cfg(feature="radix_51")]
-    pub fn to_bytes(&self) -> [u8;32] {
+    pub fn to_bytes(&self) -> [u8; 32] {
         // This reduces to the range [0,2^255), but we need [0,2^255-19)
         let mut limbs = FieldElement::reduce(self.0).0;
         // Let h = limbs[0] + limbs[1]*2^51 + ... + limbs[4]*2^204.
@@ -830,32 +829,6 @@ impl FieldElement {
         return s
     }
 
-    /// XXX clarify documentation
-    /// Determine if this field element, represented as a byte array,
-    /// is less than or equal to another field element represented as
-    /// a byte array.
-    ///
-    /// # Returns
-    ///
-    /// Returns `1u8` if `self.to_bytes() <= other.to_bytes()`, and `0u8` otherwise.
-    pub fn bytes_equal_less_than(&self, other: &[u8; 32]) -> u8 { // feBytesLess
-        // XXX cleanup
-	    let mut equal_so_far: i32 = -1i32;
-	    let mut greater:      i32 =  0i32;
-
-        let this: [u8; 32] = self.to_bytes();
-
-        for i in 32 .. 0 {
-            let x: i32 =  this[i-1] as i32;
-            let y: i32 = other[i-1] as i32;
-
-            greater      = (!equal_so_far & greater) | (equal_so_far & ((x - y) >> 31));
-            equal_so_far = equal_so_far & (((x ^ y) - 1) >> 31);
-        }
-
-        (!equal_so_far & 1 & greater) as u8
-    }
-
     /// Determine if this `FieldElement` is negative, in the sense
     /// used in the ed25519 paper: `x` is negative if the low bit is
     /// set.
@@ -907,7 +880,7 @@ impl FieldElement {
     ///
     /// If zero, return `1u8`.  Otherwise, return `0u8`.
     pub fn is_zero(&self) -> u8 {
-        return 1u8 & (!self.is_nonzero());
+        1u8 & (!self.is_nonzero())
     }
 
     /// Determine if this `FieldElement` is non-zero.
@@ -921,11 +894,11 @@ impl FieldElement {
         for b in &bytes {
             x |= *b;
         }
-        return byte_is_nonzero(x);
+        byte_is_nonzero(x)
     }
 
     #[cfg(not(feature="radix_51"))]
-    fn square_inner(&self) -> [i64;10] {
+    fn square_inner(&self) -> [i64; 10] {
         let f0     = self[0]       as i64;
         let f1     = self[1]       as i64;
         let f2     = self[2]       as i64;
@@ -964,6 +937,7 @@ impl FieldElement {
 
         h
     }
+
     #[cfg(feature="radix_51")]
     #[inline(always)]
     fn square_inner(&self) -> [u64; 5] {
@@ -1028,12 +1002,12 @@ impl FieldElement {
     /// * |h[i]| bounded by 1.1*2^25, 1.1*2^24, 1.1*2^25, 1.1*2^24, etc.
     #[cfg(not(feature="radix_51"))]
     pub fn square(&self) -> FieldElement {
-        FieldElement::reduce(&self.square_inner())
+        FieldElement::reduce(self.square_inner())
     }
     /// Compute `self^2`.
     #[cfg(feature="radix_51")]
     pub fn square(&self) -> FieldElement {
-        FieldElement::reduce( self.square_inner())
+        FieldElement::reduce(self.square_inner())
     }
 
     /// Square this field element and multiply the result by 2.
@@ -1058,7 +1032,7 @@ impl FieldElement {
         for i in 0..self.0.len() {
             coeffs[i] += coeffs[i];
         }
-        FieldElement::reduce(&coeffs)
+        FieldElement::reduce(coeffs)
     }
     /// Compute `2 * self^2`.
     #[cfg(feature="radix_51")]
@@ -1170,8 +1144,7 @@ impl FieldElement {
     /// - `(0u8, zero)`      if `v` is zero;
     /// - `(0u8, garbage)`   if `u/v` is nonsquare.
     ///
-    pub fn sqrt_ratio(u: &FieldElement, v: &FieldElement)
-                      -> (u8, FieldElement) {
+    pub fn sqrt_ratio(u: &FieldElement, v: &FieldElement) -> (u8, FieldElement) {
         // Using the same trick as in ed25519 decoding, we merge the
         // inversion, the square root, and the square test as follows.
         //
@@ -1259,28 +1232,28 @@ mod test {
     /// Random element a of GF(2^255-19), from Sage
     /// a = 1070314506888354081329385823235218444233221\
     ///     2228051251926706380353716438957572
-    pub static A_BYTES: [u8;32] =
+    pub static A_BYTES: [u8; 32] =
         [ 0x04, 0xfe, 0xdf, 0x98, 0xa7, 0xfa, 0x0a, 0x68,
           0x84, 0x92, 0xbd, 0x59, 0x08, 0x07, 0xa7, 0x03,
           0x9e, 0xd1, 0xf6, 0xf2, 0xe1, 0xd9, 0xe2, 0xa4,
           0xa4, 0x51, 0x47, 0x36, 0xf3, 0xc3, 0xa9, 0x17];
 
     /// Byte representation of a**2
-    static ASQ_BYTES: [u8;32] =
+    static ASQ_BYTES: [u8; 32] =
         [ 0x75, 0x97, 0x24, 0x9e, 0xe6, 0x06, 0xfe, 0xab,
           0x24, 0x04, 0x56, 0x68, 0x07, 0x91, 0x2d, 0x5d,
           0x0b, 0x0f, 0x3f, 0x1c, 0xb2, 0x6e, 0xf2, 0xe2,
           0x63, 0x9c, 0x12, 0xba, 0x73, 0x0b, 0xe3, 0x62];
 
     /// Byte representation of 1/a
-    static AINV_BYTES: [u8;32] =
+    static AINV_BYTES: [u8; 32] =
         [0x96, 0x1b, 0xcd, 0x8d, 0x4d, 0x5e, 0xa2, 0x3a,
          0xe9, 0x36, 0x37, 0x93, 0xdb, 0x7b, 0x4d, 0x70,
          0xb8, 0x0d, 0xc0, 0x55, 0xd0, 0x4c, 0x1d, 0x7b,
          0x90, 0x71, 0xd8, 0xe9, 0xb6, 0x18, 0xe6, 0x30];
 
     /// Byte representation of a^((p-5)/8)
-    static AP58_BYTES: [u8;32] =
+    static AP58_BYTES: [u8; 32] =
         [0x6a, 0x4f, 0x24, 0x89, 0x1f, 0x57, 0x60, 0x36,
          0xd0, 0xbe, 0x12, 0x3c, 0x8f, 0xf5, 0xb1, 0x59,
          0xe0, 0xf0, 0xb8, 0x1b, 0x20, 0xd2, 0xb5, 0x1f,
@@ -1352,7 +1325,7 @@ mod test {
 
     #[test]
     fn from_bytes_highbit_is_ignored() {
-        let mut cleared_bytes = B_BYTES.clone();
+        let mut cleared_bytes = B_BYTES;
         cleared_bytes[31] &= 127u8;
         let with_highbit_set    = FieldElement::from_bytes(&B_BYTES);
         let without_highbit_set = FieldElement::from_bytes(&cleared_bytes);
diff --git a/src/lib.rs b/src/lib.rs
index 31e7177..40efe8b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -11,10 +11,10 @@
 
 #![cfg_attr(not(feature = "std"), no_std)]
 #![cfg_attr(not(feature = "std"), feature(collections))]
-#![cfg_attr(feature = "nightly", feature(box_syntax))]
 #![cfg_attr(feature = "nightly", feature(i128_type))]
-#![allow(unused_features)]
 #![cfg_attr(feature = "bench", feature(test))]
+
+#![allow(unused_features)]
 #![deny(missing_docs)] // refuse to compile if documentation is missing
 
 //! # curve25519-dalek
diff --git a/src/scalar.rs b/src/scalar.rs
index f7808e3..34a45ba 100644
--- a/src/scalar.rs
+++ b/src/scalar.rs
@@ -29,9 +29,13 @@
 //! between two scalars, the `UnpackedScalar` struct is stored as
 //! limbs.
 
-use core::cmp::{Eq, PartialEq};
-use core::ops::{Neg, Index, IndexMut};
 use core::fmt::Debug;
+use core::ops::Neg;
+use core::ops::{Add, AddAssign};
+use core::ops::{Sub, SubAssign};
+use core::ops::{Mul, MulAssign};
+use core::ops::{Index, IndexMut};
+use core::cmp::{Eq, PartialEq};
 
 #[cfg(feature = "std")]
 use rand::Rng;
@@ -72,13 +76,7 @@ impl PartialEq for Scalar {
     ///
     /// True if they are equal, and false otherwise.
     fn eq(&self, other: &Self) -> bool {
-        let equal: u8 = arrays_equal_ct(&self.0, &other.0);
-
-        if equal == 1u8 {
-            return true;
-        } else {
-            return false;
-        }
+        arrays_equal_ct(&self.0, &other.0) == 1u8
     }
 }
 
@@ -96,28 +94,66 @@ impl CTEq for Scalar {
 impl Index<usize> for Scalar {
     type Output = u8;
 
-    fn index<'a>(&'a self, _index: usize) -> &'a u8 {
-        let ret: &'a u8 = &(self.0[_index]);
-        ret
+    fn index(&self, _index: usize) -> &u8 {
+        &(self.0[_index])
     }
 }
 
 impl IndexMut<usize> for Scalar {
-    fn index_mut<'a>(&'a mut self, _index: usize) -> &'a mut u8 {
-        let ret: &'a mut u8 = &mut(self.0[_index]);
-        ret
+    fn index_mut(&mut self, _index: usize) -> &mut u8 {
+        &mut(self.0[_index])
     }
 }
 
-impl Neg for Scalar {
-    type Output = Scalar;
-
-    /// Negate this scalar by computing (l - 1) * self - 0 (mod l).
-    fn neg(self) -> Scalar {
-        Scalar::multiply_add(&constants::lminus1, &self, &Scalar::zero())
+impl<'b> MulAssign<&'b Scalar> for Scalar {
+    fn mul_assign(&mut self, _rhs: &'b Scalar) {
+        let result = (self as &Scalar) * _rhs;
+        self.0 = result.0;
     }
 }
 
+impl<'a, 'b> Mul<&'b Scalar> for &'a Scalar {
+    type Output = Scalar;
+    fn mul(self, _rhs: &'b Scalar) -> Scalar {
+        Scalar::multiply_add(self, _rhs, &Scalar::zero())
+    }
+}
+
+impl<'b> AddAssign<&'b Scalar> for Scalar {
+    fn add_assign(&mut self, _rhs: &'b Scalar) {
+        *self = Scalar::multiply_add(&Scalar::one(), self, _rhs);
+    }
+}
+
+impl<'a, 'b> Add<&'b Scalar> for &'a Scalar {
+    type Output = Scalar;
+    fn add(self, _rhs: &'b Scalar) -> Scalar {
+        Scalar::multiply_add(&Scalar::one(), self, _rhs)
+    }
+}
+
+impl<'b> SubAssign<&'b Scalar> for Scalar {
+    fn sub_assign(&mut self, _rhs: &'b Scalar) {
+        // (l-1)*_rhs + self = self - _rhs
+        *self = Scalar::multiply_add(&constants::l_minus_1, _rhs, self);
+    }
+}
+
+impl<'a, 'b> Sub<&'b Scalar> for &'a Scalar {
+    type Output = Scalar;
+    fn sub(self, _rhs: &'b Scalar) -> Scalar {
+        // (l-1)*_rhs + self = self - _rhs
+        Scalar::multiply_add(&constants::l_minus_1, _rhs, self)
+    }
+}
+
+impl<'a> Neg for &'a Scalar {
+    type Output = Scalar;
+    fn neg(self) -> Scalar {
+        self * &constants::l_minus_1
+   }
+}
+
 impl CTAssignable for Scalar {
     /// Conditionally assign another Scalar to this one.
     ///
@@ -150,21 +186,19 @@ impl CTAssignable for Scalar {
 }
 
 impl Scalar {
-    /// Return a `Scalar` chosen uniformly at random using a CSPRNG.
-    /// Panics if the operating system's CSPRNG is unavailable.
+    /// Return a `Scalar` chosen uniformly at random using a user-provided RNG.
     ///
     /// # Inputs
     ///
-    /// * `cspring`: any cryptographically secure PRNG which
-    ///   implements the `rand::Rng` interface.
+    /// * `rng`: any RNG which implements the `rand::Rng` interface.
     ///
     /// # Returns
     ///
     /// A random scalar within ℤ/lℤ.
     #[cfg(feature = "std")]
-    pub fn random<T: Rng>(csprng: &mut T) -> Self {
+    pub fn random<T: Rng>(rng: &mut T) -> Self {
         let mut scalar_bytes = [0u8; 64];
-        csprng.fill_bytes(&mut scalar_bytes);
+        rng.fill_bytes(&mut scalar_bytes);
         Scalar::reduce(&scalar_bytes)
     }
 
@@ -173,6 +207,8 @@ impl Scalar {
     /// Takes a type parameter `D`, which is any `Digest` producing 64
     /// bytes (512 bits) of output.
     ///
+    /// Convenience wrapper around `from_hash`.
+    ///
     /// # Example
     ///
     /// ```
@@ -190,9 +226,19 @@ impl Scalar {
     /// ```
     ///
     pub fn hash_from_bytes<D>(input: &[u8]) -> Scalar
-            where D: Digest<OutputSize=U64> + Default {
+            where D: Digest<OutputSize = U64> + Default {
         let mut hash = D::default();
         hash.input(input);
+        Scalar::from_hash(hash)
+    }
+
+    /// Construct a scalar from an existing `Digest` instance.
+    ///
+    /// Use this instead of `hash_from_bytes` if it is more convenient
+    /// to stream data into the `Digest` than to pass a single byte
+    /// slice.
+    pub fn from_hash<D>(hash: D) -> Scalar
+            where D: Digest<OutputSize=U64> + Default {
         // XXX this seems clumsy
         let mut output = [0u8;64];
         output.copy_from_slice(hash.result().as_slice());
@@ -200,7 +246,7 @@ impl Scalar {
     }
 
     /// View this `Scalar` as a sequence of bytes.
-    pub fn as_bytes<'a>(&'a self) -> &'a [u8;32] {
+    pub fn as_bytes(&self) -> &[u8; 32] {
         &self.0
     }
 
@@ -215,6 +261,31 @@ impl Scalar {
                  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ])
     }
 
+    /// Construct a scalar from the given `u64`.
+    pub fn from_u64(x: u64) -> Scalar {
+        let mut s = Scalar::zero();
+        for i in 0..8 {
+            s[i] = (x >> (i*8)) as u8;
+        }
+        s
+    }
+
+    /// Compute the multiplicative inverse of this scalar.
+    pub fn invert(&self) -> Scalar {
+        self.unpack().invert().pack()
+    }
+
+    /// Get the bits of the scalar.
+    pub fn bits(&self) -> [i8;256] {
+        let mut bits = [0i8; 256];
+        for i in 0..256 {
+            // As i runs from 0..256, the bottom 3 bits index the bit,
+            // while the upper bits index the byte.
+            bits[i] = ((self.0[i>>3] >> (i&7)) & 1u8) as i8;
+        }
+        bits
+    }
+
     /// Compute a width-5 "Non-Adjacent Form" of this scalar.
     ///
     /// A width-`w` NAF of a positive integer `k` is an expression
@@ -226,14 +297,9 @@ impl Scalar {
     /// Intuitively, this is like a binary expansion, except that we
     /// allow some coefficients to grow up to `2^(w-1)` so that the
     /// nonzero coefficients are as sparse as possible.
-    pub fn non_adjacent_form(&self) -> [i8;256] {
+    pub fn non_adjacent_form(&self) -> [i8; 256] {
         // Step 1: write out bits of the scalar
-        let mut naf = [0i8; 256];
-        for i in 0..256 {
-            // As i runs from 0..256, the bottom 3 bits index the bit,
-            // while the upper bits index the byte.
-            naf[i] = ((self.0[i>>3] >> (i&7)) & 1u8) as i8;
-        }
+        let mut naf = self.bits();
 
         // Step 2: zero coefficients by carrying them upwards or downwards 
         'bits: for i in 0..256 {
@@ -270,7 +336,7 @@ impl Scalar {
     // Unpack a scalar into 12 21-bit limbs.
     fn unpack(&self) -> UnpackedScalar {
         let mask_21bits: i64 = (1 << 21) -1;
-        let mut a = UnpackedScalar([0i64;12]);
+        let mut a = UnpackedScalar([0i64; 12]);
         a[ 0]  = mask_21bits &  load3(&self.0[ 0..])      ;
         a[ 1]  = mask_21bits & (load4(&self.0[ 2..]) >> 5);
         a[ 2]  = mask_21bits & (load3(&self.0[ 5..]) >> 2);
@@ -296,7 +362,7 @@ impl Scalar {
     ///
     /// Precondition: self[31] <= 127.  This is the case whenever
     /// `self` is reduced.
-    pub fn to_radix_16(&self) -> [i8;64] {
+    pub fn to_radix_16(&self) -> [i8; 64] {
         debug_assert!(self[31] <= 127);
         let mut output = [0i8; 64];
 
@@ -339,8 +405,8 @@ impl Scalar {
     }
 
     /// Reduce a 512-bit little endian number mod l
-    pub fn reduce(input: &[u8;64]) -> Scalar {
-        let mut s = [0i64;24];
+    pub fn reduce(input: &[u8; 64]) -> Scalar {
+        let mut s = [0i64; 24];
 
         // XXX express this as two unpack_limbs
         // some issues re: masking with the top byte of the 32byte input
@@ -387,16 +453,14 @@ pub struct UnpackedScalar(pub [i64; 12]);
 impl Index<usize> for UnpackedScalar {
     type Output = i64;
 
-    fn index<'a>(&'a self, _index: usize) -> &'a i64 {
-        let ret: &'a i64 = &(self.0[_index]);
-        ret
+    fn index(&self, _index: usize) -> &i64 {
+        &(self.0[_index])
     }
 }
 
 impl IndexMut<usize> for UnpackedScalar {
-    fn index_mut<'a>(&'a mut self, _index: usize) -> &'a mut i64 {
-        let ret: &'a mut i64 = &mut(self.0[_index]);
-        ret
+    fn index_mut(&mut self, _index: usize) -> &mut i64 {
+        &mut(self.0[_index])
     }
 }
 
@@ -440,11 +504,34 @@ impl UnpackedScalar {
         s
     }
 
+    /// Return the zero scalar.
+    pub fn zero() -> UnpackedScalar {
+        UnpackedScalar([0,0,0,0,0,0,0,0,0,0,0,0])
+    }
+
+    /// Return the one scalar.
+    pub fn one() -> UnpackedScalar {
+        UnpackedScalar([1,0,0,0,0,0,0,0,0,0,0,0])
+    }
+
+    /// Compute the multiplicative inverse of this scalar.
+    pub fn invert(&self) -> UnpackedScalar {
+        let mut y = UnpackedScalar::one();
+        // Run through bits of l-2 from highest to least
+        for bit in constants::l_minus_2.bits().iter().rev() {
+            y = UnpackedScalar::multiply_add(&y, &y, &UnpackedScalar::zero());
+            if *bit == 1 {
+                y = UnpackedScalar::multiply_add(&y, self, &UnpackedScalar::zero());
+            }
+        }
+        y
+    }
+
     /// Compute `ab+c (mod l)`.
     pub fn multiply_add(a: &UnpackedScalar,
                         b: &UnpackedScalar,
                         c: &UnpackedScalar) -> UnpackedScalar {
-        let mut result = [0i64;24];
+        let mut result = [0i64; 24];
 
         // Multiply a and b, and add c
         result[0]  =         c[0] +  a[0]*b[0];
@@ -506,10 +593,10 @@ impl UnpackedScalar {
     /// limbs.  Reduction mod l amounts to eliminating all of the
     /// high limbs while carrying as appropriate to prevent
     /// overflows in the lower limbs.
-    fn reduce_limbs(mut limbs: &mut [i64;24]) -> UnpackedScalar {
+    fn reduce_limbs(mut limbs: &mut [i64; 24]) -> UnpackedScalar {
         #[inline]
         #[allow(dead_code)]
-        fn do_reduction(limbs: &mut [i64;24], i:usize) {
+        fn do_reduction(limbs: &mut [i64; 24], i:usize) {
             limbs[i - 12] += limbs[i] * 666643;
             limbs[i - 11] += limbs[i] * 470296;
             limbs[i - 10] += limbs[i] * 654183;
@@ -531,7 +618,7 @@ impl UnpackedScalar {
         #[allow(dead_code)]
         /// Carry excess from the `i`-th limb into the `(i+1)`-th limb.
         /// Postcondition: `-2^20 <= limbs[i] < 2^20`.
-        fn do_carry_centered(limbs: &mut [i64;24], i:usize) {
+        fn do_carry_centered(limbs: &mut [i64; 24], i:usize) {
             let carry: i64 = (limbs[i] + (1<<20)) >> 21;
             limbs[i+1] += carry;
             limbs[i  ] -= carry << 21;
@@ -584,8 +671,7 @@ impl UnpackedScalar {
             do_carry_uncentered(&mut limbs, i);
         }
 
-        // XXX better way to get [i64;12] from [i64;24] ?
-        UnpackedScalar(*array_ref!(limbs,0,12))
+        UnpackedScalar(*array_ref!(limbs, 0, 12))
     }
 
 }
@@ -632,7 +718,7 @@ mod test {
         0xa7, 0x58, 0xaa, 0x1b, 0x88, 0xe0, 0x40, 0xd1,
         0x58, 0x9e, 0x7b, 0x7f, 0x23, 0x76, 0xef, 0x09]);
 
-    static A_NAF: [i8;256] =
+    static A_NAF: [i8; 256] =
         [0,13,0,0,0,0,0,0,0,7,0,0,0,0,0,0,-9,0,0,0,0,-11,0,0,0,0,3,0,0,0,0,1,
          0,0,0,0,9,0,0,0,0,-5,0,0,0,0,0,0,3,0,0,0,0,11,0,0,0,0,11,0,0,0,0,0,
          -9,0,0,0,0,0,-3,0,0,0,0,9,0,0,0,0,0,1,0,0,0,0,0,0,-1,0,0,0,0,0,9,0,
@@ -650,6 +736,20 @@ mod test {
         }
     }
 
+    #[test]
+    fn from_unsigned() {
+        let val = 0xdeadbeefdeadbeef;
+        let s = Scalar::from_u64(val);
+        assert_eq!(s[7], 0xde);
+        assert_eq!(s[6], 0xad);
+        assert_eq!(s[5], 0xbe);
+        assert_eq!(s[4], 0xef);
+        assert_eq!(s[3], 0xde);
+        assert_eq!(s[2], 0xad);
+        assert_eq!(s[1], 0xbe);
+        assert_eq!(s[0], 0xef);
+    }
+
     #[test]
     fn scalar_multiply_by_one() {
         let one = Scalar::one();
@@ -661,12 +761,25 @@ mod test {
     }
 
     #[test]
-    fn scalar_multiply_only() {
-        let zero = Scalar::zero();
-        let test_scalar = Scalar::multiply_add(&X, &Y, &zero);
-        for i in 0..32 {
-            assert!(test_scalar[i] == X_TIMES_Y[i]);
-        }
+    fn impl_add() {
+        let mut two = Scalar::zero(); two[0] = 2;
+        let two = two;
+        let one = Scalar::one();
+        let should_be_two = &one + &one;
+        assert_eq!(should_be_two, two);
+    }
+
+    #[test]
+    fn impl_sub() {
+        let should_be_one = &constants::l - &constants::l_minus_1;
+        assert_eq!(should_be_one, Scalar::one());
+    }
+
+    #[allow(non_snake_case)]
+    #[test]
+    fn impl_mul() {
+        let should_be_X_times_Y = &X * &Y;
+        assert_eq!(should_be_X_times_Y, X_TIMES_Y);
     }
 
     #[test]
@@ -679,7 +792,7 @@ mod test {
 
     #[test]
     fn scalar_reduce() {
-        let mut bignum = [0u8;64];
+        let mut bignum = [0u8; 64];
         // set bignum = x + 2^256x
         for i in 0..32 {
             bignum[   i] = X[i];
@@ -697,13 +810,22 @@ mod test {
         }
     }
 
-    // Negating a scalar twice should result in the original scalar.
+    #[allow(non_snake_case)]
     #[test]
-    fn scalar_neg() {
-        let negative_x: Scalar = -X;
-        let orig: Scalar = -negative_x;
+    fn invert() {
+        let inv_X = X.invert();
+        let should_be_one = &inv_X * &X;
+        assert_eq!(should_be_one, Scalar::one());
+    }
 
-        assert!(orig == X);
+    // Negating a scalar twice should result in the original scalar.
+    #[allow(non_snake_case)]
+    #[test]
+    fn neg_twice_is_identity() {
+        let negative_X = -&X;
+        let should_be_X = -&negative_X;
+
+        assert_eq!(should_be_X, X);
     }
 }
 
@@ -727,6 +849,12 @@ mod bench {
         b.iter(|| Scalar::multiply_add(&X, &Y, &Z) );
     }
 
+    #[bench]
+    fn invert(b: &mut Bencher) {
+        let x = X.unpack();
+        b.iter(|| x.invert());
+    }
+
     #[bench]
     fn scalar_unpacked_multiply_add(b: &mut Bencher) {
         let x = X.unpack();