diff --git a/benchmarks/bench_fields_template.nim b/benchmarks/bench_fields_template.nim
index 334e7ab..d046a33 100644
--- a/benchmarks/bench_fields_template.nim
+++ b/benchmarks/bench_fields_template.nim
@@ -124,28 +124,32 @@ proc sqrtBench*(T: typedesc, iters: int) =
     discard r.sqrt_if_square()
 
 proc sqrtP3mod4Bench*(T: typedesc, iters: int) =
+  var r: T
   let x = rng.random_unsafe(T)
-  bench("SquareRoot + isSquare (p ≡ 3 (mod 4) exponentiation)", T, iters):
-    var r = x
-    discard r.sqrt_if_square_p3mod4()
+  bench("SquareRoot (p ≡ 3 (mod 4) exponentiation)", T, iters):
+    r.invsqrt_p3mod4(x)
+    r *= x
 
 proc sqrtAddChainBench*(T: typedesc, iters: int) =
+  var r: T
   let x = rng.random_unsafe(T)
-  bench("SquareRoot + isSquare (addition chain)", T, iters):
-    var r = x
-    discard r.sqrt_if_square_addchain()
+  bench("SquareRoot (addition chain)", T, iters):
+    r.invsqrt_addchain(x)
+    r *= x
 
 proc sqrtTonelliBench*(T: typedesc, iters: int) =
+  var r: T
   let x = rng.random_unsafe(T)
-  bench("SquareRoot + isSquare (constant-time Tonelli-Shanks exponentiation)", T, iters):
-    var r = x
-    discard r.sqrt_if_square_tonelli_shanks(useAddChain = false)
+  bench("SquareRoot (constant-time Tonelli-Shanks exponentiation)", T, iters):
+    r.invsqrt_tonelli_shanks(x, useAddChain = false)
+    r *= x
 
 proc sqrtTonelliAddChainBench*(T: typedesc, iters: int) =
+  var r: T
   let x = rng.random_unsafe(T)
-  bench("SquareRoot + isSquare (constant-time Tonelli-Shanks addchain)", T, iters):
-    var r = x
-    discard r.sqrt_if_square_tonelli_shanks(useAddChain = true)
+  bench("SquareRoot (constant-time Tonelli-Shanks addchain)", T, iters):
+    r.invsqrt_tonelli_shanks(x, useAddChain = true)
+    r *= x
 
 proc powBench*(T: typedesc, iters: int) =
   let x = rng.random_unsafe(T)
