diff --git a/Cargo.toml b/Cargo.toml
index 453a7be..df28fc3 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -17,6 +17,7 @@ opt-level = "z"
[workspace.dependencies]
ark-bls12-381 = "0.4.0"
+ark-bn254 = "0.4.0"
ark-ec = "0.4.2"
ark-ff = "0.4.2"
ark-std = "0.4.0"
diff --git a/[Tho09]poseidon-hash/Cargo.toml b/[Tho09]poseidon-hash/Cargo.toml
index 29080ff..dfe7b2a 100644
--- a/[Tho09]poseidon-hash/Cargo.toml
+++ b/[Tho09]poseidon-hash/Cargo.toml
@@ -5,3 +5,5 @@ version = "0.1.0"
[dependencies]
ark-ff = { workspace = true }
+ark-bn254 = { workspace = true }
+ark-std = { workspace = true }
diff --git a/[Tho09]poseidon-hash/README.md b/[Tho09]poseidon-hash/README.md
index e69de29..2d4a01e 100644
--- a/[Tho09]poseidon-hash/README.md
+++ b/[Tho09]poseidon-hash/README.md
@@ -0,0 +1,49 @@
+# Poseidon Hasher
+Poseidon Hasher is a mapping over strings of $F_p$ (for prime $p > 2^{31}$) such that it maps $F_p^* \to F_p^o$ where $o$ is the number of output elements (often chosen value is $o = 1$).
+
+Poseidon is said to be a variant of *HadesMiMC* construction however with a fixed and known key.
+
+## A primer on sponge construction
+Sponge construction looks as follows:
+![sponge_construction](sponge_construction.png)
+In this, the $I$ is maintained state that changes over time, $m_x$ are injected values to be hashed and $h_y$ are the output elements.
+
+General construction looks as follows:
+- Depending on the use case, determine the capacity element value and the input padding if needed.
+- Split the obtained input into chunks of size $r$.
+- Apply the permutation to the capacity element and the first chunk.
+- Until no more chunks are left, add them into the state and apply the permutation.
+- Output $o$ output elements out of the rate part of the state.
+If needed, iterate the permutation more times.
+
+## The HADES design strategy
+The HADES design strategy consists of:
+- First, $R_f$ rounds in the beginning, in which S-boxes
+are applied to the full state.
+- Next, $R_p$ rounds in the middle contain single S-box application. Rest of the state goes through this phase unchanged
+- Finally, $R_f$ rounds in the end, in which S-boxes
+are applied to the full state.
+
+![hades_construction](hades_construction.png)
+
+Each such round consists of the following three sub-steps:
+1. $ARC$: Add round constants
+2. $SBOX$: Application of S-Boxes
+3. $MIX$: Mix layers
+
+## Reference implementation for magic values
+Directory `reference/hadeshash` includes `generate_params_poseidong.sage` for generating values used inside poseidon hasher. We build values for BLS12-381's $F_q$ a.k.a Scalar Field. For this, use:
+
+```bash
+# For generating for BLS12-381 Fq;
+# Modulus = `52435875175126190479447740508185965837690552500527637822603658699938581184513`
+# In Hex that is `0x73EDA753299D7D483339D80809A1D80553BDA402FFFE5BFEFFFFFFFF00000001`
+# So, the following command should work:
+sage generate_params_poseidon.sage 1 0 255 3 3 128 0x73EDA753299D7D483339D80809A1D80553BDA402FFFE5BFEFFFFFFFF00000001
+```
+This generates `poseidon_params_n255_t3_alpha3_M128.txt` file. Using values in this file, we generate values ingested in the rust code.
+
+
+## References
+- https://eprint.iacr.org/2019/458.pdf
+- https://github.com/arnaucube/poseidon-ark/
diff --git a/[Tho09]poseidon-hash/hades_construction.png b/[Tho09]poseidon-hash/hades_construction.png
new file mode 100644
index 0000000..2ad9b52
Binary files /dev/null and b/[Tho09]poseidon-hash/hades_construction.png differ
diff --git a/[Tho09]poseidon-hash/reference/hadeshash/LICENSE b/[Tho09]poseidon-hash/reference/hadeshash/LICENSE
new file mode 100644
index 0000000..1522178
--- /dev/null
+++ b/[Tho09]poseidon-hash/reference/hadeshash/LICENSE
@@ -0,0 +1,19 @@
+Copyright (c) 2019 Graz University of Technology
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the ""Software""), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED *AS IS*, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
\ No newline at end of file
diff --git a/[Tho09]poseidon-hash/reference/hadeshash/README.md b/[Tho09]poseidon-hash/reference/hadeshash/README.md
new file mode 100644
index 0000000..880b457
--- /dev/null
+++ b/[Tho09]poseidon-hash/reference/hadeshash/README.md
@@ -0,0 +1,12 @@
+# Scripts and Reference Implementations of Poseidon and Starkad
+This repository contains the source code of reference implementations for various versions of Poseidon [1]. The source code is available in Sage. Moreover, scripts to calculate the round numbers, the round constants, and the MDS matrices are also included.
+
+### Update from 02/05/2023
+The script `generate_params_poseidon.sage` should be used to compute the round numbers and to generate the round constants and matrices. The scripts `calc_round_numbers.py` and `generate_parameters_grain.sage` are deprecated and should not be used anymore.
+
+### Update from 07/03/2021
+We fixed several bugs in the implementation. First, the linear layer was computed as `state = state * M` instead of `state = M * state`, and secondly the final matrix multiplication was missing. The test vectors were also changed accordingly.
+
+
+
+[1] *Poseidon: A New Hash Function for Zero-Knowledge Proof Systems*. Cryptology ePrint Archive, Report 2019/458. [https://eprint.iacr.org/2019/458](https://eprint.iacr.org/2019/458). Accepted at USENIX'21.
diff --git a/[Tho09]poseidon-hash/reference/hadeshash/code/calc_round_numbers.py b/[Tho09]poseidon-hash/reference/hadeshash/code/calc_round_numbers.py
new file mode 100644
index 0000000..f7b7881
--- /dev/null
+++ b/[Tho09]poseidon-hash/reference/hadeshash/code/calc_round_numbers.py
@@ -0,0 +1,145 @@
+from math import *
+import sys
+import Crypto.Util.number
+
+def sat_inequiv_alpha(p, t, R_F, R_P, alpha, M):
+ n = ceil(log(p, 2))
+ N = int(n * t)
+ if alpha > 0:
+ R_F_1 = 6 if M <= ((floor(log(p, 2) - ((alpha-1)/2.0))) * (t + 1)) else 10 # Statistical
+ R_F_2 = 1 + ceil(log(2, alpha) * min(M, n)) + ceil(log(t, alpha)) - R_P # Interpolation
+ #R_F_3 = ceil(min(n, M) / float(3*log(alpha, 2))) - R_P # Groebner 1
+ #R_F_3 = ((log(2, alpha) / float(2)) * min(n, M)) - R_P # Groebner 1
+ R_F_3 = 1 + (log(2, alpha) * min(M/float(3), log(p, 2)/float(2))) - R_P # Groebner 1
+ R_F_4 = t - 1 + min((log(2, alpha) * M) / float(t+1), ((log(2, alpha)*log(p, 2)) / float(2))) - R_P # Groebner 2
+ #R_F_5 = ((1.0/(2*log((alpha**alpha)/float((alpha-1)**(alpha-1)), 2))) * min(n, M) + t - 2 - R_P) / float(t - 1) # Groebner 3
+ R_F_max = max(ceil(R_F_1), ceil(R_F_2), ceil(R_F_3), ceil(R_F_4))
+ return (R_F >= R_F_max)
+ elif alpha == (-1):
+ R_F_1 = 6 if M <= ((floor(log(p, 2) - 2)) * (t + 1)) else 10 # Statistical
+ R_P_1 = 1 + ceil(0.5 * min(M, n)) + ceil(log(t, 2)) - floor(R_F * log(t, 2)) # Interpolation
+ R_P_2 = 1 + ceil(0.5 * min(M, n)) + ceil(log(t, 2)) - floor(R_F * log(t, 2))
+ R_P_3 = t - 1 + ceil(log(t, 2)) + min(ceil(M / float(t+1)), ceil(0.5*log(p, 2))) - floor(R_F * log(t, 2)) # Groebner 2
+ R_F_max = ceil(R_F_1)
+ R_P_max = max(ceil(R_P_1), ceil(R_P_2), ceil(R_P_3))
+ return (R_F >= R_F_max and R_P >= R_P_max)
+ else:
+ print("Invalid value for alpha!")
+ exit(1)
+
+def get_sbox_cost(R_F, R_P, N, t):
+ return int(t * R_F + R_P)
+
+def get_size_cost(R_F, R_P, N, t):
+ n = ceil(float(N) / t)
+ return int((N * R_F) + (n * R_P))
+
+def get_depth_cost(R_F, R_P, N, t):
+ return int(R_F + R_P)
+
+def find_FD_round_numbers(p, t, alpha, M, cost_function, security_margin):
+ n = ceil(log(p, 2))
+ N = int(n * t)
+
+ sat_inequiv = sat_inequiv_alpha
+
+ R_P = 0
+ R_F = 0
+ min_cost = float("inf")
+ max_cost_rf = 0
+ # Brute-force approach
+ for R_P_t in range(1, 500):
+ for R_F_t in range(4, 100):
+ if R_F_t % 2 == 0:
+ if (sat_inequiv(p, t, R_F_t, R_P_t, alpha, M) == True):
+ if security_margin == True:
+ R_F_t += 2
+ R_P_t = int(ceil(float(R_P_t) * 1.075))
+ cost = cost_function(R_F_t, R_P_t, N, t)
+ if (cost < min_cost) or ((cost == min_cost) and (R_F_t < max_cost_rf)):
+ R_P = ceil(R_P_t)
+ R_F = ceil(R_F_t)
+ min_cost = cost
+ max_cost_rf = R_F
+ return (int(R_F), int(R_P))
+
+def calc_final_numbers_fixed(p, t, alpha, M, security_margin):
+ # [Min. S-boxes] Find best possible for t and N
+ n = ceil(log(p, 2))
+ N = int(n * t)
+ cost_function = get_sbox_cost
+ ret_list = []
+ (R_F, R_P) = find_FD_round_numbers(p, t, alpha, M, cost_function, security_margin)
+ min_sbox_cost = cost_function(R_F, R_P, N, t)
+ ret_list.append(R_F)
+ ret_list.append(R_P)
+ ret_list.append(min_sbox_cost)
+
+ # [Min. Size] Find best possible for t and N
+ # Minimum number of S-boxes for fixed n results in minimum size also (round numbers are the same)!
+ min_size_cost = get_size_cost(R_F, R_P, N, t)
+ ret_list.append(min_size_cost)
+
+ return ret_list # [R_F, R_P, min_sbox_cost, min_size_cost]
+
+def print_latex_table_combinations(combinations, alpha, security_margin):
+ for comb in combinations:
+ N = comb[0]
+ t = comb[1]
+ M = comb[2]
+ n = int(N / t)
+ prime = Crypto.Util.number.getPrime(n)
+ ret = calc_final_numbers_fixed(prime, t, alpha, M, security_margin)
+ field_string = "\mathbb F_{p}"
+ sbox_string = "x^{" + str(alpha) + "}"
+ print("$" + str(M) + "$ & $" + str(N) + "$ & $" + str(n) + "$ & $" + str(t) + "$ & $" + str(ret[0]) + "$ & $" + str(ret[1]) + "$ & $" + field_string + "$ & $" + str(ret[2]) + "$ & $" + str(ret[3]) + "$ \\\\")
+
+# Single tests
+# print calc_final_numbers_fixed(Crypto.Util.number.getPrime(64), 24, 3, 128, True)
+# print calc_final_numbers_fixed(Crypto.Util.number.getPrime(253), 6, -1, 128, True)
+print(calc_final_numbers_fixed(Crypto.Util.number.getPrime(255), 3, 5, 128, True))
+print(calc_final_numbers_fixed(Crypto.Util.number.getPrime(255), 6, 5, 128, True))
+print(calc_final_numbers_fixed(Crypto.Util.number.getPrime(254), 3, 5, 128, True))
+print(calc_final_numbers_fixed(Crypto.Util.number.getPrime(254), 6, 5, 128, True))
+print(calc_final_numbers_fixed(Crypto.Util.number.getPrime(64), 24, 3, 128, True))
+
+# x^5 (254-bit prime number)
+#prime = 0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001
+x_5_combinations = [
+ [1536, 2, 128], [1536, 4, 128], [1536, 6, 128], [1536, 8, 128], [1536, 16, 128],
+ [1536, 2, 256], [1536, 4, 256], [1536, 6, 256], [1536, 8, 256], [1536, 16, 256]
+]
+
+# With security margin
+print("--- Table x^5 WITH security margin ---")
+print_latex_table_combinations(x_5_combinations, 5, True)
+
+# Without security margin
+print("--- Table x^5 WITHOUT security margin ---")
+print_latex_table_combinations(x_5_combinations, 5, False)
+
+x_3_combinations = [
+ [1536, 2, 128], [1536, 4, 128], [1536, 6, 128], [1536, 8, 128], [1536, 16, 128],
+ [1536, 2, 256], [1536, 4, 256], [1536, 6, 256], [1536, 8, 256], [1536, 16, 256]
+]
+
+# With security margin
+print("--- Table x^3 WITH security margin ---")
+print_latex_table_combinations(x_3_combinations, 3, True)
+
+# Without security margin
+print("--- Table x^3 WITHOUT security margin ---")
+print_latex_table_combinations(x_3_combinations, 3, False)
+
+x_inv_combinations = [
+ [1536, 2, 128], [1536, 4, 128], [1536, 6, 128], [1536, 8, 128], [1536, 16, 128],
+ [1536, 2, 256], [1536, 4, 256], [1536, 6, 256], [1536, 8, 256], [1536, 16, 256]
+]
+
+# With security margin
+print("--- Table x^(-1) WITH security margin ---")
+print_latex_table_combinations(x_inv_combinations, -1, True)
+
+# Without security margin
+print("--- Table x^(-1) WITHOUT security margin ---")
+print_latex_table_combinations(x_inv_combinations, -1, False)
diff --git a/[Tho09]poseidon-hash/reference/hadeshash/code/generate_parameters_grain.sage b/[Tho09]poseidon-hash/reference/hadeshash/code/generate_parameters_grain.sage
new file mode 100644
index 0000000..27d72f6
--- /dev/null
+++ b/[Tho09]poseidon-hash/reference/hadeshash/code/generate_parameters_grain.sage
@@ -0,0 +1,373 @@
+# Remark: This script contains functionality for GF(2^n), but currently works only over GF(p)! A few small adaptations are needed for GF(2^n).
+from sage.rings.polynomial.polynomial_gf2x import GF2X_BuildIrred_list
+
+# Note that R_P is increased to the closest multiple of t
+# GF(p), alpha=3, N = 1536, n = 64, t = 24, R_F = 8, R_P = 42: sage generate_parameters_grain.sage 1 0 64 24 8 42 0xfffffffffffffeff
+# GF(p), alpha=5, N = 1524, n = 254, t = 6, R_F = 8, R_P = 60: sage generate_parameters_grain.sage 1 0 254 6 8 60 0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001
+# GF(p), x^(-1), N = 1518, n = 253, t = 6, R_F = 8, R_P = 60: sage generate_parameters_grain.sage 1 1 253 6 8 60 0x1000000000000000000000000000000014def9dea2f79cd65812631a5cf5d3ed
+
+# GF(p), alpha=5, N = 765, n = 255, t = 3, R_F = 8, R_P = 57: sage generate_parameters_grain.sage 1 0 255 3 8 57 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
+# GF(p), alpha=5, N = 1275, n = 255, t = 5, R_F = 8, R_P = 60: sage generate_parameters_grain.sage 1 0 255 5 8 60 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
+# GF(p), alpha=5, N = 762, n = 254, t = 3, R_F = 8, R_P = 57: sage generate_parameters_grain.sage 1 0 254 3 8 57 0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001
+# GF(p), alpha=5, N = 1270, n = 254, t = 5, R_F = 8, R_P = 60: sage generate_parameters_grain.sage 1 0 254 5 8 60 0x30644e72e131a029b85045b68181585d2833e84879b9709143e1f593f0000001
+
+# GF(2^n), alpha=5, N = 768, n = 256, t = 3, R_F = 8, R_P = 57: sage generate_parameters_grain.sage 0 0 256 3 8 57 0x0
+
+if len(sys.argv) < 7:
+ print("Usage: