A fast algorithm to find reduced hyperplane unit cells and solve N-dimensional Bézout’s identities

The paper describes a method to determine a short unit cell attached to any hyperplane given by its integer vector p. Equivalently, it gives all the solutions of the N-dimensional Bézout’s identity associated with the coordinates of p.


Introduction
How to determine a unit cell attached to a plane p ¼ ðh; k; lÞ? This problem occurs for example in the crystallographic models of twinning, when the obliquity or the shear values must be calculated for many planes. It is intuitively solved for low-index planes, but the solutions are more difficult to obtain for high-index planes. In addition, if a unit cell can be found, can it be reduced to a smaller one? In dimension N, the difficulty of finding a small unit cell attached to a hyperplane of dimension N À 1 becomes even more pronounced. Let us express mathematically this 'hyperplane unit-cell problem' by the notations detailed in Appendix A. We assume that a hyperplane is known only by its Miller indices h; k; l which are coprime integers, or equivalently by its normal which is expressed as an integer vector of coordinates h; k; l in the reciprocal space. We want to determine a small unit cell such that one short integer vector of the cell points to a node of the first layer parallel to the hyperplane, and the other N À 1 short integer vectors lie in the hyperplane. The 'out-of-plane' vector and the 'in-plane' vectors are noted b 1 and b 2 ; . . . ; b j ; . . . ; b N , respectively. The vector b 1 is such that p t b 1 ¼ 1, and the N À 1 vectors b 2 ; . . . ; b j ; . . . ; b N are such that p t b j ¼ 0. The coordinates of the vector b 1 constitute a solution of the N-dimensional Bé zout's identity formed on the coordinates of p. The coordinates of any of the vectors b j , j 2 f2; . . . ; Ng, are solutions of what is called an 'integer relation' with the coordinates of p (Appendix A). For example, with N ¼ 3, the integer coordinates u; v; w of b 1 verify the equation uh þ vk þ wl = 1; the integer coordinates u; v; w of the vector b 2 (or b 3 Þ verify the equation uh þ vk þ wl ¼ 0, as illustrated in Fig. 1. Finding solutions to integer relations is not complicated. For N ¼ 3, if we know a plane p ¼ ðh; k; lÞ with let us say k 6 ¼ 0, it is not difficult to find two integer vectors b 2 and b 3 in this plane, for example, b 2 ¼ ½Àk; h; 0 and b 3 ¼ ½0; Àl; k. The difficult part of the problem is to find vectors with small coordinates by considering all the possible linear combinations of the Miller indices h; k; l. Finding the shortest solutions in dimension N is an NP-hard (non-deterministic polynomial-time) problem well known in computer science and cryptography. An algorithm called PSLQ gives short solutions to integer relations with any vector p 2 R N (Ferguson et al., 1999; see also Wikipedia, 2021a). It has permitted the discovery of numerous previously unknown identities among real numbers; one of them is the formula that allows the calculation of the nth hexadecimal digit of without computing the preceding digits (Bailey et al., 1997;Raayoni et al., 2021). The algorithm presented in the present paper gives only solutions for the vectors p 2 Z N but, as we will show, the vectors we obtain are shorter than those obtained by PSLQ. Our algorithm actually provides simultaneously a short solution to the N-dimensional Bé zout's identity and short solutions to the 'integer relations'. It gives the affine space of all the solutions of the N-dimensional Bé zout's identity. From a crystallographic point of view, it provides a small unit cell attached to a hyperplane.
Recently, Gorfman (2020) proposed a method to find some solutions to an intermediate problem that we will call the 'column-constrained unimodular matrix' (CCUM) problem in order to differentiate it from the initial 'hyperplane unit-cell problem'. The CCUM problem consists of finding a unimodular matrix M such that the first column is equal to a fixed integer vector t. We recall that a unimodular matrix has integer entries and its determinant is AE1. Note that, in Gorfman's paper, it was the last vector (and not the first one) that was imposed, but that does not change the problem. Gorfman's approach involves a series of multiplication with matrices called S containing 0, 1 and À1 in order to reduce the imposed vector t to a unit vector (a vector for which one of its coordinates is 1 and the others are 0). Gorfman showed that the same algorithm applied in the reciprocal space to a vector p gives a solution to the hyperplane unit-cell problem. Let us explain how it works with our notations. For an imposed reciprocal vector p, Gorfman's method permits one to obtain a unimodular matrix M Ã that has p for the first column vector. Then, the inverse of its transpose M ¼ ðM Ã Þ Àt is calculated. Since M Ã is a unimodular matrix, the matrix M is also unimodular, which implies that its columns are integer vectors. Let us call them b j . Since ðM Ã Þ t M is the identity matrix, its first column is a vector that has 1 as the first coordinate and 0 for all the other coordinates. This means that p t b 1 ¼ 1 and p t b j ¼ 0 for j 2 f2; . . . ; Ng, which proves that the matrix M is a solution of the hyperplane unit-cell problem. Gorfman's idea of using unimodular matrices is very interesting and his approach is innovative and inspiring, but it does not give short solutions. For example, for the plane p t ¼ ð12; 20; 225Þ, the solution determined by his algorithm in which the first imposed column vector is p is The inverse of its transpose is The reader can check that the scalar product of ð12; 20; 225Þ with the first column vector ½148; À100; 1 is 1, and that the scalar product with the last column vectors ½À595; 402; À4 and ½5; À3; 0 is 0. However, the vector ½148; À100; 1 solution of the 3D Bé zout identity and the vector ½À595; 402; À4 solution of the integer relation are large. The vector ½À595; 402; À4 is even larger than the obvious solution ½0; À4; 45. More generally, Gorfman suggests that the algorithm could be 'an alternative approach to calculate the Bezout coefficients', but we would like to show that the opposite approach is possible. The aim of the paper is to show that determination of the Bé zout's coefficients is an efficient way to find short solutions of both the CCUM problem and hyperplane unit-cell problem. The algorithm proposed in the present paper is based on Euclidean division. An algorithm to determine some short solutions to the N-dimensional Bé zout's identity is proposed in Section 2. The algorithm to solve the CCUM problem is detailed in Section 3. Sections 2 and 3 are independent. In Section 4, we explain how to combine the two algorithms to find short solutions to the hyperplane unit-cell problem. Some examples will be given and compared with the PSLQ algorithm. The method has been encoded in a Python 3.8 computer program called GeneralizedBezout. The examples given were obtained on a laptop computer equipped with an Intel Core i7-4600 CPU 2.1 GHz, 64-bit Windows system with a RAM of 8 GB. Note: the Python program General-izedBezout is freely available from the author upon request.