diff --git a/constantine.nimble b/constantine.nimble
index 6b5d7c9..abb8b8d 100644
--- a/constantine.nimble
+++ b/constantine.nimble
@@ -263,26 +263,26 @@ proc runTests(requireGMP: bool, dumpCmdFile = false, test32bit = false, testASM
         flags &= sanitizers
       test flags, td.path, dumpCmdFile
 
-proc buildAllBenches() =
+proc buildAllBenches(useAsm = true) =
   echo "\n\n------------------------------------------------------\n"
   echo "Building benchmarks to ensure they stay relevant ..."
-  buildBench("bench_fp")
-  buildBench("bench_fp_double_precision")
-  buildBench("bench_fp2")
-  buildBench("bench_fp6")
-  buildBench("bench_fp12")
-  buildBench("bench_ec_g1")
-  buildBench("bench_ec_g2")
-  buildBench("bench_pairing_bls12_377")
-  buildBench("bench_pairing_bls12_381")
-  buildBench("bench_pairing_bn254_nogami")
-  buildBench("bench_pairing_bn254_snarks")
-  buildBench("bench_summary_bls12_377")
-  buildBench("bench_summary_bls12_381")
-  buildBench("bench_summary_bn254_nogami")
-  buildBench("bench_summary_bn254_snarks")
-  buildBench("bench_sha256")
-  buildBench("bench_hash_to_curve")
+  buildBench("bench_fp", useAsm = useAsm)
+  buildBench("bench_fp_double_precision", useAsm = useAsm)
+  buildBench("bench_fp2", useAsm = useAsm)
+  buildBench("bench_fp6", useAsm = useAsm)
+  buildBench("bench_fp12", useAsm = useAsm)
+  buildBench("bench_ec_g1", useAsm = useAsm)
+  buildBench("bench_ec_g2", useAsm = useAsm)
+  buildBench("bench_pairing_bls12_377", useAsm = useAsm)
+  buildBench("bench_pairing_bls12_381", useAsm = useAsm)
+  buildBench("bench_pairing_bn254_nogami", useAsm = useAsm)
+  buildBench("bench_pairing_bn254_snarks", useAsm = useAsm)
+  buildBench("bench_summary_bls12_377", useAsm = useAsm)
+  buildBench("bench_summary_bls12_381", useAsm = useAsm)
+  buildBench("bench_summary_bn254_nogami", useAsm = useAsm)
+  buildBench("bench_summary_bn254_snarks", useAsm = useAsm)
+  buildBench("bench_sha256", useAsm = useAsm)
+  buildBench("bench_hash_to_curve", useAsm = useAsm)
   echo "All benchmarks compile successfully."
 
 # Tasks
@@ -308,7 +308,7 @@ task test_no_assembler, "Run all tests":
 
   # Ensure benchmarks stay relevant. Ignore Windows 32-bit at the moment
   if not defined(windows) or not (existsEnv"UCPU" or getEnv"UCPU" == "i686"):
-    buildAllBenches()
+    buildAllBenches(useASM = false)
 
 task test_no_gmp, "Run tests that don't require GMP":
   # -d:testingCurves is configured in a *.nim.cfg for convenience
@@ -357,7 +357,7 @@ task test_parallel_no_assembler, "Run all tests (without macro assembler) in par
   # ignore Windows 32-bit for the moment
   # Ensure benchmarks stay relevant. Ignore Windows 32-bit at the moment
   if not defined(windows) or not (existsEnv"UCPU" or getEnv"UCPU" == "i686"):
-    buildAllBenches()
+    buildAllBenches(useASM = false)
 
 task test_parallel_no_gmp, "Run all tests in parallel (via GNU parallel)":
   # -d:testingCurves is configured in a *.nim.cfg for convenience
@@ -395,7 +395,7 @@ task test_parallel_no_gmp_no_assembler, "Run all tests in parallel (via GNU para
   # ignore Windows 32-bit for the moment
   # Ensure benchmarks stay relevant. Ignore Windows 32-bit at the moment
   if not defined(windows) or not (existsEnv"UCPU" or getEnv"UCPU" == "i686"):
-    buildAllBenches()
+    buildAllBenches(useASM = false)
 
 # Finite field 𝔽p
 # ------------------------------------------
diff --git a/constantine/arithmetic/finite_fields_square_root.nim b/constantine/arithmetic/finite_fields_square_root.nim
index 2428151..ea4cf0c 100644
--- a/constantine/arithmetic/finite_fields_square_root.nim
+++ b/constantine/arithmetic/finite_fields_square_root.nim
@@ -54,8 +54,8 @@ func hasP3mod4_primeModulus(C: static Curve): static bool =
   ## Returns true iff p ≡ 3 (mod 4)
   (BaseType(C.Mod.limbs[0]) and 3) == 3
 
-func sqrt_p3mod4(a: var Fp) =
-  ## Compute the square root of ``a``
+func invsqrt_p3mod4*(r: var Fp, a: Fp) =
+  ## Compute the inverse square root of ``a``
   ##
   ## This requires ``a`` to be a square
   ## and the prime field modulus ``p``: p ≡ 3 (mod 4)
@@ -65,14 +65,6 @@ func sqrt_p3mod4(a: var Fp) =
   ## The square root, if it exist is multivalued,
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
-  static: doAssert BaseType(Fp.C.Mod.limbs[0]) mod 4 == 3
-  a.powUnsafeExponent(Fp.getPrimePlus1div4_BE())
-
-func sqrt_invsqrt_p3mod4(sqrt, invsqrt: var Fp, a: Fp) =
-  ## If ``a`` is a square, compute the square root of ``a`` in sqrt
-  ## and the inverse square root of a in invsqrt
-  ##
-  ## This assumes that the prime field modulus ``p``: p ≡ 3 (mod 4)
   # TODO: deterministic sign
   #
   # Algorithm
@@ -83,82 +75,12 @@ func sqrt_invsqrt_p3mod4(sqrt, invsqrt: var Fp, a: Fp) =
   # a^((p-3)/2))        ≡ 1/a  (mod p)
   # a^((p-3)/4))        ≡ 1/√a (mod p)      # Requires p ≡ 3 (mod 4)
   static: doAssert BaseType(Fp.C.Mod.limbs[0]) mod 4 == 3
-
-  invsqrt = a
-  invsqrt.powUnsafeExponent(Fp.getPrimeMinus3div4_BE())
-  # √a ≡ a * 1/√a ≡ a^((p+1)/4) (mod p)
-  sqrt.prod(invsqrt, a)
-
-func sqrt_invsqrt_if_square_p3mod4(sqrt, invsqrt: var Fp, a: Fp): SecretBool =
-  ## If ``a`` is a square, compute the square root of ``a`` in sqrt
-  ## and the inverse square root of a in invsqrt
-  ##
-  ## If a is not square, sqrt and invsqrt are undefined
-  ##
-  ## This assumes that the prime field modulus ``p``: p ≡ 3 (mod 4)
-  sqrt_invsqrt_p3mod4(sqrt, invsqrt, a)
-  var test {.noInit.}: Fp
-  test.square(sqrt)
-  result = test == a
-
-func sqrt_if_square_p3mod4*(a: var Fp): SecretBool =
-  ## If ``a`` is a square, compute the square root of ``a``
-  ## if not, ``a`` is unmodified.
-  ##
-  ## This assumes that the prime field modulus ``p``: p ≡ 3 (mod 4)
-  ##
-  ## The result is undefined otherwise
-  ##
-  ## The square root, if it exist is multivalued,
-  ## i.e. both x² == (-x)²
-  ## This procedure returns a deterministic result
-  var sqrt {.noInit.}, invsqrt {.noInit.}: Fp
-  result = sqrt_invsqrt_if_square_p3mod4(sqrt, invsqrt, a)
-  a.ccopy(sqrt, result)
+  r = a
+  r.powUnsafeExponent(Fp.getPrimeMinus3div4_BE())
 
 # Specialized routines for addchain-based square roots
 # ------------------------------------------------------------
 
-func sqrt_addchain(a: var Fp) =
-  ## Compute the square root of ``a``
-  ##
-  ## This requires ``a`` to be a square
-  ## The result is undefined otherwise
-  ##
-  ## The square root, if it exist is multivalued,
-  ## i.e. both x² == (-x)²
-  ## This procedure returns a deterministic result
-  var invsqrt {.noInit.}: Fp
-  invsqrt.invsqrt_addchain(a)
-  a *= invsqrt
-
-func sqrt_invsqrt_addchain(sqrt, invsqrt: var Fp, a: Fp) =
-  ## If ``a`` is a square, compute the square root of ``a`` in sqrt
-  ## and the inverse square root of a in invsqrt
-  invsqrt.invsqrt_addchain(a)
-  sqrt.prod(invsqrt, a)
-
-func sqrt_invsqrt_if_square_addchain(sqrt, invsqrt: var Fp, a: Fp): SecretBool =
-  ## If ``a`` is a square, compute the square root of ``a`` in sqrt
-  ## and the inverse square root of a in invsqrt
-  ##
-  ## If a is not square, sqrt and invsqrt are undefined
-  sqrt_invsqrt_addchain(sqrt, invsqrt, a)
-  var test {.noInit.}: Fp
-  test.square(sqrt)
-  result = test == a
-
-func sqrt_if_square_addchain*(a: var Fp): SecretBool =
-  ## If ``a`` is a square, compute the square root of ``a``
-  ## if not, ``a`` is unmodified.
-  ##
-  ## The square root, if it exist is multivalued,
-  ## i.e. both x² == (-x)²
-  ## This procedure returns a deterministic result
-  var sqrt {.noInit.}, invsqrt {.noInit.}: Fp
-  result = sqrt_invsqrt_if_square_addchain(sqrt, invsqrt, a)
-  a.ccopy(sqrt, result)
-
 {.pop.} # inline
 
 # Tonelli Shanks for any prime
@@ -174,12 +96,15 @@ func precompute_tonelli_shanks(
     a_pre_exp.powUnsafeExponent(Fp.C.tonelliShanks(exponent))
 
 func isSquare_tonelli_shanks(
-       a, a_pre_exp: Fp): SecretBool =
+       a, a_pre_exp: Fp): SecretBool {.used.} =
   ## Returns if `a` is a quadratic residue
   ## This uses common precomputation for
   ## Tonelli-Shanks based square root and inverse square root
   ##
   ## a^((p-1-2^e)/(2*2^e))
+  ##
+  ## Note: if we need to compute a candidate square root anyway
+  ##       it's faster to square it to check if we get ``a``
   const e = Fp.C.tonelliShanks(twoAdicity)
   var r {.noInit.}: Fp
   r.square(a_pre_exp)    # a^(2(q-1-2^e)/(2*2^e)) = a^((q-1)/2^e - 1)
@@ -198,10 +123,10 @@ func isSquare_tonelli_shanks(
       r.isMinusOne()
     )
 
-func sqrt_invsqrt_tonelli_shanks_pre(
-       sqrt, invsqrt: var Fp,
+func invsqrt_tonelli_shanks_pre(
+       invsqrt: var Fp,
        a, a_pre_exp: Fp) =
-  ## Compute the square_root and inverse_square_root
+  ## Compute the inverse_square_root
   ## of `a` via constant-time Tonelli-Shanks
   ##
   ## a_pre_exp is a precomputation a^((p-1-2^e)/(2*2^e))
@@ -230,12 +155,8 @@ func sqrt_invsqrt_tonelli_shanks_pre(
     t.ccopy(buf, bNotOne)
     b = t
 
-  sqrt.prod(invsqrt, a)
-
-# ----------------------------------------------
-
-func sqrt_tonelli_shanks(a: var Fp, useAddChain: static bool) =
-  ## Compute the square root of ``a``
+func invsqrt_tonelli_shanks*(r: var Fp, a: Fp, useAddChain: static bool) =
+  ## Compute the inverse square root of ``a``
   ##
   ## This requires ``a`` to be a square
   ##
@@ -245,54 +166,9 @@ func sqrt_tonelli_shanks(a: var Fp, useAddChain: static bool) =
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
   ## This procedure is constant-time
-  var a_pre_exp{.noInit.}, sqrt{.noInit.}, invsqrt{.noInit.}: Fp
-  a_pre_exp.precompute_tonelli_shanks(a, useAddChain)
-  sqrt_invsqrt_tonelli_shanks_pre(sqrt, invsqrt, a, a_pre_exp)
-  a = sqrt
-
-func sqrt_invsqrt_tonelli_shanks(sqrt, invsqrt: var Fp, a: Fp, useAddChain: static bool) =
-  ## Compute the square root and inverse square root of ``a``
-  ##
-  ## This requires ``a`` to be a square
-  ##
-  ## The result is undefined otherwise
-  ##
-  ## The square root, if it exist is multivalued,
-  ## i.e. both x² == (-x)²
-  ## This procedure returns a deterministic result
   var a_pre_exp{.noInit.}: Fp
   a_pre_exp.precompute_tonelli_shanks(a, useAddChain)
-  sqrt_invsqrt_tonelli_shanks_pre(sqrt, invsqrt, a, a_pre_exp)
-
-func sqrt_invsqrt_if_square_tonelli_shanks(sqrt, invsqrt: var Fp, a: Fp, useAddChain: static bool): SecretBool =
-  ## Compute the square root and ivnerse square root of ``a``
-  ##
-  ## This returns true if ``a`` is square and sqrt/invsqrt contains the square root/inverse square root
-  ##
-  ## The result is undefined otherwise
-  ##
-  ## The square root, if it exist is multivalued,
-  ## i.e. both x² == (-x)²
-  ## This procedure returns a deterministic result
-  var a_pre_exp{.noInit.}: Fp
-  a_pre_exp.precompute_tonelli_shanks(a, useAddChain)
-  result = isSquare_tonelli_shanks(a, a_pre_exp)
-  sqrt_invsqrt_tonelli_shanks_pre(sqrt, invsqrt, a, a_pre_exp)
-  a = sqrt
-
-func sqrt_if_square_tonelli_shanks*(a: var Fp, useAddChain: static bool): SecretBool =
-  ## If ``a`` is a square, compute the square root of ``a``
-  ## if not, ``a`` is unmodified.
-  ##
-  ## The square root, if it exist is multivalued,
-  ## i.e. both x² == (-x)²
-  ## This procedure returns a deterministic result
-  ## This procedure is constant-time
-  var a_pre_exp{.noInit.}, sqrt{.noInit.}, invsqrt{.noInit.}: Fp
-  a_pre_exp.precompute_tonelli_shanks(a, useAddChain)
-  result = isSquare_tonelli_shanks(a, a_pre_exp)
-  sqrt_invsqrt_tonelli_shanks_pre(sqrt, invsqrt, a, a_pre_exp)
-  a = sqrt
+  invsqrt_tonelli_shanks_pre(r, a, a_pre_exp)
 
 # Public routines
 # ------------------------------------------------------------
@@ -301,6 +177,24 @@ func sqrt_if_square_tonelli_shanks*(a: var Fp, useAddChain: static bool): Secret
 
 {.push inline.}
 
+func invsqrt*[C](r: var Fp[C], a: Fp[C]) =
+  ## Compute the inverse square root of ``a``
+  ##
+  ## This requires ``a`` to be a square
+  ##
+  ## The result is undefined otherwise
+  ##
+  ## The square root, if it exist is multivalued,
+  ## i.e. both x² == (-x)²
+  ## This procedure returns a deterministic result
+  ## This procedure is constant-time
+  when C.hasSqrtAddchain():
+    r.invsqrt_addchain(a)
+  elif C.hasP3mod4_primeModulus():
+    r.invsqrt_p3mod4(a)
+  else:
+    r.invsqrt_tonelli_shanks(a, useAddChain = C.hasTonelliShanksAddchain())
+
 func sqrt*[C](a: var Fp[C]) =
   ## Compute the square root of ``a``
   ##
@@ -312,12 +206,9 @@ func sqrt*[C](a: var Fp[C]) =
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
   ## This procedure is constant-time
-  when C.hasSqrtAddchain():
-    sqrt_addchain(a)
-  elif C.hasP3mod4_primeModulus():
-    sqrt_p3mod4(a)
-  else:
-    sqrt_tonelli_shanks(a, useAddChain = C.hasTonelliShanksAddchain())
+  var t {.noInit.}: Fp[C]
+  t.invsqrt(a)
+  a *= t
 
 func sqrt_invsqrt*[C](sqrt, invsqrt: var Fp[C], a: Fp[C]) =
   ## Compute the square root and inverse square root of ``a``
@@ -329,12 +220,8 @@ func sqrt_invsqrt*[C](sqrt, invsqrt: var Fp[C], a: Fp[C]) =
   ## The square root, if it exist is multivalued,
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
-  when C.hasSqrtAddchain():
-    sqrt_invsqrt_addchain(sqrt, invsqrt, a)
-  elif C.hasP3mod4_primeModulus():
-    sqrt_invsqrt_p3mod4(sqrt, invsqrt, a)
-  else:
-    sqrt_invsqrt_tonelli_shanks(sqrt, invsqrt, a, useAddChain = C.hasTonelliShanksAddchain())
+  invsqrt.invsqrt(a)
+  sqrt.prod(invsqrt, a)
 
 func sqrt_invsqrt_if_square*[C](sqrt, invsqrt: var Fp[C], a: Fp[C]): SecretBool  =
   ## Compute the square root and ivnerse square root of ``a``
@@ -346,27 +233,33 @@ func sqrt_invsqrt_if_square*[C](sqrt, invsqrt: var Fp[C], a: Fp[C]): SecretBool
   ## The square root, if it exist is multivalued,
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
-  when C.hasSqrtAddchain():
-    result = sqrt_invsqrt_if_square_addchain(sqrt, invsqrt, a)
-  elif C.hasP3mod4_primeModulus():
-    result = sqrt_invsqrt_if_square_p3mod4(sqrt, invsqrt, a)
-  else:
-    result = sqrt_invsqrt_if_square_tonelli_shanks(sqrt, invsqrt, a, useAddChain = C.hasTonelliShanksAddchain())
+  sqrt_invsqrt(sqrt, invsqrt, a)
+  var test {.noInit.}: Fp[C]
+  test.square(sqrt)
+  result = test == a
 
 func sqrt_if_square*[C](a: var Fp[C]): SecretBool =
   ## If ``a`` is a square, compute the square root of ``a``
-  ## if not, ``a`` is unmodified.
+  ## if not, ``a`` is undefined.
   ##
   ## The square root, if it exist is multivalued,
   ## i.e. both x² == (-x)²
   ## This procedure returns a deterministic result
   ## This procedure is constant-time
-  when C.hasSqrtAddchain():
-    result = sqrt_if_square_addchain(a)
-  elif C.hasP3mod4_primeModulus():
-    result = sqrt_if_square_p3mod4(a)
-  else:
-    result = sqrt_if_square_tonelli_shanks(a, useAddChain = C.hasTonelliShanksAddchain())
+  var sqrt{.noInit.}, invsqrt{.noInit.}: Fp[C]
+  result = sqrt_invsqrt_if_square(sqrt, invsqrt, a)
+  a = sqrt
+
+func invsqrt_if_square*[C](r: var Fp[C], a: Fp[C]): SecretBool =
+  ## If ``a`` is a square, compute the inverse square root of ``a``
+  ## if not, ``a`` is undefined.
+  ##
+  ## The square root, if it exist is multivalued,
+  ## i.e. both x² == (-x)²
+  ## This procedure returns a deterministic result
+  ## This procedure is constant-time
+  var sqrt{.noInit.}: Fp[C]
+  result = sqrt_invsqrt_if_square(sqrt, r, a)
 
 {.pop.} # inline
 {.pop.} # raises no exceptions
diff --git a/constantine/arithmetic/limbs_montgomery.nim b/constantine/arithmetic/limbs_montgomery.nim
index ec70cb8..b4b5608 100644
--- a/constantine/arithmetic/limbs_montgomery.nim
+++ b/constantine/arithmetic/limbs_montgomery.nim
@@ -398,10 +398,12 @@ func montySquare*[N](r: var Limbs[N], a, M: Limbs[N],
         montSquare_CIOS_asm_adx_bmi2(r, a, M, m0ninv, spareBits >= 1)
     else:
       montSquare_CIOS_asm(r, a, M, m0ninv, spareBits >= 1)
-  else:
+  elif UseASM_X86_64:
     var r2x {.noInit.}: Limbs[2*N]
     r2x.square(a)
     r.montyRedc2x(r2x, M, m0ninv, spareBits)
+  else:
+    montyMul(r, a, a, M, m0ninv, spareBits)
 
 func redc*(r: var Limbs, a, one, M: Limbs,
            m0ninv: static BaseType, spareBits: static int) =
diff --git a/constantine/curves/bls12_381_g2_hash_to_curve.nim b/constantine/curves/bls12_381_g2_hash_to_curve.nim
index 61e4345..0ffec6c 100644
--- a/constantine/curves/bls12_381_g2_hash_to_curve.nim
+++ b/constantine/curves/bls12_381_g2_hash_to_curve.nim
@@ -31,7 +31,7 @@ const BLS12_381_h2c_G2_Z* = Fp2[BLS12_381].fromHex(  # -(2 + 𝑖)
   "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaa9",
   "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaaa"
 )
-const BLS12_381_h2c_G2_minusA* = Fp2[BLS12_381].fromHex(  # -240𝑖
+const BLS12_381_h2c_G2_minus_A* = Fp2[BLS12_381].fromHex(  # -240𝑖
   "0x0",
   "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffa9bb"
 )
@@ -60,22 +60,22 @@ const BLS12_381_h2c_G2_3_isogeny_map_xnum* = [
   # The polynomial is stored as an array of coefficients ordered from k₀ to kₙ
 
   # 1
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x5c759507e8e333ebb5b7a9a47d7ed8532c52d39fd3a042a88b58423c50ae15d5c2638e343d9c71c6238aaaaaaaa97d6",
     "0x5c759507e8e333ebb5b7a9a47d7ed8532c52d39fd3a042a88b58423c50ae15d5c2638e343d9c71c6238aaaaaaaa97d6"
   ),
   # x
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x0",
     "0x11560bf17baa99bc32126fced787c88f984f87adf7ae0c7f9a208c6b4f20a4181472aaa9cb8d555526a9ffffffffc71a"
   ),
   # x^2
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x11560bf17baa99bc32126fced787c88f984f87adf7ae0c7f9a208c6b4f20a4181472aaa9cb8d555526a9ffffffffc71e",
     "0x8ab05f8bdd54cde190937e76bc3e447cc27c3d6fbd7063fcd104635a790520c0a395554e5c6aaaa9354ffffffffe38d"
   ),
   # x^3
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x171d6541fa38ccfaed6dea691f5fb614cb14b4e7f4e810aa22d6108f142b85757098e38d0f671c7188e2aaaaaaaa5ed1",
     "0x0"
   )
@@ -85,17 +85,17 @@ const BLS12_381_h2c_G2_3_isogeny_map_xden* = [
   # The polynomial is stored as an array of coefficients ordered from k₀ to kₙ
 
   # 1
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x0",
     "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaa63"
   ),
   # x
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0xc",
     "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaa9f"
   ),
   # x^2
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x1",
     "0x0"
   )
