diff --git a/src/scalar.rs b/src/scalar.rs
index 7e899d8..da555fb 100644
--- a/src/scalar.rs
+++ b/src/scalar.rs
@@ -731,6 +731,38 @@ mod test {
         assert_eq!(should_be_unpacked.0, unpacked.0);
     }
 
+
+    #[test]
+    fn montgomery_reduce_matches_reduce() {
+        let mut bignum = [0u8; 64];
+
+        // set bignum = x + 2^256x
+        for i in 0..32 {
+            bignum[   i] = X[i];
+            bignum[32+i] = X[i];
+        }
+        // x + 2^256x (mod l)
+        //         = 3958878930004874126169954872055634648693766179881526445624823978500314864344
+        let expected = Scalar([216, 154, 179, 139, 210, 121,   2,  71,
+                                69,  99, 158, 216,  23, 173,  63, 100,
+                               204,   0,  91,  50, 219, 153,  57, 249,
+                                28,  82,  31, 197, 100, 165, 192,   8]);
+        let reduced = Scalar::reduce(&bignum);
+
+        // The reduced scalar should match the expected
+        assert_eq!(reduced.0, expected.0);
+
+        //  (x + 2^256x) * R
+        let interim = UnpackedScalar::mul_internal(&UnpackedScalar::from_bytes_wide(&bignum),
+                                                   &constants::R);
+        // ((x + 2^256x) * R) / R  (mod l)
+        let montgomery_reduced = UnpackedScalar::montgomery_reduce(&interim);
+
+        // The Montgomery reduced scalar should match the reduced one, as well as the expected
+        assert_eq!(montgomery_reduced.0, reduced.unpack().0);
+        assert_eq!(montgomery_reduced.0, expected.unpack().0)
+    }
+
     #[cfg(feature = "serde")]
     use serde_cbor;
 
diff --git a/src/scalar_32bit.rs b/src/scalar_32bit.rs
index 3aa9f4c..238ee09 100644
--- a/src/scalar_32bit.rs
+++ b/src/scalar_32bit.rs
@@ -197,7 +197,7 @@ impl Scalar32 {
     ///
     /// This is implemented with a one-level refined Karatsuba decomposition
     #[inline(always)]
-    fn mul_internal(a: &Scalar32, b: &Scalar32) -> [u64; 17] {
+    pub (crate) fn mul_internal(a: &Scalar32, b: &Scalar32) -> [u64; 17] {
         let mut z = [0u64; 17];
 
         z[0] = m(a[0],b[0]);                                                             // c00
@@ -289,7 +289,7 @@ impl Scalar32 {
 
     /// Compute `limbs/R` (mod l), where R is the Montgomery modulus 2^261
     #[inline(always)]
-    fn montgomery_reduce(limbs: &[u64; 17]) -> Scalar32 {
+    pub (crate) fn montgomery_reduce(limbs: &[u64; 17]) -> Scalar32 {
 
         #[inline(always)]
         fn part1(sum: u64) -> (u64, u32) {
diff --git a/src/scalar_64bit.rs b/src/scalar_64bit.rs
index 8814732..b0ad405 100644
--- a/src/scalar_64bit.rs
+++ b/src/scalar_64bit.rs
@@ -181,18 +181,20 @@ impl Scalar64 {
 
     /// Compute `a * b`
     #[inline(always)]
-    fn mul_internal(a: &Scalar64, b: &Scalar64) -> [u128; 9] {
-        [
-            m(a[0],b[0]),
-            m(a[0],b[1]) + m(a[1],b[0]),
-            m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]),
-            m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]),
-            m(a[0],b[4]) + m(a[1],b[3]) + m(a[2],b[2]) + m(a[3],b[1]) + m(a[4],b[0]),
-                           m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]),
-                                          m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]),
-                                                         m(a[3],b[4]) + m(a[4],b[3]),
-                                                                        m(a[4],b[4])
-        ]
+    pub (crate) fn mul_internal(a: &Scalar64, b: &Scalar64) -> [u128; 9] {
+        let mut z = [0u128; 9];
+
+        z[0] = m(a[0],b[0]);
+        z[1] = m(a[0],b[1]) + m(a[1],b[0]);
+        z[2] = m(a[0],b[2]) + m(a[1],b[1]) + m(a[2],b[0]);
+        z[3] = m(a[0],b[3]) + m(a[1],b[2]) + m(a[2],b[1]) + m(a[3],b[0]);
+        z[4] = m(a[0],b[4]) + m(a[1],b[3]) + m(a[2],b[2]) + m(a[3],b[1]) + m(a[4],b[0]);
+        z[5] =                m(a[1],b[4]) + m(a[2],b[3]) + m(a[3],b[2]) + m(a[4],b[1]);
+        z[6] =                               m(a[2],b[4]) + m(a[3],b[3]) + m(a[4],b[2]);
+        z[7] =                                              m(a[3],b[4]) + m(a[4],b[3]);
+        z[8] =                                                             m(a[4],b[4]);
+
+        z
     }
 
     /// Compute `a^2`
@@ -220,7 +222,7 @@ impl Scalar64 {
 
     /// Compute `limbs/R` (mod l), where R is the Montgomery modulus 2^260
     #[inline(always)]
-    fn montgomery_reduce(limbs: &[u128; 9]) -> Scalar64 {
+    pub (crate) fn montgomery_reduce(limbs: &[u128; 9]) -> Scalar64 {
 
         #[inline(always)]
         fn part1(sum: u128) -> (u128, u64) {