Designing Efficient Dyadic Operations for Cryptographic Applications

Our Contribution

We analyze two separate aspects of dyadic operations. First, we present three different algorithms for that are aimed specifically at computing multiplication of dyadic matrices. These are, respectively, a “standard” approach that makes use of dyadic signatures, a specialized Karatsuba-like algorithm, and a procedure based on the Fast Walsh-Hadamard Transform (FWHT) [7], also called dyadic convolution. We analyze the performance of all three methods and report our timings.

As a second contribution, we describe a procedure that applies the LUP decomposition [5] to the dyadic case. The method effectively factors every quasi-dyadic matrix into a product of two triangular matrices and a permutation matrix. This leads to the possibility of a very efficient algorithm for computing the inverse of a matrix, which is particularly useful in code-based cryptography, for instance for computing the systematic form of a parity-check (or generator) matrix. According to our measurements, this improved inversion procedure is extremely fast, and provides a very large speedup during DAGS Key Generation.

Organization of the Paper

This paper is organized as follows. We start with some preliminary definitions in Section 2. We then present our main contributions: the various multiplication techniques are described in Section 3 and the improved inversion algorithm is presented in Section 4. We conclude by showing the results obtained when applying our techniques to DAGS; this is done in Section 5.

3.1 Standard Multiplication

The first algorithm we analyze is described in Algorithm 1; we refer to it as the standard multiplication. The element of C in position (i, j) is obtained as the multiplication between the i-th row of A and the j-th column of B. Since dyadic matrices are symmetric, this is equivalent to the inner product between the i-th row of A and the j-th one of B. The signature c (i.e., the first row of C) is obtained by inner products involving only a (i.e., the first row of A). Thus, we can just construct the rows of B, by permutations of the elements in b, and then compute the inner products.

The complexity of the algorithm is due to two different types of operations:

1. In order to construct the rows of B, we need the indexes of the corresponding permutations. Each index is computed as the modulo 2 sum of two binary vectors of length r, so can be obtained with a complexity of r binary operations. Thus, considering that we need to repeat this operation for 2^r − 1 rows (for the first one, no permutation is needed), the complexity of this procedure can be estimated as r ⋅ 2^r ⋅ ( 2^r − 1 ).

2. Each element of c is obtained as the inner product between two vectors of 2^r elements, assuming values in 𝔽. This operation requires 2^r multiplication and 2^r − 1 sums in 𝔽. If we denote as C_mult and C_sum the costs of, respectively, a multiplication and a sum in 𝔽, the total number of binary operations needed to compute 2^r inner products can be estimated as 2^2r ⋅ C_mult + (2^2r − 2^r)⋅ C_sum.

The complexity of a standard multiplication between two dyadic signatures can be estimated as:

3.2 Dyadic Convolution

Of particular interest to us is the case where 𝓡ⁿ is actually a field 𝔽. Dyadic matrices over 𝔽 form a commutative subring Δ(𝔽ⁿ) ⊂ 𝔽^n×n, and this property gives rise to efficient arithmetic algorithms to compute the dyadic convolution. In particular, we here consider the fast Walsh-Hadamard transform (FWHT), which is well known [7] but seldom found in a cryptographic context. We describe it here for ease of reference. We firstly recall the FWHT for the case of a field 𝔽 such that char(𝔽) ≠ 2, and then describe how this technique can be generalized to consider also the case of char(𝔽) = 2 (which, again, is the one we are interested in).

One can show by straightforward induction that Hr2=2rIr and hence Hr−1=2−rHr, which can also be expressed recursively as

Lemma 3.1 establishes that Sylvester-Hadamard matrices diagonalize all dyadic matrices. In particular, the factors in a product of dyadic matrices are thus simultaneously diagonalized, suggesting an efficient way to carry out the matrix multiplication, namely, computing Hr−1(MN)Hr=(Hr−1MHr)(Hr−1rHr) given the diagonal forms Hr−1MHr  and  Hr−1rHr of two dyadic matrices M and N requires only n multiplications of the diagonal elements.

In fact, it is not necessary to compute Hr−1MHr in full to obtain the diagonal form of M, as indicated by the following result:

The structure of Sylvester-Hadamard matrices leads to an efficient algorithm to compute aH_r for a ∈ 𝔽ⁿ, which is known as the fast Walsh-Hadamard transform. Let [a₀, a₁] be the two halves of a. Thus

This recursive algorithm, which can be easily written in purely sequential fashion (Algorithm 2), has complexity Θ(n log n), specifically, rn additions or subtractions in 𝔽. It is therefore somewhat more efficient than the fast Fourier transform, which involves multiplications by n-th roots of unity, when they are available at all (otherwise working in extension fields is unavoidable, and more expensive).

