diff --git a/src/fields/fp.rs b/src/fields/fp.rs
index 3264cba..5e17c99 100644
--- a/src/fields/fp.rs
+++ b/src/fields/fp.rs
@@ -308,33 +308,8 @@ impl Fp {
     /// Squares this element.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn square(&self) -> Fp {
-        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
-        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
-        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
-        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
-
-        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
-
-        let r7 = r6 >> 63;
-        let r6 = (r6 << 1) | (r5 >> 63);
-        let r5 = (r5 << 1) | (r4 >> 63);
-        let r4 = (r4 << 1) | (r3 >> 63);
-        let r3 = (r3 << 1) | (r2 >> 63);
-        let r2 = (r2 << 1) | (r1 >> 63);
-        let r1 = r1 << 1;
-
-        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
-        let (r1, carry) = adc(0, r1, carry);
-        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
-        let (r3, carry) = adc(0, r3, carry);
-        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
-        let (r5, carry) = adc(0, r5, carry);
-        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
-        let (r7, _) = adc(0, r7, carry);
-
-        Fp::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
+        let u = self.square_unreduced();
+        Fp::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
     }
 
     #[allow(clippy::too_many_arguments)]
@@ -388,29 +363,8 @@ impl Fp {
     /// Multiplies `rhs` by `self`, returning the result.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn mul(&self, rhs: &Self) -> Self {
-        // Schoolbook multiplication
-
-        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
-        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
-        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
-        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
-
-        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
-        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
-        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
-        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
-
-        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
-        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
-        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
-        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
-        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
-        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
-        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
-
-        Fp::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
+        let u = self.mul_unreduced(rhs);
+        Fp::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
     }
 
     /// Subtracts `rhs` from `self`, returning the result.
@@ -461,6 +415,66 @@ impl Fp {
 
         Fp([d0 & mask, d1 & mask, d2 & mask, d3 & mask])
     }
