diff --git a/src/field.rs b/src/field.rs
index 8bd65ed..87b5c4c 100644
--- a/src/field.rs
+++ b/src/field.rs
@@ -33,49 +33,6 @@ use utils::{load3, load4, load8};
 
 use constants;
 
-/// doc
-pub fn mul64(a: &[u64;5], b: &[u64;5]) -> [u64;5] {
-    #[inline(always)]
-    fn m(a: u64, b: u64) -> u128 { (a as u128) * (b as u128) }
-    // Multiply to get 128-bit coefficients of output
-    let mut c0: u128 = m(a[0],b[0]) + ( m(a[4],b[1]) +   m(a[3],b[2]) +   m(a[2],b[3]) +   m(a[1],b[4]) )*19;
-    let mut c1: u128 = m(a[1],b[0]) +   m(a[0],b[1]) + ( m(a[4],b[2]) +   m(a[3],b[3]) +   m(a[2],b[4]) )*19;
-    let mut c2: u128 = m(a[2],b[0]) +   m(a[1],b[1]) +   m(a[0],b[2]) + ( m(a[4],b[3]) +   m(a[3],b[4]) )*19;
-    let mut c3: u128 = m(a[3],b[0]) +   m(a[2],b[1]) +   m(a[1],b[2]) +   m(a[0],b[3]) + ( m(a[4],b[4]) )*19;
-    let mut c4: u128 = m(a[4],b[0]) +   m(a[3],b[1]) +   m(a[2],b[2]) +   m(a[1],b[3]) +   m(a[0],b[4]);
-    // Now c[i] < 2^2b * (1+i + (4-i)*19) < 2^(2b + lg(1+4*19)) < 2^(2b + 6.27)
-    // where b is the bitlength of the input limbs.
-
-    // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54.
-    // After the first carry pass, all c[i] fit into u64.
-
-    let low_51_bit_mask = (1u64 << 51) - 1;
-    c1 +=  (c0 >> 51) as u128;
-    let mut c0: u64 = (c0 as u64) & low_51_bit_mask;
-    c2 +=  (c1 >> 51) as u128;
-    let mut c1: u64 = (c1 as u64) & low_51_bit_mask;
-    c3 +=  (c2 >> 51) as u128;
-    let mut c2: u64 = (c2 as u64) & low_51_bit_mask;
-    c4 +=  (c3 >> 51) as u128;
-    let mut c3: u64 = (c3 as u64) & low_51_bit_mask;
-    c0 += ((c4 >> 51) as u64) * 19;
-    let mut c4: u64 = (c4 as u64) & low_51_bit_mask;
-    
-    c1 +=  c0 >> 51;
-    c0 = c0 & low_51_bit_mask;
-    c2 +=  c1 >> 51;
-    c1 = c1 & low_51_bit_mask;
-    c3 +=  c2 >> 51;
-    c2 = c2 & low_51_bit_mask;
-    c4 +=  c3 >> 51;
-    c3 = c3 & low_51_bit_mask;
-    c0 += (c4 >> 51) * 19;
-    c4 = c4 & low_51_bit_mask;
-
-    [c0,c1,c2,c3,c4]
-}
-
-
 /// With the `radix51` feature enabled, `FieldElements` are represented
 /// in radix 2^51 as five `u64`s.
 #[cfg(feature="radix_51")]
@@ -852,7 +809,54 @@ impl FieldElement {
     /// Compute `self * _rhs`.
     #[cfg(feature="radix_51")]
     pub fn multiply(&self, _rhs: &FieldElement) -> FieldElement {
-        unimplemented!();
+        /// Multiply two 64-bit integers with 128 bits of output.
+        #[inline(always)]
+        fn m(x: u64, y: u64) -> u128 { (x as u128) * (y as u128) }
+
+        // Alias self, _rhs for more readable formulas
+        let a: &[u64; 5] = &self.0;
+        let b: &[u64; 5] = &_rhs.0;
+
+        // Multiply to get 128-bit coefficients of output
+        let mut c0: u128 = m(a[0],b[0]) + ( m(a[4],b[1]) +   m(a[3],b[2]) +   m(a[2],b[3]) +   m(a[1],b[4]) )*19;
+        let mut c1: u128 = m(a[1],b[0]) +   m(a[0],b[1]) + ( m(a[4],b[2]) +   m(a[3],b[3]) +   m(a[2],b[4]) )*19;
+        let mut c2: u128 = m(a[2],b[0]) +   m(a[1],b[1]) +   m(a[0],b[2]) + ( m(a[4],b[3]) +   m(a[3],b[4]) )*19;
+        let mut c3: u128 = m(a[3],b[0]) +   m(a[2],b[1]) +   m(a[1],b[2]) +   m(a[0],b[3]) + ( m(a[4],b[4]) )*19;
+        let mut c4: u128 = m(a[4],b[0]) +   m(a[3],b[1]) +   m(a[2],b[2]) +   m(a[1],b[3]) +   m(a[0],b[4]);
+
+        // Now c[i] < 2^2b * (1+i + (4-i)*19) < 2^(2b + lg(1+4*19)) < 2^(2b + 6.27)
+        // where b is the bitlength of the input limbs.
+
+        // The carry (c[i] >> 51) fits into a u64 iff 2b+6.27 < 64+51 iff b <= 54.
+        // After the first carry pass, all c[i] fit into u64.
+
+        // The 128-bit output limbs are stored in two 64-bit registers (low/high part).
+        // By rebinding the names after carrying, we free the upper registers for reuse.
+        let low_51_bit_mask = (1u64 << 51) - 1;
+        c1 +=  (c0 >> 51) as u128;
+        let mut c0: u64 = (c0 as u64) & low_51_bit_mask;
+        c2 +=  (c1 >> 51) as u128;
+        let mut c1: u64 = (c1 as u64) & low_51_bit_mask;
+        c3 +=  (c2 >> 51) as u128;
+        let mut c2: u64 = (c2 as u64) & low_51_bit_mask;
+        c4 +=  (c3 >> 51) as u128;
+        let mut c3: u64 = (c3 as u64) & low_51_bit_mask;
+        c0 += ((c4 >> 51) as u64) * 19;
+        let mut c4: u64 = (c4 as u64) & low_51_bit_mask;
+
+        // Second carry pass to enforce 2^51 bound
+        c1 +=  c0 >> 51;
+        c0 = c0 & low_51_bit_mask;
+        c2 +=  c1 >> 51;
+        c1 = c1 & low_51_bit_mask;
+        c3 +=  c2 >> 51;
+        c2 = c2 & low_51_bit_mask;
+        c4 +=  c3 >> 51;
+        c3 = c3 & low_51_bit_mask;
+        c0 += (c4 >> 51) * 19;
+        c4 = c4 & low_51_bit_mask;
+
+        FieldElement([c0,c1,c2,c3,c4])
     }
 
     #[cfg(feature="radix_25_5")]