add: referece impl, README, matrices

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"

[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 }

[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/

[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.

[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.
+
+<br>
+
+[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.

[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)