@@ -105,22 +105,22 @@ const BLS12_381_h2c_G2_3_isogeny_map_ynum* = [
   # The polynomial is stored as an array of coefficients ordered from k₀ to kₙ
 
   # y
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x1530477c7ab4113b59a4c18b076d11930f7da5d4a07f649bf54439d87d27e500fc8c25ebf8c92f6812cfc71c71c6d706",
     "0x1530477c7ab4113b59a4c18b076d11930f7da5d4a07f649bf54439d87d27e500fc8c25ebf8c92f6812cfc71c71c6d706"
   ),
   # x*y
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x0",
     "0x5c759507e8e333ebb5b7a9a47d7ed8532c52d39fd3a042a88b58423c50ae15d5c2638e343d9c71c6238aaaaaaaa97be"
   ),
   # x^2*y
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x11560bf17baa99bc32126fced787c88f984f87adf7ae0c7f9a208c6b4f20a4181472aaa9cb8d555526a9ffffffffc71c",
     "0x8ab05f8bdd54cde190937e76bc3e447cc27c3d6fbd7063fcd104635a790520c0a395554e5c6aaaa9354ffffffffe38f"
   ),
   # x^3*y
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x124c9ad43b6cf79bfbf7043de3811ad0761b0f37a1e26286b0e977c69aa274524e79097a56dc4bd9e1b371c71c718b10",
     "0x0"
   )
