Better path-finding algorithms in LPS Ramanujan graphs

1 Introduction

LPS graphs were introduced by Lubotzky, Philipps and Sarnak [9]. LPS graphs are Ramanujan graphs, namely, they are graphs with optimal expansion properties, and they have many applications in mathematics and computer science [8]. In [2], Charles, Goren and Lauter introduced a cryptographic hash function whose security relied on the presumed computational intractability of computing short paths between arbitrary pairs of vertices in the graph. Due to the algebraic nature of LPS graphs this problem can also be stated in group-theoretic terms: given the group G=PSL⁢(2,𝔽p) and particular set of generators, and given an element of G, find a short factorization of this element as a product of the generators. This problem was later solved in [16, 11], leading to a full cryptanalysis of the cryptographic construction.

The quality of a path-finding algorithm is measured by its running time, but also by the lengths of the paths returned. In this paper we are interested in minimizing this length while keeping a polynomial run time complexity. While the algorithm in [11] returns paths that are a priori optimal up to a constant factor 16/3, we improve this algorithm and decrease this factor down to 7/3 under several plausible heuristic conditions. Similarly to [11], we first solve the case of diagonal matrices, and then we reduce the general case to this case. Our algorithm improves both the diagonal case and the reduction, and for the diagonal case, it returns factorizations with lengths very close to the expected minimal.

LPS graphs have recently been suggested for quantum circuit design. In this setting LPS generators correspond to elementary quantum gates, and they are combined to approximate arbitrary quantum gates on a single qbit. Interestingly, an algorithm similar to our diagonal algorithm has been independently developed in that context [1, 12]. It was later transposed to our setting [13], resulting in an algorithm very similar to our algorithm in the diagonal case. It seems very likely that our other improvement (reducing from the general case to the diagonal case) can also be adapted to the problem of quantum circuit design; this will be the purpose of a future paper.

2 Path-finding algorithms in LPS Ramanujan graphs

In this section we recall the definition of LPS graphs and we briefly describe previous path-finding algorithms in these graphs. In our exposition we mostly follow [16, 11], and we refer to those papers for details and proofs.

2.2 Finding paths in LPS graphs

In this paper, we are interested in algorithms to solve Problem 2.2 or, equivalently, to compute short paths between any two vertices in LPS graphs.

Let all notations be as above, and let m be a matrix in G that we want to write as a short product of elements from the set S.

The algorithm of Petit, Lauter and Quisquater [11] first decomposes the matrix m as a product

(2.1)m=λ⋅D1⋅s⋅D2⋅s⋅D3⋅s⋅D4,

where λ∈𝔽p*, the factors Di are diagonal matrices with a non-zero square determinant, and s is a particular (arbitrary) generator in the set S. As the equation is over the projective special linear group, each diagonal matrix can be normalized as Di=(100αi). Equation (2.1) then amounts to a small polynomial system with four equations and five variables. The algorithm given in [11] to solve this system picks random solutions until all αi are square, and it is heuristically expected to need 16 trials on average.

Petit, Lauter and Quisquater [11] also provide an algorithm to factor any diagonal matrix with square determinant in a short product of the generators S. This algorithm extends a previous algorithm from Tillich and Zémor [16] for the identity matrix. First the diagonal matrix is lifted into an element of Ω, then a factorization of this element as a product of the elements of Σ is computed. The factorization of the input diagonal matrix in the elements of S=ϕ⁢(Σ) follows by the group homomorphism ϕ. We now give some details on the first and second step.

Let m=(A+B⁢ι00A-B⁢ι)∈G be a diagonal matrix with det⁡m=A2+B2 a square. The lifting step consists in finding e∈ℕ and λ,w,x,y,z∈ℤ with