N-dimensional Bézout's identity
Given a set of integers fp i ; i ¼ 1; . . . ; Ng we look for another set of integers fu i ; i ¼ 1; . . . ; Ng such that P N i¼1 p i u i ¼ 1. In other words, given an integer vector p of coordinates p i , we want to get the coordinates u i of an integer vector u that is such that p t u ¼ 1. If N = 2, the fast and well known algorithm based on Euclidean division gives a solution that is also the shortest one (Capparelli, 2020;Wikipedia, 2021b). Surprisingly, we could not find in the literature algorithms in high dimensions N. We propose here two recursive algorithms. Unit cell associated with the plane p ¼ ðh; k; lÞ. The out-of-plane vector b 1 points to a node of the layer q = 1, and the in-plane vectors b 2 and b 3 lie in the layer q = 0. The vector b 1 is a solution of the Bé zout's identity p t b 1 ¼ 1, and the vectors b 2 and b 3 are solutions of the integer relations p t b 2 ¼ 0 and p t b 3 ¼ 0.
They give different solutions that are all valuable, but we will see that the second one gives shorter solutions.
Method-0. We consider p 1 and p 2 the two first coordinates of p, and we call ðu; vÞ their Bé zout numbers, i.e. up 1 þ vp 2 ¼ gcd ðp 1 ; p 2 Þ. If we note fk i ; i ¼ 2; . . . ; Ng the Bé zout numbers in dimension N À 1 associated with the set fgcd ðp 1 ; p 2 Þ; p 3 ; . . . ; p N g, a solution of the N-dimensional Bé zout's identity is thus fuk 2 ; vk 2 ; k 3 ; . . . ; k N g. This method is easy to compute by recursion until the dimension decreases to N = 2 for which the solution is given by the classical Bé zout's algorithm. The problem related to this method is that the absolute values of the Bé zout numbers u i can be quite high. One could screen all the pairs ðp i ; p j Þ in place of ðp 1 ; p 2 Þ to determine the lowest Bé zout numbers but this method would be unrealistic for high dimensions N. We could find another method for which the values are lower than those determined by method-0.
Method-1. We consider the set of integers fp i ; i ¼ 1; 2; . . . ; Ng. If 9i; jp i j ¼ 1, the solution of the Bé zout identity is immediately f0; . . . 0; p i ; 0 . . . ; 0g. If none of the p i 's has 1 as absolute value, the set fp i g is sorted in decreasing order of the absolute values of p i . The sorting permutation is kept in memory. The smaller non-null value is called p i 0 . We calculate the quotient set and the residue set {q i ; i = 1; 2; . . . ; i 0 À 1} and fr i ; i ¼ 1; 2; . . . ; i 0 À 1g with q i ¼ bp i =p i 0 c and r i ¼ p i À q i p i 0 , quotient and remainder of the Euclidean division by p i 0 . If we note fu 1 ; u 2 ; . . . ; u i 0 À1 ; 0; . . . ; 0g the Bé zout numbers associated with the set {r 1 ; r 2 ; . . . ; r i 0 À1 ; 0; . . . ; 0}, a solution of the N-dimensional Bé zout's identity is fu 1 ; u 2 ; . . . ; u i 0 À1 ; À P i 0 À1 i¼1 q i u i ; 0; . . . ; 0g. This method is easy to compute by recursion until one of the absolute values of the input vector is 1. The correct order of the Bé zout numbers associated with the initial set fp i ; i ¼ 1; 2; . . . ; Ng is restored by applying À1 . The pseudocode is given in Fig. 2.
The Bé zout numbers calculated with method-1 are smaller than those obtained by method-0. Only method-1 will be considered in the rest of the paper. With the vector p t ¼ ð12; 20; 225Þ, it gives u ¼ ½À17; À1; 1. With the vector p t ¼ ð51; 450; À102; 240; À277; 54; 450; 532Þ, it gives u = [À3, 0, 0, 0, 0, À3, 0, 1]. The calculation lasts only a few ms. Even if method-1 gives small Bé zout vectors u, it may not give systematically the smallest ones. We will see in Section 4 how 'hyperplane shearing' can give shorter Bé zout vectors u with the help of the CCUM algorithm detailed in Section 3.
3. Algorithm to solve the column-constrained unimodular matrix problem 3.1. Case where one of the coordinates of t is AE1 Now we consider the CCUM problem. There is a simple and immediate solution if the first coordinate of t is 1. In that case, any diagonal or even triangular matrix M with 1 in the diagonal and with t as the first column checks the condition det(M) = 1. If the first coordinate of t is À1, changing one 1 into À1 in the diagonal is sufficient to maintain det(M) = 1. The example used by Gorfman (2020)  Note that the result is obtained without any calculation. If one of the coordinates of t is 1 in a position i > 1, then a simple matrix of permutation P is sufficient to recalculate the matrix M. We will not give more details here because the solutions are actually included in the more general method based on Bé zout's identity explained as follows.