@@ -130,23 +130,23 @@ const BLS12_381_h2c_G2_3_isogeny_map_yden* = [
   # The polynomial is stored as an array of coefficients ordered from k₀ to kₙ
 
   # 1
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffa8fb",
     "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffa8fb"
   ),
   # x
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x0",
     "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffa9d3"
   ),
   # x^2
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x12",
     "0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaa99"
   ),
   # x^3
-  Fp2[BLS12_381].fromHex( 
+  Fp2[BLS12_381].fromHex(
     "0x1",
     "0x0"
   )
-]
\ No newline at end of file
+]
diff --git a/constantine/ec_shortweierstrass.nim b/constantine/ec_shortweierstrass.nim
new file mode 100644
index 0000000..98f09b8
--- /dev/null
+++ b/constantine/ec_shortweierstrass.nim
@@ -0,0 +1,31 @@
+# Constantine
+# Copyright (c) 2018-2019    Status Research & Development GmbH
+# Copyright (c) 2020-Present Mamy André-Ratsimbazafy
+# Licensed and distributed under either of
+#   * MIT license (license terms in the root directory or at http://opensource.org/licenses/MIT).
+#   * Apache v2 license (license terms in the root directory or at http://www.apache.org/licenses/LICENSE-2.0).
+# at your option. This file may not be copied, modified, or distributed except according to those terms.
+
+# ############################################################
+#
+#                Short Weierstrass Elliptic Curves
+#
+# ############################################################
+
+import
+  elliptic/[
+    ec_shortweierstrass_affine,
+    ec_shortweierstrass_jacobian,
+    ec_shortweierstrass_projective
+  ]
+
+export ec_shortweierstrass_affine, ec_shortweierstrass_jacobian, ec_shortweierstrass_projective
+
+func projectiveFromJacobian*[F; Tw](
+       prj: var ECP_ShortW_Prj[F, Tw],
+       jac: ECP_ShortW_Jac[F, Tw]) {.inline.} =
+  prj.x.prod(jac.x, jac.z)
+  prj.y = jac.y
+  prj.z.square(jac.z)
+  prj.z *= jac.z
+
diff --git a/constantine/elliptic/ec_shortweierstrass_affine.nim b/constantine/elliptic/ec_shortweierstrass_affine.nim
index 8379b59..51d3df0 100644
--- a/constantine/elliptic/ec_shortweierstrass_affine.nim
+++ b/constantine/elliptic/ec_shortweierstrass_affine.nim
@@ -34,6 +34,11 @@ type
 
   SexticNonResidue* = NonResidue
 
+func `==`*(P, Q: ECP_ShortW_Aff): SecretBool =
+  ## Constant-time equality check
+  result = P.x == Q.x
+  result = result and (P.y == Q.y)
+
 func isInf*(P: ECP_ShortW_Aff): SecretBool =
   ## Returns true if P is an infinity point
   ## and false otherwise
diff --git a/constantine/elliptic/ec_shortweierstrass_jacobian.nim b/constantine/elliptic/ec_shortweierstrass_jacobian.nim
index 0febadb..0d42000 100644
--- a/constantine/elliptic/ec_shortweierstrass_jacobian.nim
+++ b/constantine/elliptic/ec_shortweierstrass_jacobian.nim
@@ -134,11 +134,13 @@ func cneg*(P: var ECP_ShortW_Jac, ctl: CTBool)  {.inline.} =
   ## Negate if ``ctl`` is true
   P.y.cneg(ctl)
 
