diff --git a/tests/math_fields/t_finite_fields.nim b/tests/math_fields/t_finite_fields.nim
index 79afa4b96..d070b91f9 100644
--- a/tests/math_fields/t_finite_fields.nim
+++ b/tests/math_fields/t_finite_fields.nim
@@ -166,6 +166,240 @@ proc main() =
           # Check equality when converting back to natural domain
           9'u64 == cast[uint64](r_bytes)
 
+    test "Addition mod 2^31 - 2^27 + 1":
+      block:
+        var x, y, z: Fp[BabyBear]
+
+        x.fromUint(80'u32)
+        y.fromUint(10'u32)
+        z.fromUint(90'u32)
+
+        x += y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 90'u32
+
+      block:
+        var x, y, z: Fp[BabyBear]
+
+        x.fromUint(0x78000000'u32) # p-1
+        y.fromUint(1'u32)
+        z.fromUint(0'u32)
+
+        x += y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 0'u32
+
+    test "Substraction mod 2^31 - 2^27 + 1":
+      block:
+        var x, y, z: Fp[BabyBear]
+
+        x.fromUint(80'u32)
+        y.fromUint(10'u32)
+        z.fromUint(70'u32)
+
+        x -= y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 70'u32
+
+      block:
+        var x, y, z: Fp[BabyBear]
+
+        x.fromUint(0'u32)
+        y.fromUint(1'u32)
+        z.fromUint(0x78000000'u32) # p-1
+
+        x -= y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 0x78000000'u32
+
+    test "Multiplication mod 2^31 - 2^27 + 1":
+      block:
+        var x, y, z, r: Fp[BabyBear]
+
+        x.fromUint(10'u32)
+        y.fromUint(10'u32)
+        z.fromUint(100'u32)
+
+        r.prod(x, y)
+
+        var r_bytes: array[4, byte]
+        r_bytes.marshal(r, cpuEndian)
+        let new_r = cast[uint32](r_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == r)
+          # Check equality when converting back to natural domain
+          new_r == 100'u32
+
+      block:
+        var x, y, z, r: Fp[BabyBear]
+
+        x.fromUint(0x10000'u32)
+        y.fromUint(0x10000'u32)
+        z.fromUint(uint32((0x10000'u64 * 0x10000'u64) mod 0x78000001'u64))
+
+        r.prod(x, y)
+
+        var r_bytes: array[4, byte]
+        r_bytes.marshal(r, cpuEndian)
+        let new_r = cast[uint32](r_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == r)
+          # Check equality when converting back to natural domain
+          new_r == uint32((0x10000'u64 * 0x10000'u64) mod 0x78000001'u64)
+
+    test "Addition mod 2^31 - 2^24 + 1":
+      block:
+        var x, y, z: Fp[KoalaBear]
+
+        x.fromUint(50'u32)
+        y.fromUint(20'u32)
+        z.fromUint(70'u32)
+
+        x += y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 70'u32
+
+      block:
+        var x, y, z: Fp[KoalaBear]
+
+        x.fromUint(0x7f000000'u32) # p-1
+        y.fromUint(1'u32)
+        z.fromUint(0'u32)
+
+        x += y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 0'u32
+
+    test "Substraction mod 2^31 - 2^24 + 1":
+      block:
+        var x, y, z: Fp[KoalaBear]
+
+        x.fromUint(70'u32)
+        y.fromUint(20'u32)
+        z.fromUint(50'u32)
+
+        x -= y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 50'u32
+
+      block:
+        var x, y, z: Fp[KoalaBear]
+
+        x.fromUint(0'u32)
+        y.fromUint(1'u32)
+        z.fromUint(0x7f000000'u32) # p-1
+
+        x -= y
+
+        var x_bytes: array[4, byte]
+        x_bytes.marshal(x, cpuEndian)
+        let new_x = cast[uint32](x_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == x)
+          # Check equality when converting back to natural domain
+          new_x == 0x7f000000'u32
+
+    test "Multiplication mod 2^31 - 2^24 + 1":
+      block:
+        var x, y, z, r: Fp[KoalaBear]
+
+        x.fromUint(12'u32)
+        y.fromUint(12'u32)
+        z.fromUint(144'u32)
+
+        r.prod(x, y)
+
+        var r_bytes: array[4, byte]
+        r_bytes.marshal(r, cpuEndian)
+        let new_r = cast[uint32](r_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == r)
+          # Check equality when converting back to natural domain
+          new_r == 144'u32
+
+      block:
+        var x, y, z, r: Fp[KoalaBear]
+
+        x.fromUint(0x10000'u32)
+        y.fromUint(0x20000'u32)
+        z.fromUint(uint32((0x10000'u64 * 0x20000'u64) mod 0x7f000001'u64))
+
+        r.prod(x, y)
+
+        var r_bytes: array[4, byte]
+        r_bytes.marshal(r, cpuEndian)
+        let new_r = cast[uint32](r_bytes)
+
+        check:
+          # Check equality in the Montgomery domain
+          bool(z == r)
+          # Check equality when converting back to natural domain
+          new_r == uint32((0x10000'u64 * 0x20000'u64) mod 0x7f000001'u64)
+
     test "Addition mod 2^61 - 1":
       block:
         var x, y, z: Fp[Mersenne61]
@@ -302,7 +536,6 @@ proc main() =
           # Check equality when converting back to natural domain
           new_r == 2'u64
 
-
 main()
 
 proc largeField() =
diff --git a/tests/math_fields/t_finite_fields_mulsquare.nim b/tests/math_fields/t_finite_fields_mulsquare.nim
index f33b33201..c306df2ff 100644
--- a/tests/math_fields/t_finite_fields_mulsquare.nim
+++ b/tests/math_fields/t_finite_fields_mulsquare.nim
@@ -81,6 +81,8 @@ proc mainSanity() =
   suite "Modular squaring is consistent with multiplication on special elements" & " [" & $WordBitWidth & "-bit words]":
     sanity Fake101
     sanity Mersenne61
+    sanity BabyBear
+    sanity KoalaBear
     sanity Mersenne127
     sanity P224         # P224 uses the fast-path with 64-bit words and the slow path with 32-bit words
     sanity P256
diff --git a/tests/math_fields/t_io_fields.nim b/tests/math_fields/t_io_fields.nim
index c22a05be6..642c7c300 100644
--- a/tests/math_fields/t_io_fields.nim
+++ b/tests/math_fields/t_io_fields.nim
@@ -46,6 +46,96 @@ proc main() =
         marshal(r_bytes, f, littleEndian)
         check: x_bytes == r_bytes
 
+      # BabyBear ---------------------------------
+      block:
+        # "Little-endian" - 0
+        let x = 0'u64
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[BabyBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      block:
+        # "Little-endian" - 1
+        let x = 1'u64
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[BabyBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      block:
+        # "Little-endian" - 2^15
+        let x = 1'u64 shl 15
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[BabyBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      block:
+        # "Little-endian" - p-1
+        let x = 2013265921'u64 - 1
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[BabyBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      # KoalaBear ---------------------------------
+      block:
+        # "Little-endian" - 0
+        let x = 0'u64
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[KoalaBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      block:
+        # "Little-endian" - 1
+        let x = 1'u64
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[KoalaBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      block:
+        # "Little-endian" - 2^15
+        let x = 1'u64 shl 15
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[KoalaBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
+      block:
+        # "Little-endian" - p-1
+        let x = 2130706433'u64 - 1
+        let x_bytes = cast[array[8, byte]](x)
+        var f: Fp[KoalaBear]
+        f.fromUint(x)
+
+        var r_bytes: array[8, byte]
+        marshal(r_bytes, f, littleEndian)
+        check: x_bytes == r_bytes
+
       # Mersenne 61 ---------------------------------
       block:
         # "Little-endian" - 0