The product of two dyadic matrices Δ(a) and Δ(b), or equivalently the dyadic convolution a△b, can thus be efficiently computed as described in Algorithm 3. The total cost is 3rn additions or subtractions and 2n multiplications (half of these by the constant 2^−r = 1/n) in 𝔽, with an overall complexity Θ(n log n). Notice that this is also the complexity of computing det Δ(a).

The fast Walsh-Hadamard transform itself is not immediately possible on fields of characteristic 2, since it depends on Sylvester-Hadamard matrices which must contain a primitive square root of unity. Yet the FWHT algorithm can be lifted to characteristic 0, namely, from 𝔽₂ = ℤ/2ℤ to ℤ, or more generally from 𝔽_2^N = (ℤ/2ℤ)[x]/P(x) (for some irreducible P(x) of degree N) to ℤ[x]. Algorithm 3 can then be applied, and its output mapped back to the relevant binary field by modular reduction. This incurs a space expansion by a logarithmic factor, though. Each bit from 𝔽₂ is mapped to intermediate values that can occupy as much as 3r + 1 bits; correspondingly, each element from 𝔽_2^N is mapped to intermediate values that can occupy as much as (3r + 1)N bits. Thus the component-wise multiplication in Algorithm 3 becomes more complicated to implement for large N. However, the method remains very efficient for the binary case as long as each expanded integer component fits a computer word. For a typical word size of 32 bits and each binary component being expanded by a factor of 3r + 1, this means that blocks as large as 1024 × 1024 can be tackled efficiently. On more restricted platforms where the maximum available word size is 16 bits, dyadic blocks of size 32 × 32 can still be handled with relative ease.

3.3 Karatsuba Multiplication

In this section we propose a method which is inspired by Karatsuba’s algorithm for the multiplication of two integers [9]. Let us denote by a₀ and a₁, respectively, the first and second halves of a, i.e.:

The same notation is used for b₀ and b₁ and c₀ and c₁, corresponding to the halves of B and C. Some straightforward computations show that the following relations hold:

The iterative application of equation (5) allows to compute multiplications between dyadic matrices of any size. Let us denote as Cmul(2z)  and  Csum(2z) the complexities of a multiplication and a sum between two signatures of length 2^z. For the sum of two dyadic signatures of size 2^z we have:

where C_sum again denotes the complexity of a sum in the finite field.

The complexity of this algorithm can thus be estimated as:

Taking into account the well known sum of a geometric series, we have:

Considering this result, equation (7) leads to:

4.1 Dyadic Matrices

The inverse of a dyadic matrix (i.e. d₁ = d₂ = 1) can be efficiently computed, using only the signature, as described by the following Lemma.

As we mentioned before, the above Lemma yields a very simple way for computing the inverse of a dyadic matrix: given a signature a, we just need to compute its diagonalized form as aH_r, compute the reciprocals of its elements and put it in a vector b̃. Finally, the inverse of Δ(a) can be obtained as 12rb~Hr. This property also leads to a very simple way to check the singularity of Δ(a): if its diagonalized form contains some null elements, then it is singular.

We now focus on the case of dyadic matrices over a field 𝔽 with characteristic 2. One can show by induction that in such a case a dyadic matrix Δ(a) of dimension n satisfies Δ(a)² = (∑_ia_i)²I, and hence its inverse, when it exists, is Δ(a)⁻¹ = (∑_ia_i)⁻²Δ(a), which can be computed in O(n) steps since it is entirely determined by its first row. It is equally clear that det Δ(a) = (∑_ia_i)ⁿ, which can be computed with the same complexity (notice that raising to the power of n = 2^r only involves r squarings). Basically, verifying whether a dyadic matrix has full rank or not can be easily done by checking whether the sum of the elements of the signature equals 0.

4.2 Quasi-Dyadic Matrices

Consider a quasi-dyadic matrix M. Since the matrix has to be square, we have d₁ = d₂ = d, and the matrix has dimension dn × dn. Such a matrix can be compactly represented just by the signatures of the dyadic blocks. To simplify notation, we can denote the signature of the dyadic-block in position (i, j) as m̂_i,j, and store all such vectors in a matrix M̂ ∈ 𝓡^d×dn:

We focus again on the special case of quasi-dyadic matrices over a field 𝔽 with characteristic 2.

The LUP decomposition is a method which factorizes a matrix M as LUP, where L and U are lower triangular and upper triangular matrices, respectively, and P is a permutation.

Exploiting this factorization, the inverse of M can thus be expressed as:

The advantage of this method is that the inverses appearing in (10) can be easily computed, because of their particular structures. In fact, the inverse of an upper (lower) triangular matrix is obtained via a simple backward (forward) substitution procedure, while the inverse of P is its transpose.