+
+    /// Multiplies `rhs` by `self`, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn mul_unreduced(&self, rhs: &Self) -> [u64; 8] {
+        // Schoolbook multiplication
+
+        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
+        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
+        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
+
+        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
+        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
+        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
+        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
+
+        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
+        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
+        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
+        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
+        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
+        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
+        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
+
+    /// Squares this element, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn square_unreduced(&self) -> [u64; 8] {
+        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
+        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
+        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
+
+        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
+
+        let r7 = r6 >> 63;
+        let r6 = (r6 << 1) | (r5 >> 63);
+        let r5 = (r5 << 1) | (r4 >> 63);
+        let r4 = (r4 << 1) | (r3 >> 63);
+        let r3 = (r3 << 1) | (r2 >> 63);
+        let r2 = (r2 << 1) | (r1 >> 63);
+        let r1 = r1 << 1;
+
+        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
+        let (r1, carry) = adc(0, r1, carry);
+        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
+        let (r3, carry) = adc(0, r3, carry);
+        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
+        let (r5, carry) = adc(0, r5, carry);
+        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
+        let (r7, _) = adc(0, r7, carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
 }
 
 impl From<Fp> for [u8; 32] {
diff --git a/src/fields/fq.rs b/src/fields/fq.rs
index 8177fa4..acb6914 100644
--- a/src/fields/fq.rs
+++ b/src/fields/fq.rs
@@ -308,33 +308,8 @@ impl Fq {
     /// Squares this element.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn square(&self) -> Fq {
-        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
-        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
-        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
-        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
-
-        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
-
-        let r7 = r6 >> 63;
-        let r6 = (r6 << 1) | (r5 >> 63);
-        let r5 = (r5 << 1) | (r4 >> 63);
-        let r4 = (r4 << 1) | (r3 >> 63);
-        let r3 = (r3 << 1) | (r2 >> 63);
-        let r2 = (r2 << 1) | (r1 >> 63);
-        let r1 = r1 << 1;
-
-        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
-        let (r1, carry) = adc(0, r1, carry);
-        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
-        let (r3, carry) = adc(0, r3, carry);
-        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
-        let (r5, carry) = adc(0, r5, carry);
-        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
-        let (r7, _) = adc(0, r7, carry);
-
-        Fq::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
+        let u = self.square_unreduced();
+        Fq::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
     }
 
     #[allow(clippy::too_many_arguments)]
@@ -388,29 +363,8 @@ impl Fq {
     /// Multiplies `rhs` by `self`, returning the result.
     #[cfg_attr(not(feature = "uninline-portable"), inline)]
     pub const fn mul(&self, rhs: &Self) -> Self {
-        // Schoolbook multiplication
-
-        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
-        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
-        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
-        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
-
-        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
-        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
-        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
-        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
-
-        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
-        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
-        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
-        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
-
-        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
-        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
-        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
-        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
-
-        Fq::montgomery_reduce(r0, r1, r2, r3, r4, r5, r6, r7)
+        let u = self.mul_unreduced(rhs);
+        Fq::montgomery_reduce(u[0], u[1], u[2], u[3], u[4], u[5], u[6], u[7])
     }
 
     /// Subtracts `rhs` from `self`, returning the result.
@@ -461,6 +415,66 @@ impl Fq {
 
         Fq([d0 & mask, d1 & mask, d2 & mask, d3 & mask])
     }
+
+    /// Multiplies `rhs` by `self`, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn mul_unreduced(&self, rhs: &Self) -> [u64; 8] {
+        // Schoolbook multiplication
+
+        let (r0, carry) = mac(0, self.0[0], rhs.0[0], 0);
+        let (r1, carry) = mac(0, self.0[0], rhs.0[1], carry);
+        let (r2, carry) = mac(0, self.0[0], rhs.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], rhs.0[3], carry);
+
+        let (r1, carry) = mac(r1, self.0[1], rhs.0[0], 0);
+        let (r2, carry) = mac(r2, self.0[1], rhs.0[1], carry);
+        let (r3, carry) = mac(r3, self.0[1], rhs.0[2], carry);
+        let (r4, r5) = mac(r4, self.0[1], rhs.0[3], carry);
+
+        let (r2, carry) = mac(r2, self.0[2], rhs.0[0], 0);
+        let (r3, carry) = mac(r3, self.0[2], rhs.0[1], carry);
+        let (r4, carry) = mac(r4, self.0[2], rhs.0[2], carry);
+        let (r5, r6) = mac(r5, self.0[2], rhs.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[3], rhs.0[0], 0);
+        let (r4, carry) = mac(r4, self.0[3], rhs.0[1], carry);
+        let (r5, carry) = mac(r5, self.0[3], rhs.0[2], carry);
+        let (r6, r7) = mac(r6, self.0[3], rhs.0[3], carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
+
+    /// Squares this element, returning the unreduced 512-bit product.
+    #[cfg_attr(not(feature = "uninline-portable"), inline)]
+    pub(crate) const fn square_unreduced(&self) -> [u64; 8] {
+        let (r1, carry) = mac(0, self.0[0], self.0[1], 0);
+        let (r2, carry) = mac(0, self.0[0], self.0[2], carry);
+        let (r3, r4) = mac(0, self.0[0], self.0[3], carry);
+
+        let (r3, carry) = mac(r3, self.0[1], self.0[2], 0);
+        let (r4, r5) = mac(r4, self.0[1], self.0[3], carry);
+
+        let (r5, r6) = mac(r5, self.0[2], self.0[3], 0);
+
+        let r7 = r6 >> 63;
+        let r6 = (r6 << 1) | (r5 >> 63);
+        let r5 = (r5 << 1) | (r4 >> 63);
+        let r4 = (r4 << 1) | (r3 >> 63);
+        let r3 = (r3 << 1) | (r2 >> 63);
+        let r2 = (r2 << 1) | (r1 >> 63);
+        let r1 = r1 << 1;
+
+        let (r0, carry) = mac(0, self.0[0], self.0[0], 0);
+        let (r1, carry) = adc(0, r1, carry);
+        let (r2, carry) = mac(r2, self.0[1], self.0[1], carry);
+        let (r3, carry) = adc(0, r3, carry);
+        let (r4, carry) = mac(r4, self.0[2], self.0[2], carry);
+        let (r5, carry) = adc(0, r5, carry);
+        let (r6, carry) = mac(r6, self.0[3], self.0[3], carry);
+        let (r7, _) = adc(0, r7, carry);
+
+        [r0, r1, r2, r3, r4, r5, r6, r7]
+    }
 }
 
 impl From<Fq> for [u8; 32] {