diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..d915edc
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,21 @@
+name: CI
+
+on: [push, pull_request]
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        nim: [2.2.4]
+    steps:
+    - uses: actions/checkout@v2
+    - uses: iffy/install-nim@v5
+      with:
+        version: ${{ matrix.nim }}
+    - name: Install dependencies
+      run: nimble install -d -y
+    - name: Test
+      run: nimble test -y
+    - name: Build binary
+      run: nimble build -y
diff --git a/cli/cli_main.nim b/cli/cli_main.nim
index b6a5a0e..b126503 100644
--- a/cli/cli_main.nim
+++ b/cli/cli_main.nim
@@ -9,6 +9,8 @@ import std/times
 import std/options
 # import strformat
 
+import taskpools
+
 import groth16/prover
 import groth16/verifier
 import groth16/files/witness
@@ -203,11 +205,12 @@ proc cliMain(cfg: Config) =
     else:
       echo("generating proof...")
       let print_timings = cfg.measure_time and cfg.verbose
+      var pool = Taskpool.new(cfg.nthreads)
       withMeasureTime(cfg.measure_time,"proving"):
         if cfg.no_masking:
-          proof = generateProofWithTrivialMask(cfg.nthreads, print_timings, zkey, wtns)
+          proof = generateProofWithTrivialMask(zkey, wtns, pool, print_timings)
         else:
-          proof = generateProof(cfg.nthreads, print_timings, zkey, wtns)
+          proof = generateProof(zkey, wtns, pool, print_timings)
     
       if not (cfg.output_file == ""):
         echo("exporting the proof to " & quoted(cfg.output_file))
diff --git a/groth16.nimble b/groth16.nimble
index a6566ea..168d302 100644
--- a/groth16.nimble
+++ b/groth16.nimble
@@ -8,4 +8,4 @@ binDir      = "build"
 namedBin    = {"cli/cli_main": "nim-groth16"}.toTable()
 
 requires "https://github.com/status-im/nim-taskpools"
-requires "https://github.com/mratsim/constantine#5f7ba18f2ed351260015397c9eae079a6decaee1"
\ No newline at end of file
+requires "https://github.com/mratsim/constantine >= 0.2.0"
\ No newline at end of file
diff --git a/groth16/bn128/curves.nim b/groth16/bn128/curves.nim
index 42611b2..f8bb9f1 100644
--- a/groth16/bn128/curves.nim
+++ b/groth16/bn128/curves.nim
@@ -13,13 +13,12 @@
 #import constantine/platforms/abstractions
 #import constantine/math/isogenies/frobenius
 
-import constantine/math/arithmetic    except Fp, Fr
-import constantine/math/io/io_fields  except Fp, Fr
 import constantine/math/io/io_bigints
-import constantine/math/config/curves  
+import constantine/math/arithmetic
+import constantine/math/io/io_fields
 
-import constantine/math/config/type_ff          as tff except Fp, Fr
-import constantine/math/extension_fields/towers as ext except Fp, Fp2, Fp12, Fr
+import constantine/named/properties_fields      as tff
+import constantine/math/extension_fields/towers as ext
 
 import constantine/math/elliptic/ec_shortweierstrass_affine     as aff 
 import constantine/math/elliptic/ec_shortweierstrass_projective as prj 
@@ -30,19 +29,19 @@ import groth16/bn128/fields
 
 #-------------------------------------------------------------------------------
 
-type G1*   = aff.ECP_ShortW_Aff[Fp , aff.G1]
-type G2*   = aff.ECP_ShortW_Aff[Fp2, aff.G2]
+type G1*   = aff.EC_ShortW_Aff[Fp[BN254_Snarks] , aff.G1]
+type G2*   = aff.EC_ShortW_Aff[Fp2[BN254_Snarks], aff.G2]
 
-type ProjG1*  = prj.ECP_ShortW_Prj[Fp , prj.G1]
-type ProjG2*  = prj.ECP_ShortW_Prj[Fp2, prj.G2]
+type ProjG1*  = prj.EC_ShortW_Prj[Fp[BN254_Snarks] , prj.G1]
+type ProjG2*  = prj.EC_ShortW_Prj[Fp2[BN254_Snarks], prj.G2]
 
 #-------------------------------------------------------------------------------
 
-func unsafeMkG1* ( X, Y: Fp ) : G1 =
-  return aff.ECP_ShortW_Aff[Fp, aff.G1](x: X, y: Y)
+func unsafeMkG1* ( X, Y: Fp[BN254_Snarks] ) : G1 =
+  return aff.EC_ShortW_Aff[Fp[BN254_Snarks], aff.G1](x: X, y: Y)
 
-func unsafeMkG2* ( X, Y: Fp2 ) : G2 =
-  return aff.ECP_ShortW_Aff[Fp2, aff.G2](x: X, y: Y)
+func unsafeMkG2* ( X, Y: Fp2[BN254_Snarks] ) : G2 =
+  return aff.EC_ShortW_Aff[Fp2[BN254_Snarks], aff.G2](x: X, y: Y)
 
 #-------------------------------------------------------------------------------
 
@@ -51,15 +50,15 @@ const infG2*   : G2  = unsafeMkG2( zeroFp2 , zeroFp2 )
 
 #-------------------------------------------------------------------------------
 
-func checkCurveEqG1*( x, y: Fp ) : bool =
+func checkCurveEqG1*( x, y: Fp[BN254_Snarks] ) : bool =
   if bool(isZero(x)) and bool(isZero(y)):
     # the point at infinity is on the curve by definition
     return true
   else:
-    var x2 : Fp = squareFp(x)
-    var y2 : Fp = squareFp(y)
-    var x3 : Fp = x2 * x
-    var eq : Fp
+    var x2 = squareFp(x)
+    var y2 = squareFp(y)
+    var x3 = x2 * x
+    var eq : Fp[BN254_Snarks]
     eq =  x3
     eq += intToFp(3)
     eq -= y2
@@ -72,19 +71,19 @@ func checkCurveEqG1*( x, y: Fp ) : bool =
 # B = b1 + bu*u
 # b1 = 19485874751759354771024239261021720505790618469301721065564631296452457478373
 # b2 = 266929791119991161246907387137283842545076965332900288569378510910307636690
-const twistCoeffB_1 : Fp  = fromHex(Fp, "0x2b149d40ceb8aaae81be18991be06ac3b5b4c5e559dbefa33267e6dc24a138e5")
-const twistCoeffB_u : Fp  = fromHex(Fp, "0x009713b03af0fed4cd2cafadeed8fdf4a74fa084e52d1852e4a2bd0685c315d2")
-const twistCoeffB   : Fp2 = mkFp2( twistCoeffB_1 , twistCoeffB_u )
+const twistCoeffB_1 = fromHex(Fp[BN254_Snarks], "0x2b149d40ceb8aaae81be18991be06ac3b5b4c5e559dbefa33267e6dc24a138e5")
+const twistCoeffB_u = fromHex(Fp[BN254_Snarks], "0x009713b03af0fed4cd2cafadeed8fdf4a74fa084e52d1852e4a2bd0685c315d2")
+const twistCoeffB   = mkFp2( twistCoeffB_1 , twistCoeffB_u )
 
-func checkCurveEqG2*( x, y: Fp2 ) : bool =
+func checkCurveEqG2*( x, y: Fp2[BN254_Snarks] ) : bool =
   if isZeroFp2(x) and isZeroFp2(y):
     # the point at infinity is on the curve by definition
     return true
   else:
-    var x2 : Fp2 = squareFp2(x)
-    var y2 : Fp2 = squareFp2(y)
-    var x3 : Fp2 = x2 * x;
-    var eq : Fp2
+    var x2 = squareFp2(x)
+    var y2 = squareFp2(y)
+    var x3 = x2 * x
+    var eq : Fp2[BN254_Snarks]
     eq =  x3
     eq += twistCoeffB
     eq -= y2
@@ -92,14 +91,14 @@ func checkCurveEqG2*( x, y: Fp2 ) : bool =
 
 #-------------------------------------------------------------------------------
 
-func mkG1*( x, y: Fp ) : G1 =
+func mkG1*( x, y: Fp[BN254_Snarks] ) : G1 =
   if isZeroFp(x) and isZeroFp(y):
     return infG1
   else:
     assert( checkCurveEqG1(x,y) , "mkG1: not a G1 curve point" )
     return unsafeMkG1(x,y)
 
-func mkG2*( x, y: Fp2 ) : G2 =
+func mkG2*( x, y: Fp2[BN254_Snarks] ) : G2 =
   if isZeroFp2(x) and isZeroFp2(y):
     return infG2
   else:
@@ -109,16 +108,16 @@ func mkG2*( x, y: Fp2 ) : G2 =
 #-------------------------------------------------------------------------------
 # group generators
 
-const gen1_x  : Fp = fromHex(Fp, "0x01")
-const gen1_y  : Fp = fromHex(Fp, "0x02")
+const gen1_x  = fromHex(Fp[BN254_Snarks], "0x01")
+const gen1_y  = fromHex(Fp[BN254_Snarks], "0x02")
 
-const gen2_xi : Fp = fromHex(Fp, "0x1adcd0ed10df9cb87040f46655e3808f98aa68a570acf5b0bde23fab1f149701")
-const gen2_xu : Fp = fromHex(Fp, "0x09e847e9f05a6082c3cd2a1d0a3a82e6fbfbe620f7f31269fa15d21c1c13b23b")
-const gen2_yi : Fp = fromHex(Fp, "0x056c01168a5319461f7ca7aa19d4fcfd1c7cdf52dbfc4cbee6f915250b7f6fc8")
-const gen2_yu : Fp = fromHex(Fp, "0x0efe500a2d02dd77f5f401329f30895df553b878fc3c0dadaaa86456a623235c")
+const gen2_xi  = fromHex(Fp[BN254_Snarks], "0x1adcd0ed10df9cb87040f46655e3808f98aa68a570acf5b0bde23fab1f149701")
+const gen2_xu  = fromHex(Fp[BN254_Snarks], "0x09e847e9f05a6082c3cd2a1d0a3a82e6fbfbe620f7f31269fa15d21c1c13b23b")
+const gen2_yi  = fromHex(Fp[BN254_Snarks], "0x056c01168a5319461f7ca7aa19d4fcfd1c7cdf52dbfc4cbee6f915250b7f6fc8")
+const gen2_yu  = fromHex(Fp[BN254_Snarks], "0x0efe500a2d02dd77f5f401329f30895df553b878fc3c0dadaaa86456a623235c")
 
-const gen2_x  : Fp2 = mkFp2( gen2_xi, gen2_xu )
-const gen2_y  : Fp2 = mkFp2( gen2_yi, gen2_yu )
+const gen2_x   = mkFp2( gen2_xi, gen2_xu )
+const gen2_y   = mkFp2( gen2_yi, gen2_yu )
 
 const gen1* : G1 = unsafeMkG1( gen1_x, gen1_y )
 const gen2* : G2 = unsafeMkG2( gen2_x, gen2_y )
@@ -215,9 +214,9 @@ func `**`*( coeff: BigInt , point: G2 ) : G2 =
 
 #-------------------------------------------------------------------------------
 
-func pairing* (p: G1, q: G2) : Fp12 =
-  var t : Fp12
-  ate.pairing_bn[BN254Snarks]( t, p, q )
+func pairing* (p: G1, q: G2) : Fp12[BN254_Snarks] =
+  var t : Fp12[BN254_Snarks]
+  ate.pairing_bn( t, p, q )
   return t
 
 #-------------------------------------------------------------------------------
