diff --git a/bellman/src/domain.rs b/bellman/src/domain.rs
index a1e4a4c..6e8a6a7 100644
--- a/bellman/src/domain.rs
+++ b/bellman/src/domain.rs
@@ -63,7 +63,7 @@ impl<E: ScalarEngine, G: Group<E>> EvaluationDomain<E, G> {
         // Compute omega, the 2^exp primitive root of unity
         let mut omega = E::Fr::root_of_unity();
         for _ in exp..E::Fr::S {
-            omega.square();
+            omega = omega.square();
         }
 
         // Extend the coeffs vector with zeroes if necessary
diff --git a/bellman/src/gadgets/num.rs b/bellman/src/gadgets/num.rs
index 08a175b..65bee25 100644
--- a/bellman/src/gadgets/num.rs
+++ b/bellman/src/gadgets/num.rs
@@ -254,8 +254,7 @@ impl<E: ScalarEngine> AllocatedNum<E> {
         let var = cs.alloc(
             || "squared num",
             || {
-                let mut tmp = *self.value.get()?;
-                tmp.square();
+                let tmp = self.value.get()?.square();
 
                 value = Some(tmp);
 
diff --git a/bellman/src/groth16/tests/dummy_engine.rs b/bellman/src/groth16/tests/dummy_engine.rs
index 224f84e..0e722f2 100644
--- a/bellman/src/groth16/tests/dummy_engine.rs
+++ b/bellman/src/groth16/tests/dummy_engine.rs
@@ -140,8 +140,8 @@ impl Field for Fr {
         (self.0).0 == 0
     }
 
-    fn square(&mut self) {
-        self.0 = (self.0 * self.0) % MODULUS_R;
+    fn square(&self) -> Self {
+        Fr((self.0 * self.0) % MODULUS_R)
     }
 
     fn double(&self) -> Self {
@@ -191,22 +191,21 @@ impl SqrtField for Fr {
                 while t != <Fr as Field>::one() {
                     let mut i = 1;
                     {
-                        let mut t2i = t;
-                        t2i.square();
+                        let mut t2i = t.square();
                         loop {
                             if t2i == <Fr as Field>::one() {
                                 break;
                             }
-                            t2i.square();
+                            t2i = t2i.square();
                             i += 1;
                         }
                     }
 
                     for _ in 0..(m - i - 1) {
-                        c.square();
+                        c = c.square();
                     }
                     MulAssign::mul_assign(&mut r, &c);
-                    c.square();
+                    c = c.square();
                     MulAssign::mul_assign(&mut t, &c);
                     m = i;
                 }
diff --git a/bellman/tests/mimc.rs b/bellman/tests/mimc.rs
index 0792af5..a1de0f1 100644
--- a/bellman/tests/mimc.rs
+++ b/bellman/tests/mimc.rs
@@ -41,8 +41,7 @@ fn mimc<E: Engine>(mut xl: E::Fr, mut xr: E::Fr, constants: &[E::Fr]) -> E::Fr {
     for i in 0..MIMC_ROUNDS {
         let mut tmp1 = xl;
         tmp1.add_assign(&constants[i]);
-        let mut tmp2 = tmp1;
-        tmp2.square();
+        let mut tmp2 = tmp1.square();
         tmp2.mul_assign(&tmp1);
         tmp2.add_assign(&xr);
         xr = xl;
@@ -88,8 +87,7 @@ impl<'a, E: Engine> Circuit<E> for MiMCDemo<'a, E> {
             // tmp = (xL + Ci)^2
             let tmp_value = xl_value.map(|mut e| {
                 e.add_assign(&self.constants[i]);
-                e.square();
-                e
+                e.square()
             });
             let tmp = cs.alloc(
                 || "tmp",
diff --git a/ff/ff_derive/src/lib.rs b/ff/ff_derive/src/lib.rs
index 88146b3..73b560b 100644
--- a/ff/ff_derive/src/lib.rs
+++ b/ff/ff_derive/src/lib.rs
@@ -447,8 +447,7 @@ fn prime_field_constants_and_sqrt(
 
                         let mut a1 = self.pow(#mod_minus_3_over_4);
 
-                        let mut a0 = a1;
-                        a0.square();
+                        let mut a0 = a1.square();
                         a0.mul_assign(self);
 
                         if a0.0 == #repr(#rneg) {
@@ -484,22 +483,21 @@ fn prime_field_constants_and_sqrt(
                                 while t != Self::one() {
                                     let mut i = 1;
                                     {
-                                        let mut t2i = t;
-                                        t2i.square();
+                                        let mut t2i = t.square();
                                         loop {
                                             if t2i == Self::one() {
                                                 break;
                                             }
-                                            t2i.square();
+                                            t2i = t2i.square();
                                             i += 1;
                                         }
                                     }
 
                                     for _ in 0..(m - i - 1) {
-                                        c.square();
+                                        c = c.square();
                                     }
                                     r.mul_assign(&c);
-                                    c.square();
+                                    c = c.square();
                                     t.mul_assign(&c);
                                     m = i;
                                 }
@@ -715,7 +713,9 @@ fn prime_field_impl(
         );
 
         gen.extend(quote! {
-            self.mont_reduce(#mont_calling);
+            let mut ret = *self;
+            ret.mont_reduce(#mont_calling);
+            ret
         });
 
         gen
@@ -1113,7 +1113,7 @@ fn prime_field_impl(
             }
 
             #[inline]
-            fn square(&mut self)
+            fn square(&self) -> Self
             {
                 #squaring_impl
             }
diff --git a/ff/src/lib.rs b/ff/src/lib.rs
index fd8d45d..10975a3 100644
--- a/ff/src/lib.rs
+++ b/ff/src/lib.rs
@@ -51,7 +51,8 @@ pub trait Field:
     fn is_zero(&self) -> bool;
 
     /// Squares this element.
-    fn square(&mut self);
+    #[must_use]
+    fn square(&self) -> Self;
 
     /// Doubles this element.
     #[must_use]
@@ -73,7 +74,7 @@ pub trait Field:
 
         for i in BitIterator::new(exp) {
             if found_one {
-                res.square();
+                res = res.square();
             } else {
                 found_one = i;
             }
diff --git a/pairing/benches/bls12_381/fq.rs b/pairing/benches/bls12_381/fq.rs
index 046491e..bd4ed15 100644
--- a/pairing/benches/bls12_381/fq.rs
+++ b/pairing/benches/bls12_381/fq.rs
@@ -210,8 +210,7 @@ fn bench_fq_square(b: &mut ::test::Bencher) {
 
     let mut count = 0;
     b.iter(|| {
-        let mut tmp = v[count];
-        tmp.square();
+        let tmp = v[count].square();
         count = (count + 1) % SAMPLES;
         tmp
     });
@@ -264,11 +263,7 @@ fn bench_fq_sqrt(b: &mut ::test::Bencher) {
     ]);
 
     let v: Vec<Fq> = (0..SAMPLES)
-        .map(|_| {
-            let mut tmp = Fq::random(&mut rng);
-            tmp.square();
-            tmp
-        })
+        .map(|_| Fq::random(&mut rng).square())
         .collect();
 
     let mut count = 0;
diff --git a/pairing/benches/bls12_381/fq12.rs b/pairing/benches/bls12_381/fq12.rs
index 0eab387..b79bf5c 100644
--- a/pairing/benches/bls12_381/fq12.rs
+++ b/pairing/benches/bls12_381/fq12.rs
@@ -84,8 +84,7 @@ fn bench_fq12_squaring(b: &mut ::test::Bencher) {
 
     let mut count = 0;
     b.iter(|| {
-        let mut tmp = v[count];
-        tmp.square();
+        let tmp = v[count].square();
         count = (count + 1) % SAMPLES;
         tmp
     });
diff --git a/pairing/benches/bls12_381/fq2.rs b/pairing/benches/bls12_381/fq2.rs
index 7862ea3..ace0171 100644
--- a/pairing/benches/bls12_381/fq2.rs
+++ b/pairing/benches/bls12_381/fq2.rs
@@ -84,8 +84,7 @@ fn bench_fq2_squaring(b: &mut ::test::Bencher) {
 
     let mut count = 0;
     b.iter(|| {
-        let mut tmp = v[count];
-        tmp.square();
+        let tmp = v[count].square();
         count = (count + 1) % SAMPLES;
         tmp
     });
diff --git a/pairing/benches/bls12_381/fr.rs b/pairing/benches/bls12_381/fr.rs
index 16a5c08..e4d07d2 100644
--- a/pairing/benches/bls12_381/fr.rs
+++ b/pairing/benches/bls12_381/fr.rs
@@ -210,8 +210,7 @@ fn bench_fr_square(b: &mut ::test::Bencher) {
 
     let mut count = 0;
     b.iter(|| {
-        let mut tmp = v[count];
-        tmp.square();
+        let tmp = v[count].square();
         count = (count + 1) % SAMPLES;
         tmp
     });
@@ -264,11 +263,7 @@ fn bench_fr_sqrt(b: &mut ::test::Bencher) {
     ]);
 
     let v: Vec<Fr> = (0..SAMPLES)
-        .map(|_| {
-            let mut tmp = Fr::random(&mut rng);
-            tmp.square();
-            tmp
-        })
+        .map(|_| Fr::random(&mut rng).square())
         .collect();
 
     let mut count = 0;
diff --git a/pairing/src/bls12_381/ec.rs b/pairing/src/bls12_381/ec.rs
index 6f45592..abf6873 100644
--- a/pairing/src/bls12_381/ec.rs
+++ b/pairing/src/bls12_381/ec.rs
@@ -54,10 +54,8 @@ macro_rules! curve_impl {
                 // are equal when (X * Z^2) = (X' * Z'^2)
                 // and (Y * Z^3) = (Y' * Z'^3).
 
-                let mut z1 = self.z;
-                z1.square();
-                let mut z2 = other.z;
-                z2.square();
+                let mut z1 = self.z.square();
+                let mut z2 = other.z.square();
 
                 let mut tmp1 = self.x;
                 tmp1.mul_assign(&z2);
@@ -101,8 +99,7 @@ macro_rules! curve_impl {
             /// largest y-coordinate be selected.
             fn get_point_from_x(x: $basefield, greatest: bool) -> Option<$affine> {
                 // Compute x^3 + b
-                let mut x3b = x;
-                x3b.square();
+                let mut x3b = x.square();
                 x3b.mul_assign(&x);
                 x3b.add_assign(&$affine::get_coeff_b());
 
@@ -122,11 +119,9 @@ macro_rules! curve_impl {
                     true
                 } else {
                     // Check that the point is on the curve
-                    let mut y2 = self.y;
-                    y2.square();
+                    let y2 = self.y.square();
 
-                    let mut x3b = self.x;
-                    x3b.square();
+                    let mut x3b = self.x.square();
                     x3b.mul_assign(&self.x);
                     x3b.add_assign(&Self::get_coeff_b());
 
@@ -283,8 +278,7 @@ macro_rules! curve_impl {
 
                 // Perform affine transformations
                 for g in v.iter_mut().filter(|g| !g.is_normalized()) {
-                    let mut z = g.z; // 1/z
-                    z.square(); // 1/z^2
+                    let mut z = g.z.square(); // 1/z^2
                     g.x.mul_assign(&z); // x/z^2
                     z.mul_assign(&g.z); // 1/z^3
                     g.y.mul_assign(&z); // y/z^3
@@ -305,21 +299,18 @@ macro_rules! curve_impl {
                 // http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
 
                 // A = X1^2
-                let mut a = self.x;
-                a.square();
+                let a = self.x.square();
 
                 // B = Y1^2
-                let mut b = self.y;
-                b.square();
+                let b = self.y.square();
 
                 // C = B^2
-                let mut c = b;
-                c.square();
+                let mut c = b.square();
 
                 // D = 2*((X1+B)2-A-C)
                 let mut d = self.x;
                 d.add_assign(&b);
-                d.square();
+                d = d.square();
                 d.sub_assign(&a);
                 d.sub_assign(&c);
                 d = d.double();
@@ -329,8 +320,7 @@ macro_rules! curve_impl {
                 e.add_assign(&a);
 
                 // F = E^2
-                let mut f = e;
-                f.square();
+                let f = e.square();
 
                 // Z3 = 2*Y1*Z1
                 self.z.mul_assign(&self.y);
@@ -362,12 +352,10 @@ macro_rules! curve_impl {
                 // http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl
 
                 // Z1Z1 = Z1^2
-                let mut z1z1 = self.z;
-                z1z1.square();
+                let z1z1 = self.z.square();
 
                 // Z2Z2 = Z2^2
-                let mut z2z2 = other.z;
-                z2z2.square();
+                let z2z2 = other.z.square();
 
                 // U1 = X1*Z2Z2
                 let mut u1 = self.x;
@@ -398,8 +386,7 @@ macro_rules! curve_impl {
                     h.sub_assign(&u1);
 
                     // I = (2*H)^2
-                    let mut i = h.double();
-                    i.square();
+                    let i = h.double().square();
 
                     // J = H*I
                     let mut j = h;
@@ -415,8 +402,7 @@ macro_rules! curve_impl {
                     v.mul_assign(&i);
 
                     // X3 = r^2 - J - 2*V
-                    self.x = r;
-                    self.x.square();
+                    self.x = r.square();
                     self.x.sub_assign(&j);
                     self.x.sub_assign(&v);
                     self.x.sub_assign(&v);
@@ -431,7 +417,7 @@ macro_rules! curve_impl {
 
                     // Z3 = ((Z1+Z2)^2 - Z1Z1 - Z2Z2)*H
                     self.z.add_assign(&other.z);
-                    self.z.square();
+                    self.z = self.z.square();
                     self.z.sub_assign(&z1z1);
                     self.z.sub_assign(&z2z2);
                     self.z.mul_assign(&h);
@@ -453,8 +439,7 @@ macro_rules! curve_impl {
                 // http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-madd-2007-bl
 
                 // Z1Z1 = Z1^2
-                let mut z1z1 = self.z;
-                z1z1.square();
+                let z1z1 = self.z.square();
 
                 // U2 = X2*Z1Z1
                 let mut u2 = other.x;
@@ -476,8 +461,7 @@ macro_rules! curve_impl {
                     h.sub_assign(&self.x);
 
                     // HH = H^2
-                    let mut hh = h;
-                    hh.square();
+                    let hh = h.square();
 
                     // I = 4*HH
                     let i = hh.double().double();
@@ -496,8 +480,7 @@ macro_rules! curve_impl {
                     v.mul_assign(&i);
 
                     // X3 = r^2 - J - 2*V
-                    self.x = r;
-                    self.x.square();
+                    self.x = r.square();
                     self.x.sub_assign(&j);
                     self.x.sub_assign(&v);
                     self.x.sub_assign(&v);
@@ -512,7 +495,7 @@ macro_rules! curve_impl {
 
                     // Z3 = (Z1+H)^2-Z1Z1-HH
                     self.z.add_assign(&h);
-                    self.z.square();
+                    self.z = self.z.square();
                     self.z.sub_assign(&z1z1);
                     self.z.sub_assign(&hh);
                 }
@@ -589,8 +572,7 @@ macro_rules! curve_impl {
                 } else {
                     // Z is nonzero, so it must have an inverse in a field.
                     let zinv = p.z.inverse().unwrap();
-                    let mut zinv_powered = zinv;
-                    zinv_powered.square();
+                    let mut zinv_powered = zinv.square();
 
                     // X/Z^2
                     let mut x = p.x;
@@ -933,8 +915,7 @@ pub mod g1 {
         let mut i = 0;
         loop {
             // y^2 = x^3 + b
-            let mut rhs = x;
-            rhs.square();
+            let mut rhs = x.square();
             rhs.mul_assign(&x);
             rhs.add_assign(&G1Affine::get_coeff_b());
 
@@ -1638,8 +1619,7 @@ pub mod g2 {
         let mut i = 0;
         loop {
             // y^2 = x^3 + b
-            let mut rhs = x;
-            rhs.square();
+            let mut rhs = x.square();
             rhs.mul_assign(&x);
             rhs.add_assign(&G2Affine::get_coeff_b());
 
diff --git a/pairing/src/bls12_381/fq.rs b/pairing/src/bls12_381/fq.rs
index 77c85ea..28326f4 100644
--- a/pairing/src/bls12_381/fq.rs
+++ b/pairing/src/bls12_381/fq.rs
@@ -1930,7 +1930,7 @@ fn test_fq_mul_assign() {
 
 #[test]
 fn test_fq_squaring() {
-    let mut a = Fq(FqRepr([
+    let a = Fq(FqRepr([
         0xffffffffffffffff,
         0xffffffffffffffff,
         0xffffffffffffffff,
@@ -1939,9 +1939,8 @@ fn test_fq_squaring() {
         0x19ffffffffffffff,
     ]));
     assert!(a.is_valid());
-    a.square();
     assert_eq!(
-        a,
+        a.square(),
         Fq::from_repr(FqRepr([
             0x1cfb28fe7dfbbb86,
             0x24cbe1731577a59,
@@ -1961,14 +1960,7 @@ fn test_fq_squaring() {
     for _ in 0..1000000 {
         // Ensure that (a * a) = a^2
         let a = Fq::random(&mut rng);
-
-        let mut tmp = a;
-        tmp.square();
-
-        let mut tmp2 = a;
-        tmp2.mul_assign(&a);
-
-        assert_eq!(tmp, tmp2);
+        assert_eq!(a.square(), a * a);
     }
 }
 
@@ -2071,8 +2063,7 @@ fn test_fq_sqrt() {
         // Ensure sqrt(a^2) = a or -a
         let a = Fq::random(&mut rng);
         let nega = a.neg();
-        let mut b = a;
-        b.square();
+        let b = a.square();
 
         let b = b.sqrt().unwrap();
 
@@ -2083,10 +2074,8 @@ fn test_fq_sqrt() {
         // Ensure sqrt(a)^2 = a for random a
         let a = Fq::random(&mut rng);
 
-        if let Some(mut tmp) = a.sqrt() {
-            tmp.square();
-
-            assert_eq!(a, tmp);
+        if let Some(tmp) = a.sqrt() {
+            assert_eq!(a, tmp.square());
         }
     }
 }
diff --git a/pairing/src/bls12_381/fq12.rs b/pairing/src/bls12_381/fq12.rs
index d072725..77670e4 100644
--- a/pairing/src/bls12_381/fq12.rs
+++ b/pairing/src/bls12_381/fq12.rs
@@ -199,7 +199,7 @@ impl Field for Fq12 {
         self.c1.c2.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
     }
 
-    fn square(&mut self) {
+    fn square(&self) -> Self {
         let mut ab = self.c0;
         ab.mul_assign(&self.c1);
         let mut c0c1 = self.c0;
@@ -209,18 +209,16 @@ impl Field for Fq12 {
         c0.add_assign(&self.c0);
         c0.mul_assign(&c0c1);
         c0.sub_assign(&ab);
-        self.c1 = ab;
-        self.c1.add_assign(&ab);
+        let mut c1 = ab;
+        c1.add_assign(&ab);
         ab.mul_by_nonresidue();
         c0.sub_assign(&ab);
-        self.c0 = c0;
+        Fq12 { c0, c1 }
     }
 
     fn inverse(&self) -> Option<Self> {
-        let mut c0s = self.c0;
-        c0s.square();
-        let mut c1s = self.c1;
-        c1s.square();
+        let mut c0s = self.c0.square();
+        let mut c1s = self.c1.square();
         c1s.mul_by_nonresidue();
         c0s.sub_assign(&c1s);
 
diff --git a/pairing/src/bls12_381/fq2.rs b/pairing/src/bls12_381/fq2.rs
index 137c1ba..f271913 100644
--- a/pairing/src/bls12_381/fq2.rs
+++ b/pairing/src/bls12_381/fq2.rs
@@ -46,10 +46,8 @@ impl Fq2 {
 
     /// Norm of Fq2 as extension field in i over Fq
     pub fn norm(&self) -> Fq {
-        let mut t0 = self.c0;
-        let mut t1 = self.c1;
-        t0.square();
-        t1.square();
+        let t0 = self.c0.square();
+        let mut t1 = self.c1.square();
         t1.add_assign(&t0);
 
         t1
@@ -198,7 +196,7 @@ impl Field for Fq2 {
         self.c0.is_zero() && self.c1.is_zero()
     }
 
-    fn square(&mut self) {
+    fn square(&self) -> Self {
         let mut ab = self.c0;
         ab.mul_assign(&self.c1);
         let mut c0c1 = self.c0;
@@ -207,10 +205,10 @@ impl Field for Fq2 {
         c0.add_assign(&self.c0);
         c0.mul_assign(&c0c1);
         c0.sub_assign(&ab);
-        self.c1 = ab;
-        self.c1.add_assign(&ab);
+        let mut c1 = ab;
+        c1.add_assign(&ab);
         c0.add_assign(&ab);
-        self.c0 = c0;
+        Fq2 { c0, c1 }
     }
 
     fn double(&self) -> Self {
@@ -221,10 +219,8 @@ impl Field for Fq2 {
     }
 
     fn inverse(&self) -> Option<Self> {
-        let mut t1 = self.c1;
-        t1.square();
-        let mut t0 = self.c0;
-        t0.square();
+        let t1 = self.c1.square();
+        let mut t0 = self.c0.square();
         t0.add_assign(&t1);
         t0.inverse().map(|t| Fq2 {
             c0: self.c0.mul(&t),
@@ -257,8 +253,7 @@ impl SqrtField for Fq2 {
                 0x92c6e9ed90d2eb35,
                 0x680447a8e5ff9a6,
             ]);
-            let mut alpha = a1;
-            alpha.square();
+            let mut alpha = a1.square();
             alpha.mul_assign(self);
             let mut a0 = alpha;
             a0.frobenius_map(1);
@@ -353,32 +348,30 @@ fn test_fq2_squaring() {
     use super::fq::FqRepr;
     use ff::PrimeField;
 
-    let mut a = Fq2 {
+    let a = Fq2 {
         c0: Fq::one(),
         c1: Fq::one(),
     }; // u + 1
-    a.square();
     assert_eq!(
-        a,
+        a.square(),
         Fq2 {
             c0: Fq::zero(),
             c1: Fq::from_repr(FqRepr::from(2)).unwrap(),
         }
     ); // 2u
 
-    let mut a = Fq2 {
+    let a = Fq2 {
         c0: Fq::zero(),
         c1: Fq::one(),
     }; // u
-    a.square();
-    assert_eq!(a, {
+    assert_eq!(a.square(), {
         Fq2 {
             c0: Fq::one().neg(),
             c1: Fq::zero(),
         }
     }); // -1
 
-    let mut a = Fq2 {
+    let a = Fq2 {
         c0: Fq::from_repr(FqRepr([
             0x9c2c6309bbf8b598,
             0x4eef5c946536f602,
@@ -398,9 +391,8 @@ fn test_fq2_squaring() {
         ]))
         .unwrap(),
     };
-    a.square();
     assert_eq!(
-        a,
+        a.square(),
         Fq2 {
             c0: Fq::from_repr(FqRepr([
                 0xf262c28c538bcf68,
diff --git a/pairing/src/bls12_381/fq6.rs b/pairing/src/bls12_381/fq6.rs
index 8e5b394..f8cb619 100644
--- a/pairing/src/bls12_381/fq6.rs
+++ b/pairing/src/bls12_381/fq6.rs
@@ -303,35 +303,35 @@ impl Field for Fq6 {
         self.c2.mul_assign(&FROBENIUS_COEFF_FQ6_C2[power % 6]);
     }
 
-    fn square(&mut self) {
-        let mut s0 = self.c0;
-        s0.square();
+    fn square(&self) -> Self {
+        let s0 = self.c0.square();
         let mut ab = self.c0;
         ab.mul_assign(&self.c1);
         let s1 = ab.double();
         let mut s2 = self.c0;
         s2.sub_assign(&self.c1);
         s2.add_assign(&self.c2);
-        s2.square();
+        s2 = s2.square();
         let mut bc = self.c1;
         bc.mul_assign(&self.c2);
         let s3 = bc.double();
-        let mut s4 = self.c2;
-        s4.square();
+        let s4 = self.c2.square();
 
-        self.c0 = s3;
-        self.c0.mul_by_nonresidue();
-        self.c0.add_assign(&s0);
+        let mut c0 = s3;
+        c0.mul_by_nonresidue();
+        c0.add_assign(&s0);
 
-        self.c1 = s4;
-        self.c1.mul_by_nonresidue();
-        self.c1.add_assign(&s1);
+        let mut c1 = s4;
+        c1.mul_by_nonresidue();
+        c1.add_assign(&s1);
 
-        self.c2 = s1;
-        self.c2.add_assign(&s2);
-        self.c2.add_assign(&s3);
-        self.c2.sub_assign(&s0);
-        self.c2.sub_assign(&s4);
+        let mut c2 = s1;
+        c2.add_assign(&s2);
+        c2.add_assign(&s3);
+        c2.sub_assign(&s0);
+        c2.sub_assign(&s4);
+
+        Fq6 { c0, c1, c2 }
     }
 
     fn inverse(&self) -> Option<Self> {
@@ -340,20 +340,17 @@ impl Field for Fq6 {
         c0.mul_assign(&self.c1);
         c0 = c0.neg();
         {
-            let mut c0s = self.c0;
-            c0s.square();
+            let c0s = self.c0.square();
             c0.add_assign(&c0s);
         }
-        let mut c1 = self.c2;
-        c1.square();
+        let mut c1 = self.c2.square();
         c1.mul_by_nonresidue();
         {
             let mut c01 = self.c0;
             c01.mul_assign(&self.c1);
             c1.sub_assign(&c01);
         }
-        let mut c2 = self.c1;
-        c2.square();
+        let mut c2 = self.c1.square();
         {
             let mut c02 = self.c0;
             c02.mul_assign(&self.c2);
diff --git a/pairing/src/bls12_381/fr.rs b/pairing/src/bls12_381/fr.rs
index 21832de..777aa87 100644
--- a/pairing/src/bls12_381/fr.rs
+++ b/pairing/src/bls12_381/fr.rs
@@ -693,16 +693,15 @@ fn test_fr_mul_assign() {
 
 #[test]
 fn test_fr_squaring() {
-    let mut a = Fr(FrRepr([
+    let a = Fr(FrRepr([
         0xffffffffffffffff,
         0xffffffffffffffff,
         0xffffffffffffffff,
         0x73eda753299d7d47,
     ]));
     assert!(a.is_valid());
-    a.square();
     assert_eq!(
-        a,
+        a.square(),
         Fr::from_repr(FrRepr([
             0xc0d698e7bde077b8,
             0xb79a310579e76ec2,
@@ -720,14 +719,7 @@ fn test_fr_squaring() {
     for _ in 0..1000000 {
         // Ensure that (a * a) = a^2
         let a = Fr::random(&mut rng);
-
-        let mut tmp = a;
-        tmp.square();
-
-        let mut tmp2 = a;
-        tmp2.mul_assign(&a);
-
-        assert_eq!(tmp, tmp2);
+        assert_eq!(a.square(), a * a);
     }
 }
 
@@ -830,8 +822,7 @@ fn test_fr_sqrt() {
         // Ensure sqrt(a^2) = a or -a
         let a = Fr::random(&mut rng);
         let nega = a.neg();
-        let mut b = a;
-        b.square();
+        let b = a.square();
 
         let b = b.sqrt().unwrap();
 
@@ -842,10 +833,8 @@ fn test_fr_sqrt() {
         // Ensure sqrt(a)^2 = a for random a
         let a = Fr::random(&mut rng);
 
-        if let Some(mut tmp) = a.sqrt() {
-            tmp.square();
-
-            assert_eq!(a, tmp);
+        if let Some(tmp) = a.sqrt() {
+            assert_eq!(a, tmp.square());
         }
     }
 }
diff --git a/pairing/src/bls12_381/mod.rs b/pairing/src/bls12_381/mod.rs
index 0a3813d..5843e94 100644
--- a/pairing/src/bls12_381/mod.rs
+++ b/pairing/src/bls12_381/mod.rs
@@ -97,7 +97,7 @@ impl Engine for Bls12 {
                 }
             }
 
-            f.square();
+            f = f.square();
         }
 
         for &mut (p, ref mut coeffs) in &mut pairs {
@@ -131,8 +131,7 @@ impl Engine for Bls12 {
                 }
 
                 let mut x = BLS_X;
-                let mut y0 = r;
-                y0.square();
+                let y0 = r.square();
                 let mut y1 = y0;
                 exp_by_x(&mut y1, x);
                 x >>= 1;
@@ -185,18 +184,15 @@ impl G2Prepared {
 
         fn doubling_step(r: &mut G2) -> (Fq2, Fq2, Fq2) {
             // Adaptation of Algorithm 26, https://eprint.iacr.org/2010/354.pdf
-            let mut tmp0 = r.x;
-            tmp0.square();
+            let mut tmp0 = r.x.square();
 
-            let mut tmp1 = r.y;
-            tmp1.square();
+            let mut tmp1 = r.y.square();
 
-            let mut tmp2 = tmp1;
-            tmp2.square();
+            let mut tmp2 = tmp1.square();
 
             let mut tmp3 = tmp1;
             tmp3.add_assign(&r.x);
-            tmp3.square();
+            tmp3 = tmp3.square();
             tmp3.sub_assign(&tmp0);
             tmp3.sub_assign(&tmp2);
             tmp3 = tmp3.double();
@@ -207,18 +203,16 @@ impl G2Prepared {
             let mut tmp6 = r.x;
             tmp6.add_assign(&tmp4);
 
-            let mut tmp5 = tmp4;
-            tmp5.square();
+            let tmp5 = tmp4.square();
 
-            let mut zsquared = r.z;
-            zsquared.square();
+            let zsquared = r.z.square();
 
             r.x = tmp5;
             r.x.sub_assign(&tmp3);
             r.x.sub_assign(&tmp3);
 
             r.z.add_assign(&r.y);
-            r.z.square();
+            r.z = r.z.square();
             r.z.sub_assign(&tmp1);
             r.z.sub_assign(&zsquared);
 
@@ -234,7 +228,7 @@ impl G2Prepared {
             tmp3.mul_assign(&zsquared);
             tmp3 = tmp3.double().neg();
 
-            tmp6.square();
+            tmp6 = tmp6.square();
             tmp6.sub_assign(&tmp0);
             tmp6.sub_assign(&tmp5);
 
@@ -251,18 +245,16 @@ impl G2Prepared {
 
         fn addition_step(r: &mut G2, q: &G2Affine) -> (Fq2, Fq2, Fq2) {
             // Adaptation of Algorithm 27, https://eprint.iacr.org/2010/354.pdf
-            let mut zsquared = r.z;
-            zsquared.square();
+            let zsquared = r.z.square();
 
-            let mut ysquared = q.y;
-            ysquared.square();
+            let ysquared = q.y.square();
 
             let mut t0 = zsquared;
             t0.mul_assign(&q.x);
 
             let mut t1 = q.y;
             t1.add_assign(&r.z);
-            t1.square();
+            t1 = t1.square();
             t1.sub_assign(&ysquared);
             t1.sub_assign(&zsquared);
             t1.mul_assign(&zsquared);
@@ -270,8 +262,7 @@ impl G2Prepared {
             let mut t2 = t0;
             t2.sub_assign(&r.x);
 
-            let mut t3 = t2;
-            t3.square();
+            let t3 = t2.square();
 
             let t4 = t3.double().double();
 
@@ -288,14 +279,13 @@ impl G2Prepared {
             let mut t7 = t4;
             t7.mul_assign(&r.x);
 
-            r.x = t6;
-            r.x.square();
+            r.x = t6.square();
             r.x.sub_assign(&t5);
             r.x.sub_assign(&t7);
             r.x.sub_assign(&t7);
 
             r.z.add_assign(&t2);
-            r.z.square();
+            r.z = r.z.square();
             r.z.sub_assign(&zsquared);
             r.z.sub_assign(&t3);
 
@@ -313,11 +303,10 @@ impl G2Prepared {
             r.y = t8;
             r.y.sub_assign(&t0);
 
-            t10.square();
+            t10 = t10.square();
             t10.sub_assign(&ysquared);
 
-            let mut ztsquared = r.z;
-            ztsquared.square();
+            let ztsquared = r.z.square();
 
             t10.sub_assign(&ztsquared);
 
diff --git a/pairing/src/bls12_381/tests/mod.rs b/pairing/src/bls12_381/tests/mod.rs
index 636fcfe..2fd4deb 100644
--- a/pairing/src/bls12_381/tests/mod.rs
+++ b/pairing/src/bls12_381/tests/mod.rs
@@ -189,8 +189,7 @@ fn test_g1_uncompressed_invalid_vectors() {
         let mut x = Fq::one();
 
         loop {
-            let mut x3b = x;
-            x3b.square();
+            let mut x3b = x.square();
             x3b.mul_assign(&x);
             x3b.add_assign(&Fq::from_repr(FqRepr::from(4)).unwrap()); // TODO: perhaps expose coeff_b through API?
 
@@ -326,8 +325,7 @@ fn test_g2_uncompressed_invalid_vectors() {
         let mut x = Fq2::one();
 
         loop {
-            let mut x3b = x;
-            x3b.square();
+            let mut x3b = x.square();
             x3b.mul_assign(&x);
             x3b.add_assign(&Fq2 {
                 c0: Fq::from_repr(FqRepr::from(4)).unwrap(),
@@ -422,8 +420,7 @@ fn test_g1_compressed_invalid_vectors() {
         let mut x = Fq::one();
 
         loop {
-            let mut x3b = x;
-            x3b.square();
+            let mut x3b = x.square();
             x3b.mul_assign(&x);
             x3b.add_assign(&Fq::from_repr(FqRepr::from(4)).unwrap()); // TODO: perhaps expose coeff_b through API?
 
@@ -447,8 +444,7 @@ fn test_g1_compressed_invalid_vectors() {
         let mut x = Fq::one();
 
         loop {
-            let mut x3b = x;
-            x3b.square();
+            let mut x3b = x.square();
             x3b.mul_assign(&x);
             x3b.add_assign(&Fq::from_repr(FqRepr::from(4)).unwrap()); // TODO: perhaps expose coeff_b through API?
 
@@ -553,8 +549,7 @@ fn test_g2_compressed_invalid_vectors() {
         };
 
         loop {
-            let mut x3b = x;
-            x3b.square();
+            let mut x3b = x.square();
             x3b.mul_assign(&x);
             x3b.add_assign(&Fq2 {
                 c0: Fq::from_repr(FqRepr::from(4)).unwrap(),
@@ -585,8 +580,7 @@ fn test_g2_compressed_invalid_vectors() {
         };
 
         loop {
-            let mut x3b = x;
-            x3b.square();
+            let mut x3b = x.square();
             x3b.mul_assign(&x);
             x3b.add_assign(&Fq2 {
                 c0: Fq::from_repr(FqRepr::from(4)).unwrap(),
diff --git a/pairing/src/tests/field.rs b/pairing/src/tests/field.rs
index 304e142..eaba476 100644
--- a/pairing/src/tests/field.rs
+++ b/pairing/src/tests/field.rs
@@ -31,8 +31,7 @@ pub fn random_sqrt_tests<F: SqrtField>() {
 
     for _ in 0..10000 {
         let a = F::random(&mut rng);
-        let mut b = a;
-        b.square();
+        let b = a.square();
         assert_eq!(b.legendre(), LegendreSymbol::QuadraticResidue);
 
         let b = b.sqrt().unwrap();
@@ -43,8 +42,7 @@ pub fn random_sqrt_tests<F: SqrtField>() {
 
     let mut c = F::one();
     for _ in 0..10000 {
-        let mut b = c;
-        b.square();
+        let mut b = c.square();
         assert_eq!(b.legendre(), LegendreSymbol::QuadraticResidue);
 
         b = b.sqrt().unwrap();
@@ -218,12 +216,8 @@ fn random_doubling_tests<F: Field, R: RngCore>(rng: &mut R) {
 
 fn random_squaring_tests<F: Field, R: RngCore>(rng: &mut R) {
     for _ in 0..10000 {
-        let mut a = F::random(rng);
-        let mut b = a;
-        a.mul_assign(&b);
-        b.square();
-
-        assert_eq!(a, b);
+        let a = F::random(rng);
+        assert_eq!(a * a, a.square());
     }
 }
 
diff --git a/zcash_primitives/src/jubjub/edwards.rs b/zcash_primitives/src/jubjub/edwards.rs
index 45975cb..4ae2442 100644
--- a/zcash_primitives/src/jubjub/edwards.rs
+++ b/zcash_primitives/src/jubjub/edwards.rs
@@ -107,8 +107,7 @@ impl<E: JubjubEngine> Point<E, Unknown> {
         // as dy^2 + 1 = 0 has no solution in Fr.
 
         // tmp1 = y^2
-        let mut tmp1 = y;
-        tmp1.square();
+        let mut tmp1 = y.square();
 
         // tmp2 = (y^2 * d) + 1
         let mut tmp2 = tmp1;
@@ -335,17 +334,13 @@ impl<E: JubjubEngine, Subgroup> Point<E, Subgroup> {
         //     http://hyperelliptic.org/EFD/g1p/auto-twisted-extended.html#doubling-dbl-2008-hwcd
 
         // A = X1^2
-        let mut a = self.x;
-        a.square();
+        let a = self.x.square();
 
         // B = Y1^2
-        let mut b = self.y;
-        b.square();
+        let b = self.y.square();
 
         // C = 2*Z1^2
-        let mut c = self.z;
-        c.square();
-        c = c.double();
+        let c = self.z.square().double();
 
         // D = a*A
         //   = -A
@@ -354,7 +349,7 @@ impl<E: JubjubEngine, Subgroup> Point<E, Subgroup> {
         // E = (X1+Y1)^2 - A - B
         let mut e = self.x;
         e.add_assign(&self.y);
-        e.square();
+        e = e.square();
         e.add_assign(&d); // -A = D
         e.sub_assign(&b);
 
diff --git a/zcash_primitives/src/jubjub/fs.rs b/zcash_primitives/src/jubjub/fs.rs
index aa18590..da9bb3e 100644
--- a/zcash_primitives/src/jubjub/fs.rs
+++ b/zcash_primitives/src/jubjub/fs.rs
@@ -575,7 +575,7 @@ impl Field for Fs {
     }
 
     #[inline]
-    fn square(&mut self) {
+    fn square(&self) -> Self {
         let mut carry = 0;
         let r1 = mac_with_carry(0, (self.0).0[0], (self.0).0[1], &mut carry);
         let r2 = mac_with_carry(0, (self.0).0[0], (self.0).0[2], &mut carry);
@@ -606,7 +606,10 @@ impl Field for Fs {
         let r5 = adc(r5, 0, &mut carry);
         let r6 = mac_with_carry(r6, (self.0).0[3], (self.0).0[3], &mut carry);
         let r7 = adc(r7, 0, &mut carry);
-        self.mont_reduce(r0, r1, r2, r3, r4, r5, r6, r7);
+
+        let mut ret = *self;
+        ret.mont_reduce(r0, r1, r2, r3, r4, r5, r6, r7);
+        ret
     }
 }
 
@@ -736,8 +739,7 @@ impl SqrtField for Fs {
             0x4199cec0404d0ec0,
             0x39f6d3a994cebea,
         ]);
-        let mut a0 = a1;
-        a0.square();
+        let mut a0 = a1.square();
         a0.mul_assign(self);
 
         if a0 == NEGATIVE_ONE {
@@ -1403,16 +1405,15 @@ fn test_fs_mul_assign() {
 
 #[test]
 fn test_fr_squaring() {
-    let mut a = Fs(FsRepr([
+    let a = Fs(FsRepr([
         0xffffffffffffffff,
         0xffffffffffffffff,
         0xffffffffffffffff,
         0xe7db4ea6533afa8,
     ]));
     assert!(a.is_valid());
-    a.square();
     assert_eq!(
-        a,
+        a.square(),
         Fs::from_repr(FsRepr([
             0x12c7f55cbc52fbaa,
             0xdedc98a0b5e6ce9e,
@@ -1431,8 +1432,7 @@ fn test_fr_squaring() {
         // Ensure that (a * a) = a^2
         let a = Fs::random(&mut rng);
 
-        let mut tmp = a;
-        tmp.square();
+        let tmp = a.square();
 
         let mut tmp2 = a;
         tmp2.mul_assign(&a);
@@ -1538,8 +1538,7 @@ fn test_fs_sqrt() {
         // Ensure sqrt(a^2) = a or -a
         let a = Fs::random(&mut rng);
         let nega = a.neg();
-        let mut b = a;
-        b.square();
+        let b = a.square();
 
         let b = b.sqrt().unwrap();
 
@@ -1550,10 +1549,8 @@ fn test_fs_sqrt() {
         // Ensure sqrt(a)^2 = a for random a
         let a = Fs::random(&mut rng);
 
-        if let Some(mut tmp) = a.sqrt() {
-            tmp.square();
-
-            assert_eq!(a, tmp);
+        if let Some(tmp) = a.sqrt() {
+            assert_eq!(a, tmp.square());
         }
     }
 }
diff --git a/zcash_primitives/src/jubjub/montgomery.rs b/zcash_primitives/src/jubjub/montgomery.rs
index 5bb052f..a708d02 100644
--- a/zcash_primitives/src/jubjub/montgomery.rs
+++ b/zcash_primitives/src/jubjub/montgomery.rs
@@ -50,8 +50,7 @@ impl<E: JubjubEngine> Point<E, Unknown> {
     pub fn get_for_x(x: E::Fr, sign: bool, params: &E::Params) -> Option<Self> {
         // Given an x on the curve, y = sqrt(x^3 + A*x^2 + x)
 
-        let mut x2 = x;
-        x2.square();
+        let mut x2 = x.square();
 
         let mut rhs = x2;
         rhs.mul_assign(params.montgomery_a());
@@ -220,8 +219,7 @@ impl<E: JubjubEngine, Subgroup> Point<E, Subgroup> {
             delta.add_assign(&tmp);
         }
         {
-            let mut tmp = self.x;
-            tmp.square();
+            let mut tmp = self.x.square();
             delta.add_assign(&tmp);
             tmp = tmp.double();
             delta.add_assign(&tmp);
@@ -231,8 +229,7 @@ impl<E: JubjubEngine, Subgroup> Point<E, Subgroup> {
             delta.mul_assign(&tmp.inverse().expect("y is nonzero so this must be nonzero"));
         }
 
-        let mut x3 = delta;
-        x3.square();
+        let mut x3 = delta.square();
         x3.sub_assign(params.montgomery_a());
         x3.sub_assign(&self.x);
         x3.sub_assign(&self.x);
@@ -281,8 +278,7 @@ impl<E: JubjubEngine, Subgroup> Point<E, Subgroup> {
                         );
                     }
 
-                    let mut x3 = delta;
-                    x3.square();
+                    let mut x3 = delta.square();
                     x3.sub_assign(params.montgomery_a());
                     x3.sub_assign(&self.x);
                     x3.sub_assign(&other.x);
diff --git a/zcash_primitives/src/jubjub/tests.rs b/zcash_primitives/src/jubjub/tests.rs
index 2469d58..f0bb464 100644
--- a/zcash_primitives/src/jubjub/tests.rs
+++ b/zcash_primitives/src/jubjub/tests.rs
@@ -20,11 +20,9 @@ pub fn test_suite<E: JubjubEngine>(params: &E::Params) {
 }
 
 fn is_on_mont_curve<E: JubjubEngine, P: JubjubParams<E>>(x: E::Fr, y: E::Fr, params: &P) -> bool {
-    let mut lhs = y;
-    lhs.square();
+    let lhs = y.square();
 
-    let mut x2 = x;
-    x2.square();
+    let x2 = x.square();
 
     let mut x3 = x2;
     x3.mul_assign(&x);
@@ -42,11 +40,9 @@ fn is_on_twisted_edwards_curve<E: JubjubEngine, P: JubjubParams<E>>(
     y: E::Fr,
     params: &P,
 ) -> bool {
-    let mut x2 = x;
-    x2.square();
+    let x2 = x.square();
 
-    let mut y2 = y;
-    y2.square();
+    let y2 = y.square();
 
     // -x^2 + y^2
     let mut lhs = y2;
@@ -346,8 +342,7 @@ fn test_jubjub_params<E: JubjubEngine>(params: &E::Params) {
 
     {
         // Check that A^2 - 4 is nonsquare:
-        let mut tmp = params.montgomery_a().clone();
-        tmp.square();
+        let mut tmp = params.montgomery_a().square();
         tmp.sub_assign(&E::Fr::from_str("4").unwrap());
         assert!(tmp.legendre() == LegendreSymbol::QuadraticNonResidue);
     }
diff --git a/zcash_proofs/src/circuit/ecc.rs b/zcash_proofs/src/circuit/ecc.rs
index 9889a50..29c14a4 100644
--- a/zcash_proofs/src/circuit/ecc.rs
+++ b/zcash_proofs/src/circuit/ecc.rs
@@ -323,8 +323,7 @@ impl<E: JubjubEngine> EdwardsPoint<E> {
 
         // Compute C = d*A*A
         let c = AllocatedNum::alloc(cs.namespace(|| "C"), || {
-            let mut t0 = *a.get_value().get()?;
-            t0.square();
+            let mut t0 = a.get_value().get()?.square();
             t0.mul_assign(params.edwards_d());
 
             Ok(t0)
@@ -612,8 +611,7 @@ impl<E: JubjubEngine> MontgomeryPoint<E> {
 
         // Compute x'' = lambda^2 - A - x - x'
         let xprime = AllocatedNum::alloc(cs.namespace(|| "xprime"), || {
-            let mut t0 = *lambda.get_value().get()?;
-            t0.square();
+            let mut t0 = lambda.get_value().get()?.square();
             t0.sub_assign(params.montgomery_a());
             t0.sub_assign(self.x.get_value().get()?);
             t0.sub_assign(other.x.get_value().get()?);