Case where t has at least one pair of coprime coordinates
In the case N = 2, the general solution to the CCUM problem is given by the classical 2D Bé zout's identity. We note There is a solution if and only if the integers t 1 and t 2 are coprime, and the solution is simply where u, v are the Bé zout numbers associated with t 1 and t 2 , i.e. solutions of the equation ut 1 þ vt 2 ¼ 1. If t 1 and t 2 are not coprime, the determinant of any matrix M with t in the first column would be a multiple of gcd(t 1 , t 2 ), the greatest common divisor of t 1 and t 2 , and thus cannot be equal to AE1. The resulting vector Pseudocode to find Bé zout numbers associated with the coordinates of a vector p. Now, we consider the case where N > 2 and the vector t has its two first coordinates t 1 and t 2 that are coprime numbers. We consider the matrix M made of two blocks, the top left one is where u, v are the Bé zout numbers associated with t 1 and t 2 , and the bottom right one is the ðN À 2Þ Â ðN À 2Þ identity matrix. Then, the first column of M is replaced by t (t 1 and t 2 are not changed, and the zeros in M i,1 are replaced by t i ; i > 2Þ. The matrix M is the solution of the CCUM problem. When the two coprime coordinates of vector t, t 1 and t 2 , are not the first ones, the permutation matrices Pði; 1Þ and Pðj; 2Þ are used to return to the previous case. We recall that a permutation matrix Pði; jÞ is a N Â N identity matrix, except for the line i for which 1 is written in the column j, and for the column j where 1 is written in the line i. Permutation matrices are unimodular matrices and are equal to their inverse. The unimodular matrix P = Pði; 1ÞPðj; 2Þ is such that the vector P Á t has for first coordinates the coprime numbers t 1 and t 2 . We thus return to the previous case. If we call M the two-block solution of that case, the solution of the problem is given by the matrix P À1 M. Note that P À1 ¼ Pðj; 2ÞPði; 1Þ 6 ¼ P. With t of coordinates ½12; 20; 225, the algorithm gives Now, let us consider the rarer cases in which none of the pairs ðt i ; t j Þ of coordinates of t are coprime despite the fact that the set of coordinates of t are coprime (as mentioned previously, if they are not, there is no solution to the problem). The set of integers ft i ; i ¼ 1; . . . ; Ng is said to be 'coprime but not pairwise coprime'. A classical example is {6, 10, 15}. Let us recall that in large dimensions N, the probability that a set of integers ft i g is coprime but not pairwise coprime is very small because the probability that two randomly chosen integers are coprime is quite high: it is equal to 1=½ð2Þ ¼ 6= 2 ' 61%, where refers to the Riemann zeta function (Wikipedia, 2021c). The exact calculation of the probability for a set of N integers ft i g to be coprime but not pairwise coprime as a function of N is not straightforward and is beyond the scope of the present study. Even if rare, these cases can be solved as follows. We consider the two first coordinates t 1 and t 2 of the vector t (any pair of coordinates would also work). As t 1 and t 2 are not coprime, they can be written t 1 ¼ xy and t 2 ¼ yz, where x, y, z are three integers and y ¼ gcd ðt 1 ; t 2 Þ > 1. It is important to note here that there is at least another coordinate t i with i > 2 that cannot be divided by y because if it were not so the set ft i g would not be coprime. We call ðu; vÞ the Bé zout numbers associated with ðt 1 ; t 2 Þ, ut 1 þ vt 2 ¼ y. The pair ðu; vÞ are also the Bé zout numbers associated with ðx; zÞ, ux þ vz ¼ 1, i.e. ðu; vÞ are also coprime. We call ð; Þ the Bé zout numbers associated with ðu; vÞ. We consider the matrix u v À ; its determinant is 1, and with k 2 Z. We build the N Â N matrix K from the 2 Â 2 block B and from the ðN À 2Þ Â ðN À 2Þ identity matrix. The first coordinate ðK Á tÞ ð1Þ of the new vector K Á t is coprime with at least one of the coordinates ðK Á tÞ ðiÞ with i > 2. It means that the method described in Section 3.2 can be applied to calculate a matrix L such that detðLÞ ¼ 1, and its last column is the vector K Á t. The matrix M ¼ K À1 L is then such that its determinant is also 1 and its last column is t. As the determinant of K is 1, K À1 is the adjugate (transpose of the cofactor matrix) of K, and is thus an integer matrix. Consequently, M is also an integer matrix, solution of the problem. The algorithm is effective and fast, whatever the dimension N of the vector t. The pseudocode is shown in Fig. 3.
We give an example with the classical set of coprime but not coprime coordinates [6,10,15]. The algorithm gives immediately a solution (the vectors are written in columns):

