This repository contains a golang implementation of the PSLQ Algorithm as described in the original 1992 PSLQ paper, with modifications geared to cryptanalysis. PSLQ is an algorithm that can find a small but not-all-zero integer-only solution m₁,m₂,...,m_n of the equation

x₁m₁+x₂m₂+...+x_nm_n=0

where the x_i are real numbers.

The purpose of this repository is to demonstrate the usefulness of PSLQ in cryptanalysis, for example against a type of cryptography based on something called LWE.

"LWE" stands for Learning with Errors, and the "LWE problem" refers to a problem that a cryptanalyst must solve in order to break the class of crypto-systems bearing the same name. The LWE problem can be reduced to solving Am = 0 where

• A is a pxn public matrix with entries from ℤ_q.
• m is a short n-long vector with integer entries.

Right-multiplying A by b^{^} = [1, b, b², ..., b^p-1] transforms A into a 1xn matrix, appropriate for input to PSLQ. Solutions of <b^{^}A, m> = 0 are likely solutions of Am = 0 for sufficiently large b.

There are lots of issues to work out, including the fact that the entries of A are elements of ℤ_q, not ℤ. There are also answers to these issues, some better than others. For example, extending A with a certain p x p matrix containing 0s and non-zero integers of the form b^kq, gets around the ℤ_q, vs. ℤ issue. Anyway, the cryptanalyst only needs to win now and then to succeed.

Rather than run down all the details of the transformation of the LWE problem into a PSLQ problem, the point here is to raise the possibility that PSLQ is a powerful, if under-appreciated, tool for cryptanalysis. To be a credible tool, PSLQ should be able to solve <x,m> = 0 with the smallest possible m when x has hundreds of large integer entries. This repository begins to tackle this problem.

Without some changes, PSLQ is not a useful tool for cryptanalysis. PSLQ was originally designed to handle real (or complex or quaternion), non-integer input. Given integer input, PSLQ as defined in the original 1992 paper quickly finds a bad (high-norm) solution, m, and terminates.

Read on to see in detail how this problem can be fixed. For now, let it suffice to say that the modification of PSLQ proposed here delays the termination of the algorithm until the solution is optimal.

PSLQ, as originally defined, is both a framework and a strategy. The framework consists of a matrix equation and a set of allowed operations on this equation that change its components until a solution of <x,m> = 0 is found. The strategy specifies what operations to perform, and when to perform them.

The framework cannot change; the strategy can. In fact, the strategy has been modified in the literature, most notably in this paper, where one of the authors of the original 1992 paper proposes a strategy that can take advantage of parallel processing. The reason the framework cannot change is that its set of allowed operations is what guarantees some invariants needed to solve <x,m> = 0.

The adjective, "classic", in both this README and in the code, refers to the original strategy suggested in the 1992 paper. The classic strategy is geared towards proving bounds on the performance of PSLQ, so it can be considered a polynomial time algorithm that finds solutions within a certain factor of the optimal solution.

In this repository, the classic strategy is modified greatly. The purpose of the modifications is to delay termination and optimize the solution.

The PSLQ algorithm performs iterations that maintain a matrix equation,

xBAH_xQ = 0 (equation 1)

while updating the factors B, A and Q.

Here,

• x is the sequence of real numbers for which we are trying to find a relation. In the context of matrix equations like equation 1, consider x to be a 1 x n matrix of real numbers
• B and A are n x n matrices with integer entries and determinant 1 or -1. They are identity matrices at initialization. After each iteration, B = A^-1.
• H_x is an n x n-1 matrix with real entries, non-zero diagonal entries, and zeroes above the diagonal
• Q is a rotation matrix that keeps 0s above the diagonal of AH_xQ

PSLQ stores AH_xQ, which we will call H (as opposed to H_x) in this section. In later sections, when the value of H at a particular iteration k matters more than it does here, the notation H_k will be used. Substituting H for AH_xQ in equation 1,

xBH = 0 (equation 2)

Both x and H_x remain fixed at their initial values, whereas B, A and H are updated while PSLQ runs. Q itself is not stored directly, but is implicitly stored by storing H. H, like H_x, is non-zero on the diagonal and zero above the diagonal.

Because H is zero above the diagonal, the last column of H -- column n-1 -- is mostly 0. It contains non-zero values in at most its last two entries. It would be a shame if something happened to one of those entries, like H_n-1,n-1 becoming 0. That would leave just one non-zero entry in column n-1, namely H_n,n-1. Then the fact that (xBH)_n-1 = 0 (along with the other entries of xBH) would mean that (xB)_n = 0.

Of course this wouldn't really be a shame. According to equation 2, H_n-1,n-1 = 0 and H_n,n-1 ≠ 0 would mean that the last entry of xB -- the coefficient of H_n,n-1 in equation 2 -- has to be 0. In other words,

<x,last column B> = 0

This means that -- once H_n-1,n-1 is forced to be zero -- setting m to the last column of B makes m a solution with integer entries of <x,m> = 0.

So the idea of PSLQ is to force a 0 into the last diagonal element of H while maintaining 0s above the diagonal of H and the validity of equation 2. When that happens, PSLQ reaps a solution of <x,m> = 0 from column n of B!

The way B, A and H are modified is through row operations on A, and their inverses as column operations on B. These row operations put non-zero values above the diagonal of H, which the last factor in equation 1, Q, zeroes out with a rotation. Zeroing out entries of H above the diagonal is called "corner removal", because it removes the non-zero upper-right corner entry from a 2x2 sub-matrix of H. It is essential that Q be a rotation so the diagonal entries of H track the norms of potential solutions of <x,m> = 0.

The row operations on A are designed not only to force a zero into H_n-1,n-1, but to nudge the large diagonal entries of H towards the right, for reasons explained later. These goals have random effects on the entries of H below the diagonal. Left unchecked, these effects would cause the sub-diagonal of H to grow without bound. To offset this, the original 1992 PSLQ paper specifies the use of what it calls Hermite reduction, which aggressively minimizes the entries below the diagonal of H.

Up to the mention of Hermite reduction, everything in the description above can be considered the framework of PSLQ. Hermite reduction, in lieu of other kinds of reduction of H below the diagonal, can be considered part of the strategy for putting a zero in H_n-1,n-1 to force a solution of <x,m> = 0 into a column of B.