@@ -225,7 +224,8 @@ func pairing* (p: G1, q: G2) : Fp12 =
 proc sanityCheckGroupGen*() =
   echo( "gen1 on the curve  = ", checkCurveEqG1(gen1.x,gen1.y) )
   echo( "gen2 on the curve  = ", checkCurveEqG2(gen2.x,gen2.y) )
-  echo( "order of gen1 is R = ", (not bool(isInf(gen1))) and bool(isInf(primeR ** gen1)) )
-  echo( "order of gen2 is R = ", (not bool(isInf(gen2))) and bool(isInf(primeR ** gen2)) )
+  # TODO: fix compilation error with Constantine 0.2.0:
+  # echo( "order of gen1 is R = ", (not bool(isNeutral(gen1))) and bool(isNeutral(primeR ** gen1)) )
+  # echo( "order of gen2 is R = ", (not bool(isNeutral(gen2))) and bool(isNeutral(primeR ** gen2)) )
 
 #-------------------------------------------------------------------------------
diff --git a/groth16/bn128/debug.nim b/groth16/bn128/debug.nim
index ab7403c..7795294 100644
--- a/groth16/bn128/debug.nim
+++ b/groth16/bn128/debug.nim
@@ -9,23 +9,26 @@
 # equation: y^2 = x^3 + 3
 #
 
+import constantine/named/properties_fields
+import constantine/math/extension_fields/towers
+
 import groth16/bn128/fields
 import groth16/bn128/curves
 import groth16/bn128/io
 
 #-------------------------------------------------------------------------------
 
-proc debugPrintFp*(prefix: string, x: Fp) =
+proc debugPrintFp*(prefix: string, x: Fp[BN254_Snarks]) =
   echo(prefix & toDecimalFp(x))
 
-proc debugPrintFp2*(prefix: string, z: Fp2) =
+proc debugPrintFp2*(prefix: string, z: Fp2[BN254_Snarks]) =
   echo(prefix & " 1 ~> " & toDecimalFp(z.coords[0]))
   echo(prefix & " u ~> " & toDecimalFp(z.coords[1]))
 
-proc debugPrintFr*(prefix: string, x: Fr) =
+proc debugPrintFr*(prefix: string, x: Fr[BN254_Snarks]) =
   echo(prefix & toDecimalFr(x))
 
-proc debugPrintFrSeq*(msg: string, xs: seq[Fr]) =
+proc debugPrintFrSeq*(msg: string, xs: seq[Fr[BN254_Snarks]]) =
   echo "---------------------"
   echo msg
   for x in xs:
diff --git a/groth16/bn128/fields.nim b/groth16/bn128/fields.nim
index 5ffa8a1..2f52c50 100644
--- a/groth16/bn128/fields.nim
+++ b/groth16/bn128/fields.nim
@@ -14,22 +14,16 @@ import std/sequtils
 import constantine/math/arithmetic
 import constantine/math/io/io_fields
 import constantine/math/io/io_bigints
-import constantine/math/config/curves
-import constantine/math/config/type_ff          as tff
+import constantine/named/properties_fields      as tff
 import constantine/math/extension_fields/towers as ext
 
 #-------------------------------------------------------------------------------
 
 type B*    = BigInt[256]
-type Fr*   = tff.Fr[BN254Snarks]
-type Fp*   = tff.Fp[BN254Snarks]
 
-type Fp2*  = ext.QuadraticExt[Fp]
-type Fp12* = ext.Fp12[BN254Snarks]
-
-func mkFp2* (i: Fp, u: Fp) : Fp2 =
-  let c : array[2, Fp] = [i,u]
-  return ext.QuadraticExt[Fp]( coords: c )
+func mkFp2* (i: Fp[BN254_Snarks], u: Fp[BN254_Snarks]) : Fp2[BN254_Snarks] =
+  let c : array[2, Fp[BN254_Snarks]] = [i,u]
+  return ext.QuadraticExt[Fp[BN254_Snarks]]( coords: c )
 
 #-------------------------------------------------------------------------------
 