Hyperplane unit cell by oblique projection
Let us recall the hyperplane unit-cell problem. We are looking for a set of N vectors fb 1 ; . . . ; b j ; . . . ; b N g such that the first out-of-plane vector b 1 is such that p t b 1 ¼ 1 (pointing to a node of the layer q ¼ 1), and the N À 1 in-plane vectors b 2 ; . . . ; b j ; . . . ; b N are such that p t b j ¼ 0 (lying in the layer q ¼ 0). The solution b 1 of Bezout's identity is placed in the first position to be coherent with the notations we used in a separate paper dedicated to the lattice reduction (Cayron, 2021). The method we propose uses a short solution to the Ndimensional Bé zout identity (Section 2) and a solution of the CCUM problem (Section 3). We start from the input vector p. The coordinates of the vector b 1 pointing to a node of the layer q ¼ 1 are the solution of the Bé zout's identity associated with p. They are obtained by the algorithm detailed in Section 2 (method-1). Now, how to determine the N À 1 vectors in the layer q ¼ 0? We consider the unimodular matrix M that is such that the first column is the vector b 1 . The N À 1 next column vectors of the matrix M are called v j for j 2 f2; . . . Ng. Each of these vectors belongs to the lattice; thus, they verify p t v j ¼ q j 2 Z. The vectors b j ¼ v j À q j b 1 verify p t b j ¼ 0 for j 2 f2; . . . Ng; i.e. they are in-plane vectors lying in the layer q = 0. Geometrically, the vectors b j are obtained by oblique projection of the vectors v j along b 1 onto the plane p, as illustrated for N ¼ 3 in Fig. 4. Now, we have a cell U ¼ ðb 1 ; . . . ; b j ; . . . ; b N Þ attached to the plane p such that detðUÞ ¼ 1, p t b 1 ¼ 1 and p t b j ¼ 0 for i 2 f2; . . . Ng. It is thus the unit cell we were looking for. As the vectors used for the projection are short, the unit cell is not large. It can be reduced even more. There are different methods to find a reduced unit cell U 0 ¼ ðb 0 1 ; . . . ; b 0 j ; . . . ; b 0 N Þ, with b 0 j that have the same properties as b j with the vector p, but with shorter lengths and with angles between them closer to orthogonality. One could apply for example the LLL algorithm well known in computer science (Lenstra et al., 1982). We realized however that the algorithm developed in Section 4 can also be used to define the operation of 'hyperplane shearing' which consists of shearing the unit cell such that the vector b 1 becomes b 0 1 nearly normal to the plane p as illustrated in Fig. 4, and that this operation can be coupled with other lattice reduction methods to rival LLL implemented in Mathematica. The hyperplane reduction and its application to lattice reduction are detailed in a separate paper (Cayron, 2021). The pseudocode of the set of operations Bé zout-CCUM-Projection-Hyperplane reduction is shown in Fig. 5.
The program written in Python 3.8 called General-izedBezout incorporates the lattice reduction operation described by Cayron (2021). Let us give some examples we obtained: (i) With p t ¼ ð12; 20; 225Þ. The Bé zout vector associated with the plane p given by method-1 described in Section 2 is b 1 ¼ ½À17; À1; 1. After determining a first unit cell by projections along b 1 , and after lattice reduction, this vector becomes b 0 1 ¼ ½À7; À7; 1. The final reduced unit cell is given by the matrix Projections of the vectors v 2 ; v 3 along b 1 onto the plane p ¼ ðh; k; lÞ. The vector b 1 is a short solution of the Bé zout's identity p t 1 b 1 ¼ 1. The vectors v 2 ; v 3 are the solutions of the b 1 -CCUM problem. The vector b 1 is a short vector that can be further shortened into a vector b 0 1 by 'hyperplane shearing' (Cayron, 2021).
where the vectors are written in columns. The first vector is the short solution of the 3D Bé zout's identity and the other vectors are short solutions of the integer relations with the coordinates of p.
The matrix U 0 is interpreted crystallographically/geometrically as the unit cell attached to the hyperplane p. From an algebraic point of view, equivalently be understood as the infinite set of solutions of the N-dimensional Bé zout's identity, where b 0 1 is a short solution of the equation p t b 0 1 ¼ 1, and the other vectors are short solutions of the integer relation p t b 0 j ¼ 0, j 2 f2; . . . Ng. The set of solutions of Bé zout's identity is thus b 0 1 þ fZ: b 0 j g with j 2 f2; . . . Ng, where {Z:} means all the linear combinations with integer coefficients. This N À 1-dimensional affine space represents all the solutions of Bé zout's identity made on the coordinates of p.