-func sum*[F; Tw: static Twisted](
+template sumImpl[F; Tw: static Twisted](
        r: var ECP_ShortW_Jac[F, Tw],
-       P, Q: ECP_ShortW_Jac[F, Tw]
+       P, Q: ECP_ShortW_Jac[F, Tw],
+       CoefA: untyped
      ) =
   ## Elliptic curve point addition for Short Weierstrass curves in Jacobian coordinates
+  ## with the curve ``a`` being a parameter for summing on isogenous curves.
   ##
   ##   R = P + Q
   ##
@@ -148,6 +150,8 @@ func sum*[F; Tw: static Twisted](
   ##   y² = x³ + a x + b
   ##
   ## ``r`` is initialized/overwritten with the sum
+  ## ``CoefA`` allows fast path for curve with a == 0 or a == -3
+  ##            and also allows summing on curve isogenies.
   ##
   ## Implementation is constant-time, in particular it will not expose
   ## that P == Q or P == -Q or P or Q are the infinity points
@@ -218,7 +222,16 @@ func sum*[F; Tw: static Twisted](
   V_or_S *= HH_or_YY       # V = U₁*HH (add) or S = X₁*YY (dbl)
 
   block: # Compute M for doubling
-    when F.C.getCoefA() == 0:
+    # "when" static evaluation doesn't shortcut booleans :/
+    # which causes issues when CoefA isn't an int but Fp or Fp2
+    when CoefA is int:
+      const CoefA_eq_zero = CoefA == 0
+      const CoefA_eq_minus3 = CoefA == -3
+    else:
+      const CoefA_eq_zero = false
+      const CoefA_eq_minus3 = false
+
+    when CoefA_eq_zero:
       var a = H
       var b = HH_or_YY
       a.ccopy(P.x, isDbl)           # H or X₁
@@ -230,7 +243,7 @@ func sum*[F; Tw: static Twisted](
       M += HHH_or_Mpre              # 3X₁²/2
       R_or_M.ccopy(M, isDbl)
 
-    elif F.C.getCoefA() == -3:
+    elif CoefA_eq_minus3:
       var a{.noInit.}, b{.noInit.}: F
       a.sum(P.x, Z1Z1)
       b.diff(P.z, Z1Z1)
@@ -247,7 +260,7 @@ func sum*[F; Tw: static Twisted](
       # TODO: Costly `a` coefficients can be computed
       # by merging their computation with Z₃ = Z₁*Z₂*H (add) or Z₃ = Y₁*Z₁ (dbl)
       var a = H
-      var b = HH
+      var b = HH_or_YY
       a.ccopy(P.x, isDbl)
       b.ccopy(P.x, isDbl)
       HHH_or_Mpre.prod(a, b)        # HHH or X₁²
@@ -256,7 +269,8 @@ func sum*[F; Tw: static Twisted](
       a.square(HHH_or_Mpre)
       a *= HHH_or_Mpre              # a = 3X₁²
       b.square(Z1Z1)
-      b *= F.C.getCoefA()           # b = αZZ, with α the "a" coefficient of the curve
+      # b.mulCheckSparse(CoefA)     # TODO: broken static compile-time type inference
+      b *= CoefA                    # b = αZZ, with α the "a" coefficient of the curve
 
       a += b
       a.div2()
@@ -292,6 +306,52 @@ func sum*[F; Tw: static Twisted](
     r.ccopy(Q, P.isInf())
     r.ccopy(P, Q.isInf())
 
+func sum*[F; Tw: static Twisted](
+       r: var ECP_ShortW_Jac[F, Tw],
+       P, Q: ECP_ShortW_Jac[F, Tw],
+       CoefA: static F
+     ) =
+  ## Elliptic curve point addition for Short Weierstrass curves in Jacobian coordinates
+  ## with the curve ``a`` being a parameter for summing on isogenous curves.
+  ##
+  ##   R = P + Q
+  ##
+  ## Short Weierstrass curves have the following equation in Jacobian coordinates
+  ##   Y² = X³ + aXZ⁴ + bZ⁶
+  ## from the affine equation
+  ##   y² = x³ + a x + b
+  ##
+  ## ``r`` is initialized/overwritten with the sum
+  ## ``CoefA`` allows fast path for curve with a == 0 or a == -3
+  ##            and also allows summing on curve isogenies.
+  ##
+  ## Implementation is constant-time, in particular it will not expose
+  ## that P == Q or P == -Q or P or Q are the infinity points
+  ## to simple side-channel attacks (SCA)
+  ## This is done by using a "complete" or "exception-free" addition law.
+  r.sumImpl(P, Q, CoefA)
+
+func sum*[F; Tw: static Twisted](
+       r: var ECP_ShortW_Jac[F, Tw],
+       P, Q: ECP_ShortW_Jac[F, Tw]
+     ) =
+  ## Elliptic curve point addition for Short Weierstrass curves in Jacobian coordinates
+  ##
+  ##   R = P + Q
+  ##
+  ## Short Weierstrass curves have the following equation in Jacobian coordinates
+  ##   Y² = X³ + aXZ⁴ + bZ⁶
+  ## from the affine equation
+  ##   y² = x³ + a x + b
+  ##
+  ## ``r`` is initialized/overwritten with the sum
+  ##
+  ## Implementation is constant-time, in particular it will not expose
+  ## that P == Q or P == -Q or P or Q are the infinity points
+  ## to simple side-channel attacks (SCA)
+  ## This is done by using a "complete" or "exception-free" addition law.
+  r.sumImpl(P, Q, F.C.getCoefA())
+
 func madd*[F; Tw: static Twisted](
        r: var ECP_ShortW_Jac[F, Tw],
        P: ECP_ShortW_Jac[F, Tw],
diff --git a/constantine/hash_to_curve/cofactors.nim b/constantine/hash_to_curve/cofactors.nim
index eb5cce1..5f8a452 100644
--- a/constantine/hash_to_curve/cofactors.nim
+++ b/constantine/hash_to_curve/cofactors.nim
@@ -14,7 +14,8 @@ import
   ../towers,
   ../config/curves,
   ../io/io_bigints,
-  ../elliptic/[ec_shortweierstrass_projective, ec_scalar_mul]
+  ../elliptic/[ec_shortweierstrass_projective, ec_scalar_mul],
+  ../isogeny/frobenius
 
 # ############################################################
 #
@@ -96,3 +97,97 @@ func clearCofactorReference*(P: var ECP_ShortW_Prj[Fp[BW6_761], OnTwist]) {.inli
   ## Clear the cofactor of BW6_761 G2
   # Endomorphism acceleration cannot be used if cofactor is not cleared
   P.scalarMulGeneric(Cofactor_Eff_BW6_761_G2)
+
+# ############################################################
+#
+#                Clear Cofactor - Optimized
+#
+# ############################################################
+
+# BLS12 G2
+# ------------------------------------------------------------
+# From any point on the elliptic curve E2 of a BLS12 curve
+# Obtain a point in the G2 prime-order subgroup
+#
+# Described in https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#appendix-G.4
+#
+# Implementations, multiple implementations are possible in increasing order of speed:
+#
+# - The default, canonical, implementation is h_eff * P
+# - Scott et al, "Fast Hashing to G2 on Pairing-Friendly Curves", https://doi.org/10.1007/978-3-642-03298-1_8
+# - Fuentes-Castaneda et al, "Fast Hashing to G2 on Pairing-Friendly Curves", https://doi.org/10.1007/978-3-642-28496-0_25
+# - Budroni et al, "Hashing to G2 on BLS pairing-friendly curves", https://doi.org/10.1145/3313880.3313884
+# - Wahby et al "Fast and simple constant-time hashing to the BLS12-381 elliptic curve", https://eprint.iacr.org/2019/403
+# - IETF "Hashing to Elliptic Curves", https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#appendix-G.4
+#
+# In summary, the elliptic curve point multiplication is very expensive,
+# the fast methods uses endomorphism acceleration instead.
+#
+# The method described in Wahby et al is implemented by Riad Wahby
+# in C at: https://github.com/kwantam/bls12-381_hash/blob/23c1930039f58606138459557677668fabc8ce39/src/curve2/ops2.c#L106-L204
+# following Budroni et al, "Efficient hash maps to G2 on BLS curves"
+# https://eprint.iacr.org/2017/419
+#
+# "P -> [x² - x - 1] P + [x - 1] ψ(P) + ψ(ψ([2]P))"
+#
+# with Psi (ψ) - untwist-Frobenius-Twist function
+# and x the curve BLS parameter
+
+func double_repeated*[EC](P: var EC, num: int) {.inline.} =
+  ## Repeated doublings
+  for _ in 0 ..< num:
+    P.double()
+
+func pow_x(
+       r{.noalias.}: var ECP_ShortW_Prj[Fp2[BLS12_381], OnTwist],
+       P{.noalias.}: ECP_ShortW_Prj[Fp2[BLS12_381], OnTwist],
+     ) =
+  ## Does the scalar multiplication [x]P
+  ## with x the BLS12 curve parameter
+  ## For BLS12_381 [-0xd201000000010000]P
+  ## Requires r and P to not alias
+
+  # In binary
+  # 0b11
+  r.double(P)
+  r += P
+  # 0b1101
+  r.double_repeated(2)
+  r += P
+  # 0b1101001
+  r.double_repeated(3)
+  r += P
+  # 0b1101001000000001
+  r.double_repeated(9)
+  r += P
+  # 0b110100100000000100000000000000000000000000000001
+  r.double_repeated(32)
+  r += P
+  # 0b1101001000000001000000000000000000000000000000010000000000000000
+  r.double_repeated(16)
+
+  # Negative, x = -0xd201000000010000
+  r.neg(r)
+
+
+func clearCofactorFast*(P: var ECP_ShortW_Prj[Fp2[BLS12_381], OnTwist]) =
+  ## Clear the cofactor of BLS12_381 G2
+  ## Optimized using endomorphisms
+  ## P -> [x²-x-1]P + [x-1] ψ(P) + ψ²([2]P)
+
+  var xP{.noInit.}, x2P{.noInit.}: typeof(P)
+
+  xP.pow_x(P)             # 1. xP = [x]P
+  x2P.pow_x(xP)           # 2. x2P = [x²]P
+
+  x2P.diff(x2P, xP)       # 3. x2P = [x²-x]P
+  x2P.diff(x2P, P)        # 4. x2P = [x²-x-1]P
+
+  xP.diff(xP, P)          # 5. xP = [x-1]P
+  xP.frobenius_psi(xP)    # 6. xP = ψ([x-1]P) = [x-1] ψ(P)
+
+  P.double(P)             # 7. P = [2]P
+  P.frobenius_psi(P, k=2) # 8. P = ψ²([2]P)
+
+  P.sum(P, x2P)           # 9. P = [x²-x-1]P + ψ²([2]P)
+  P.sum(P, xP)            # 10. P = [x²-x-1]P + [x-1] ψ(P) + ψ²([2]P)
diff --git a/constantine/hash_to_curve/h2c_hash_to_field.nim b/constantine/hash_to_curve/h2c_hash_to_field.nim
index 03ef25c..aa53539 100644
--- a/constantine/hash_to_curve/h2c_hash_to_field.nim
+++ b/constantine/hash_to_curve/h2c_hash_to_field.nim
@@ -30,12 +30,11 @@ func ceilDiv(a, b: uint): uint =
   ## ceil(a / b)
   (a + b - 1) div b
 
-proc copyFrom[N](output: var openarray[byte], bi: array[N, byte], cur: var uint) =
-  var b_index = 0'u
-  while b_index < bi.len.uint and cur < output.len.uint:
-    output[cur] = bi[b_index]
-    inc cur
-    inc b_index
+proc copyFrom[M, N: static int](output: var array[M, byte], bi: array[N, byte], cur: var uint) =
+  static: doAssert M mod N == 0
+  for i in 0'u ..< N:
+    output[cur+i] = bi[i]
+  cur += N.uint
 
 template strxor(b_i: var array, b0: array): untyped =
   for i in 0 ..< b_i.len:
@@ -57,9 +56,9 @@ func shortDomainSepTag[DigestSize: static int, B: byte|char](
   ctx.update oversizedDST
   ctx.finish(output)
 
-func expandMessageXMD*[B1, B2, B3: byte|char](
+func expandMessageXMD*[B1, B2, B3: byte|char, len_in_bytes: static int](
        H: type CryptoHash,
-       output: var openarray[byte],
+       output: var array[len_in_bytes, byte],
        augmentation: openarray[B1],
        message: openarray[B2],
        domainSepTag: openarray[B3]
@@ -110,23 +109,23 @@ func expandMessageXMD*[B1, B2, B3: byte|char](
   const DigestSize = Hash.digestSize()
   const BlockSize = Hash.internalBlockSize()
 
-  assert output.len mod 8 == 0
+  static:
+    doAssert output.len mod 8 == 0  # By spec
+    doAssert output.len mod 32 == 0 # Assumed by copy optimization
 
   let ell = ceilDiv(output.len.uint, DigestSize.uint)
   const zPad = default(array[BlockSize, byte])
-  let l_i_b_str = output.len.uint16.toBytesBE()
+  var l_i_b_str0 {.noInit.}: array[3, byte]
+  l_i_b_str0.asBytesBE(output.len.uint16, pos = 0)
+  l_i_b_str0[2] = 0
 
   var b0 {.noinit, align: DigestSize.}: array[DigestSize, byte]
-  func ctZpad(): Hash =
-    # Compile-time precompute
-    # TODO upstream: `toOpenArray` throws "cannot generate code for: mSlice"
-    result.init()
-    result.update zPad
-  var ctx = ctZpad() # static(ctZpad())
+  var ctx {.noInit.}: Hash
+  ctx.initZeroPadded()
   ctx.update augmentation
   ctx.update message
-  ctx.update l_i_b_str
-  ctx.update [byte 0]
+  ctx.update l_i_b_str0
+  # ctx.update [byte 0] # already appended to l_i_b_str
   ctx.update domainSepTag
   ctx.update [byte domainSepTag.len] # DST_prime
   ctx.finish(b0)
diff --git a/constantine/hash_to_curve/h2c_map_to_isocurve_swu.nim b/constantine/hash_to_curve/h2c_map_to_isocurve_swu.nim
index 159c166..dedb940 100644
--- a/constantine/hash_to_curve/h2c_map_to_isocurve_swu.nim
+++ b/constantine/hash_to_curve/h2c_map_to_isocurve_swu.nim
@@ -109,8 +109,7 @@ func invsqrt_if_square[C: static Curve](
     t2 *= NonResidue
     t1 -= t2
 
-  # TODO: implement invsqrt alone
-  result = sqrt_invsqrt_if_square(sqrt = r.c1, invsqrt = t3, t1) # 1/sqrt(a0² - β a1²)
+  result = t3.invsqrt_if_square(t1)    # 1/sqrt(a0² - β a1²)
 
   # If input is not a square in Fp2, multiply by 1/Z³
   inp.prod(a, h2cConst(C, G2, inv_Z3)) # inp = a / Z³
@@ -131,8 +130,7 @@ func invsqrt_if_square[C: static Curve](
   t1.ccopy(t2, t1.isZero())
   t1.div2()    # (a0 ± sqrt(a0² - βa1²))/2
 
-  # TODO: implement invsqrt alone
-  sqrt_invsqrt(sqrt = r.c1, invsqrt = r.c0, t1)
+  r.c0.invsqrt(t1)
 
   r.c1 = inp.c1
   r.c1.div2()
@@ -152,19 +150,22 @@ func invsqrt_if_square[C: static Curve](
 
 func mapToIsoCurve_sswuG2_opt9mod16*[C: static Curve](
        xn, xd, yn: var Fp2[C],
-       u: Fp2[C]) =
+       u: Fp2[C], xd3: var Fp2[C]) =
   ## Given G2, the target prime order subgroup of E2 we want to hash to,
   ## this function maps any field element of Fp2 to E'2
   ## a curve isogenous to E2 using the Simplified Shallue-van de Woestijne method.
   ##
   ## This requires p² ≡ 9 (mod 16).
-  #
-  # Input:
-  # - u, an Fp2 element
-  # Output:
-  # - (xn, xd, yn, yd) such that (x', y') = (xn/xd, yn/yd)
-  #   is a point of E'2
-  # - yd is implied to be 1
+  ##
+  ## Input:
+  ## - u, an Fp2 element
+  ## Output:
+  ## - (xn, xd, yn, yd) such that (x', y') = (xn/xd, yn/yd)
+  ##   is a point of E'2
+  ## - yd is implied to be 1
+  ## Scratchspace:
+  ## - xd3 is temporary scratchspace that will hold xd³
+  ##   after execution (which might be useful for Jacobian coordinate conversion)
   #
   # Paper: https://eprint.iacr.org/2019/403
   # Spec: https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#appendix-G.2.3
@@ -175,7 +176,6 @@ func mapToIsoCurve_sswuG2_opt9mod16*[C: static Curve](
   var
     uu {.noInit.}, tv2 {.noInit.}: Fp2[C]
     tv4 {.noInit.}, x2n {.noInit.}, gx1 {.noInit.}: Fp2[C]
-    gxd {.noInit.}: Fp2[C]
     y2 {.noInit.}: Fp2[C]
     e1, e2: SecretBool
 
@@ -184,6 +184,7 @@ func mapToIsoCurve_sswuG2_opt9mod16*[C: static Curve](
   template x1n: untyped = xn
   template y1: untyped = yn
   template Zuu: untyped = x2n
+  template gxd: untyped = xd3
 
   # x numerators
   uu.square(u)                      # uu = u²
@@ -203,7 +204,7 @@ func mapToIsoCurve_sswuG2_opt9mod16*[C: static Curve](
   # y numerators
   tv2.square(xd)
   gxd.prod(xd, tv2)                         # gxd = xd³
-  tv2 *= h2CConst(C, G2, Aprime_E2)
+  tv2.mulCheckSparse(h2CConst(C, G2, Aprime_E2))
   gx1.square(x1n)
   gx1 += tv2                                # x1n² + A * xd²
   gx1 *= x1n                                # x1n³ + A * x1n * xd²
diff --git a/constantine/hash_to_curve/hash_to_curve.nim b/constantine/hash_to_curve/hash_to_curve.nim
index 76a967b..0fd6c88 100644
--- a/constantine/hash_to_curve/hash_to_curve.nim
+++ b/constantine/hash_to_curve/hash_to_curve.nim
@@ -11,7 +11,7 @@ import
   ../config/[common, curves],
   ../primitives, ../arithmetic, ../towers,
   ../curves/zoo_hash_to_curve,
-  ../elliptic/ec_shortweierstrass_projective,
+  ../ec_shortweierstrass,
   ./h2c_hash_to_field,
   ./h2c_map_to_isocurve_swu,
   ./cofactors,
@@ -34,8 +34,28 @@ import
 # Map to curve
 # ----------------------------------------------------------------
 
+func mapToIsoCurve_sswuG2_opt9mod16[F; Tw: static Twisted](
+       r: var ECP_ShortW_Jac[F, Tw],
+       u: F) =
+  var
+    xn{.noInit.}, xd{.noInit.}: F
+    yn{.noInit.}: F
+    xd3{.noInit.}: F
+
+  mapToIsoCurve_sswuG2_opt9mod16(
+    xn, xd,
+    yn,
+    u, xd3
+  )
+
+  # Convert to Jacobian
+  r.z = xd          # Z = xd
+  r.x.prod(xn, xd)  # X = xZ² = xn/xd * xd² = xn*xd
+  r.y.prod(yn, xd3) # Y = yZ³ = yn * xd³
+
 func mapToCurve[F; Tw: static Twisted](
-       r: var ECP_ShortW_Prj[F, Tw], u: F) =
+       r: var (ECP_ShortW_Prj[F, Tw] or ECP_ShortW_Jac[F, Tw]),
+       u: F) =
   ## Map an element of the
   ## finite or extension field F
   ## to an elliptic curve E
@@ -48,11 +68,12 @@ func mapToCurve[F; Tw: static Twisted](
     var
       xn{.noInit.}, xd{.noInit.}: F
       yn{.noInit.}: F
+      xd3{.noInit.}: F
 
     mapToIsoCurve_sswuG2_opt9mod16(
       xn, xd,
       yn,
-      u
+      u, xd3
     )
 
     # 2. Map from E'2 to E2
@@ -64,6 +85,39 @@ func mapToCurve[F; Tw: static Twisted](
   else:
     {.error: "Not implemented".}
 
+func mapToCurve_fusedAdd[F; Tw: static Twisted](
+       r: var ECP_ShortW_Jac[F, Tw],
+       u0, u1: F) =
+  ## Map 2 elements of the
+  ## finite or extension field F
+  ## to an elliptic curve E
+  ## and add them
+  # Optimization suggested in https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-6.6.3
+  #   Note that iso_map is a group homomorphism, meaning that point
+  #   addition commutes with iso_map.  Thus, when using this mapping in the
+  #   hash_to_curve construction of Section 3, one can effect a small
+  #   optimization by first mapping u0 and u1 to E', adding the resulting
+  #   points on E', and then applying iso_map to the sum.  This gives the
+  #   same result while requiring only one evaluation of iso_map.
+
+  # Jacobian formulae are independent of the curve equation B'
+  # y² = x³ + A'*x + B'
+  # unlike the complete projective formulae which heavily depends on it
+  # So we use jacobian coordinates for computation on isogenies.
+
+  var P0{.noInit.}, P1{.noInit.}: ECP_ShortW_Jac[F, Tw]
+  when F.C == BLS12_381 and F is Fp2:
+    # 1. Map to E'2 isogenous to E2
+    P0.mapToIsoCurve_sswuG2_opt9mod16(u0)
+    P1.mapToIsoCurve_sswuG2_opt9mod16(u1)
+
+    P0.sum(P0, P1, h2CConst(F.C, G2, Aprime_E2))
+
+    # 2. Map from E'2 to E2
+    r.h2c_isogeny_map(P0, isodegree = 3)
+  else:
+    {.error: "Not implemented".}
+
 # Hash to curve
 # ----------------------------------------------------------------
 
@@ -102,9 +156,14 @@ func hashToCurve*[
   var u{.noInit.}: array[2, F]
   H.hashToField(k, u, augmentation, message, domainSepTag)
 
-  var P{.noInit.}: array[2, ECP_ShortW_Prj[F, Tw]]
-  P[0].mapToCurve(u[0])
-  P[1].mapToCurve(u[1])
+  when false:
+    var P{.noInit.}: array[2, ECP_ShortW_Prj[F, Tw]]
+    P[0].mapToCurve(u[0])
+    P[1].mapToCurve(u[1])
+    output.sum(P[0], P[1])
+  else:
+    var Pjac{.noInit.}: ECP_ShortW_Jac[F, Tw]
+    Pjac.mapToCurve_fusedAdd(u[0], u[1])
+    output.projectiveFromJacobian(Pjac)
 
-  output.sum(P[0], P[1])
-  output.clearCofactorReference() # TODO - fast cofactor clear
+  output.clearCofactorFast()
diff --git a/constantine/hashes/h_sha256.nim b/constantine/hashes/h_sha256.nim
index 63c6250..873a4fd 100644
--- a/constantine/hashes/h_sha256.nim
+++ b/constantine/hashes/h_sha256.nim
@@ -315,14 +315,37 @@ func init*(ctx: var Sha256Context) =
   ctx.buf.setZero()
   ctx.bufIdx = 0
 
-  ctx.H[0] = 0x6a09e667'u32;
-  ctx.H[1] = 0xbb67ae85'u32;
-  ctx.H[2] = 0x3c6ef372'u32;
-  ctx.H[3] = 0xa54ff53a'u32;
-  ctx.H[4] = 0x510e527f'u32;
-  ctx.H[5] = 0x9b05688c'u32;
-  ctx.H[6] = 0x1f83d9ab'u32;
-  ctx.H[7] = 0x5be0cd19'u32;
+  ctx.H[0] = 0x6a09e667'u32
+  ctx.H[1] = 0xbb67ae85'u32
+  ctx.H[2] = 0x3c6ef372'u32
+  ctx.H[3] = 0xa54ff53a'u32
+  ctx.H[4] = 0x510e527f'u32
+  ctx.H[5] = 0x9b05688c'u32
+  ctx.H[6] = 0x1f83d9ab'u32
+  ctx.H[7] = 0x5be0cd19'u32
+
+func initZeroPadded*(ctx: var Sha256Context) =
+  ## Initialize a Sha256 context
+  ## with the result of
+  ## ctx.init()
+  ## ctx.update default(array[BlockSize, byte])
+  #
+  # This work arounds `toOpenArray`
+  # not working in the Nim VM, preventing `sha256.update`
+  # at compile-time
+
+  ctx.msgLen = 64
+  ctx.buf.setZero()
+  ctx.bufIdx = 0
+
+  ctx.H[0] = 0xda5698be'u32
+  ctx.H[1] = 0x17b9b469'u32
+  ctx.H[2] = 0x62335799'u32
+  ctx.H[3] = 0x779fbeca'u32
+  ctx.H[4] = 0x8ce5d491'u32
+  ctx.H[5] = 0xc0d26243'u32
+  ctx.H[6] = 0xbafef9ea'u32
+  ctx.H[7] = 0x1837a9d8'u32
 
 func update*[T: char|byte](ctx: var Sha256Context, message: openarray[T]) =
   ## Append a message to a SHA256 context
diff --git a/constantine/io/endians.nim b/constantine/io/endians.nim
index 44003cc..beedd40 100644
--- a/constantine/io/endians.nim
+++ b/constantine/io/endians.nim
@@ -77,8 +77,13 @@ func dumpRawInt*[T: byte|char](
     for i in 0'u ..< L:
       dst[cursor+i] = toByte(src shr ((L-i-1) * 8))
 
-func toBytesBE*(num: SomeUnsignedInt): array[sizeof(num), byte] {.inline.}=
-  ## Convert an integer to an array of bytes
+func asBytesBE*[N: static int](
+       r: var array[N, byte],
+       num: SomeUnsignedInt,
+       pos: static int) {.inline.}=
+  ## Store an integer into an array of bytes
+  ## in big endian representation
   const L = sizeof(num)
+  static: doAssert N - (pos + L) >= 0
   for i in 0 ..< L:
-    result[i] = toByte(num shr ((L-1-i) * 8))
+    r[i] = toByte(num shr ((L-1-i) * 8))
diff --git a/constantine/isogeny/h2c_isogeny_maps.nim b/constantine/isogeny/h2c_isogeny_maps.nim
index ec5b3b7..4e97bea 100644
--- a/constantine/isogeny/h2c_isogeny_maps.nim
+++ b/constantine/isogeny/h2c_isogeny_maps.nim
@@ -10,7 +10,10 @@ import
   # Internals
   ../primitives, ../arithmetic, ../towers,
   ../curves/zoo_hash_to_curve,
-  ../elliptic/ec_shortweierstrass_projective
+  ../elliptic/[
+    ec_shortweierstrass_projective,
+    ec_shortweierstrass_jacobian,
+  ]
 
 # ############################################################
 #
@@ -78,6 +81,51 @@ func poly_eval_horner_scaled[F; D, N: static int](
     # Missing scaling factor
     r *= xd_pow[isodegree - poly_degree - 1]
 
+func h2c_isogeny_map[F](
+       rxn, rxd, ryn, ryd: var F,
+       xn, xd, yn: F, isodegree: static int) =
+  ## Given G2, the target prime order subgroup of E2,
+  ## this function maps an element of
+  ## E'2 a curve isogenous to E2
+  ## to E2.
+  ##
+  ## The E'2 input is represented as
+  ## (x', y') with x' = xn/xd and y' = yn/yd
+  ##
+  ## yd is assumed to be 1 hence y' == yn
+  ##
+  ## The E2 output is represented as
+  ## (rx, ry) with rx = rxn/rxd and ry = ryn/ryd
+
+  # xd^e with e in [1, N], for example [xd, xd², xd³]
+  static: doAssert isodegree >= 2
+  var xd_pow{.noInit.}: array[isodegree, F]
+  xd_pow[0] = xd
+  xd_pow[1].square(xd_pow[0])
+  for i in 2 ..< xd_pow.len:
+    xd_pow[i].prod(xd_pow[i-1], xd_pow[0])
+
+  rxn.poly_eval_horner_scaled(
+    xn, xd_pow,
+    h2cIsomapPoly(F.C, G2, isodegree, xnum)
+  )
+  rxd.poly_eval_horner_scaled(
+    xn, xd_pow,
+    h2cIsomapPoly(F.C, G2, isodegree, xden)
+  )
+
+  ryn.poly_eval_horner_scaled(
+    xn, xd_pow,
+    h2cIsomapPoly(F.C, G2, isodegree, ynum)
+  )
+  ryd.poly_eval_horner_scaled(
+    xn, xd_pow,
+    h2cIsomapPoly(F.C, G2, isodegree, yden)
+  )
+
+  # y coordinate is y' * poly_yNum(x)
+  ryn *= yn
+
 func h2c_isogeny_map*[F; Tw: static Twisted](
        r: var ECP_ShortW_Prj[F, Tw],
        xn, xd, yn: F, isodegree: static int) =
@@ -95,33 +143,15 @@ func h2c_isogeny_map*[F; Tw: static Twisted](
   ## - https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-6.6.3
   ## - https://github.com/cfrg/draft-irtf-cfrg-hash-to-curve
 
-  # xd^e with e in [1, N], for example [xd, xd², xd³]
-  var xd_pow: array[isodegree, F]
-  xd_pow[0] = xd
-  for i in 1 ..< xd_pow.len:
-    xd_pow[i].prod(xd_pow[i-1], xd)
+  var t{.noInit.}: F
 
-  r.x.poly_eval_horner_scaled(
-    xn, xd_pow,
-    h2cIsomapPoly(F.C, G2, isodegree, xnum)
+  h2c_isogeny_map(
+    rxn = r.x,
+    rxd = r.z,
+    ryn = r.y,
+    ryd = t,
+    xn, xd, yn, isodegree
   )
-  r.z.poly_eval_horner_scaled(
-    xn, xd_pow,
-    h2cIsomapPoly(F.C, G2, isodegree, xden)
-  )
-
-  r.y.poly_eval_horner_scaled(
-    xn, xd_pow,
-    h2cIsomapPoly(F.C, G2, isodegree, ynum)
-  )
-  var t {.noInit.}: F
-  t.poly_eval_horner_scaled(
-    xn, xd_pow,
-    h2cIsomapPoly(F.C, G2, isodegree, yden)
-  )
-
-  # y coordinate is y' * poly_yNum(x)
-  r.y *= yn
 
   # Now convert to projective coordinates
   # (x, y) => (xnum/xden, ynum/yden)
@@ -129,3 +159,103 @@ func h2c_isogeny_map*[F; Tw: static Twisted](
   r.y *= r.z
   r.x *= t
   r.z *= t
+
+func h2c_isogeny_map*[F; Tw: static Twisted](
+       r: var ECP_ShortW_Jac[F, Tw],
+       xn, xd, yn: F, isodegree: static int) =
+  ## Given G2, the target prime order subgroup of E2,
+  ## this function maps an element of
+  ## E'2 a curve isogenous to E2
+  ## to E2.
+  ##
+  ## The E'2 input is represented as
+  ## (x', y') with x' = xn/xd and y' = yn/yd
+  ##
+  ## yd is assumed to be 1 hence y' == yn
+  ##
+  ## Reference:
+  ## - https://datatracker.ietf.org/doc/html/draft-irtf-cfrg-hash-to-curve-11#section-6.6.3
+  ## - https://github.com/cfrg/draft-irtf-cfrg-hash-to-curve
+
+  var rxn{.noInit.}, rxd{.noInit.}: F
+  var ryn{.noInit.}, ryd{.noInit.}: F
+
+  h2c_isogeny_map(
+    rxn, rxd,
+    ryn, ryd,
+    xn, xd, yn, isodegree
+  )
+
+  # Now convert to jacobian coordinates
+  # (x, y) => (xnum/xden, ynum/yden)
+  #       <=> (xZ², yZ³, Z)
+  #       <=> (xnum*xden*yden², ynum*yden²*xden³, xden*yden)
+  r.z.prod(rxd, ryd) # Z = xn * xd
+  r.x.prod(rxn, ryd) # X = xn * yd
+  r.x *= r.z         # X = xn * xd * yd²
+  r.y.square(r.z)    # Y = xd² * yd²
+  r.y *= rdx         # Y = yd² * xd³
+  r.y *= ryn         # Y = yn * yd² * xd³
+
+func h2c_isogeny_map*[F; Tw: static Twisted](
+       r: var ECP_ShortW_Jac[F, Tw],
+       P: ECP_ShortW_Jac[F, Tw],
+       isodegree: static int) =
+  ## Map P in isogenous curve E'2
+  ## to r in E2
+  ##
+  ## r and P are NOT on the same curve.
+  #
+  # We have in affine <=> jacobian representation
+  # (x, y) <=> (xn/xd, yn/yd)
+  #        <=> (xZ², yZ³, Z)
+  #
+  # We scale the isogeny coefficients by powers
+  # of Z²
+
+  var xn{.noInit.}, xd{.noInit.}: F
+  var yn{.noInit.}, yd{.noInit.}: F
+
+  # Z²^e with e in [1, N], for example [Z², Z⁴, Z⁶]
+  static: doAssert isodegree >= 2
+  var ZZpow{.noInit.}: array[isodegree, F]
+  ZZpow[0].square(P.z)
+  ZZpow[1].square(ZZpow[0])
+  for i in 2 ..< ZZpow.len:
+    ZZpow[i].prod(ZZpow[i-1], ZZpow[0])
+
+  xn.poly_eval_horner_scaled(
+    P.x, ZZpow,
+    h2cIsomapPoly(F.C, G2, isodegree, xnum)
+  )
+  xd.poly_eval_horner_scaled(
+    P.x, ZZpow,
+    h2cIsomapPoly(F.C, G2, isodegree, xden)
+  )
+
+  yn.poly_eval_horner_scaled(
+    P.x, ZZpow,
+    h2cIsomapPoly(F.C, G2, isodegree, ynum)
+  )
+  yd.poly_eval_horner_scaled(
+    P.x, ZZpow,
+    h2cIsomapPoly(F.C, G2, isodegree, yden)
+  )
+
+  # yn = y' * poly_yNum(x) = yZ³ * poly_yNum(x)
+  yn *= P.y
+
+  # Scale yd by Z³
+  yd *= P.z
+  yd *= ZZpow[0]
+
+  # Now convert to jacobian coordinates
+  # (x, y) => (xnum/xden, ynum/yden)
+  #       <=> (xZ², yZ³, Z)
+  #       <=> (xnum*xden*yden², ynum*yden²*xden³, xden*yden)
+  r.z.prod(xd, yd) # Z = xn * xd
+  r.x.prod(xn, yd) # X = xn * yd
+  r.x *= r.z       # X = xn * xd * yd²
+  r.y.square(r.z)  # Y = xd² * yd²
+  r.y *= xd        # Y = yd² * xd³
+  r.y *= yn        # Y = yn * yd² * xd³
diff --git a/constantine/tower_field_extensions/square_root_fp2.nim b/constantine/tower_field_extensions/square_root_fp2.nim
index 4bb1767..d67c35d 100644
--- a/constantine/tower_field_extensions/square_root_fp2.nim
+++ b/constantine/tower_field_extensions/square_root_fp2.nim
@@ -133,10 +133,8 @@ func sqrt_if_square_opt(a: var Fp2): SecretBool =
   t1.ccopy(t2, t1.isZero())
   t1.div2()                           # (a0 ± sqrt(a0² - β a1²))/2
 
-  # TODO: implement invsqrt alone
   # t1 being an actual sqrt will be tested in sqrt_rotate_extension
-  # 1/sqrt((a0 ± sqrt(a0² - β b²))/2)
-  sqrt_invsqrt(sqrt = cand.c1, invsqrt = cand.c0, t1) # we only want invsqrt = cand.c0
+  cand.c0.invsqrt(t1)                 # 1/sqrt((a0 ± sqrt(a0² - β b²))/2)
 
   cand.c1 = a.c1
   cand.c1.div2()
diff --git a/sage/derive_hash_to_curve.sage b/sage/derive_hash_to_curve.sage
index b3589fd..1a11a07 100644
--- a/sage/derive_hash_to_curve.sage
+++ b/sage/derive_hash_to_curve.sage
@@ -181,7 +181,7 @@ def genBLS12381G2_H2C_constants(curve_config):
   buf += field_to_nim(Z, 'Fp2', curve_name, comment_right = "-(2 + 𝑖)")
   buf += '\n'
 
-  buf += f'const {curve_name}_h2c_G2_minusA* = '
+  buf += f'const {curve_name}_h2c_G2_minus_A* = '
   buf += field_to_nim(minus_A, 'Fp2', curve_name, comment_right = "-240𝑖")
   buf += '\n'
 
diff --git a/tests/t_finite_fields_sqrt.nim b/tests/t_finite_fields_sqrt.nim
index 8664df5..d536125 100644
--- a/tests/t_finite_fields_sqrt.nim
+++ b/tests/t_finite_fields_sqrt.nim
@@ -80,7 +80,6 @@ proc exhaustiveCheck(C: static Curve, modulus: static int) =
         check:
           bool not a.isSquare()
           bool not a.sqrt_if_square()
-          bool (a == a2) # a shouldn't be modified
 
 template testImpl(a: untyped): untyped {.dirty.} =
   var na{.noInit.}: typeof(a)
diff --git a/tests/t_fp_tower_frobenius_template.nim b/tests/t_fp_tower_frobenius_template.nim
index 5e9f4ba..4e4e986 100644
--- a/tests/t_fp_tower_frobenius_template.nim
+++ b/tests/t_fp_tower_frobenius_template.nim
@@ -61,7 +61,7 @@ proc runFrobeniusTowerTests*[N](
     ) =
   # Random seed for reproducibility
   var rng: RngState
-  let seed = 0 # uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
+  let seed = uint32(getTime().toUnix() and (1'i64 shl 32 - 1)) # unixTime mod 2^32
   rng.seed(seed)
   echo moduleName, " xoshiro512** seed: ", seed