As noted earlier, there are many reasonable strategies. The primary example is the classic strategy from the original 1992 paper: Swap rows according to a criterion governed by a parameter, γ. When j < n-1 and γ^j |H_j,j| ≥ γⁱ |H_i,i| for all i, it is rows j and j+1 that are swapped. If no j satisfies this criterion, rows n-1 and n are swapped. After swapping, update Q to remove the corner just created (unless swapping rows n-1 and n), and perform Hermite reduction.

The classic strategy is good if there is just one independent solution of <x,m> = 0, but not good if there are many. In cryptanalytic use cases, there are many independent solutions of <x,m> = 0, so this README develops alternative strategies to find good solutions among the many.

The diagonal of H is crucial. It is the key to "accuracy", which is the term used here for the Euclidean length ("norm") of the solution m of <x,m> = 0 that PSLQ calculates. The smaller the norm, the more accurate the output.

As shown in the section "A Sharper Lower Bound on the Smallest Solution While PSLQ is Running" below, the larger the diagonal elements, the smaller the norm of the relation m that PSLQ finds (details are deferred to that section).

Of greatest importance is the last diagonal element, H_n-1,n-1. Lemma 10 in a 1999 paper analyzing PSLQ states that the norm of the solution is the value of 1/|H_n-1,n-1| at the iteration of PSLQ before H_n-1,n-1 becomes 0. So it would improve accuracy to keep |H_n-1,n-1| as large as possible while PSLQ is running.

Lemma 10 is crucial to the strategies employed in this repository, so we need to delve into some of the details of its proof. This section is a guide to that proof, not a full proof. It fills in some details, and fleshes out the assumptions behind the lemma. That way, we can be assured that adding row operations to the toolkit used in classical PSLQ does not violate the assumptions.

Lemma 10 relies on the assumptions listed below, which are not broken by any row operation with determinant 1 or -1 ("unit determinant"), or any rotation to remove zeroes above the diagonal of H. These assumptions refer to

• A matrix P_x, defined in the statement of lemma 2 in the same paper
• m, the solution the PSLQ algorithm is about to output, when rows n-1 and n are swapped.
• A, the matrix with integer entries and unit determinant that is the product of all previous row operations.
• B=A^-1, following the notation used here, though in the proof of lemma 10 "A^-1" is the notation for what we call B throughout this README.

The assumptions are:

• AP_x = TDQ^tH_x^t is a decomposition of AP_x into into the product of a lower trapezoidal matrix T with diagonal 1s, an invertible diagonal matrix D with the same diagonal as H, and an n−1×n matrix Q^tH_x^t with orthonormal rows. This, by the way, is copied from the proof of theorem 1, not lemma 10. The proof of lemma 10 implicitly refers back to that proof.
• At the point where a zero appears in H_n,n-1, the (n−1)-st column of B is m.

Based on the second assumption,

Am^t =<B^-1, column n-1 of B> = e_n−1, the (n−1)-st standard basis vector.

The second of the two "=" signs is true for any B^-1 and B. The left and right quantities are equated in the proof of lemma 10 without connecting this equality to the second assumption. So Am^t=e_n-1 can appear to be a separate assumption but it's not.

The first assumption depends only on the initial setup of PSLQ and the fact that A has integer entries and unit determinant. So no row operation with unit determinant falsifies the first assumption.

The second assumption relies only on the fact that H_n,n-1 = 0, leaving H_n-1,n-1 as the lone non-zero entry in column n-1 of H. Equation 1, xBAH_xQ = 0, tells us that in particular, if you just focus on the last coordinate of xBAH_xQ, that coordinate is zero. From this and H_n-1,n-1 being the only non-zero entry in column n-1, we get:

0 = (xBAH_xQ)_n-1 = (xBH)_n-1 = (xB)_n-1 (H_n-1,n-1) ⇒ (xB)_n-1 = 0

The rightmost equality above just says that column n-1 of B is a solution we can call "m" of <x,m> = 0.

Each iteration of PSLQ performs a pair of row operations on H -- one to tame the diagonal of H, another to reduce the entries below the diagonal. In this section, "row operation" refers to the former, designed to tame the diagonal of H. Any row operation with determinant 1 or -1 is acceptable. But the original 1992 PSLQ paper and the 1999 analysis of PSLQ consider only swaps of adjacent rows.

The reason for re-implementing PSLQ here, rather than using an existing implementation, is to replace the classic strategy in the original 1992 PSLQ paper with

• Row operations other than swaps of adjacent rows
• Swaps of adjacent rows chosen with criteria other than the ones specified in the original 1992 PSLQ paper
• Delaying termination until best possible solution becomes available.

Using these three extensions to "improve" the diagonal of H trades provable performance for empirically verified accuracy. Proofs of both accuracy and speed are presented in the original 1992 PSLQ paper and in the 1999 paper analyzing PSLQ. But the accuracy bounds these proofs promise are poor, as Table 1 below shows in its results for the classic strategy, which the 1992 (and 1999) papers propose.

Because the accuracy guarantees in the 1992 and 1999 PSLQ papers do not meet the needs of cryptographic use cases, the extensions in this repository sacrifice them, along with speed guarantees, in exchange for empirically demonstrated, albeit not mathematically proven, improvements in accuracy. Empirical results in Table 1 show that the improved accuracy comes at the cost of speed. It would not be surprising if the speed remains polynomial, but with an increase of 1 in the degree of the polynomial.

Table 1 contains a column called strategy. As noted earlier, a "strategy" is a set of rules for choosing row operations to imrpove the diagonal of H. The strategies compared in Table 1 are:

• Classic: Swap rows to improve the diagonal of H as recommended in the original PSLQ paper, until a zero-valued entry is swapped into the last diagonal element of H; terminate when that happens.
• IDASIF: "Improve diagonal after solution is found". Use the Classic strategy until a zero is about to be swapped into the last diagonal entry of H. Then instead of swapping in that zero and terminating, use row operations to improve the last three columns of the table below, until there are no row operations left to perform that improve the diagonal. IDASIF is an early version of the as-yet untested "Swap, Reduce, Solve" strategy. See the section below, "The Swap, Reduce, Solve Strategy", for details about this strategy.