(2.2){(A⁢λ+w⁢p)2+(B⁢λ+x⁢p)2+4⁢p2⁢(y2+z2)=ℓe,A⁢λ+w⁢p=1mod2,B⁢λ+x⁢p=0mod2.

After fixing e large enough, Petit, Lauter and Quisquater solve the norm equation modulo p; as A2+B2 is a square, this gives two possible values for λ and one is picked randomly. Next, the norm equation is considered modulo p2; this gives a bilinear equation in w and x, and a random solution is selected. At this point, the norm equation is considered over the integers, and it gives

4⁢(y2+z2)=n,

where

(2.3)n:=(ℓe-(A⁢λ+w⁢p)2-(B⁢λ+x⁢p)2)/p2∈ℤ

for λ,x,z chosen as before. This equation has a solution whenever n/4 is an integer, and all prime factors of n congruent to 3 modulo 4 appear an even number of times in the factorization of n. To avoid a costly factorization step, Petit, Lauter and Quisquater suggest to pick random solutions for (w,x) until n is 4 times a prime congruent to 1 modulo 4. Suitable values for (y,z) are then computed with Cornacchia’s algorithm [4]. Taking e larger than log⁡8⁢p4≈4⁢log⁡p ensures that n is positive in this algorithm.

Once the integers e,λ,w,x,y,z are found, they give a matrix

m~=((A+w⁢p)+(B+x⁢p)⁢i2⁢y⁢p+2⁢z⁢p⁢i-2⁢y⁢p+2⁢z⁢p⁢i(A+w⁢p)-(B+x⁢p)⁢i)∈Ω,

reducing to m modulo p. Remember that by Lemma 2.1 factorization is essentially unique in Ω. Tillich and Zémor showed in [16] how to recover all the factors successively: for any s∈Σ, we have m~=m~′⁢s with m~′∈Ω if and only if m~⁢s-1∈Ω.

Algorithm 1 (Factorization in Ω.).

The correctness and running time of this algorithm are analyzed in the following lemma.

Lemma 2.3.

For any m~∈Ω, Algorithm 1 successively recovers the factorization given by Lemma 2.1, in a time O⁢(logℓ2⁢det⁡m~).

Proof.

Let m~∈Ω and det⁡(m~)≠1. Then there is a unique s∈Σ such that m~=m~′⁢s with m~′∈Σ, by Lemma 2.1. Therefore, there is an equally unique s′=s¯∈Σ such that s⁢s′=ℓ⁢I and m~′⁢s⁢s′ is divisible by ℓ. This way, by iterating over Σ, the algorithm always finds the matrix s if det⁡(m~)≠1, adding it to the beginning of the list ℒ, and finishes when det⁡(m~)=1. In the end, ℒ holds the complete factorization in the right order for m~ on elements of Σ.

As for the running time, one can write all coefficients of m~ in base ℓ in O⁢(logℓ⁢det⁡m~) bit operations, after which step 3 becomes trivial. The while loop is executed logℓ⁢det⁡m~ times and each iteration requires O⁢(ℓ2)=O⁢(1) multiplications of a O⁢(logℓ⁢det⁡m~) length element by a O⁢(ℓ)=O⁢(1) length element, and each of them requires O⁢(logℓ⁢det⁡m~) bit operations. ∎

The algorithm of this paper will incorporate two crucial improvements over Petit, Lauter and Quisquater’s algorithm to reduce factorization lengths from 16⁢logℓ⁡p to (7+o⁢(1))⁢logℓ⁡p. In Section 4 we will reduce the factorization lengths for diagonal matrices from 4⁢logℓ⁡p to (3+o⁢(1))⁢logℓ⁡p. In Section 5 we show how to factor arbitrary matrices with two diagonal matrices plus a third matrix with factorization length logℓ⁡p.

We will provide an alternative method to reduce the general factorization problem to the factorization of only two diagonal matrices instead of four.

3 Algorithm for short paths

Our first algorithm computes the shortest path between any two matrices when the length of this shortest path is short enough. This is useful in a cryptanalysis of the Charles–Goren–Lauter hash function assuming short messages are hashed.

The algorithm relies on the following simple observation: when the products are short, then no reduction modulo p occurs and the lifting step required in previous algorithms [11, 16] becomes trivial.

More generally, let m1,m2∈G such that there exists a product m of at most ⌊2⁢logl⁡(p/2)⌋ matrices in ϕ⁢(Σ) such that m2=m1⁢m. For every possible length e up to ⌊2⁢logl⁡(p/2)⌋, Algorithm 2 normalizes the matrix to ensure its determinant is ℓe. Then it writes the result in the form (a+b⁢ιc+d⁢ι-c+d⁢ιa-b⁢ι) and it checks whether a,b,c,d are small enough to be lifted directly to Ω. When the correct factorization length is identified the algorithm then proceeds as in the Tillich–Zémor algorithm to compute all the factors.

Algorithm 2 (Short path algorithm.).

(3.1)

We remark that there is at most one value of e≤⌊logl⁡(p)⌋ so that the algorithm succeeds: indeed, the girth of LPS graphs is larger than 2⁢logℓ⁡p when ℓ is a quadratic residue modulo p, see [9]. Lemma 3.2 below gives the complexity of this algorithm.

4 Improved algorithm for diagonal matrices

We now turn to diagonal matrices, and improve the algorithm of Petit, Lauter and Quisquater [11] to return factorizations of length (3+o⁢(1))⁢log⁡p instead of 4⁢log⁡p. A similar algorithm was independently suggested in [13].

Let m=(A+B⁢ι00A-B⁢ι)∈G be a diagonal matrix with determinant (A2+B2) a quadratic residue modulo p. If B=0 then m is the identity matrix and it admits a trivial factorization. So we now assume B≠0.

Our algorithm follows the lines of Petit, Lauter and Quisquater’s algorithm sketched in Section 2.2. Recall that in this algorithm we need e≈4⁢log⁡p to ensure that the value n defined by equation (2.3) is positive. In their algorithm (w,x) is chosen as a random solution to the norm equation modulo p2 such that n is of a form that makes Cornacchia’s algorithm efficient. Our main idea for the algorithm of this section is to choose (w,x) not randomly but

1. to minimize the value (A⁢λ+w⁢p)2+(B⁢λ+x⁢p)2,

2. while keeping n=(ℓe-(A⁢λ+w⁢p)2-(B⁢λ+x⁢p)2)/p2 of a form that makes Cornacchia’s algorithm efficient.

We now show how to model the minimization problem as a variant of a closest vector problem in a two-dimensional lattice, which can then be solved using classical techniques such as lattice reduction, the Euclidean algorithm, or equivalently continued fractions.

Let e be a positive integer (to be fixed later) and let λ∈ℤ such that λ2⁢(A2+B2)=ℓemodp. Let w0=1 if A⁢λ is even and w0=0 otherwise. Similarly, let x0=0 if B⁢λ is even and x0=1 otherwise. For any couple (x,y) that is part of a solution to system (2.2), there corresponds a couple of integers (w′,x′) such that w=w0+2⁢w′ and x=x0+2⁢x′. Moreover, we have

or equivalently

where we defined

We observe that the solutions to equation (4.1) are obtained by translation of a certain lattice.

We then re-write the norm we want to minimize as the distance between a particular element in the Euclidean plane and a lattice element:

where v:=(w′,x′)-(0,β) is a lattice point and t:=(-A⁢λ+w0⁢p2⁢p,-B⁢λ+x0⁢p2⁢p-β). Minimizing this norm is an instance of the closest vector problem for which standard algorithmic solutions exist [10]. We adapt these solutions to incorporate the additional requirement that n:=ℓe-4⁢p2⁢∥v-t∥2 is of a form that makes Cornacchia’s algorithm efficient.

Our strategy is to first compute a minimal basis for the lattice, namely, two vectors that generate the lattice and have minimal norms. This can be done using Gauss reduction, or equivalently the Euclidean algorithm. We then write t as a linear combination of the two short vectors, and we round the coefficients to obtain a lattice vector v that is close to t. Finally, we add small lattice vectors to v and compute the corresponding value of n, until it is of a suitable form. Once we find n of a suitable form, we compute the values of w and x and we apply Cornacchia’s algorithm to compute y and z. In our algorithm we start with an a priori sufficient value e=3⁢logℓ⁡p+log⁡logℓ⁡p, determined by our complexity analysis below, and then we increase e during the course of the algorithm if no n value of a suitable form can be found. This strategy is described in Algorithm 3.

Algorithm 3 (Diagonal matrix algorithm.).

(4.2)

While the algorithm of Petit, Lauter and Quisquater returns factorizations of length 4⁢logℓ⁡p, we expect our algorithm to return factorizations of length (3+o⁢(1))⁢logℓ⁡p. Our analysis in the proof of the following lemma relies on the probability that the number n generated by the algorithm is a prime. We approximate this by the probability that random numbers of the same size are prime. The analysis also relies on the so-called “Gaussian heuristic”: the two vectors in a minimal basis are expected to be of norm roughly det⁡ℒ. This will be true for most inputs, but for exceptional inputs one of these vectors may be significantly larger than the other one, in which case the best achievable bound is (4+o⁢(1))⁢logℓ⁡p. We refer to [13] for a more thorough analysis on a very similar algorithm, and to [14] for some progress on the heuristic assumptions involved in this analysis.

Note that paths of length about 3⁢logℓ⁡p are expected to be necessary and sufficient for random matrices in G (since G is of size roughly p3/2). We therefore expect that the factorizations we compute are close to optimal ones unless diagonal matrices are particularly close to the identity in LPS graphs, which there is no reason to believe a priori.

5 Improved algorithm in the general case

We now consider the case of arbitrary matrices in G. In [11] Petit, Lauter and Quisquater showed how to factor general matrices with four diagonal matrices, for a total factorization length of 16⁢logℓ⁡p. Here we factor general matrices with just two diagonal matrices and one matrix in Ω, for a total expected length of (7+o⁢(1))⁢log⁡p.

Let J:=(01-10). We will use the following lemma.

These properties are of course preserved by reduction modulo p through the homomorphism ϕ. For any diagonal matrix D=(a+b⁢i00a-b⁢i), we define its conjugate D¯=(a-b⁢i00a+b⁢i), which is equal to D-1 up to a scalar factor.

Let m∈G. Our strategy is to find a matrix m2∈Ω and two diagonal matrices m1,m3∈G such that

and then to use the factorization algorithms for Ω elements and diagonal matrices to compute a factorization of m. (Note that since the inverse of any element in Σ is also in Σ a factorization of m1, it will also give a factorization of m1-1.)

Following Lemma 5.1, we write m=D+E⁢J, with D,E diagonal matrices in G, and similarly we write m2=X2+Y2⁢J, with X2,Y2 matrices of the form (w+x⁢i00w-x⁢i). We also define X1:=m1 and X3:=m3 to highlight that they are diagonal matrices. Equation (5.1) and the requirement that m2 is in Ω then lead to the following system:

Taking the determinants of both equations and dividing one by another, we deduce

Since we require m2∈Ω, we also have

for some positive integer e. These last two equations together fix the determinants of X2 and Y2 up to a multiple of p, for each value of e:

Given the determinant d of a matrix of the form (w+x⁢i00w-x⁢i), we can recover suitable values for w and x by solving the Diophantine equation w2+x2=d; this can be done efficiently as long as a factorization of d is known and all prime factors congruent to 3 modulo 4 appear an even number of times.

In our algorithm (Algorithm 4 below), we will fix e a priori to ⌈logℓ⁡p⌉, then try every possible value of det⁡X2<ℓe consistent with its value modulo p, and deduce det⁡Y2=ℓe-det⁡X2, until both det⁡X2 and det⁡Y2 are a power of 2 times a prime congruent to 1 modulo 4. As in previous algorithms, this restriction may lead to suboptimal factorization lengths but it allows to avoid a costly factorization step. Once the condition is satisfied, we apply Cornacchia’s algorithm on both det⁡X2 and det⁡Y2 to obtain X2 and Y2, and we deduce the matrix m2. If no suitable values for det⁡X2 and det⁡Y2 are found, we simply increase e and try again.

Let us now assume that X2 and Y2 are fixed as above. From system (5.2), we deduce

As diagonal matrices commute, each of these two equations can be rewritten in the form

for i=1,3 and some known diagonal matrices Ki. Moreover, by the way we constructed X2 and Y2, we have det⁡Ki=1.

Algorithm 4 (General factorization algorithm.).

(5.5)

If Ki is the identity, then we can take the identity for Xi. Otherwise, we write Xi=(w+x⁢ι00w-x⁢ι) and Ki=(w0+x0⁢ι00w0-x0⁢ι) for unknown w,x and known w0,x0, and we deduce the following system:

Moreover, since det⁡(w0-1-x0x0w0+1)=w02+x02-1=det⁡Ki-1=0, this system is equivalent to the single equation

where w0-1 and x0 are both non-zero, as otherwise Ki would be the identity. We can therefore choose w=x0 and x=w0-1, in other words,

In order to call our diagonal matrix factorization algorithm on X1 and X3, we also need their determinants to be quadratic residues modulo p. By assumption, det⁡m and det⁡m2=ℓe are quadratic residues, so either both det⁡X1 and det⁡X3 are quadratic residues or none of them is. In our algorithm when they are non-quadratic residues we just search for new values of det⁡X2 and det⁡Y2.

In addition to the heuristic involved in Lemma 4.2, Algorithm 4 requires other heuristics related to the probability to find suitable determinants for X2 and Y2 that are not too large, and the probability that both diagonal matrices have square determinant. Under these heuristics, the runtime complexity and the output of this algorithm are given by the following lemma.

Note that paths of length about 3⁢logℓ⁡p are expected to be necessary and sufficient for random matrices in G. The factorizations produced by our algorithm are therefore expected to be a factor 7/3 larger than optimal ones, instead of 16/3 larger than optimal ones for the algorithm of Petit, Lauter and Quisquater [11]. Moreover, the best diameter bound on LPS graphs is 6⁢log⁡p, so in that sense we can expect that our factorizations are only 16% larger than optimal ones in the worst case.

6 Experimental validation

Figure 1

Factorization lengths of our algorithm for general matrices.

Because of the heuristic nature of our analysis we have experimentally checked our results using the computer algebra package Magma [6]. In our experiments we fixed ℓ=5 and we generated various p such that p=1mod4 and (ℓp)=1. For each such set of parameters, we generated various matrices of square determinant, and we ran Algorithm 4 on them. The resulting factorization lengths are displayed in Figure 1 as functions of the parameter logℓ⁡p. The best near square approximation line for the data points has a slope equal to 7.03, confirming our analysis that the lengths are (7+o⁢(1))⁢logℓ⁡p.

7 Conclusion

We presented a new algorithm that computes short paths between arbitrary pair of vertices in LPS Ramanujan graphs. Subject to several plausible heuristics, the algorithm is polynomial time and for LPS graphs with parameters p and ℓ, it returns paths of length (7+o⁢(1))⁢logℓ⁡p.