Module utils 

Equivalent matrix decomposition for efficient partial rounds.

Overview

This module implements the sparse matrix optimization described in Appendix B of the Poseidon1 paper (Grassi et al., USENIX Security 2021). It transforms the RP partial rounds from their textbook form (dense MDS multiply per round, O(t^2) each) into an equivalent form using sparse matrices (O(t) each).

Background: Textbook Partial Rounds

In the unoptimized Poseidon1, each partial round applies:

  state <- M * SBox(state + rc)

where M is the dense txt MDS matrix, SBox applies x^D only to state[0], and rc is a full t-vector of round constants. The cost per round is dominated by the dense matrix multiply: O(t^2).

The Sparse Factorization (Poseidon1 Paper, Appendix B, Eq. 5)

The key insight is that M can be factored as:

  M = M' * M''

where:

  M' = ┌───┬───┐       M'' = ┌─────────┬─────┐
       │ 1 │ 0 │             │ M[0][0] │  v  │
       ├───┼───┤             ├─────────┼─────┤
       │ 0 │ M̂ │             │   ŵ     │  I  │
       └───┴───┘             └─────────┴─────┘

Here M̂ is the (t-1)x(t-1) submatrix M[1..t, 1..t], v is the first row of M (excluding M[0][0]), ŵ = M̂^{-1} * w where w is the first column of M (excluding M[0][0]), and I is the (t-1)x(t-1) identity.

Since the partial S-box only touches state[0], the dense M’ factor can be “absorbed” into the next round’s M, leaving only the sparse M’’ to be applied per round via an O(t) matrix-vector product.

After iterating this factorization across all RP rounds, we obtain:

• One dense transition matrix m_i (applied once before the loop)
• RP sparse matrices S_r, each parameterized by vectors v_r and ŵ_r of length t-1

Round Constant Compression

In parallel, the round constants are compressed via backward substitution. Since M^{-1} is linear, we can “push” round constants backward through the inverse matrix, accumulating them into the first partial round. After this transformation:

• The first partial round uses a full t-vector of constants.
• Each subsequent partial round uses a single scalar constant (for state[0]).
• The last partial round has no additive constant at all.

Implementation Note

This implementation follows the HorizenLabs reference (plain_implementations/src/poseidon/poseidon_params.rs) which works on the transposed MDS matrix internally and reverses the sparse matrix ordering before returning.

Functions

circulant_to_dense

Expand a circulant matrix from its first column into a dense NxN matrix.

forward_constant_substitution

Forward constant substitution for the textbook partial round path.