diff --git a/constantine/math/bigints_checked.nim b/constantine/math/bigints_checked.nim
index a7f92a7..bb434c8 100644
--- a/constantine/math/bigints_checked.nim
+++ b/constantine/math/bigints_checked.nim
@@ -158,7 +158,12 @@ func montyMul*[mBits](r: var BigInt[mBits], a, b, M: BigInt[mBits], negInvModWor
   ## to avoid duplicating with Nim zero-init policy
   montyMul(r.view, a.view, b.view, M.view, Word(negInvModWord))
 
-import stew/byteutils
+func montySquare*[mBits](r: var BigInt[mBits], a, M: BigInt[mBits], negInvModWord: static BaseType) =
+  ## Compute r <- a^2 (mod M) in the Montgomery domain
+  ##
+  ## This resets r to zero before processing. Use {.noInit.}
+  ## to avoid duplicating with Nim zero-init policy
+  montySquare(r.view, a.view, M.view, Word(negInvModWord))
 
 func montyPow*[mBits, eBits: static int](
        a: var BigInt[mBits], exponent: BigInt[eBits],
diff --git a/constantine/math/bigints_raw.nim b/constantine/math/bigints_raw.nim
index 76745e9..4672bf6 100644
--- a/constantine/math/bigints_raw.nim
+++ b/constantine/math/bigints_raw.nim
@@ -550,6 +550,17 @@ func montyResidue*(
 
   montyMul(r, a, r2ModN, N, negInvModWord)
 
+func montySquare(
+       r: BigIntViewMut, a: BigIntViewAny,
+       M: BigIntViewConst, negInvModWord: Word) {.inline.} =
+  ## Compute r <- a^2 (mod M) in the Montgomery domain
+  ## `negInvModWord` = -1/M (mod Word). Our words are 2^31 or 2^63
+
+  # TODO: specialized implementation when optimizing for speed
+  #       and montyMul when optimizing for size
+  # - https://www.intel.com/content/dam/www/public/us/en/documents/white-papers/ia-large-integer-arithmetic-paper.pdf
+  montyMul(r, a, a, M, negInvModWord)
+
 # Montgomery Modular Exponentiation
 # ------------------------------------------
 # We use fixed-window based exponentiation
@@ -664,7 +675,7 @@ func montyPowSquarings(
 
   # We have k bits and can do k squaring
   for i in 0 ..< k:
-    tmp.montyMul(a, a, M, negInvModWord)
+    tmp.montySquare(a, M, negInvModWord)
     copyMem(pointer a, pointer tmp, bigIntSize)
 
   return (k, bits)
diff --git a/constantine/math/finite_fields.nim b/constantine/math/finite_fields.nim
index 9fcf273..2ca17de 100644
--- a/constantine/math/finite_fields.nim
+++ b/constantine/math/finite_fields.nim
@@ -129,6 +129,15 @@ func `*`*(a, b: Fq): Fq {.noInit.} =
   result.mres.setInternalBitLength()
   result.mres.montyMul(a.mres, b.mres, Fq.C.Mod.mres, Fq.C.getNegInvModWord())
 
+func square*(a: Fq): Fq {.noInit.} =
+  ## Squaring over Fq
+  ##
+  ## It is recommended to assign with {.noInit.}
+  ## as Fq elements are usually large and this
+  ## routine will zero init internally the result.
+  result.mres.setInternalBitLength()
+  result.mres.montySquare(a.mres, Fq.C.Mod.mres, Fq.C.getNegInvModWord())
+
 func pow*(a: var Fq, exponent: BigInt) =
   ## Exponentiation over Fq
   ## ``a``: a field element to be exponentiated