In some cases, applying a block-wise LUP decomposition might lead to some complexity reduction; for instance, see [13] for the inversion of a sparse matrix. Here, we consider the case of a quasi-dyadic matrix; the corresponding procedure is shown in Algorithm 4.

Our proposed procedure consists in using a block decomposition, which works directly on the signatures, in order to exploit the simple and efficient algebra of dyadic matrices. The operations in Algorithm 4 only refer to the signatures in M̂: for instance, the expression m̂_{i, j} m̂_{i, l} means the product between the dyadics having as signatures m̂_{i, j} and m̂_{i, l}. This choice may result in some abuse of notation, but is useful to emphasize the fact that, as we have explained in the previous sections, operations with dyadics can be efficiently computed just by taking into account their signatures. It can be easily shown that, for a quasi-dyadic matrix, its factors L, U and P are in quasi-dyadic form as well: as we have done for the matrix M, we refer to their compact representations as L̂, Û and P̂, respectively.

The algorithm takes as input a matrix M̂, as in (9), and computes its LUP factorization; outputs of the algorithm are the modified matrix M̂, having as elements the ones of its factors L̂ and Û, and the permutation P̂. As in (9), we denote as m̂_i,j the signature in position (i, j) in the output matrix M̂. The matrices L̂ and Û can then be expressed as:

where 1̂ and 0̂ denote, respectively, the signature of the identity matrix and the one of the null matrix (i.e. the length-k vectors [1, 0, ⋯, 0] and [0, 0, ⋯, 0]).

The matrix P̂ is represented through a length-d vector [p₀, p₁, ⋯, p_d−1], containing a permutation of the integers [0, 1, ⋯, d − 1]; the rows of M̂ get permuted according to the elements of P̂. In particular, the elements of P̂ are obtained through a block pivoting procedure, which is described in Algorithm 5.

This function takes as input M̂, P̂ and an integer j, and searches for a pivot (i.e., a non singular signature) in the j-th column of M̂, starting from m̂_j,j, and places it in position (j, j). As the procedure goes on, every time a singular signature is tested, the rows of M̂ get permuted; the elements of P̂ are accordingly modified. If the j-th column contains all singular blocks, this means that the matrix M̂ is singular; in such a case, this event is notified by setting u = 0.

We point out that, for the matrices we are considering, we expect Algorithm 4 to be particularly efficient. First of all, as we have already said, this is due to the possibility of efficiently performing operations involving dyadic matrices; in addition, the dyadic structure should also speed-up the pivoting procedure. In fact, we can consider a signature in M̂ as a collection of k random elements picked from GF(2^N): thus, their sum can be assumed to be a random variable with uniform distribution among the elements of the field GF(2^N). So, the probability of it being equal to 0, which corresponds to the probability of the corresponding signature to be singular, equals 2^−m. This probability gets lower as m increases: this fact means that the expected number of operations performed by Algorithm 5 should be particularly low. Basically, most of the times the function will just compute the sum of the elements in m̂_j,j and verify whether it is null or not.

Once the factorization of M̂ has been obtained, we just need to perform the computation of M⁻¹ through (10). Since the inverse of a triangular matrix maintains the original triangular structure, the computation of the inverses L̂⁻¹ and Û⁻¹ can be efficiently performed. A possible way for computing these matrices is to store the elements of both matrices in just one output matrix T̂. We do this in Algorithm 6.

The matrix Î_d is the compact representation of a dn × dn identity matrix, and so is composed of signatures δ_i,j1̂, where δ_i,j denotes the Kronecker delta.

If we denote as t̂_i,j the signature in position (i, j), we have:

Implementation Details

We developed a code in “C'' to implement our procedures. In all cases, we use no optimizations apart from the optimization from the GCC compiler (“-O3”). The GCC version used was 7.3.1 20180406, the code was compiled for the processor Intel(R) Core(TM) i5-5300U CPU @ 2.30GHz with 16GB of memory and operating system Arch linux version 2018.05.01 with Kernel 4.16.5. We ran 100 times each piece of code and computed the average of all measurements; to obtain the number of cycles, we used the file “cpucycles.h” from supercop^[1].

Fast Multiplication

To compare our methods, we fix a dyadic order r and measure the cost of a multiplication of two matrices of size n = 2^r. Relevant dyadic orders for DAGS are for instance r = 4 and r = 5. We do this over different fields to highlight the difference in performance when changing fields: we tested 𝔽_2⁵ and 𝔽_2⁶ which are the fields currently used by DAGS.

Efficient Inversion

We report here the results of the improved inversion procedure (Algorithms 4, 5 and 6). We compared our procedure with the equivalent portion of the DAGS implementation that we extrapolated from the publicly available source code [2]. In particular, we measured the piece of code that begins with the creation of the Cauchy matrix and ends with the generation of the systematic matrix. Table 2 shows the comparison, measured in cpu cycles.