Conclusion
The problem treated in the present paper called the 'hyperplane unit-cell problem' consists of finding, for any hyperplane p of N dimensions, one short vector b 1 that is such that p t b 1 ¼ 1 and N À 1 short integer vectors b 2 ; . . . ; b j ; . . . ; b N that are such that p t b j ¼ 0. The short outof-plane vector b 1 is the solution of Bezout's identity with p, and the short in-plane vectors b j , j ! 2, are solutions of the integer relation with p. These vectors constitute a unit cell attached to the hyperplane p. The algorithm to find a short solution to the N-dimensional Bé zout's identity is presented in Section 2. The algorithm to find a solution to a connected problem called the column-constrained unimodular matrix (CCUM) is detailed in Section 3. Both algorithms are then combined with the help of an oblique projection to determine a small unit cell attached to any hyperplane p (Section 4). The vectors b 1 ; . . . ; b N are short and can be further shortened by lattice reduction. We have shown in some examples that the solutions of the integer relation are even shorter than those determined by the PSLQ algorithm computed with Mathematica.

APPENDIX A Notations, Bézout's identity and integer relations
We note u i the ith coordinate of a vector u. Sometimes, the notation u ðiÞ will be equivalently used. It should not be confused with u i which is the ith vector in a set of vectors fu i g. The coordinates of a vector u are written in columns and those of a vector u t are in rows. From a crystallographic point of view, column and row vectors belong to direct and reciprocal spaces, respectively. The matrix multiplication notation is adopted. It means that even a 'simple' scalar (inner) product p Á u ¼ P i p i u i is written p t u where p t means transpose of p.
Bé zout's identity in 2D is an arithmetic theorem that states that for a and b which are integers with d for greatest common Pseudocode of the sequence of operations to determine a short unit cell associated with a hyperplane p. The unit cell is made of N short vectors U 0 ¼ ðb 0 1 ; . . . ; b 0 j ; . . . ; b 0 N Þ, such that p t b 0 1 ¼ 1 and p t b 0 j ¼ 0.
divisor, there exist integers u and v such that au þ bv ¼ d.
More generally, the linear combinations of the form au þ bv are exactly the multiples of d. It can be generalized to any dimension N and written as follows. For any integer vector p, calling d the greatest common divisor of the coordinates of p, there exist integer vectors u such that p t u ¼ P i p i u i ¼ d. In the paper, we suppose that the coordinates of p are coprime, i.e. d = 1.
An integer relation between a real vector x exists if and only if there is an integer vector u such that x t u ¼ P i x i u i ¼ 0. There are different algorithms to determine integer relations, such as the PSLQ (Ferguson et al., 1999). Searching for an integer relation between a set of powers of x {1, x, x 2 , . . . , x n } permits one to determine whether a given real number x is likely to be algebraic. Integer relations are also searched between some mathematical constants such as e, and ln(2) in order to establish new arithmetic conjectures. In the paper, only the cases where x is an integer vector (called p) are studied. (The equivalences between the mathematical and crystallographic terms used in the paper are given in Table 1.)