diff --git a/tests/poseidon2/fuzzing.nim b/tests/poseidon2/fuzzing.nim
new file mode 100644
index 0000000..9affa0c
--- /dev/null
+++ b/tests/poseidon2/fuzzing.nim
@@ -0,0 +1,5 @@
+import std/random
+
+randomize()
+let seed*: int64 = int64.rand()
+randomize(seed)
diff --git a/tests/poseidon2/reference.nim b/tests/poseidon2/reference.nim
new file mode 100644
index 0000000..0e4376b
--- /dev/null
+++ b/tests/poseidon2/reference.nim
@@ -0,0 +1,50 @@
+import std/sequtils
+import constantine/math/arithmetic
+import constantine/math/io/io_fields
+import poseidon2/types
+import poseidon2/io
+import poseidon2/compress
+
+const KeyNone = F.fromHex("0x0")
+const KeyBottomLayer = F.fromHex("0x1")
+const KeyOdd = F.fromHex("0x2")
+const KeyOddAndBottomLayer = F.fromhex("0x3")
+
+# Reference implementation of merkle root algorithm
+# Only used in tests
+func merkleRoot(xs: openArray[F], isBottomLayer: static bool) : F =
+  let a = low(xs)
+  let b = high(xs)
+  let m = b-a+1
+
+  when isBottomLayer:
+    assert m > 0, "merkle root of empty sequence is not defined"
+
+  when not isBottomLayer:
+    if m==1:
+      return xs[a]
+
+  let halfn  : int  = m div 2
+  let n      : int  = 2*halfn
+  let isOdd : bool = (n != m)
+
+  var ys : seq[F]
+  if not isOdd:
+    ys = newSeq[F](halfn)
+  else:
+    ys = newSeq[F](halfn+1)
+
+  for i in 0..<halfn:
+    const key = when isBottomLayer: KeyBottomLayer else: KeyNone
+    ys[i] = compress( xs[a+2*i], xs[a+2*i+1], key = key )
+  if isOdd:
+    const key = when isBottomLayer: KeyOddAndBottomLayer else: KeyOdd
+    ys[halfn] = compress( xs[n], zero, key = key )
+
+  return merkleRoot(ys, isBottomLayer = false)
+
+func merkleRoot*(xs: openArray[F]) : F =
+  merkleRoot(xs, isBottomLayer = true)
+
+func merkleRoot*(bytes: openArray[byte]): F =
+  merkleRoot(toSeq bytes.elements(F))
diff --git a/tests/poseidon2/testMerkle.nim b/tests/poseidon2/testMerkle.nim
index 84bfd66..da8d2cf 100644
--- a/tests/poseidon2/testMerkle.nim
+++ b/tests/poseidon2/testMerkle.nim
@@ -2,6 +2,7 @@ import std/unittest
 import std/math
 import std/sequtils
 import std/sugar
+import std/random
 
 import constantine/math/arithmetic
 import constantine/math/io/io_fields
@@ -13,6 +14,9 @@ import poseidon2/io
 import poseidon2/compress
 import poseidon2/merkle
 
+import ./fuzzing
+import ./reference
+
 suite "merkle root":
 
   const isBottomLayer = 1
@@ -219,3 +223,22 @@ suite "merkle root test vectors":
       let input = collect(newSeq, (for i in 1..n: byte(i)))
       let root = Merkle.digest(input)
       check root.toDecimal == expected[n]
+
+suite "merkle root fuzzing (seed: " & $fuzzing.seed & ")":
+
+  test "merkle root algorithm matches reference implementation":
+
+    proc checkInput(input: seq[byte]) =
+      let expected = reference.merkleRoot(input)
+      check bool(Merkle.digest(input) == expected)
+
+    # a couple of tests with small input
+    for _ in 0..<1000:
+      let len = rand(1024)
+      let input = newSeqwith(len, byte.rand())
+      checkInput(input)
+
+    # one test with larger input
+    let len = rand(1024 * 1024)
+    let input = newSeqwith(len, byte.rand())
+    checkInput(input)