It is understood that, just based on the description above, IDASIF is not a well-defined strategy. To learn the details, search improveDiagonalWhenAboutToTerminate in strategy/strategyv1.go.

Entries in Table 1 were copied from the output of the test, TestGetRImprovingDiagonal. In that test, the input to PSLQ is an n-long challenge vector (x, in the notation above) with a known small solution m₀, i.e. <x,m₀> = 0. Each entry of x is chosen from the uniform distribution on [-maxX/2,maxX/2], where maxX is chosen so that the chance of at least one random vector of norm |m₀| or less being perpendicular to x is deemed to be about 0.001. The effort that went into this probability calculation is minimal compared to an exact calculation (it's not an easy calculation). But maxX is in the ballpark of having the desired property.

The metric |largest diagonal element / last| refers to the diagonal just before PSLQ terminates. | output of PSLQ | is the norm of the solution PSLQ computes, and | m₀ | is the norm of the causal vector PSLQ should ideally find. Note that two rows with the same | m₀ | are likely to have the same input and causal relation, especially for large n.

Table 1 - Test results comparing Classic and IDASIF strategies

An advanced version of "IDASIF" from the table above can be implemented using this repository. This strategy, called "Swap, Reduce, Solve" (SRS), has two phases. Phase 1 refers to the time before the maximum number of possible solutions of <x,m>=0 is found; and phase 2 refines those solutions.

A summary of phases 1 and 2 is:

• Phase 1: H contains non-zero elements in its last row, row n. Under the right conditions, row combinations of non-zero elements of row n and diagonal can be used to reduce the absolute values of diagonal elements of H.
• Phase 2: H no longer contains non-zero elements in row n. For any column with zero in row n (which is all columns in phase 2), the smallest index, i, of a non-zero entry, is the index of a column in B with a solution of <x,m>=0.

In both phases, the priority is to "swap" values of diagonal elements to make them increase towards the bottom right. In phase 1, when no further swaps can be made on adjacent rows that improve the diagonal of H, diagonal elements are reduced using a row operation involving the diagonal element and the last row of H. In phase 2, the equivalent situation terminates the entire algorithm.

In phase 1 only, when there are no swaps to perform, the SRS strategy (reluctantly) reduces a diagonal element with a row operation involving row n. Thus the second priority is to "reduce", and sometimes reducing creates a solution in B. Hence the name, "Swap, Reduce, Solve". As noted earlier, phase 2 has no reductions, so the only option is to swap or terminate. Still, "Swap, Reduce, Solve" is an accurate overall description.

For now, we will focus on the first phase of SRS. Sub-section headings below indicate which phase the section refers to.

Once the right-most column is fully reduced, putting a zero in H_n,n-1, the same procedure that did that for column n-1 works for column n-2, then n-3, etc. This procedure only works in columns to the right of which row n is zero. Though the zeroes do appear in row n eventually, they only appear when it improves the diagonal of H.

The reason that reducing (the absolute value of) diagonal elements of H makes progress is that it isolates H_n-1,n-1 as an increasingly large diagonal element, compared to the others, which are being reduced. Remember, a corollary of lemma 10 in the 1999 paper analyzing PSLQ is that the solution is optimal when H_n-1,n-1 is the largest diagonal element.

The technique for reducing |H_n-p,n-p| involves a continued fraction approximation of H_n-p,n-p / H_n,n-p for p=2,3,.... It replaces H_n-p,n-p and H_n,n-p with errors from successive iterations of this approximation. For example, if H_n-p,n-p=.5 and H_n,n-p=.3, a row operation would replace H_n-p,n-p by .5-.3=.2; and a second one would replace H_n,n-p with .3-.2=.1, etc. If these row operations terminate with a zero in H_n-p,n-p, a row swap puts that zero in H_n,n-p and keeps H_n-p,n-p non-zero.

If and when |H_n-p,n-p| is small compared to its neighbor to the left, |H_n-p,n-p-1|, the reduction can stop, because what makes a swap of rows n-p-1 and n-p reduce the upper-left diagonal element, |H_n-p-1,n-p-1|, is a small enough Euclidean length of the vector, (H_n-p,n-p-1, H_n-p,n-p). When to stop reducing is a parameter of the strategy that governs the choice of row operations. Later, we will specify the details of the "Swap, Reduce, Solve" strategy, which say to reduce down to near the precision of the numbers in H. But that is just one strategy, and others may be recommended that stop reducing well above the numerical precision.

If H_n,n-p starts off at zero, no reduction can occur in column n-p, but reduction can be attempted in column n-p-1. This is because the condition that makes the reductions of H_n-2,n-2, H_n-3,n-3 ... possible is that there are only zeroes to the right of these entries and their counterparts in row n of the same column. These zeroes enable arbitrary integer row operations involving rows n-p and n for p=2,3,... until for some p, a zero cannot be made to appear in H_n,n-p. For that p, |H_n-p,n-p| is reduced but no reduction of the same kind is possible in columns to the left of column p.

Once a zero appears in H_n,n-1, it makes reduction work for H_n-2,n-2 and H_n,n-2. The zero in (xB)_n-1 assures that row operations can put a zero in H_n,n-2 while reducing H_n-2,n-2, provided x contains only integers (more on this below). There is no guarantee that reducing H_n-3,n-3 can put a zero in H_n,n-3, but if it can, that zero makes reduction possible in column n-4, etc.

Let's pause here to note what happens to H in large dimensions, like 50 and above. In spite of the best efforts to move large diagonal elements towards the bottom right using adjacent-row swaps, the largest diagonal elements end up in the upper left. Even small numbers in the sub-diagonal, one below the main diagonal -- typically a hundredth to a tenth the size of the main diagonal elements -- prevent swaps of larger diagonal elements from left to right. Diagonal improvement via standard row swaps comes to a halt.

All this changes, once small numbers appear in the diagonal close to the right-hand side of H. Row swaps are unstuck, as they can readily move these small diagonal elements to the upper left. After that is done, new large diagonal elements appear in H_n-2,n-2, H_n-3,n-3 ... and the cycle of diagonal reduction and standard row-swaps repeats.

When a zero appears in H_n,n-1, every entry of column n-1 of H is zero except H_n-1,n-1. This means that, since xBH=0, (xBH)_n-1=0, and therefore (xB)_n-1=0 -- making column n-1 of B an integer-valued solution of <x,m>=0. All this has been stated above in the section "How PSLQ Works".

But the same that goes for column n-1 of H goes for the other columns. Once a zero appears in H_n,n-p, there is

• An integer matrix D with determinant 1 or -1 and inverse E, and
• A non-integer matrix N with small entries off the diagonal

for which DHN has zeroes in column n-p, except (DHN)_n-p,n-p≠0.

On the final step of PSLQ, A is replaced by DA, H by DHN and B by BE. The fact that (xB)_n-p=0 means that column n-p of B is a solution of <x,m> = 0. After these replacements, in every column n-p with H_n,n-p=0, column n-p of B is a solution of <x,m> = 0.

It is expected, but not proven here, that when PSLQ terminates this way, an analog of lemma 10 from the 1997 analysis of PSLQ extends to |H_n-p,n-p|: the solution obtained by putting a zero in H_n,n-p has norm 1/|H_n-p,n-p|; but in this case, that norm is distorted to the extent that N is not a rotation. If so, it is important that D be a full Hermite reduction of H as described in the original 1992 PSLQ paper, making the interior of H below the diagonal small so that N can have small elements off its diagonal. That way, the actual norm of column n-p of B is as close as possible to 1/|H_n-p,n-p|.

It is now possible to give the details of phase 1. These use "gentle Hermite reduction". The code in this repository supports three kinds of gentle Hermite reduction, along with full Hermite reduction as in the classic PSLQ algorithm: See GetInt64D for details. Only the first kind of gentle reduction listed below is under active consideration at the moment, for reasons indicated below.

• Reduce the first n-1 rows of H, but not the last row. In the code, this is indicated by the constant, ReductionAllButLastRow.
• Do not reduce any rows of H, except in the sub-diagonal, and when doing so is necessary to keep the interior of H from blowing up in absolute value. In the code, this is indicated by the constant, ReductionGentle. There are signs from experimentation that this kind of reduction requires the use of BigNumbers when computing the Hermite reduction matrix, D. Currently, D is computed as an int64 matrix.
• Do not reduce any rows of H, except in the sub-diagonal, period. This is indicated by the constant, ReductionSubDiagonal. In the first experiment on live data using this mode, the PSLQ algorithm entered an infinite loop with no changes to H.

During phase 1, the SRS strategy runs in three modes, Swap, Reduce and Solve. Upon exiting Solve mode, a termination check is performed, but Termination Check is not considered a mode, because termination happens just once and it would mess up the name of this strategy. In the description below, "zero" refers to 0 up to the precision of the numbers in H. The modes and transitions between them are as follows.

• Swap: Swap rows (or perform a general integer-valued, unit-determinant row operation) when doing so improves the diagonal of H. Improving the diagonal is defined as starting with |H_j,j| > |H_j+1,j+1| for some j with 1 ≤ j < n-1, performing a row operation and corner removal, and thereby reducing |H_j,j|. The row operation to perform on each iteration is the one that reduces |H_j,j,| by the greatest factor over all choices of j and all operations on rows j and j+1. At the end of each iteration in Swap mode, perform gentle Hermite reduction on the modified rows.
• Reduce: At some point, no swap or other row operation involving adjacent rows among rows 1 through n-1 is left that improves the diagonal of H. Let n-p be the number of the rightmost column for which H_n,n-p is not zero. Reduce H_n-p,n-p against H_p,n-p with continued fraction approximations, until one of the pair reaches a pre-determined, non-zero threshold and Swap mode can reduce |H_n-p-1,n-p-1|; or one of the pair reaches zero. If zero = H_n-p,n-p, or zero < |H_n,n-p| < |H_n-p,n-p|, swap rows n-p and n to keep the diagonal non-zero while minimizing it. Perform gentle Hermite reduction on rows n-p and n.
• Solve: If H_n-p,n-p or H_p,n-p reaches zero in the Reduce stage, a solution of <x,m> = 0 has become available, as described in the section, "A Zero in Row n Generates a Solution". First, perform the Termination Check based on this new solution. Failing that, if the new value of H_n-p,n-p has enabled a row operation to reduce |H_n-p-1,n-p-1|, use that fact to return to the Swap mode. Otherwise, proceed to Reduce mode using H_n-p-1,n-p-1 and H_n,n-p-1.
• Termination Check: If the maximum diagonal element is in a column j for which H_n,j is zero, or if all of row n is zero, terminate the PSLQ algorithm and use the procedure described in the section, "A Zero in Row n Generates a Solution" to generate the available solutions.

Note that after Solve mode puts a zero in row n, Swap mode may overwrite it with corner removal. This happens if, in column n-p with a zero in row n, the diagonal element H_n-p,n-p is small enough that Swap mode performs a row operation on rows n-p-1 and n-p, then removes the corner in H_n-p-1,n-p. In the course of that corner removal, the zero in H_n,n-p is replaced by a non-zero.

The reason to delay the Reduce mode until no row operations are available to Swap mode is that when Reduce mode leads to Solve mode, zeroing out H_n,n-p, the initial value of |H_n-p,n-p| should be as large as possible. Starting with a large |H_n-p,n-p| gives the greatest chance of ending with a large |H_n-p,n-p|. Running in Swap mode increases the diagonal elements in columns like n-p with solutions, as they are on the right side of H where Swap mode is putting large diagonal elements.

The reason to stop reducing H_n-p,n-p against H_n,n-p when one of the pair reaches a pre-determined minimum is to avoid underflow, while introducing small elements into the diagonal of H so the "solved" columns are large by comparison.

Waiting for the maximum diagonal element to appear in a solution column (i.e. one with a zero in row n) seems to offer the best chance of solving the shortest vector problem. But it is not a guarantee of solving it, because only column n-1 contains a diagonal element with the undistorted reciprocal of a solution norm.

As promised, here is an explanation of why both H_n-2,n-2 and H_n-3,n-3 can be reduced given integer-valued input x and a non-zero H_n,n-2 and H_n,n-3. As mentioned above, this is because a zero appears in H_n,n-2 when reducing |H_n-2,n-2|. The key to why that zero appears is that

0 = <((xB)_n-2, (xB)_n-1, (xB)_n), (H_n-2,n-2, H_n-1,n-2, H_n,n-2)>

    =<((xB)_n-2, 0, (xB)_n), (H_n-2,n-2, H_n-1,n-2, H_n,n-2)>

    =<((xB)_n-2, (xB)_n), (H_n-2,n-2, H_n,n-2):

is an integer relation between H_n-2,n-2 and H_n,n-2. This guarantees that H_n-2,n-2 / H_n,n-2 is rational. The row operations that mirror the continued fraction approximation of this ratio put an error of zero in H_n,n-2 (or H_n-2,n-2) on the last of finitely many steps. If the zero appears in H_n-2,n-2, you would just swap rows n-2 and n to put the zero in H_n,n-2.

Since we are contemplating the placement of very small entries in the diagonal of H, it may take several rounds of row swaps, corner removals and reduction of sub-diagonal elements in the same 2x2 sub-matrix before this sub-matrix has its best ordering of diagonal elements. Consider, for example, the 2x2 sub-matrix M_k =

After swapping rows and zeroing the corner, M_k+1 =

A new round of row reduction, swap and corner removal finally yields the best form M can take in isolation from the rest of H: M_k+2 =

But there is a way to collapse the two rounds of swap, reduction and corner removal into one "general" (non-swap) row operation: RM_kQ =

To save operations when putting small elements in the diagonal, it could be worth the while to look for general row operations -- even though it does cost a bit to check for them. The reason small diagonal elements create the need for more than one round of row swaps will become apparent below.

Suppose a small ε₁ has just been placed in H_j+1,j+1. The introduction of ε₁ changes the balance of the diagonal elements in the 2x2 sub-matrix, M_k, containing t=H_j,j, u=H_j+1,j and ε₁=H_j+1,j+1. Using continued fraction reduction, find relatively prime, non-zero integers a and b such that at+bu=ε₀ is small. Let R be the 2x2 matrix with rows [a, b] and [-w, v] with minimum |v| and determinant 1. R improves the diagonal of M_k: RM_k=

After zeroing out the upper right corner, the diagonal elements of M_k+1 have absolute values

|H_j,j| ← |ε₂| := sqrt(ε₀² + b² ε₁²) and

|H_j+1,j+1| ← |ε₃| := |t ε₁ / ε₂| (since |det(M_k+1)| = |det(M_k)| =|t ε₁|)

Let's compare ε₂ to what |H_j,j| gets after a row swap and corner removal, which yeilds:

|H_j,j| ← ε₄ := sqrt(u² + ε₁²)

|H_j+1,j+1| ← ε₅ := |t ε₁ / ε₄|

The row operation, R, performs better than a swap if it reduces |H_j,j| more, i.e.

|ε₂| < |ε₄| ⇔ sqrt(ε₀² + b² ε₁²) < sqrt(u² + ε₁²)

    ⇔ ε₀² + b² ε₁² < u² + ε₁²

    ⇔ (b² - 1) ε₁² < u² - ε₀²

From this we can conclude that

• The general row operation, R, performs best when b and ε₀ are small. This is possible when -a/b is a good approximation of t/u with a small denominator -- the kind you get from continued fractions.

• b is not a unit, because that would be incompatible with the fact that after Hermite reduction, or reduction of just the sub-diagonal of H, |u| ≤ |t|/2. If you set b to a unit, you get a contradiction as follows. Since R is not a row swap, a and b are not zero. Therefore,

       b² = 1 ⇒ |a| ≥ |b| = 1

          ⇒ |ε₀| ≥ min(|at+u|, |at-u|) ≥ min(|t+u|, |t-u|) ≥ |t|/2 (since |u| ≤ |t|/2)

          ⇒ u² - ε₀² ≤ u² - t²/4 (since |ε₀| ≥ |t|/2)

          ⇒ (b²-1) ε₁² < u² - ε₀² ≤ u² - t²/4 (since R is better than a row swap)

          ⇒ 0 < u² - ε₀² ≤ u² - t²/4 (since b² = 1)

          ⇒ t²/4 < u² (using the two ends of the inequality on the previous line)

          ⇒ |u| > |t|/2,

       which contradicts the premise that |u| ≤ |t|/2

• Since |b| ≥ 2 as just shown,

       (b² - 1) ε₁² < u² ⇒ 3 ε₁² < u²,

       so only row swaps make sense unless |ε₁| < |u| / sqrt(3).

• Since |b| ≥ 2 and the condition (b² - 1) ε₁² < u² - ε₀² is what makes R better than a row swap, the range of values of b to consider as the upper-right entry of R satisfies

       4 ≤ b² ≤ 1 + (u² - ε₀²) / ε₁² ≤ 1 + u² / ε₁²

          ⇒ 2 ≤ |b| ≤ sqrt(1 + u² / ε₁²),

       which gives a stopping condition for the reduction of t and u.

The last bullet puts us in a position to circle back to a claim above. The claim was that introducing small elements into the diagonal of H creates the need for multiple rounds of swaps, corner removal and reduction, unless general row operations are used. It is believed that each such round corresponds to one round of continued fraction approximation of t and u. Though this is not proven here, the last bullet shows that when |ε₁| is small (compared to |u|), an R with a large upper-right entry b still has the potential to out-perform a row swap. Large entries comparable in absolute value to sqrt(1 + u² / ε₁²) only appear in R after multiple rounds of continued fraction approximation.

To take advantage of the foregoing analysis, a strategy, like SRS, that places very small elements in the diagonal of H should save time using general row operations. After placing ε₁ in H_j+1,j+1, the strategy would reduce t=H_j,j against u=H_j+1,j, storing the reduced entry in H_j,j and the version of R that puts it there at each stage. If the smallest absolute value of the stored H_j,j entries is smaller than sqrt(u² + ε₁²), the strategy would use the corresponding R instead of a row swap in the next PSLQ iteration.

Phase 2 is, thankfully, much easier to explain than phase 1. As noted earlier, phase 2 begins when all columns of H represent solutions of <x,m>=0, and no swapping of adjacent diagonal elements can improve the diagonal of H. In phase 2, any two rows j₀ < j₁ < n for which the Euclidean length L₁ of (H_j1,j0, H_j1,j0+1, ... H_j1,j1) is less than L₀ = |H_j0,j0| can be swapped. After performing this swap and rotating to remove the non-zero elements to the right of H_j0,j0, |H_j0,j0| = L₁ < L₀ = |H_j1,j1|. This improves the diagonal of H because

• Before swapping: |H_j1,j1| ≤ L₁ < L₀ = H_j0,j0. Diagonal elements H_j0,j0 and H_j1,j1 are out of order.
• After swapping: |H_j1,j1| = L₀ > L₁ = H_j0,j0. Diagonal elements H_j0,j0 and H_j1,j1 are now in order.

In phase 2, each possible pair, (j₀, j₁), is ranked by how large L₀/L₁ is -- the larger the better. The pair with the largest L₀/L₁ is swapped. If no pair has L₀/L₁ > 1, phase 2 and the entire algorithm terminates.

Phase 2 works surprisingly well, albeit slowly. The reason it works well is the same reason it works slowly: There is no shortage of non-adjacent row swaps to perform, even though each requires not just one sub-diagonal element to be accounted for when computing L₁, but j₁ - j₀ of them. The abundance of possible non-adjacent row swaps makes for slow, steady progress. The time it would take to terminate phase 2 may still be polynomial, as phase 1 probably is, but the polynomial would have a degree one higher than that of phase 1 because it operates on arbitrary pairs of rows, not single rows and their immediate successors.

Since the inputs for high-dimension problems are exquisitely precise, so too are the elements of H during phase 1. This enables the algorithm to recognize when a solution is found by the fact that an entry in the last row of H is essentially zero.

After phase 1, however, there is no longer a need to recognize solutions. The solutions, which already lie in the columns of B, are just being combined by integer column operations (whose inverses are row operations in H). H is no more than a rough guide, relatively speaking, for which rows to swap.

Because H no longer needs to be kept with high precision, the bignumber package used in phase 1 is no longer needed. H can be kept in float64 or even float32 throughout phase 2. However, the repository does not yet take advantage of this insight.

Initial experiments with SRS have shown that it finds solutions in higher degrees than "IDASIF". The building blocks of phase 1 can be found in the BottomRightOfH and RowOpGenerator in the code. Phase 2 uses HPairStatistics. Details on the performance of SRS, and actual code for SRS, can be obtained under a separate, individually negotiated, license.

To run PSLQ using this repository as a library, you can emulate the way the tests in pslqops/run_test.go and strategy/strategy_test.go use pslqops.OneIteration. See below the overview of the pslqops and strategy packages, and the bignumber and bigmatrix packages they depend on.

The bignumber package enables computation of sums, differences, products, quotients and square roots of numbers with arbitrary precision. The precision must be set at most once (it has a default of 1000 bits) using the function, bignumber.Init. bignumber also implements operations between bignumber and int64 instances, so PSLQ can run a bit faster before it switches over to using bignumber for almost all operations, deep into most runs.

bignumber is similar to the native golang big.Float, except that

• Overall minimum precision is set for all instances of bignumber, whereas the precision of big.Float is per instance.
• There is hereby an explicit guarantee that bignumber instances initialized with int64 instances, and any combinations of these, are integer-valued with no round-off error. The underlying big.Int incorporated into bignumber makes this self-evident by reading the code. Though this is most likely true of big.Float, it is not guaranteed as far as we know.

Though bignumber has a few constructors, the one the pslqops package uses to take input is the one that parses base 10 (decimal) string input into a bignumber.

The bigmatrix package enables multiplication, addition and subtraction of matrices with entries that are bignumber instances. It also implements operations between bignumber instances and int64 instances, so PSLQ can run a bit faster before it switches over to using bignumber for almost all operations, deep into most runs.

Though bigmatrix has a few constructors, the one the pslqops package uses to take input is the one that parses arrays of base 10 (decimal) string inputs into a bigmatrix.

The pslqops package contains all the building blocks of PSLQ, except for non-standard strategies for choosing row operations (see the strategy package for those).

This package includes

• New for constructing a pslqops.State, which keeps track of the matrices x, B, A and H. New takes an array of strings representing the input to PSLQ in decimal form.
• pslqops.State.OneIteration, the top-level function that performs one iteration of PSLQ.

OneIteration takes as an argument a function, getR, that examines H and returns R, a row operation that performs a row operation on H. In the classic PSLQ from the original 1992 PSLQ paper, R is a swap of adjacent rows j and j+1 for which a certain quantity is maximized (see pslqops.GetRClassic or the original 1992 PSLQ paper). Other rules for choosing R are implemented in the strategy package. One of these strategies is what Table 1 shows results for in rows labeled IDASIF.

A point of confusion could be that getR does not return anything called "R". It returns a RowOperation type saying what rows to operate on and what matrix, or permutation, to apply. That's OK, it's still an R matrix -- in the sense that the original 1992 PSLQ paper uses that notation -- in the form that OneIteration accepts.

PSLQ maintains invariants like equation 2, xBH = 0, which you can verify with GetObservedRoundOffError. Another invariant verifier is CheckInvariants, which verifies that B = A^-1 and that the upper right of H contains zeroes, up to round-off error.

The strategy package is where all the fun ideas for improving the empirical performance of PSLQ are defined. These are the functions passed to pslqopa.OneIteration as parameter getR. This package is in flux as new ideas are tried. In order to avoid making Table 1 out of date, only the IDASIF strategy (constant improveDiagonalWhenAboutToTerminate in strategyv1.go) will necessarily be retained as-is.

All the mathematical variable names in the original PSLQ paper are incorporated into the golang variable and/or function names in the source code. The main examples are

• the matrices A, B, D, E, G, H and R have variable and/or function names to match.
• the vectors s and x have variable and/or function names to match

Though there are exceptions, the rule of thumb is, if a variable is mentioned in the paper, it is used under the same name in the code; but not vice-versa.

The original PSLQ paper and 1999 paper analyzing PSLQ cover a variety of properties of the H matrix. But it was left to later research to discover a key fact about the diagonal of H: The reciprocal of the largest diagonal element in H is a bound on the size of any solution.

This section fills in gaps like this in the two papers, using other references and original results.

Here we present one geometric view of PSLQ. There is at least one other geometric view, not covered in detail here: PSLQ finds an integer matrix with determinant 1 whose columns approximate the solution plane,

S = {m : <x,m> = 0}

So what is presented in detail here is not the geometric view of PSLQ, but it is one very appealing interpretation.

PSLQ computes a matrix, A_k, at every iteration k. This "A" is the same A as in the PSLQ paper, in the invariants above and in the section below, "A Sharper Lower Bound on the Smallest Solution While PSLQ is Running" -- the "Sharper Bound" section for short. Here as in that section, the subscript k is useful to track A through different iterations k=1,2,3,... of PSLQ.

Successive A_k get closer and closer to a change of basis, followed by a rotation and rounding, when applied to S. To see why, let

• (H_x)_p be column p of H_x for p=1,...,n-1
• m = ∑_p y_p (H_x)_p be an arbitrary element of S

Then A_km = H_k(Q_kH_x^tm) (equation 9 from the "Sharper Bound" section)

The change of basis comes from the product H_x^tm in the right-hand side of equation 9. In the context of A_km, m is expressed in terms of the basis (e₁, ..., e_n). But H_x^tm gives m in terms of the basis ((H_x)₁, ..., (H_x)_n-1). In other words,

(H_x^tm)_i = y_i (equation 3)

The reason for this is that

H_x^tH_x = I_n-1,

as noted in section 3 of the 1992 PSLQ paper. The following calculation uses this identity to prove equation 3:

(H_x^tm)_i = (H_x^t (∑_p y_p (H_x)_p)_i

     = (∑_p y_p H_x^t (H_x)_p)_i

     = y_i

In equation 9, once H_x^t applies a change of basis to m, the result is rotated by Q_k, then dilated with error by H_k. The fact that Q_k is a rotation matrix is explained in the "Sharper Bound" section. The dilation comes from the diagonal elements of H_k, and the error comes from the off-diagonal elements -- including all of row n (so there are really just n-1 meaningful entries in A_km).

Some error is necessary because the left-hand side of equation 9 is an integer matrix. This means that the error in the off-diagonal elements of H_k can be considered to be a rounding error. But this rounding error decreases with each iteration of the PSLQ algorithm, because H_k tends towards a diagonal matrix as k increases.

In summary, A_km is m written as a combination of the columns of H_x, rotated and dilated with rounding error.

There is a sharper bound than 1/|H| (from the original 1992 PSLQ paper) on the size of a solution while the algorithm is still running (i.e., when H_k has no 0s in its diagonal -- a fact used below). This bound,

1/max(|H_1,1|, |H_2,2|, ..., |H_n-1,n-1|) ≤ |m| for any solution m of <x, m> = 0,

is found, among other places, on pages 97-99 of linear Algebra in Situ, CAAM 335, Fall 2016 by Steven J. Cox. This is a textbook that covers many topics, including QR decomposition. QR decomposition is the same as LQ decomposition, used in PSLQ, except every matrix is transposed. Because the overall topic is QR decomposition in this work, every matrix in the PSLQ algorithm is transposed there; and many are renamed. In what follows, the argument in "Linear Algebra in Situ" is repeated here, but in the LQ context, using similar names to those in the original PSLQ paper and in the source code of this repository.

The notation used below follows the original PSLQ paper, except many matrices are indexed by an iteration number denoted k. Initial matrices are:

• x, the input to PSLQ. It is a unit vector of real numbers, none of which is 0.
• H_x is the initial value of the n x n-1 matrix H.
• P = H_xH_x^t

Below is notation for a specific iteration k of the PSLQ algorithm as presented in the original PSLQ paper. k starts at 1 (as opposed to 0). If k = 1, H_k-1 = H_x.

• Step 1
  • H_k is the n x n-1 matrix H after iteration k.
  • D_k is the n x n integer matrix used to update H_k-1 in step 1 of iteration k.
• Step 2
  • j is the integer selected in step 2 of iteration k.
• Step 3
  • R_k is the n x n permutation matrix such that R_kM swaps rows j and j+1 of M. The PSLQ paper names this R_j, after the starting index j of the row swap. In what follows we need to track R over multiple iterations, so the k subscript is necessary.
  • G_k is the n-1 x n-1 orhtogonal matrix that the PSLQ paper calls G_j. The same comment about subscript j vs. k applies to G that applied to R.

Using this notation, iteration k can be interpreted as:

1. H <- D_kH_k-1. H is an intermediate value, not quite H_k yet.
2. Choose j so R_k and G_k are defined.
3. H_k <- R_kHG_k = R_kD_kH_k-1G_k

After iteration k,

H_k = R_kD_kH_k-1G_k

    = R_kD_kR_k-1D_k-1H_k-2G_k-1G_k

    = ...

    = R_kD_kR_k-1D_k-1...R₁D₁ H_x G₁...G_k-1G_k (equation 4)

Let A_k and Q_k^-1 be what lie to the left and right of H_x, respectively, in equation 4:

• A_k = R_kD_kR_k-1D_k-1...R₁D₁
• Q_k = (G₁...G_k-1G_k)^-1

A_k is the same "A" as in the original PSLQ paper.

Computation of the bound mentioned at the beginning of this section,

1/max(|H_1,1|, |H_2,2|, ..., |H_n-1,n-1|) ≤ |m| for any solution m of <x, m> = 0 (equation 5),

begins with the LQ decomposition of A_kH_x:

H_k = A_kH_xQ^-1

A_kH_x = H_kQ (equation 6)

Equation 6 is an LQ decomposition of non-singular A_kH_x, because

• A_k is an n x n integer matrix with determinant 1, like all of the R_i and D_i in the original PSLQ paper.
• Q_k is orthonormal, like all of the G_i in the original PSLQ paper.

As noted earlier, the PSLQ paper defines a matrix P = H_xH_x^t. P fixes any m for which <x,m> = 0. In other words,

xm = 0 ⇒ Pm = m (equation 7)

From equation 7 comes the following proposition: If (A_km)_p,1 = 0 for p < i, then

(A_km)_i,1 = (H_k)_i,i (Q_kH_x^tm)_i,1 (equation 8)

Substituting from equation 7 in the first line and equation 6 in the fourth line,

A_km = A_kPm

    = A_k(H_xH_x^t)m

    = (A_kH_x)(H_x^tm)

    = (H_kQ)(H_x^tm)

    = H_k(Q_kH_x^tm) (equation 9)

Using equation 9, we will now calculate (A_km)_i,1, starting with i = 1, until (A_km)_i,1 ≠ 0. The index p in the summations below ranges from 1 to i, after which (H_k)_i,p = 0.

If i = 1, then using equation 9 in the second line below,

(A_km)_i,1 = (A_km)_1,1

     = ∑_p (H_k)_1,p (Q_kH_x^tm)_p,1

     = (H_k)_1,1 (Q_kH_x^tm)_1,1 (equation 10)

If (A_km)_1,1 = 0, we continue with i = 2.

0 = (A_km)_1,1 = (H_k)_1,1 (Q_kH_x^tm)_1,1

Since H_k has no 0s on its diagonal,

(Q_kH_x^tm)_1,1 = 0 (equation 11)

(A_km)_i,1 = (A_km)_2,1

     = ∑_p (H_k)_2,p (Q_kH_x^tm)_p,1

     = ((H_k)_2,1)(0) + (H_k)_2,2 (Q_kH_x^tm)_2,1

     = (H_k)_2,2 (Q_kH_x^tm)_2,1 (equation 12)

If (A_km)_1,1 = (A_km)_2,1 = 0, we continue with i = 3.

0 = (A_km)_2,1 = (H_k)_2,2 (Q_kH_x^tm)_2,1

From equation 11 and since H_k has no 0s on its diagonal,

(Q_kH_x^tm)_1,1 = (Q_kH_x^tm)_2,1 = 0 (equation 13)

(A_km)_i,1 = (A_km)_3,1

     = ∑_p (H_k)_3,p (Q_kH_x^tm)_p,1

     = ((H_k)_3,1)(0) + ((H_k)_3,2)(0) + (H_k)_3,3 (Q_kH_x^tm)_3,1

     = (H_k)_3,3 (Q_kH_x^tm)_3,1

This reasoning continues until the first i for which (A_km)_i,1 ≠ 0. The formula for (A_km)_i,1 is

(A_km)_i,1 = (H_k)_i,i (Q_kH_x^tm)_i,1 (proving equation 8)

Recall that the bound to prove is

1/max(H_1,1, H_2,2, ..., H_n-1,n-1) ≤ |m| for any solution m of <x, m> = 0 (repeating equation 5)

Let

• i be the smallest index for which (A_km)_i,1 ≠ 0
• (Q_kH_x^t)_i denote row i of Q_kH_x^t

Note that

• A_k and m are non-zero integer matrices and A_k is non-singular, which makes the first line work in the calculation below
• Equation 8 from the section, "A Formula for (A_km)_i,1", permits the replacement of (A_km)_i,1 in the second line below.
• Q_k is a product of the inverses of matrices G_k, defined in equations 10 through 15 of the original PSLQ paper. These equations define G_k as a rotation matrix. This makes Q_k a rotation matrix, which is one of two facts used in the fourth line below to conclude that the norm of a row in Q_kH_x^t is 1.
• H_x^tH_x = I_n-1, which is the second fact needed to conclude that the norm of a row in Q_kH_x^t is 1.

1 ≤ |(A_km)_i,1|

     = |(H_k)_i,i (Q_kH_x^tm)_i,1|

     ≤ |(H_k)_i,i| |(Q_kH_x^t)_i| |m|

     = |(H_k)_i,i| |m|

     ≤ max(H_1,1, H_2,2, ..., H_n-1,n-1) |m| (proving equation 5)

This section explains the ideas behind the alternative row operations for improving the diagonal of H. As seen in the previous section, "A Sharper Lower Bound on the Smallest Solution While PSLQ is Running", reducing the maximum diagonal element of H_k sharpens the bound on the smallest solution m to the integer relation problem, <x,m> = 0. This raises the question, when would a given row swap reduce the maximum diagonal element of H_k and thereby reduce this lower bound?

Following the notation in a 1999 paper analyzing PSLQ, by the same authors as the original PSLQ paper, the two rows and columns involved in both the row swap and corner steps of an iteration of PSLQ are

The 1999 paper also defines δ = sqrt(β² + λ²).

The 1999 paper analyzing PSLQ derives the formula

for the result of the row swap and cornering. Up to absolute value, Λ₁ is obtained by left-multiplying Λ₀ by

The row swap and corner steps can be considered to reduce the maximum diagonal element if

max(|δ|, |α λ / δ|) < max(|α|, |λ|) (equation 14)

The row swap and corner steps reduce the maximum diagonal element if and only if

|α| > |δ| > |λ| (equation 15)

To prove this, first assume equation 14 and argue for equation 13. Equation 14 precludes the possibility that |α| < |λ|, since

|α| < |λ| ⇒ |λ| = max(|α|, |λ|) > max(|δ|, |α λ / δ|) ≥ |δ| = sqrt(β² + λ²) ≥ |λ|, a contradiction.

Therefore, |α| ≥ |λ|. Using equation 14,

|α| = max(|α|, |λ|) > max(|δ|, |α λ / δ|)

        ⇔ |α| > |δ| and |α| > |αλ / δ|

        ⇔ |α| > |δ| and 1 > |λ / δ|

        ⇔ |α| > |δ| and |δ| > |λ| (a restatement of equation 15)

For the reverse direction, assume equation 15. Then since |α| > |δ|,

|α λ / δ| > |δ λ / δ| = |λ|

Therefore,

max(|δ|, |α λ / δ|) ≥ max(|δ|, |λ|)

     = max(sqrt(β² + λ²), |λ|)

     = sqrt(β² + λ²)

     = |δ| (equation 16)

Equation 16 selects which of max(|δ|, |α λ / δ|) is the maximum, namely

max(|δ|, |α λ / δ|) = |α λ / δ| (equation 17)

Using equation 17 and the fact from the premise, equation 15, that |λ| < |δ|,

max(|δ|, |α λ / δ|) = |α λ / δ| < |α| ≤ max(|α|, |λ|) (equation 18)

Equation 18 proves equation 14.