@@ -38,16 +32,16 @@ const primeR* : B = fromHex( B, "0x30644e72e131a029b85045b68181585d2833e84879b97
 
 #-------------------------------------------------------------------------------
 
-const zeroFp*  : Fp = fromHex( Fp, "0x00" )
-const zeroFr*  : Fr = fromHex( Fr, "0x00" )
-const oneFp*   : Fp = fromHex( Fp, "0x01" )
-const oneFr*   : Fr = fromHex( Fr, "0x01" )
+const zeroFp*   = fromHex( Fp[BN254_Snarks], "0x00" )
+const zeroFr*   = fromHex( Fr[BN254_Snarks], "0x00" )
+const oneFp*    = fromHex( Fp[BN254_Snarks], "0x01" )
+const oneFr*    = fromHex( Fr[BN254_Snarks], "0x01" )
 
-const zeroFp2* : Fp2 = mkFp2( zeroFp, zeroFp )
-const oneFp2*  : Fp2 = mkFp2( oneFp , zeroFp )
+const zeroFp2* = mkFp2( zeroFp, zeroFp )
+const oneFp2*  = mkFp2( oneFp , zeroFp )
 
-const minusOneFp* : Fp = fromHex( Fp, "0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd46" )
-const minusOneFr* : Fr = fromHex( Fr, "0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000000" )
+const minusOneFp* = fromHex( Fp[BN254_Snarks], "0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd46" )
+const minusOneFr* = fromHex( Fr[BN254_Snarks], "0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000000" )
 
 #-------------------------------------------------------------------------------
 
@@ -56,13 +50,13 @@ func intToB*(a: uint): B =
   y.setUint(a)
   return y
 
-func intToFp*(a: int): Fp =
-  var y : Fp
+func intToFp*(a: int): Fp[BN254_Snarks] =
+  var y : Fp[BN254_Snarks]
   y.fromInt(a)
   return y
 
-func intToFr*(a: int): Fr =
-  var y : Fr
+func intToFr*(a: int): Fr[BN254_Snarks] =
+  var y : Fr[BN254_Snarks]
   y.fromInt(a)
   return y
 
@@ -154,27 +148,27 @@ func smallPowFr*(base: Fr, expo: int): Fr =
 
 #-------------------------------------------------------------------------------
 
-func deltaFr*(i, j: int) : Fr =
+func deltaFr*[T](i, j: int) : Fr[T] =
   return (if (i == j): oneFr else: zeroFr)
 
 #-------------------------------------------------------------------------------
 
 # Montgomery batch inversion
-func batchInverseFr*( xs: seq[Fr] ) : seq[Fr] =
+func batchInverseFr*( xs: seq[Fr[BN254_Snarks]] ) : seq[Fr[BN254_Snarks]] =
   let n = xs.len
   assert(n>0)
-  var us : seq[Fr] = newSeq[Fr](n+1)
+  var us : seq[Fr[BN254_Snarks]] = newSeq[Fr[BN254_Snarks]](n+1)
   var a = xs[0]
   us[0] = oneFr
   us[1] = a
   for i in 1..<n: ( a *= xs[i] ; us[i+1] = a )
-  var vs : seq[Fr] = newSeq[Fr](n)
+  var vs : seq[Fr[BN254_Snarks]] = newSeq[Fr[BN254_Snarks]](n)
   vs[n-1] = invFr( us[n] )
   for i in countdown(n-2,0): vs[i] = vs[i+1] * xs[i+1]
   return collect( newSeq, (for i in 0..<n: us[i]*vs[i] ) )
 
 proc sanityCheckBatchInverseFr*() =
-  let xs : seq[Fr] = map( toSeq(101..137) , intToFr )
+  let xs = map( toSeq(101..137) , intToFr )
   let ys = batchInverseFr( xs )
   let zs = collect( newSeq, (for x in xs: invFr(x)) )
   let n = xs.len
diff --git a/groth16/bn128/io.nim b/groth16/bn128/io.nim
index 085898a..6426142 100644
--- a/groth16/bn128/io.nim
+++ b/groth16/bn128/io.nim
@@ -3,11 +3,11 @@
 import std/strutils
 import std/streams
 
-import constantine/math/arithmetic    except Fp, Fp2, Fr
-import constantine/math/io/io_fields  except Fp, Fp2, Fp
+import constantine/math/arithmetic
+import constantine/math/io/io_fields
 import constantine/math/io/io_bigints
-import constantine/math/config/curves
-import constantine/math/config/type_ff as tff except Fp, Fp2, Fr
+import constantine/named/properties_fields
+import constantine/math/extension_fields/towers
 
 import groth16/bn128/fields
 import groth16/bn128/curves
@@ -25,13 +25,13 @@ func toDecimalBig*[n](a : BigInt[n]): string =
   if s.len == 0: s="0"
   return s
 
-func toDecimalFp*(a : Fp): string =
+func toDecimalFp*(a : Fp[BN254_Snarks]): string =
   var s : string = toDecimal(a)
   s = s.strip( leading=true, trailing=false, chars={'0'} )
   if s.len == 0: s="0"
   return s
 
-func toDecimalFr*(a : Fr): string =
+func toDecimalFr*(a : Fr[BN254_Snarks]): string =
   var s : string = toDecimal(a)
   s = s.strip( leading=true, trailing=false, chars={'0'} )
   if s.len == 0: s="0"
@@ -41,13 +41,13 @@ func toDecimalFr*(a : Fr): string =
 
 const k65536 : BigInt[254] = fromHex( BigInt[254], "0x10000", bigEndian )
 
-func signedToDecimalFp*(a : Fp): string =
+func signedToDecimalFp*(a : Fp[BN254_Snarks]): string =
   if bool( a.toBig() > primeP_254 - k65536 ):
     return "-" & toDecimalFp(negFp(a))
   else:
     return toDecimalFp(a)
 
-func signedToDecimalFr*(a : Fr): string =
+func signedToDecimalFr*(a : Fr[BN254_Snarks]): string =
   if bool( a.toBig() > primeR_254 - k65536 ):
     return "-" & toDecimalFr(negFr(a))
   else:
@@ -58,38 +58,38 @@ func signedToDecimalFr*(a : Fr): string =
 #
 
 # R=2^256; this computes 2^256 mod Fp
-func calcFpMontR*() : Fp =
-  var x : Fp = intToFp(2)
+func calcFpMontR*() : Fp[BN254_Snarks] =
+  var x : Fp[BN254_Snarks] = intToFp(2)
   for i in 1..8:
     square(x)
   return x
 
 # R=2^256; this computes the inverse of (2^256 mod Fp)
-func calcFpInvMontR*() : Fp =
-  var x : Fp = calcFpMontR()
+func calcFpInvMontR*() : Fp[BN254_Snarks] =
+  var x : Fp[BN254_Snarks] = calcFpMontR()
   inv(x)
   return x
 
 # R=2^256; this computes 2^256 mod Fr
-func calcFrMontR*() : Fr =
-  var x : Fr = intToFr(2)
+func calcFrMontR*() : Fr[BN254_Snarks] =
+  var x : Fr[BN254_Snarks] = intToFr(2)
   for i in 1..8:
     square(x)
   return x
 
 # R=2^256; this computes the inverse of (2^256 mod Fp)
-func calcFrInvMontR*() : Fr =
-  var x : Fr = calcFrMontR()
+func calcFrInvMontR*() : Fr[BN254_Snarks] =
+  var x : Fr[BN254_Snarks] = calcFrMontR()
   inv(x)
   return x
 
 # apparently we cannot compute these in compile time for some reason or other... (maybe because `intToFp()`?)
-const fpMontR*    : Fp = fromHex( Fp, "0x0e0a77c19a07df2f666ea36f7879462c0a78eb28f5c70b3dd35d438dc58f0d9d" )
-const fpInvMontR* : Fp = fromHex( Fp, "0x2e67157159e5c639cf63e9cfb74492d9eb2022850278edf8ed84884a014afa37" )
+const fpMontR*    = fromHex( Fp[BN254_Snarks], "0x0e0a77c19a07df2f666ea36f7879462c0a78eb28f5c70b3dd35d438dc58f0d9d" )
+const fpInvMontR* = fromHex( Fp[BN254_Snarks], "0x2e67157159e5c639cf63e9cfb74492d9eb2022850278edf8ed84884a014afa37" )
 
 # apparently we cannot compute these in compile time for some reason or other... (maybe because `intToFp()`?)
-const frMontR*    : Fr = fromHex( Fr, "0x0e0a77c19a07df2f666ea36f7879462e36fc76959f60cd29ac96341c4ffffffb" )
-const frInvMontR* : Fr = fromHex( Fr, "0x15ebf95182c5551cc8260de4aeb85d5d090ef5a9e111ec87dc5ba0056db1194e" )
+const frMontR*    = fromHex( Fr[BN254_Snarks], "0x0e0a77c19a07df2f666ea36f7879462e36fc76959f60cd29ac96341c4ffffffb" )
+const frInvMontR* = fromHex( Fr[BN254_Snarks], "0x15ebf95182c5551cc8260de4aeb85d5d090ef5a9e111ec87dc5ba0056db1194e" )
 
 proc checkMontgomeryConstants*() =
   assert( bool( fpMontR    == calcFpMontR()    ) )
@@ -103,18 +103,18 @@ proc checkMontgomeryConstants*() =
 # the binary file `.zkey` used by the `circom` ecosystem uses little-endian 
 # Montgomery representation. So when we unmarshal with Constantine, it will 
 # give the wrong result. Calling this function on the result fixes that.
-func fromMontgomeryFp*(x : Fp) : Fp =
-  var y : Fp = x;
+func fromMontgomeryFp*(x : Fp[BN254_Snarks]) : Fp[BN254_Snarks] =
+  var y : Fp[BN254_Snarks] = x;
   y *= fpInvMontR
   return y
 
-func fromMontgomeryFr*(x : Fr) : Fr =
-  var y : Fr = x;
+func fromMontgomeryFr*(x : Fr[BN254_Snarks]) : Fr[BN254_Snarks] =
+  var y = x;
   y *= frInvMontR
   return y
 
-func toMontgomeryFr*(x : Fr) : Fr =
-  var y : Fr = x;
+func toMontgomeryFr*(x : Fr[BN254_Snarks]) : Fr[BN254_Snarks] =
+  var y = x;
   y *= frMontR
   return y
 
@@ -123,47 +123,47 @@ func toMontgomeryFr*(x : Fr) : Fr =
 # Note: in the `.zkey` coefficients, e apparently DOUBLE Montgomery encoding is used ?!?
 #
 
-func unmarshalFpMont* ( bs: array[32,byte] ) : Fp =
+func unmarshalFpMont* ( bs: array[32,byte] ) : Fp[BN254_Snarks] =
   var big : BigInt[254]
   unmarshal( big, bs, littleEndian );
-  var x : Fp
+  var x : Fp[BN254_Snarks]
   x.fromBig( big )
   return fromMontgomeryFp(x)
 
 # WTF Jordi, go home you are drunk
-func unmarshalFrWTF* ( bs: array[32,byte] ) : Fr =
+func unmarshalFrWTF* ( bs: array[32,byte] ) : Fr[BN254_Snarks] =
   var big : BigInt[254]
   unmarshal( big, bs, littleEndian );
-  var x : Fr
+  var x : Fr[BN254_Snarks]
   x.fromBig( big )
   return fromMontgomeryFr(fromMontgomeryFr(x))
 
-func unmarshalFrStd* ( bs: array[32,byte] ) : Fr =
+func unmarshalFrStd* ( bs: array[32,byte] ) : Fr[BN254_Snarks] =
   var big : BigInt[254]
   unmarshal( big, bs, littleEndian );
-  var x : Fr
+  var x : Fr[BN254_Snarks]
   x.fromBig( big )
   return x
 
-func unmarshalFrMont* ( bs: array[32,byte] ) : Fr =
+func unmarshalFrMont* ( bs: array[32,byte] ) : Fr[BN254_Snarks] =
   var big : BigInt[254]
   unmarshal( big, bs, littleEndian );
-  var x : Fr
+  var x : Fr[BN254_Snarks]
   x.fromBig( big )
   return fromMontgomeryFr(x)
 
 #-------------------------------------------------------------------------------
 
-func unmarshalFpMontSeq* ( len: int,  bs: openArray[byte] ) : seq[Fp] =
-  var vals  : seq[Fp] = newSeq[Fp]( len )
+func unmarshalFpMontSeq* ( len: int,  bs: openArray[byte] ) : seq[Fp[BN254_Snarks]] =
+  var vals  : seq[Fp[BN254_Snarks]] = newSeq[Fp[BN254_Snarks]]( len )
   var bytes : array[32,byte]
   for i in 0..<len:
     copyMem( addr(bytes) , unsafeAddr(bs[32*i]) , 32 )
     vals[i] = unmarshalFpMont( bytes )
   return vals
 
-func unmarshalFrMontSeq* ( len: int,  bs: openArray[byte] ) : seq[Fr] =
-  var vals  : seq[Fr] = newSeq[Fr]( len )
+func unmarshalFrMontSeq* ( len: int,  bs: openArray[byte] ) : seq[Fr[BN254_Snarks]] =
+  var vals  : seq[Fr[BN254_Snarks]] = newSeq[Fr[BN254_Snarks]]( len )
   var bytes : array[32,byte]
   for i in 0..<len:
     copyMem( addr(bytes) , unsafeAddr(bs[32*i]) , 32 )
@@ -172,7 +172,7 @@ func unmarshalFrMontSeq* ( len: int,  bs: openArray[byte] ) : seq[Fr] =
 
 #-------------------------------------------------------------------------------
 
-proc loadValueFrWTF*( stream: Stream ) : Fr =
+proc loadValueFrWTF*( stream: Stream ) : Fr[BN254_Snarks] =
   var bytes : array[32,byte]
   let n = stream.readData( addr(bytes), 32 )
   # for i in 0..<32: stdout.write(" " & toHex(bytes[i]))
@@ -180,45 +180,45 @@ proc loadValueFrWTF*( stream: Stream ) : Fr =
   assert( n == 32 )
   return unmarshalFrWTF(bytes)
 
-proc loadValueFrStd*( stream: Stream ) : Fr =
+proc loadValueFrStd*( stream: Stream ) : Fr[BN254_Snarks] =
   var bytes : array[32,byte]
   let n = stream.readData( addr(bytes), 32 )
   assert( n == 32 )
   return unmarshalFrStd(bytes)
 
-proc loadValueFrMont*( stream: Stream ) : Fr =
+proc loadValueFrMont*( stream: Stream ) : Fr[BN254_Snarks] =
   var bytes : array[32,byte]
   let n = stream.readData( addr(bytes), 32 )
   assert( n == 32 )
   return unmarshalFrMont(bytes)
 
-proc loadValueFpMont*( stream: Stream ) : Fp =
+proc loadValueFpMont*( stream: Stream ) : Fp[BN254_Snarks] =
   var bytes : array[32,byte]
   let n = stream.readData( addr(bytes), 32 )
   assert( n == 32 )
   return unmarshalFpMont(bytes)
 
-proc loadValueFp2Mont*( stream: Stream ) : Fp2 =
+proc loadValueFp2Mont*( stream: Stream ) : Fp2[BN254_Snarks] =
   let i = loadValueFpMont( stream )
   let u = loadValueFpMont( stream )
   return mkFp2(i,u)
 
 #---------------------------------------
 
-proc loadValuesFrStd*( len: int, stream: Stream ) : seq[Fr] =
-  var values : seq[Fr]
+proc loadValuesFrStd*( len: int, stream: Stream ) : seq[Fr[BN254_Snarks]] =
+  var values: seq[Fr[BN254_Snarks]]
   for i in 1..len:
     values.add( loadValueFrStd(stream) )
   return values
 
-proc loadValuesFpMont*( len: int, stream: Stream ) : seq[Fp] =
-  var values : seq[Fp]
+proc loadValuesFpMont*( len: int, stream: Stream ) : seq[Fp[BN254_Snarks]] =
+  var values : seq[Fp[BN254_Snarks]]
   for i in 1..len:
     values.add( loadValueFpMont(stream) )
   return values
 
-proc loadValuesFrMont*( len: int, stream: Stream ) : seq[Fr] =
-  var values : seq[Fr]
+proc loadValuesFrMont*( len: int, stream: Stream ) : seq[Fr[BN254_Snarks]] =
+  var values: seq[Fr[BN254_Snarks]]
   for i in 1..len:
     values.add( loadValueFrMont(stream) )
   return values
diff --git a/groth16/bn128/msm.nim b/groth16/bn128/msm.nim
index 71f51e9..b4b4005 100644
--- a/groth16/bn128/msm.nim
+++ b/groth16/bn128/msm.nim
@@ -10,15 +10,14 @@ import taskpools
 # import constantine/curves_primitives except Fp, Fp2, Fr
  
 import constantine/platforms/abstractions   except Subgroup
-import constantine/math/isogenies/frobenius except Subgroup
+import constantine/math/endomorphisms/frobenius except Subgroup
 
-import constantine/math/arithmetic     except Fp, Fp2, Fr
-import constantine/math/io/io_fields   except Fp, Fp2, Fr
 import constantine/math/io/io_bigints
-import constantine/math/config/curves  except G1, G2, Subgroup
-import constantine/math/config/type_ff except Fp, Fr, Subgroup
+import constantine/named/properties_fields except Subgroup
+import constantine/math/arithmetic
+import constantine/math/io/io_fields
 
-import constantine/math/extension_fields/towers                 as ext except Fp, Fp2, Fp12, Fr
+import constantine/math/extension_fields/towers                 as ext
 import constantine/math/elliptic/ec_shortweierstrass_affine     as aff except Subgroup
 import constantine/math/elliptic/ec_shortweierstrass_projective as prj except Subgroup
 import constantine/math/elliptic/ec_scalar_mul_vartime          as scl except Subgroup
@@ -32,7 +31,7 @@ import std/times
 
 #-------------------------------------------------------------------------------
 
-proc msmConstantineG1*( coeffs: openArray[Fr] , points: openArray[G1] ): G1 =
+proc msmConstantineG1*( coeffs: openArray[Fr[BN254_Snarks]] , points: openArray[G1] ): G1 =
 
   # let start = cpuTime()
 
@@ -60,7 +59,7 @@ proc msmConstantineG1*( coeffs: openArray[Fr] , points: openArray[G1] ): G1 =
 
 #---------------------------------------
 
-func msmConstantineG2*( coeffs: openArray[Fr] , points: openArray[G2] ): G2 =
+func msmConstantineG2*( coeffs: openArray[Fr[BN254_Snarks]] , points: openArray[G2] ): G2 =
 
   let N = coeffs.len
   assert( N == points.len, "incompatible sequence lengths" )
@@ -86,7 +85,7 @@ func msmConstantineG2*( coeffs: openArray[Fr] , points: openArray[G2] ): G2 =
 
 const task_multiplier : int = 1
 
-proc msmMultiThreadedG1*( nthreads_hint: int, coeffs: seq[Fr] , points: seq[G1] ): G1 =
+proc msmMultiThreadedG1*( coeffs: seq[Fr[BN254_Snarks]] , points: seq[G1], pool: Taskpool ): G1 =
 
   # for N <= 255 , we use 1 thread
   # for N == 256 , we use 2 threads
@@ -95,11 +94,10 @@ proc msmMultiThreadedG1*( nthreads_hint: int, coeffs: seq[Fr] , points: seq[G1]
 
   let N = coeffs.len
   assert( N == points.len, "incompatible sequence lengths" )
-  let nthreads_target = if (nthreads_hint<=0): countProcessors() else: min( nthreads_hint, 256 )
+  let nthreads_target = min( pool.numThreads, 256 )
   let nthreads = max( 1 , min( N div 128 , nthreads_target ) )
   let ntasks   = if nthreads>1: (nthreads*task_multiplier) else: 1
 
-  var pool = Taskpool.new(num_threads = nthreads)
   var pending : seq[FlowVar[mycurves.G1]] = newSeq[FlowVar[mycurves.G1]](ntasks)
 
   var a : int = 0
@@ -118,22 +116,18 @@ proc msmMultiThreadedG1*( nthreads_hint: int, coeffs: seq[Fr] , points: seq[G1]
   for k in 0..<ntasks:
     res += sync pending[k]
 
-  pool.syncAll()    
-  pool.shutdown()
-
   return res
 
 #---------------------------------------
 
-proc msmMultiThreadedG2*( nthreads_hint: int, coeffs: seq[Fr] , points: seq[G2] ): G2 =
+proc msmMultiThreadedG2*( coeffs: seq[Fr[BN254_Snarks]] , points: seq[G2], pool: Taskpool ): G2 =
 
   let N = coeffs.len
   assert( N == points.len, "incompatible sequence lengths" )
-  let nthreads_target = if (nthreads_hint<=0): countProcessors() else: min( nthreads_hint, 256 )
+  let nthreads_target = min( pool.numThreads, 256 )
   let nthreads = max( 1 , min( N div 128 , nthreads_target ) )
   let ntasks   = if nthreads>1: (nthreads*task_multiplier) else: 1
 
-  var pool = Taskpool.new(num_threads = nthreads)
   var pending : seq[FlowVar[mycurves.G2]] = newSeq[FlowVar[mycurves.G2]](ntasks)
 
   var a : int = 0
@@ -152,19 +146,16 @@ proc msmMultiThreadedG2*( nthreads_hint: int, coeffs: seq[Fr] , points: seq[G2]
   for k in 0..<ntasks:
     res += sync pending[k]
 
-  pool.syncAll()    
-  pool.shutdown()
-
   return res
 
 #-------------------------------------------------------------------------------
 
-func msmNaiveG1*( coeffs: seq[Fr] , points: seq[G1] ): G1 =
+func msmNaiveG1*( coeffs: seq[Fr[BN254_Snarks]] , points: seq[G1] ): G1 =
   let N = coeffs.len
   assert( N == points.len, "incompatible sequence lengths" )
 
   var s : ProjG1
-  s.setInf()
+  s.setNeutral()
 
   for i in 0..<N:
     var t : ProjG1
@@ -179,12 +170,12 @@ func msmNaiveG1*( coeffs: seq[Fr] , points: seq[G1] ): G1 =
 
 #---------------------------------------
 
-func msmNaiveG2*( coeffs: seq[Fr] , points: seq[G2] ): G2 =
+func msmNaiveG2*( coeffs: seq[Fr[BN254_Snarks]] , points: seq[G2] ): G2 =
   let N = coeffs.len
   assert( N == points.len, "incompatible sequence lengths" )
 
   var s : ProjG2
-  s.setInf()
+  s.setNeutral()
 
   for i in 0..<N:
     var t : ProjG2
@@ -199,8 +190,8 @@ func msmNaiveG2*( coeffs: seq[Fr] , points: seq[G2] ): G2 =
 
 #-------------------------------------------------------------------------------
 
-proc msmG1*( coeffs: seq[Fr] , points: seq[G1] ): G1 = msmConstantineG1(coeffs, points)
-proc msmG2*( coeffs: seq[Fr] , points: seq[G2] ): G2 = msmConstantineG2(coeffs, points)
+proc msmG1*( coeffs: seq[Fr[BN254_Snarks]] , points: seq[G1] ): G1 = msmConstantineG1(coeffs, points)
+proc msmG2*( coeffs: seq[Fr[BN254_Snarks]] , points: seq[G2] ): G2 = msmConstantineG2(coeffs, points)
 
 #-------------------------------------------------------------------------------
 
diff --git a/groth16/bn128/rnd.nim b/groth16/bn128/rnd.nim
index 10200fe..f2f3607 100644
--- a/groth16/bn128/rnd.nim
+++ b/groth16/bn128/rnd.nim
@@ -3,9 +3,10 @@ import std/random
 
 # import constantine/platforms/abstractions
 
-import constantine/math/arithmetic       except Fp, Fp2, Fr
-import constantine/math/io/io_fields     except Fp, Fp2, Fr
+import constantine/math/arithmetic
+import constantine/math/io/io_fields
 import constantine/math/io/io_bigints
+import constantine/named/properties_fields
 
 import groth16/bn128/fields
 
@@ -59,9 +60,9 @@ proc randBig*[bits: static int](): BigInt[bits] =
 
   return d
 
-proc randFr*(): Fr =
+proc randFr*(): Fr[BN254_Snarks] =
   let b : BigInt[254] = randBig[254]()
-  var y : Fr
+  var y : Fr[BN254_Snarks]
   y.fromBig( b )
   return y
 
diff --git a/groth16/fake_setup.nim b/groth16/fake_setup.nim
index b75dcf4..f98bd20 100644
--- a/groth16/fake_setup.nim
+++ b/groth16/fake_setup.nim
@@ -9,7 +9,8 @@
 import sugar
 import std/tables
 
-import constantine/math/arithmetic except Fp, Fr
+import constantine/math/arithmetic
+import constantine/named/properties_fields
 
 import groth16/bn128
 import groth16/math/domain
@@ -22,11 +23,11 @@ import groth16/misc
 
 type 
   ToxicWaste* = object
-    alpha*: Fr
-    beta*:  Fr
-    gamma*: Fr
-    delta*: Fr
-    tau*:   Fr
+    alpha*: Fr[BN254_Snarks]
+    beta*:  Fr[BN254_Snarks]
+    gamma*: Fr[BN254_Snarks]
+    delta*: Fr[BN254_Snarks]
+    tau*:   Fr[BN254_Snarks]
 
 proc randomToxicWaste*(): ToxicWaste = 
   let a = randFr()
@@ -73,9 +74,9 @@ type DenseMatrix*[T] = seq[DenseColumn[T]]
 
 type 
   DenseMatrices* = object
-    A* : DenseMatrix[Fr]
-    B* : DenseMatrix[Fr]
-    C* : DenseMatrix[Fr]
+    A* : DenseMatrix[Fr[BN254_Snarks]]
+    B* : DenseMatrix[Fr[BN254_Snarks]]
+    C* : DenseMatrix[Fr[BN254_Snarks]]
 
 #[
 
@@ -137,13 +138,13 @@ func denseMatricesToCoeffs*(matrices: DenseMatrices): seq[Coeff] =
 
 type SparseColumn*[T] = Table[int,T]
 
-proc columnInsertWithAddFr( col: var SparseColumn[Fr] , i: int,  y: Fr ) =
+proc columnInsertWithAddFr( col: var SparseColumn[Fr[BN254_Snarks]] , i: int,  y: Fr[BN254_Snarks] ) =
   var x = getOrDefault( col, i, zeroFr )
   x += y
   col[i] = x
 
-proc sparseDenseDotProdFr( U: SparseColumn[Fr], V: DenseColumn[Fr] ): Fr =
-  var acc : Fr = zeroFr
+proc sparseDenseDotProdFr( U: SparseColumn[Fr[BN254_Snarks]], V: DenseColumn[Fr[BN254_Snarks]] ): Fr[BN254_Snarks] =
+  var acc : Fr[BN254_Snarks] = zeroFr
   for i,x in U.pairs:
     acc += x * V[i]
   return acc
@@ -152,9 +153,9 @@ type SparseMatrix*[T] = seq[SparseColumn[T]]
 
 type 
   SparseMatrices* = object
-    A* : SparseMatrix[Fr]
-    B* : SparseMatrix[Fr]
-    C* : SparseMatrix[Fr]
+    A* : SparseMatrix[Fr[BN254_Snarks]]
+    B* : SparseMatrix[Fr[BN254_Snarks]]
+    C* : SparseMatrix[Fr[BN254_Snarks]]
 
 func r1csToSparseMatrices*(r1cs: R1CS): SparseMatrices =
   let n = r1cs.constraints.len
@@ -164,11 +165,11 @@ func r1csToSparseMatrices*(r1cs: R1CS): SparseMatrices =
   let logDomSize = ceilingLog2(n+p+1)
   let domSize    = 1 shl logDomSize
 
-  var matA, matB, matC: SparseMatrix[Fr]
+  var matA, matB, matC: SparseMatrix[Fr[BN254_Snarks]]
   for i in 0..<m:
-    var colA : SparseColumn[Fr] = initTable[int,Fr]()
-    var colB : SparseColumn[Fr] = initTable[int,Fr]()
-    var colC : SparseColumn[Fr] = initTable[int,Fr]()
+    var colA : SparseColumn[Fr[BN254_Snarks]] = initTable[int,Fr[BN254_Snarks]]()
+    var colB : SparseColumn[Fr[BN254_Snarks]] = initTable[int,Fr[BN254_Snarks]]()
+    var colC : SparseColumn[Fr[BN254_Snarks]] = initTable[int,Fr[BN254_Snarks]]()
     matA.add( colA )
     matB.add( colB )
     matC.add( colC )
@@ -188,10 +189,10 @@ func r1csToSparseMatrices*(r1cs: R1CS): SparseMatrices =
 
 #-------------------------------------------------------------------------------
 
-func dotProdFr(xs, ys: seq[Fr]): Fr = 
+func dotProdFr(xs, ys: seq[Fr[BN254_Snarks]]): Fr[BN254_Snarks] =
   let n = xs.len
   assert( n == ys.len, "dotProdFr: incompatible vector lengths" )
-  var s : Fr = zeroFr
+  var s : Fr[BN254_Snarks] = zeroFr
   for i in 0..<n:
     s += xs[i] * ys[i]
   return s
@@ -252,8 +253,8 @@ func fakeCircuitSetup*(r1cs: R1CS, toxic: ToxicWaste, flavour=Snarkjs): ZKey =
 
   # the Lagrange polynomials L_k(x) evaluated at x=tau
   # we can then simply take the dot product of these with the column vectors to compute the points A,B1,B2,C
-  let lagrangeTaus : seq[Fr] = collect( newSeq, (for k in 0..<domSize: evalLagrangePolyAt(D, k, toxic.tau) ))
-  
+  let lagrangeTaus : seq[Fr[BN254_Snarks]] = collect( newSeq, (for k in 0..<domSize: evalLagrangePolyAt(D, k, toxic.tau) ))
+
 #[
   # dense matrices use way too much memory
   let columnTausA  : seq[Fr] = collect( newSeq, (for col in matrices.A: dotProdFr(col,lagrangeTaus) ))
@@ -261,22 +262,22 @@ func fakeCircuitSetup*(r1cs: R1CS, toxic: ToxicWaste, flavour=Snarkjs): ZKey =
   let columnTausC  : seq[Fr] = collect( newSeq, (for col in matrices.C: dotProdFr(col,lagrangeTaus) ))
 ]#
 
-  let columnTausA  : seq[Fr] = collect( newSeq, (for col in matrices.A: sparseDenseDotProdFr(col,lagrangeTaus) ))
-  let columnTausB  : seq[Fr] = collect( newSeq, (for col in matrices.B: sparseDenseDotProdFr(col,lagrangeTaus) ))
-  let columnTausC  : seq[Fr] = collect( newSeq, (for col in matrices.C: sparseDenseDotProdFr(col,lagrangeTaus) ))
+  let columnTausA  : seq[Fr[BN254_Snarks]] = collect( newSeq, (for col in matrices.A: sparseDenseDotProdFr(col,lagrangeTaus) ))
+  let columnTausB  : seq[Fr[BN254_Snarks]] = collect( newSeq, (for col in matrices.B: sparseDenseDotProdFr(col,lagrangeTaus) ))
+  let columnTausC  : seq[Fr[BN254_Snarks]] = collect( newSeq, (for col in matrices.C: sparseDenseDotProdFr(col,lagrangeTaus) ))
 
   let pointsA  : seq[G1] = collect( newSeq , (for y in columnTausA: (y ** gen1) ))
   let pointsB1 : seq[G1] = collect( newSeq , (for y in columnTausB: (y ** gen1) ))
   let pointsB2 : seq[G2] = collect( newSeq , (for y in columnTausB: (y ** gen2) ))
   let pointsC  : seq[G1] = collect( newSeq , (for y in columnTausC: (y ** gen1) ))
 
-  let gammaInv : Fr = invFr(toxic.gamma)
-  let deltaInv : Fr = invFr(toxic.delta)
+  let gammaInv : Fr[BN254_Snarks] = invFr(toxic.gamma)
+  let deltaInv : Fr[BN254_Snarks] = invFr(toxic.delta)
 
-  let pointsL  : seq[G1] = collect( newSeq , (for j in 0..npub: 
+  let pointsL  : seq[G1] = collect( newSeq , (for j in 0..npub:
         gammaInv ** ( toxic.beta ** pointsA[j] + toxic.alpha ** pointsB1[j] + pointsC[j] ) ))
 
-  let pointsK  : seq[G1] = collect( newSeq , (for j in npub+1..nvars-1: 
+  let pointsK  : seq[G1] = collect( newSeq , (for j in npub+1..nvars-1:
         deltaInv ** ( toxic.beta ** pointsA[j] + toxic.alpha ** pointsB1[j] + pointsC[j] ) ))
  
   let polyZ    = vanishingPoly(D)
diff --git a/groth16/files/container.nim b/groth16/files/container.nim
index ac75f3c..28e2a0d 100644
--- a/groth16/files/container.nim
+++ b/groth16/files/container.nim
@@ -25,7 +25,7 @@ import std/streams
 
 import sugar
 
-import constantine/math/arithmetic except Fp, Fr
+import constantine/math/arithmetic
 import constantine/math/io/io_bigints
 
 #-------------------------------------------------------------------------------
@@ -61,7 +61,8 @@ proc readSection[T] ( expectedMagic: string
                     , stream: Stream
                     , user: var T
                     , callback: SectionCallback[T] 
-                    , filt: (int) -> bool ) =
+                    , filt: proc(_: int): bool {.raises: [IOError, OSError].}
+                    ) {.raises: [IOError, OSError].} =
 
   let sectId  = int( stream.readUint32() )
   let sectLen = int( stream.readUint64() )
@@ -77,7 +78,8 @@ proc parseContainer*[T] ( expectedMagic: string
                         , fname: string
                         , user: var T
                         , callback: SectionCallback[T] 
-                        , filt: (int) -> bool ) =
+                        , filt: proc(_: int): bool {.raises: [IOError, OSError].}
+                        ) {.raises: [IOError, OSError].} =
 
   let stream = newFileStream(fname, mode = fmRead)
   defer: stream.close()
diff --git a/groth16/files/export_json.nim b/groth16/files/export_json.nim
index 53b2828..a4d991f 100644
--- a/groth16/files/export_json.nim
+++ b/groth16/files/export_json.nim
@@ -3,19 +3,20 @@
 # export proof and public input in `circom`-compatible JSON files
 #
 
-import constantine/math/arithmetic   except Fp, Fr
-#import constantine/math/io/io_fields except Fp, Fr
+import constantine/math/arithmetic
+import constantine/named/properties_fields
+import constantine/math/extension_fields/towers
 
 import groth16/bn128
 from groth16/prover import Proof
 
 #-------------------------------------------------------------------------------
 
-func toQuotedDecimalFp(x: Fp): string = 
+func toQuotedDecimalFp(x: Fp[BN254_Snarks]): string =
   let s : string = toDecimalFp(x)
   return ("\"" & s & "\"")
 
-func toQuotedDecimalFr(x: Fr): string = 
+func toQuotedDecimalFr(x: Fr[BN254_Snarks]): string =
   let s : string = toDecimalFr(x)
   return ("\"" & s & "\"")
 
diff --git a/groth16/files/r1cs.nim b/groth16/files/r1cs.nim
index cb65712..ad3d4fa 100644
--- a/groth16/files/r1cs.nim
+++ b/groth16/files/r1cs.nim
@@ -51,8 +51,9 @@
 
 import std/streams
 
-import constantine/math/arithmetic except Fp, Fr
+import constantine/math/arithmetic
 import constantine/math/io/io_bigints
+import constantine/named/properties_fields
 
 import groth16/bn128
 import groth16/files/container
@@ -67,8 +68,8 @@ type
     nPubIn*  : int           # number of public inputs
     nPrivIn* : int           # number of private inputs
     nLabels* : int           # number of labels
-  
-  Term*       = tuple[ wireIdx: int, value: Fr ]
+
+  Term*       = tuple[ wireIdx: int, value: Fr[BN254_Snarks] ]
   LinComb*    = seq[Term]
   Constraint* = tuple[ A: LinComb, B: LinComb, C: LinComb ]
 
diff --git a/groth16/files/witness.nim b/groth16/files/witness.nim
index f849eb3..381e85c 100644
--- a/groth16/files/witness.nim
+++ b/groth16/files/witness.nim
@@ -16,8 +16,9 @@
 
 import std/streams
 
-import constantine/math/arithmetic except Fp, Fr
+import constantine/math/arithmetic
 import constantine/math/io/io_bigints
+import constantine/named/properties_fields
 
 import groth16/bn128
 import groth16/files/container
@@ -29,7 +30,7 @@ type
     curve*  : string
     r*      : BigInt[256] 
     nvars*  : int
-    values* : seq[Fr]
+    values* : seq[Fr[BN254_Snarks]]
 
 #-------------------------------------------------------------------------------
 
diff --git a/groth16/files/zkey.nim b/groth16/files/zkey.nim
index 73ce1e5..8ee4129 100644
--- a/groth16/files/zkey.nim
+++ b/groth16/files/zkey.nim
@@ -94,7 +94,7 @@
 
 import std/streams
 
-import constantine/math/arithmetic except Fp, Fr
+import constantine/math/arithmetic
 #import constantine/math/io/io_bigints
  
 import groth16/bn128
diff --git a/groth16/math/domain.nim b/groth16/math/domain.nim
index 2185814..160b418 100644
--- a/groth16/math/domain.nim
+++ b/groth16/math/domain.nim
@@ -3,9 +3,10 @@
 # power-of-two sized multiplicative FFT domains in the scalar field
 #
 
-import constantine/math/arithmetic   except Fp,Fr
-import constantine/math/io/io_fields except Fp,Fr
+import constantine/math/arithmetic
+import constantine/math/io/io_fields
 #import constantine/math/io/io_bigints
+import constantine/named/properties_fields
 
 import groth16/bn128
 import groth16/misc
@@ -14,27 +15,27 @@ import groth16/misc
 
 type 
   Domain* = object
-    domainSize*    : int         # `N = 2^n`
-    logDomainSize* : int         # `n = log2(N)`
-    domainGen*     : Fr          # `g`
-    invDomainGen*  : Fr          # `g^-1`
-    invDomainSize* : Fr          # `1/n`
+    domainSize*    : int              # `N = 2^n`
+    logDomainSize* : int              # `n = log2(N)`
+    domainGen*     : Fr[BN254_Snarks] # `g`
+    invDomainGen*  : Fr[BN254_Snarks] # `g^-1`
+    invDomainSize* : Fr[BN254_Snarks] # `1/n`
 
 #-------------------------------------------------------------------------------
 
 # the generator of the multiplicative subgroup with size `2^28`
-const gen28 : Fr = fromHex( Fr, "0x2a3c09f0a58a7e8500e0a7eb8ef62abc402d111e41112ed49bd61b6e725b19f0" )
+const gen28 = fromHex( Fr[BN254_Snarks], "0x2a3c09f0a58a7e8500e0a7eb8ef62abc402d111e41112ed49bd61b6e725b19f0" )
 
 func createDomain*(size: int): Domain = 
   let log2 = ceilingLog2(size)
   assert( (1 shl log2) == size , "domain must have a power-of-two size" )
 
   let expo : uint = 1'u shl (28 - log2)
-  let gen  : Fr   = smallPowFr(gen28, expo)
+  let gen  = smallPowFr(gen28, expo)
 
   let halfSize = size div 2
-  let a : Fr = smallPowFr(gen, size    )
-  let b : Fr = smallPowFr(gen, halfSize)
+  let a = smallPowFr(gen, size    )
+  let b = smallPowFr(gen, halfSize)
   assert(     bool(a == oneFr) , "domain generator sanity check /A" )
   assert( not bool(b == oneFr) , "domain generator sanity check /B" )
 
@@ -47,9 +48,9 @@ func createDomain*(size: int): Domain =
 
 #-------------------------------------------------------------------------------
 
-func enumerateDomain*(D: Domain): seq[Fr] =
-  var xs : seq[Fr] = newSeq[Fr](D.domainSize)
-  var g  : Fr = oneFr
+func enumerateDomain*(D: Domain): seq[Fr[BN254_Snarks]] =
+  var xs = newSeq[Fr[BN254_Snarks]](D.domainSize)
+  var g = oneFr
   for i in 0..<D.domainSize:
     xs[i] = g
     g *= D.domainGen
diff --git a/groth16/math/ntt.nim b/groth16/math/ntt.nim
index cea1281..b69a8ab 100644
--- a/groth16/math/ntt.nim
+++ b/groth16/math/ntt.nim
@@ -6,8 +6,9 @@
 
 #-------------------------------------------------------------------------------
 
-import constantine/math/arithmetic except Fp,Fr
+import constantine/math/arithmetic
 import constantine/math/io/io_fields
+import constantine/named/properties_fields
 
 import groth16/bn128
 import groth16/math/domain
@@ -16,10 +17,10 @@ import groth16/math/domain
 
 func forwardNTT_worker( m: int
                       , srcStride: int
-                      , gpows: seq[Fr]
-                      , src: seq[Fr]     , srcOfs: int
-                      , buf: var seq[Fr] , bufOfs: int
-                      , tgt: var seq[Fr] , tgtOfs: int ) =
+                      , gpows: seq[Fr[BN254_Snarks]]
+                      , src: seq[Fr[BN254_Snarks]]     , srcOfs: int
+                      , buf: var seq[Fr[BN254_Snarks]] , bufOfs: int
+                      , tgt: var seq[Fr[BN254_Snarks]] , tgtOfs: int ) =
   case m 
 
     of 0: 
@@ -45,25 +46,25 @@ func forwardNTT_worker( m: int
                        , buf , bufOfs + N
                        , buf , bufOfs + halfN )
       for j in 0..<halfN:
-        let y : Fr = gpows[j*srcStride] * buf[bufOfs+j+halfN] 
+        let y = gpows[j*srcStride] * buf[bufOfs+j+halFN]
         tgt[tgtOfs+j      ] = buf[bufOfs+j] + y
         tgt[tgtOfs+j+halfN] = buf[bufOfs+j] - y
 
 #---------------------------------------
 
 # forward number-theoretical transform (corresponds to polynomial evaluation)
-func forwardNTT*(src: seq[Fr], D: Domain): seq[Fr] =
+func forwardNTT*(src: seq[Fr[BN254_Snarks]], D: Domain): seq[Fr[BN254_Snarks]] =
   assert( D.domainSize == (1 shl D.logDomainSize) , "domain must have a power-of-two size" )
   assert( D.domainSize == src.len , "input must have the same size as the domain" )
-  var buf : seq[Fr] = newSeq[Fr]( 2 * D.domainSize )
-  var tgt : seq[Fr] = newSeq[Fr](     D.domainSize )
+  var buf = newSeq[Fr[BN254_Snarks] ]( 2 * D.domainSize )
+  var tgt = newSeq[Fr[BN254_Snarks]](     D.domainSize )
 
   # precalc powers of gen
   let N     = D.domainSize
   let halFN = N div 2
-  var gpows : seq[Fr] = newSeq[Fr]( halfN )
-  var x   : Fr = oneFr
-  let gen : Fr = D.domainGen
+  var gpows = newSeq[Fr[BN254_Snarks]]( halFN )
+  var x     = oneFr
+  let gen   = D.domainGen
   for i in 0..<halfN:
     gpows[i] = x
     x *= gen
@@ -78,28 +79,28 @@ func forwardNTT*(src: seq[Fr], D: Domain): seq[Fr] =
 
 # pads the input with zeros to get a pwoer of two size
 # TODO: optimize the FFT so that it doesn't do the multiplications with zeros 
-func extendAndForwardNTT*(src: seq[Fr], D: Domain): seq[Fr] =
+func extendAndForwardNTT*(src: seq[Fr[BN254_Snarks]], D: Domain): seq[Fr[BN254_Snarks]] =
   let n = src.len
   let N = D.domainSize 
   assert( n <= N )
   if n == N:
     return forwardNTT(src, D)
   else:
-    var padded : seq[Fr] = newSeq[Fr]( N )
+    var padded = newSeq[Fr[BN254_Snarks]]( N )
     for i in 0..<n: padded[i] = src[i]
     # for i in n..<N: padded[i] = zeroFr 
     return forwardNTT(padded, D)
 
 #-------------------------------------------------------------------------------
 
-const oneHalfFr* : Fr = fromHex(Fr, "0x183227397098d014dc2822db40c0ac2e9419f4243cdcb848a1f0fac9f8000001")
+const oneHalfFr* = fromHex(Fr[BN254_Snarks], "0x183227397098d014dc2822db40c0ac2e9419f4243cdcb848a1f0fac9f8000001")
 
 func inverseNTT_worker( m: int
                       , tgtStride: int
-                      , gpows: seq[Fr]
-                      , src: seq[Fr]     , srcOfs: int
-                      , buf: var seq[Fr] , bufOfs: int
-                      , tgt: var seq[Fr] , tgtOfs: int ) =
+                      , gpows: seq[Fr[BN254_Snarks]]
+                      , src: seq[Fr[BN254_Snarks]]     , srcOfs: int
+                      , buf: var seq[Fr[BN254_Snarks]] , bufOfs: int
+                      , tgt: var seq[Fr[BN254_Snarks]] , tgtOfs: int ) =
   case m 
 
     of 0: 
@@ -136,18 +137,18 @@ func inverseNTT_worker( m: int
 #---------------------------------------
 
 # inverse number-theoretical transform (corresponds to polynomial interpolation)
-func inverseNTT*(src: seq[Fr], D: Domain): seq[Fr] =
+func inverseNTT*(src: seq[Fr[BN254_Snarks]], D: Domain): seq[Fr[BN254_Snarks]] =
   assert( D.domainSize == (1 shl D.logDomainSize) , "domain must have a power-of-two size" )
   assert( D.domainSize == src.len , "input must have the same size as the domain" )
-  var buf : seq[Fr] = newSeq[Fr]( 2 * D.domainSize )
-  var tgt : seq[Fr] = newSeq[Fr](     D.domainSize )
+  var buf = newSeq[Fr[BN254_Snarks]]( 2 * D.domainSize )
+  var tgt = newSeq[Fr[BN254_Snarks]](     D.domainSize )
 
   # precalc 1/2 times powers of gen^-1
   let N     = D.domainSize
   let halFN = N div 2
-  var gpows : seq[Fr] = newSeq[Fr]( halfN )
-  var x    : Fr = oneHalfFr
-  let ginv : Fr = invFr( D.domainGen )
+  var gpows = newSeq[Fr[BN254_Snarks]]( halFN )
+  var x     = oneHalfFr
+  let ginv  = invFr( D.domainGen )
   for i in 0..<halfN:
     gpows[i] = x
     x *= ginv
diff --git a/groth16/math/poly.nim b/groth16/math/poly.nim
index f9e82d3..ff20e94 100644
--- a/groth16/math/poly.nim
+++ b/groth16/math/poly.nim
@@ -9,8 +9,9 @@
 import std/sequtils
 import std/sugar
 
-import constantine/math/arithmetic except Fp,Fr
+import constantine/math/arithmetic
 #import constantine/math/io/io_fields
+import constantine/named/properties_fields
 
 import groth16/bn128
 import groth16/math/domain
@@ -21,7 +22,7 @@ import groth16/misc
 
 type 
   Poly* = object
-    coeffs* : seq[Fr]
+    coeffs* : seq[Fr[BN254_Snarks]]
 
 #-------------------------------------------------------------------------------
 
@@ -41,8 +42,8 @@ func polyIsZero*(P: Poly) : bool =
   return b
 
 func polyIsEqual*(P, Q: Poly) : bool = 
-  let xs : seq[Fr] = P.coeffs ; let n = xs.len
-  let ys : seq[Fr] = Q.coeffs ; let m = ys.len
+  let xs = P.coeffs ; let n = xs.len
+  let ys = Q.coeffs ; let m = ys.len
   var b = true
   if n >= m:
     for i in 0..<m: ( if not isEqualFr(xs[i], ys[i]): ( b = false ; break ) )
@@ -54,10 +55,10 @@ func polyIsEqual*(P, Q: Poly) : bool =
 
 #-------------------------------------------------------------------------------
 
-func polyEvalAt*(P: Poly, x0: Fr): Fr =
+func polyEvalAt*(P: Poly, x0: Fr[BN254_Snarks]): Fr[BN254_Snarks] =
   let cs = P.coeffs ; let n = cs.len
-  var y : Fr = zeroFr
-  var r : Fr = oneFr
+  var y : Fr[BN254_Snarks] = zeroFr
+  var r : Fr[BN254_Snarks] = oneFr
   if n > 0: y = cs[0]
   for i in 1..<n:
     r *= x0
@@ -67,13 +68,13 @@ func polyEvalAt*(P: Poly, x0: Fr): Fr =
 #-------------------------------------------------------------------------------
 
 func polyNeg*(P: Poly) : Poly =
-  let zs : seq[Fr] = map( P.coeffs , negFr )
+  let zs = map( P.coeffs , negFr )
   return Poly(coeffs: zs)
 
 func polyAdd*(P, Q: Poly) : Poly =
   let xs = P.coeffs ; let n = xs.len
   let ys = Q.coeffs ; let m = ys.len
-  var zs : seq[Fr] = newSeq[Fr](max(n,m))
+  var zs = newSeq[Fr[BN254_Snarks]](max(n,m))
   if n >= m:
     for i in 0..<m: zs[i] = ( xs[i] + ys[i] )
     for i in m..<n: zs[i] = ( xs[i]         )
@@ -85,7 +86,7 @@ func polyAdd*(P, Q: Poly) : Poly =
 func polySub*(P, Q: Poly) : Poly =
   let xs = P.coeffs ; let n = xs.len
   let ys = Q.coeffs ; let m = ys.len
-  var zs : seq[Fr] = newSeq[Fr](max(n,m))
+  var zs = newSeq[Fr[BN254_Snarks]](max(n,m))
   if n >= m:
     for i in 0..<m: zs[i] = ( xs[i]  - ys[i] )
     for i in m..<n: zs[i] = ( xs[i]          )
@@ -97,7 +98,7 @@ func polySub*(P, Q: Poly) : Poly =
 #-------------------------------------------------------------------------------
 
 func polyScale*(s: Fr, P: Poly): Poly =
-  let zs : seq[Fr] = map( P.coeffs , proc (x: Fr): Fr = s*x )
+  let zs = map( P.coeffs , proc (x: Fr[BN254_Snarks]): Fr[BN254_Snarks] = s*x )
   return Poly(coeffs: zs)
 
 #-------------------------------------------------------------------------------
@@ -106,7 +107,7 @@ func polyMulNaive*(P, Q : Poly): Poly =
   let xs = P.coeffs ; let n1 = xs.len
   let ys = Q.coeffs ; let n2 = ys.len
   let N  = n1 + n2 - 1
-  var zs : seq[Fr] = newSeq[Fr](N) 
+  var zs = newSeq[Fr[BN254_Snarks]](N)
   for k in 0..<N:
     # 0 <= i <= min(k , n1-1)
     # 0 <= j <= min(k , n2-1)
@@ -132,10 +133,10 @@ func polyMulFFT*(P, Q: Poly): Poly =
   let N    : int    = (1 shl log2)
   let D    : Domain = createDomain( N )
 
-  let us : seq[Fr] = extendAndForwardNTT( P.coeffs, D )
-  let vs : seq[Fr] = extendAndForwardNTT( Q.coeffs, D )
-  let zs : seq[Fr] = collect( newSeq, (for i in 0..<N: us[i]*vs[i] ))
-  let ws : seq[Fr] = inverseNTT( zs, D ) 
+  let us = extendAndForwardNTT( P.coeffs, D )
+  let vs = extendAndForwardNTT( Q.coeffs, D )
+  let zs = collect( newSeq, (for i in 0..<N: us[i]*vs[i] ))
+  let ws = inverseNTT( zs, D )
 
   return Poly(coeffs: ws)
 
@@ -160,9 +161,9 @@ func `*`*(P: Poly, s: Fr  ): Poly  = return polyScale(s, P)
 #-------------------------------------------------------------------------------
 
 # the generalized vanishing polynomial `(a*x^N - b)`
-func generalizedVanishingPoly*(N: int, a: Fr, b: Fr): Poly = 
+func generalizedVanishingPoly*(N: int, a: Fr[BN254_Snarks], b: Fr[BN254_Snarks]): Poly =
   assert( N>=1 )
-  var cs : seq[Fr] = newSeq[Fr]( N+1 )
+  var cs = newSeq[Fr[BN254_Snarks]]( N+1 )
   cs[0] = negFr(b)
   cs[N] = a
   return Poly(coeffs: cs)
@@ -186,9 +187,9 @@ type
 func polyQuotRemByVanishing*(P: Poly, N: int): QuotRem[Poly] = 
   assert( N>=1 )
   let deg  : int     = polyDegree(P)
-  let src  : seq[Fr] = P.coeffs 
-  var quot : seq[Fr] = newSeq[Fr]( max(1, deg - N + 1) )
-  var rem  : seq[Fr] = newSeq[Fr]( N )
+  let src  = P.coeffs
+  var quot = newSeq[Fr[BN254_Snarks]]( max(1, deg - N + 1) )
+  var rem  = newSeq[Fr[BN254_Snarks]]( N )
 
   if deg < N:
     rem = src
@@ -222,14 +223,14 @@ func polyDivideByVanishing*(P: Poly, N: int): Poly =
 # Lagrange basis polynomials
 func lagrangePoly*(D: Domain, k: int): Poly =
   let N             = D.domainSize
-  let omMinusK : Fr = smallPowFr( D.invDomainGen , k )
-  let invN     : Fr = invFr(intToFr(N))
+  let omMinusK = smallPowFr( D.invDomainGen , k )
+  let invN     = invFr(intToFr(N))
 
-  var cs : seq[Fr]  = newSeq[Fr]( N )
+  var cs  = newSeq[Fr[BN254_Snarks]]( N )
   if k == 0:
     for i in 0..<N: cs[i] = invN
   else:
-    var s : Fr = invN
+    var s = invN
     for i in 0..<N: 
       cs[i] = s
       s *= omMinusK
@@ -239,7 +240,7 @@ func lagrangePoly*(D: Domain, k: int): Poly =
 #---------------------------------------
 
 # evaluate a Lagrange basis polynomial at a given point `zeta` (outside the domain)
-func evalLagrangePolyAt*(D: Domain, k: int, zeta: Fr): Fr =
+func evalLagrangePolyAt*(D: Domain, k: int, zeta: Fr[BN254_Snarks]): Fr[BN254_Snarks] =
   let omegaK = smallPowFr(D.domainGen, k)
   let denom  = (zeta - omegaK)
   if bool(isZero(denom)):
@@ -252,16 +253,16 @@ func evalLagrangePolyAt*(D: Domain, k: int, zeta: Fr): Fr =
 #-------------------------------------------------------------------------------
 
 # evaluates a polynomial on an FFT domain
-func polyForwardNTT*(P: Poly, D: Domain): seq[Fr] =
+func polyForwardNTT*(P: Poly, D: Domain): seq[Fr[BN254_Snarks]] =
   let n = P.coeffs.len
   assert( n <= D.domainSize , "the domain must be as least as big as the polynomial" )
-  let src : seq[Fr] = P.coeffs
+  let src = P.coeffs
   return forwardNTT(src, D)
 
 #---------------------------------------
 
 # interpolates a polynomial on an FFT domain
-func polyInverseNTT*(ys: seq[Fr], D: Domain): Poly =
+func polyInverseNTT*(ys: seq[Fr[BN254_Snarks]], D: Domain): Poly =
   let n = ys.len
   assert( n == D.domainSize , "the domain must be same size as the input" )
   let tgt = inverseNTT(ys, D)
diff --git a/groth16/prover.nim b/groth16/prover.nim
index 70fdcd1..8f1bb89 100644
--- a/groth16/prover.nim
+++ b/groth16/prover.nim
@@ -7,6 +7,8 @@
 # See <https://geometry.xyz/notebook/the-hidden-little-secret-in-snarkjs>
 #
 
+{.push raises:[].}
+
 #[
 import sugar
 import constantine/math/config/curves 
@@ -20,10 +22,11 @@ import std/times
 import std/cpuinfo
 import system
 import taskpools
+import constantine/math/arithmetic
+import constantine/named/properties_fields
+import constantine/math/extension_fields/towers
 
-import constantine/math/arithmetic except Fp, Fr
 #import constantine/math/io/io_extfields except Fp12
-#import constantine/math/extension_fields/towers except Fp2, Fp12  
 
 import groth16/bn128
 import groth16/math/domain
@@ -36,7 +39,7 @@ import groth16/misc
 
 type
   Proof* = object
-    publicIO* : seq[Fr]
+    publicIO* : seq[Fr[BN254_Snarks]]
     pi_a*     : G1
     pi_b*     : G2
     pi_c*     : G1
@@ -48,17 +51,17 @@ type
 
 type
   ABC = object
-    valuesAz : seq[Fr]
-    valuesBz : seq[Fr]
-    valuesCz : seq[Fr]
+    valuesAz : seq[Fr[BN254_Snarks]]
+    valuesBz : seq[Fr[BN254_Snarks]]
+    valuesCz : seq[Fr[BN254_Snarks]]
 
 # computes the vectors A*z, B*z, C*z where z is the witness
-func buildABC( zkey: ZKey, witness: seq[Fr] ): ABC = 
+func buildABC( zkey: ZKey, witness: seq[Fr[BN254_Snarks]] ): ABC =
   let hdr: GrothHeader = zkey.header
   let domSize = hdr.domainSize
 
-  var valuesAz : seq[Fr] = newSeq[Fr](domSize)
-  var valuesBz : seq[Fr] = newSeq[Fr](domSize)
+  var valuesAz = newSeq[Fr[BN254_Snarks]](domSize)
+  var valuesBz = newSeq[Fr[BN254_Snarks]](domSize)
 
   for entry in zkey.coeffs:
     case entry.matrix 
@@ -66,7 +69,7 @@ func buildABC( zkey: ZKey, witness: seq[Fr] ): ABC =
       of MatrixB: valuesBz[entry.row] += entry.coeff * witness[entry.col]
       else: raise newException(AssertionDefect, "fatal error")
 
-  var valuesCz : seq[Fr] = newSeq[Fr](domSize)
+  var valuesCz = newSeq[Fr[BN254_Snarks]](domSize)
   for i in 0..<domSize:
     valuesCz[i] = valuesAz[i] * valuesBz[i]
 
@@ -93,29 +96,43 @@ func computeQuotientNaive( abc: ABC ): Poly=
 #---------------------------------------
 
 # returns [ eta^i * xs[i] | i<-[0..n-1] ]
-func multiplyByPowers( xs: seq[Fr], eta: Fr ): seq[Fr] = 
+func multiplyByPowers( xs: seq[Fr[BN254_Snarks]], eta: Fr[BN254_Snarks] ): seq[Fr[BN254_Snarks]] =
   let n = xs.len
   assert(n >= 1)
-  var ys : seq[Fr] = newSeq[Fr](n)
+  var ys = newSeq[Fr[BN254_Snarks]](n)
   ys[0] = xs[0]
   if n >= 1: ys[1] = eta * xs[1]
-  var spow : Fr = eta
+  var spow = eta
   for i in 2..<n: 
     spow *= eta
     ys[i] = spow * xs[i]
   return ys
 
 # interpolates a polynomial, shift the variable by `eta`, and compute the shifted values
-func shiftEvalDomain( values: seq[Fr], D: Domain, eta: Fr ): seq[Fr] =
+func shiftEvalDomain(
+  values: seq[Fr[BN254_Snarks]],
+  D: Domain,
+  eta: Fr[BN254_Snarks] ): seq[Fr[BN254_Snarks]] =
   let poly : Poly = polyInverseNTT( values , D )
-  let cs : seq[Fr] = poly.coeffs
-  var ds : seq[Fr] = multiplyByPowers( cs, eta )
+  let cs : seq[Fr[BN254_Snarks]] = poly.coeffs
+  var ds : seq[Fr[BN254_Snarks]] = multiplyByPowers( cs, eta )
   return polyForwardNTT( Poly(coeffs:ds), D )
 
+# Wraps shiftEvalDomain such that it can be called by Taskpool.spawn. The result
+# is written to the output parameter. Has an unused return type because
+# Taskpool.spawn cannot handle a void return type.
+func shiftEvalDomainTask(
+  values: seq[Fr[BN254_Snarks]],
+  D: Domain,
+  eta: Fr[BN254_Snarks],
+  output: ptr Isolated[seq[Fr[BN254_Snarks]]]): bool =
+
+  output[] = isolate shiftEvalDomain(values, D, eta)
+
 # computes the quotient polynomial Q = (A*B - C) / Z
 # by computing the values on a shifted domain, and interpolating the result
 # remark: Q has degree `n-2`, so it's enough to use a domain of size n
-proc computeQuotientPointwise( nthreads: int, abc: ABC ): Poly =
+proc computeQuotientPointwise( abc: ABC, pool: TaskPool ): Poly =
   let n    = abc.valuesAz.len
   assert( abc.valuesBz.len == n )
   assert( abc.valuesCz.len == n )
@@ -127,24 +144,25 @@ proc computeQuotientPointwise( nthreads: int, abc: ABC ): Poly =
   let eta   = createDomain(2*n).domainGen
   let invZ1 = invFr( smallPowFr(eta,n) - oneFr )
 
-  var pool = Taskpool.new(num_threads = nthreads)
+  var outputA1, outputB1, outputC1: Isolated[seq[Fr[BN254_Snarks]]]
 
-  var A1fv : FlowVar[seq[Fr]] = pool.spawn shiftEvalDomain( abc.valuesAz, D, eta )
-  var B1fv : FlowVar[seq[Fr]] = pool.spawn shiftEvalDomain( abc.valuesBz, D, eta )
-  var C1fv : FlowVar[seq[Fr]] = pool.spawn shiftEvalDomain( abc.valuesCz, D, eta )
+  var taskA1 = pool.spawn shiftEvalDomainTask( abc.valuesAz, D, eta, addr outputA1 )
+  var taskB1 = pool.spawn shiftEvalDomainTask( abc.valuesBz, D, eta, addr outputB1 )
+  var taskC1 = pool.spawn shiftEvalDomainTask( abc.valuesCz, D, eta, addr outputC1 )
 
-  let A1 = sync A1fv
-  let B1 = sync B1fv
-  let C1 = sync C1fv
+  discard sync taskA1
+  discard sync taskB1
+  discard sync taskC1
 
-  var ys : seq[Fr] = newSeq[Fr]( n )
+  let A1 = outputA1.extract()
+  let B1 = outputB1.extract()
+  let C1 = outputC1.extract()
+
+  var ys : seq[Fr[BN254_Snarks]] = newSeq[Fr[BN254_Snarks]]( n )
   for j in 0..<n: ys[j] = ( A1[j]*B1[j] - C1[j] ) * invZ1
   let Q1 = polyInverseNTT( ys, D )
   let cs = multiplyByPowers( Q1.coeffs, invFr(eta) )
 
-  pool.syncAll() 
-  pool.shutdown()
-
   return Poly(coeffs: cs)
 
 #---------------------------------------
@@ -155,29 +173,30 @@ proc computeQuotientPointwise( nthreads: int, abc: ABC ): Poly =
 # (shifted) Lagrange bases.
 # see <https://geometry.xyz/notebook/the-hidden-little-secret-in-snarkjs>
 #
-proc computeSnarkjsScalarCoeffs( nthreads: int, abc: ABC): seq[Fr] =
+proc computeSnarkjsScalarCoeffs( abc: ABC, pool: TaskPool ): seq[Fr[BN254_Snarks]] =
   let n    = abc.valuesAz.len
   assert( abc.valuesBz.len == n )
   assert( abc.valuesCz.len == n )
   let D    = createDomain(n)
   let eta  = createDomain(2*n).domainGen
 
-  var pool = Taskpool.new(num_threads = nthreads)
+  var outputA1, outputB1, outputC1: Isolated[seq[Fr[BN254_Snarks]]]
 
-  var A1fv : FlowVar[seq[Fr]] = pool.spawn shiftEvalDomain( abc.valuesAz, D, eta )
-  var B1fv : FlowVar[seq[Fr]] = pool.spawn shiftEvalDomain( abc.valuesBz, D, eta )
-  var C1fv : FlowVar[seq[Fr]] = pool.spawn shiftEvalDomain( abc.valuesCz, D, eta )
+  var taskA1 = pool.spawn shiftEvalDomainTask( abc.valuesAz, D, eta, addr outputA1 )
+  var taskB1 = pool.spawn shiftEvalDomainTask( abc.valuesBz, D, eta, addr outputB1 )
+  var taskC1 = pool.spawn shiftEvalDomainTask( abc.valuesCz, D, eta, addr outputC1 )
 
-  let A1 = sync A1fv
-  let B1 = sync B1fv
-  let C1 = sync C1fv
+  discard sync taskA1
+  discard sync taskB1
+  discard sync taskC1
 
-  var ys : seq[Fr] = newSeq[Fr]( n )
+  let A1 = outputA1.extract()
+  let B1 = outputB1.extract()
+  let C1 = outputC1.extract()
+
+  var ys : seq[Fr[BN254_Snarks]] = newSeq[Fr[BN254_Snarks]]( n )
   for j in 0..<n: ys[j] = ( A1[j] * B1[j] - C1[j] ) 
 
-  pool.syncAll() 
-  pool.shutdown()
-
   return ys
 
 #[
@@ -209,10 +228,10 @@ proc computeSnarkjsScalarCoeffs( nthreads: int, abc: ABC ): seq[Fr] =
 
 type
   Mask* = object
-    r*: Fr              # masking coefficients 
-    s*: Fr              # for zero knowledge
+    r*: Fr[BN254_Snarks]              # masking coefficients
+    s*: Fr[BN254_Snarks]              # for zero knowledge
 
-proc generateProofWithMask*( nthreads: int, printTimings: bool, zkey: ZKey, wtns: Witness, mask: Mask ): Proof =
+proc generateProofWithMask*( zkey: ZKey, wtns: Witness, mask: Mask, pool: Taskpool, printTimings: bool): Proof =
 
   when not (defined(gcArc) or defined(gcOrc) or defined(gcAtomicArc)):
     {.fatal: "Compile with arc/orc!".}
@@ -236,7 +255,7 @@ proc generateProofWithMask*( nthreads: int, printTimings: bool, zkey: ZKey, wtns
   assert( nvars == witness.len , "wrong witness length" )
 
   # remark: with the special variable "1" we actuall have (npub+1) public IO variables
-  var pubIO : seq[Fr] = newSeq[Fr]( npubs + 1)
+  var pubIO = newSeq[Fr[BN254_Snarks]]( npubs + 1)
   for i in 0..npubs: pubIO[i] = witness[i]             
 
   start = cpuTime()
@@ -245,21 +264,21 @@ proc generateProofWithMask*( nthreads: int, printTimings: bool, zkey: ZKey, wtns
     abc = buildABC( zkey, witness )
 
   start = cpuTime()
-  var qs : seq[Fr]
+  var qs : seq[Fr[BN254_Snarks]]
   withMeasureTime(printTimings,"computing the quotient (FFTs)"):
     case zkey.header.flavour
 
       # the points H are [delta^-1 * tau^i * Z(tau)]
       of JensGroth:
-        let polyQ = computeQuotientPointwise( nthreads, abc )
+        let polyQ = computeQuotientPointwise( abc, pool )
         qs = polyQ.coeffs
   
       # the points H are `[delta^-1 * L_{2i+1}(tau)]_1`
       # where L_i are Lagrange basis polynomials on the double-sized domain
       of Snarkjs:
-        qs = computeSnarkjsScalarCoeffs( nthreads, abc )
+        qs = computeSnarkjsScalarCoeffs( abc, pool )
 
-  var zs : seq[Fr] = newSeq[Fr]( nvars - npubs - 1 )
+  var zs = newSeq[Fr[BN254_Snarks]]( nvars - npubs - 1 )
   for j in npubs+1..<nvars:
     zs[j-npubs-1] = witness[j]
 
@@ -279,41 +298,41 @@ proc generateProofWithMask*( nthreads: int, printTimings: bool, zkey: ZKey, wtns
   withMeasureTime(printTimings,"computing pi_A (G1 MSM)"):
     pi_a =  spec.alpha1
     pi_a += r ** spec.delta1
-    pi_a += msmMultiThreadedG1( nthreads , witness , pts.pointsA1 )
+    pi_a += msmMultiThreadedG1( witness , pts.pointsA1, pool )
 
   var rho : G1 
   withMeasureTime(printTimings,"computing rho (G1 MSM)"):
     rho =  spec.beta1
     rho += s ** spec.delta1
-    rho += msmMultiThreadedG1( nthreads , witness , pts.pointsB1 )
+    rho += msmMultiThreadedG1(  witness , pts.pointsB1, pool )
 
   var pi_b : G2
   withMeasureTime(printTimings,"computing pi_B (G2 MSM)"):
     pi_b =  spec.beta2
     pi_b += s ** spec.delta2
-    pi_b += msmMultiThreadedG2( nthreads , witness , pts.pointsB2 )
+    pi_b += msmMultiThreadedG2( witness , pts.pointsB2, pool )
 
   var pi_c : G1
   withMeasureTime(printTimings,"computing pi_C (2x G1 MSM)"):
     pi_c =  s ** pi_a
     pi_c += r ** rho
     pi_c += negFr(r*s) ** spec.delta1
-    pi_c += msmMultiThreadedG1( nthreads, qs , pts.pointsH1 )
-    pi_c += msmMultiThreadedG1( nthreads, zs , pts.pointsC1 )
+    pi_c += msmMultiThreadedG1( qs , pts.pointsH1, pool )
+    pi_c += msmMultiThreadedG1( zs , pts.pointsC1, pool )
 
   return Proof( curve:"bn128", publicIO:pubIO, pi_a:pi_a, pi_b:pi_b, pi_c:pi_c )
 
 #-------------------------------------------------------------------------------
 
-proc generateProofWithTrivialMask*( nthreads: int, printTimings: bool, zkey: ZKey, wtns: Witness ): Proof =
+proc generateProofWithTrivialMask*( zkey: ZKey, wtns: Witness, pool: Taskpool, printTimings: bool ): Proof =
   let mask = Mask( r: zeroFr , s: zeroFr )
-  return generateProofWithMask( nthreads, printTimings, zkey, wtns, mask )
+  return generateProofWithMask( zkey, wtns, mask, pool, printTimings )
 
-proc generateProof*( nthreads: int, printTimings: bool, zkey: ZKey, wtns: Witness ): Proof =
+proc generateProof*( zkey: ZKey, wtns: Witness, pool: Taskpool, printTimings = false ): Proof =
 
   # masking coeffs
-  let r : Fr = randFr()
-  let s : Fr = randFr()
+  let r = randFr()
+  let s = randFr()
   let mask = Mask(r: r, s: s)
 
-  return generateProofWithMask( nthreads, printTimings, zkey, wtns, mask )
+  return generateProofWithMask( zkey, wtns, mask, pool, printTimings )
diff --git a/groth16/verifier.nim b/groth16/verifier.nim
index 18bf5aa..d60be47 100644
--- a/groth16/verifier.nim
+++ b/groth16/verifier.nim
@@ -16,8 +16,9 @@ import ./zkey
 ]#
 
 # import constantine/math/arithmetic except Fp, Fr
-import constantine/math/io/io_extfields except Fp12
-import constantine/math/extension_fields/towers except Fp2, Fp12  
+import constantine/math/io/io_extfields
+import constantine/named/properties_fields
+import constantine/math/extension_fields/towers
 
 import groth16/bn128
 import groth16/zkey_types
@@ -38,12 +39,12 @@ proc verifyProof* (vkey: VKey, prf: Proof): bool =
 
   var pubG1 : G1 = msmG1( prf.publicIO , vkey.vpoints.pointsIC )  
 
-  let lhs   : Fp12 = pairing( negG1(prf.pi_a) , prf.pi_b )          # < -pi_a   , pi_b  >
-  let rhs1  : Fp12 = vkey.spec.alphaBeta                            # < alpha  , beta  >
-  let rhs2  : Fp12 = pairing( prf.pi_c , vkey.spec.delta2 )         # < pi_c   , delta >
-  let rhs3  : Fp12 = pairing( pubG1    , vkey.spec.gamma2 )         # < sum... , gamma >
+  let lhs   = pairing( negG1(prf.pi_a) , prf.pi_b )          # < -pi_a   , pi_b  >
+  let rhs1  = vkey.spec.alphaBeta                            # < alpha  , beta  >
+  let rhs2  = pairing( prf.pi_c , vkey.spec.delta2 )         # < pi_c   , delta >
+  let rhs3  = pairing( pubG1    , vkey.spec.gamma2 )         # < sum... , gamma >
 
-  var eq : Fp12
+  var eq : Fp12[BN254_Snarks]
   eq =  lhs  
   eq *= rhs1
   eq *= rhs2
diff --git a/groth16/zkey_types.nim b/groth16/zkey_types.nim
index 16b3969..1ce1d0b 100644
--- a/groth16/zkey_types.nim
+++ b/groth16/zkey_types.nim
@@ -1,5 +1,7 @@
 
-import constantine/math/arithmetic except Fp, Fr
+import constantine/math/arithmetic
+import constantine/named/properties_fields
+import constantine/math/extension_fields/towers
 
 import groth16/bn128
 
@@ -28,7 +30,7 @@ type
     gamma2*      : G2                # = gamma * g2
     delta1*      : G1                # = delta * g1
     delta2*      : G2                # = delta * g2       
-    alphaBeta*   : Fp12              # = <alpha1 , beta2>
+    alphaBeta*   : Fp12[BN254_Snarks] # = <alpha1 , beta2>
 
   VerifierPoints* = object
     pointsIC*    : seq[G1]           # the points `delta^-1 * ( beta*A_j(tau) + alpha*B_j(tau) + C_j(tau) ) * g1` (for j <= npub)
@@ -49,7 +51,7 @@ type
     matrix* : MatrixSel
     row*    : int
     col*    : int
-    coeff*  : Fr
+    coeff*  : Fr[BN254_Snarks]
 
   ZKey* = object
     # sectionMask* : uint32
diff --git a/tests/groth16/testProver.nim b/tests/groth16/testProver.nim
index cdf8a5b..3bb3eda 100644
--- a/tests/groth16/testProver.nim
+++ b/tests/groth16/testProver.nim
@@ -2,6 +2,8 @@
 import std/unittest
 import std/sequtils
 
+import taskpools
+
 import groth16/prover
 import groth16/verifier
 import groth16/fake_setup
@@ -44,7 +46,7 @@ const myR1CS =
       )
 
 # the equation we want prove is `7*11*13 + 1022 == 2023`
-let myWitnessValues : seq[Fr] = map( @[ 1, 2023, 1022, 7, 11, 13, 7*11, 7*11*13 ] , intToFr )
+let myWitnessValues = map( @[ 1, 2023, 1022, 7, 11, 13, 7*11, 7*11*13 ] , intToFr )
 # wire indices:         ^^^^^^^         0      1    2  3   4   5    6      7
 
 let myWitness = 
@@ -57,9 +59,11 @@ let myWitness =
 #-------------------------------------------------------------------------------
 
 proc testProof(zkey: ZKey, witness: Witness): bool = 
-  let proof = generateProof( zkey, witness )
+  var pool = Taskpool.new()
+  let proof = generateProof( zkey, witness, pool )
   let vkey  = extractVKey( zkey)
   let ok    = verifyProof( vkey, proof )
+  pool.shutdown()
   return ok
 
 suite "prover":