Solving p-adic polynomial systems via iterative eigenvector algorithms

In this article, we describe an implementation of a polynomial system solver to compute the approximate solutions of a 0-dimensional polynomial system with finite precision p-adic arithmetic. We also describe an improvement to an algorithm of Caruso, Roe, and Vaccon for calculating the eigenvalues and eigenvectors of a p-adic matrix.


INTRODUCTION
Let k be a field and let f 1 , . . ., f m ∈ k[x 1 , . . ., x n ] be polynomials such that the variety defined by f 1 , . . ., f m has dimension 0. It is often of interest to compute the solutions of such a system, either exactly or approximately.A significant bottleneck in computing the exact solutions of such a polynomial system is the complexity of the field extension required to write down all of the solutions.A popular method to compute the solutions exactly is to compute a triangular decomposition for the ideal f 1 , . . ., f m and solve for the coordinates via backsubstitution.A difficulty with this approach is the explosion of the size of the coefficients, which in conjunction with taking complicated field extensions, usually renders a computer algebra system unresponsive.
Computations in an exact field (such as Q) can often be approximated by computations in an inexact field (such as R, C), where by inexact, we mean that the computer representation of an element is only an approximation up to some precision.It is often far faster to compute an approximation to the solutions, though the trade-off is that the solutions are not known perfectly.In numerous applications, a real or complex approximation to the solutions is sufficient.
In Linear Algebra, p-adic methods go back to at least Dixon [Dix82].One advantage pointed out by Dixon is that a linear system can be ill-conditioned with respect to an archimedean norm, but well-conditioned p-adically; general polynomial systems exhibit this behavior as well.Recently, new developments in number theory and tropical geometry highlight that non-archimedean data is not just a route to integral solutions, but is significant in and of itself.
Our polynomial system solver is based on the truncated normal form solver of [vBMT18].To transplant their implementation to the p-adic setting, we rely on finite precision linear algebra over the p-adic field, which in analogy to the case over R we refer to as pnumerical linear algebra.It is essential for us to be able to compute the eigenvalues and eigenvectors of an approximate p-adic matrix quickly and stably.The naive algorithm to compute the eigenvalues of a matrix relies on the calculation and factorization of the characteristic polynomial.In the numerical setting (i.e, over R) it is desirable to avoid this step due to the instability of solving for the roots of a polynomial.There is a similar loss of precision in solving for the roots of a p-adic polynomial when the roots of the polynomial are p-adically close together.Unlike the archimedean setting, the accurate calculation of the characteristic polynomial is more difficult since a division free algorithm must be used [CRV17,Introduction].Our idea to improve on the algorithm of [CRV17] is to adapt the ideas from classical numerical linear algebra to the pnumerical setting; specifically to use an iterative scheme to compute the eigenvectors or eigenvalues.To our knowledge, this is the first appearance of a p-adic numerical algorithm based on iterating matrix multiplication to solve for eigenvalues and eigenvectors.
The significant deviation from the classical numerical setting is the preconditioning strategy.A general matrix over R tends to have eigenvalues of distinct absolute values, so aspects of the iteration are dominated by the unique largest eigenvalue.In the p-adic setting, it tends to be the case that all of the eigenvalues have identical absolute value, and this is true for most shifts of the form A → A − µI as well.Our strategy is to first approximate the eigenvalues by computing the characteristic polynomial over the residue field, then refine the initial approximation via iteration.The advantage of this approach is that the computation of the characteristic polynomial of a matrix in M n (F q ) is significantly faster than computing the characteristic polynomial over Q p [CRV17,Sto01].The basic Hensel-lifting based strategy applied to the system of equations Av = λv, v i = 1, with v i some fixed entry of v, would require the computation of the inverse of the Jacobian matrix of the map F : (v, λ) → (Av − λv, v i − 1) at each step in the iteration.Our cost per iteration step is much smaller.Unfortunately, matrices of the form µI + N , with N ∈ M n (Z p ) topologically nilpotent, do not allow us to apply our deflation strategy as the characteristic polynomial over the residue field is unhelpful; we hope that future improvements to our strategy can be made.
The structure of the article is as follows.In Section 2 we discuss some background and terminology regarding the subject of finite precision computation over Q p .We review some specific results on linear algebra over the p-adic field, especially the p-adic analogues of condition-stable algorithms based on matrix factorizations.We describe our iterative strategy for the computation of the eigenvectors of a matrix over Q p and compare it to the existing alternatives.In Section 3, we describe our adaptations of the algorithm of [vBMT18] to solve polynomial systems over Q p .Our algorithms are available as Julia packages, and can be obtained from: https://github.com/a-kulkarn/pAdicSolverhttps://github.com/a-kulkarn/Dory(For technical reasons related to the Julia package system, we choose to divide our implementation between two packages.)2. LINEAR ALGEBRA 2.1.Precision, norms, and condition numbers.We state some basic definitions for our discourse.We direct the reader to [CRV15,Car17,Ked10] for more details.For a matrix A ∈ M n (Z p ), we will denote by χ A the characteristic polynomial and denote χ A,p := χ A (mod p).We denote the standard basis of Q n p (or Z n p ) by {e 1 , . . ., e n }.
Definition 2.2.Let A ∈ Q n×m p be a matrix.The operator norm of A with respect to • p is The condition number for an invertible square matrix is κ(A) := A p • A −1 p .The condition number for a singular matrix is ∞.
We can identify a subgroup of GL n (Q p ) where every matrix is well-conditioned, serving the analogous role to O n (R) in the real setting.We now come to the discussion of p-adic precision.There are many possible ways to represent a p-adic element a ∈ Q p in a computer system [Car17].In the Nemo/Hecke systems [FHHJ17], an element is represented by a series a = a −r p −r + . . .
where the O(p N ) is the p-adic ball representing the uncertainty of the remaining digits.The relative precision of a is the quantity N +r, and the absolute precision is the number N .In the terminology of [Car17], we consider a system with the zealous (i.e, interval) implementation of arithmetic.The operations −, + preserve the minimum of the absolute precision of the operands, and ×, ÷ preserve the minimum relative precision of the operands.If u ∈ Z × p , a ∈ Z p , and N ≤ N ′ , then we have that Multiplication by p preserves the relative precision and increases the absolute precision by 1.The worst operation when it comes to absolute p-adic precision is dividing a small number by p.For example, the expression begins with 3 numbers with an absolute and relative precision of at least 100, and ends with a result where not even the constant term is known.
In Nemo/Hecke, elements of a matrix store their own precision.In the terminology of [CRV17], the associated precision structure is the jagged precision lattice.However, we will not keep track of this finer precision here.Instead, in the terminology of [CRV17], we will look at the flat precision of the matrix: Definition 2.4.Let A, B ∈ M n (Z p ) be matrices such that a i,j = b i,j + O(p Ni,j ).Then we write A = B + O(p N ), where N := min i,j N i,j .
To refer to a matrix A ∈ M n (Q p ) whose elements are known at an absolute precision at least N , we will simply write A + O(p N ).
2.1.1.Design choices for the implementation.We discuss some of the design choices of the computer algebra package accompanying this article.The user we had in mind while designing our system is one who is interested in data that is inherently p-adic.More specifically, we have assumed the following hypotheses: Assumptions.
• Input and output are inherently approximate.
• The desired output precision is roughly the input precision.
• Element-wise gains in absolute precision after a → pa or a → a p are usually lost in a subsequent step.
We believe that our linear algebra package will be used in the middle of an ongoing p-adic computation, such as solving polynomial systems of equations.In such a circumstance, it does not make sense for our algorithm to "request" extra digits of precision from the input data since the input itself is approximate.We also assume that if a user inputs data at precision N , they expect output at precision roughly N , or at least N minus the condition number for the problem.Computations using relaxed p-adic numbers [Car17] allow the desired output precision to be specified.Unfortunately, relaxed arithmetic does not seem well-suited for use in an iterative scheme due to the extra memory costs (see [Car17, Section 2.4]).We also assume for a dense matrix that the entries are all given at the same flat absolute precision.We believe that this situation is fairly common in practical applications, but of course not ubiquitous.That said, the implementation of p-adic matrices in [FHHJ17] uses an elementwise precision model, so it is possible for a computation using our software to validly return a result with more accuracy than we indicate here.
In our analysis, we do not account for the extra precision on an element a ∈ Q p gained after a → pa or a → a p .Experimentally, this seems to be a reasonable assumption for randomly selected dense matrices over a p-adic field, with p ≥ 40.For a more structured scenario this assumption might not be reasonable.A nice feature of capped precision is that we do not have to account for the growth of arithmetic costs due to increased precision.
With these assumptions made, it is often useful to think of the computer as performing arithmetic in a Z/p N Z-module for a computation with no divisions by p.It is useful to consider the image of a matrix A ∈ M n (Z p ) in the ring M n (Z/p N Z).It is an essential difficulty of finite precision p-adic linear algebra that Z/p N Z is a local principal ideal ring, but not a domain when N > 1.
Remark 2.5.Our eigenvector algorithms are focused on finding only the eigenvectors defined over Q p .We made this choice for the following two reasons.First, in the last section of the article, we focus on finding the Q p -solutions of a 0-dimensional system of polynomial equations over Q p .Secondly, our present implementation is written using components of the OSCAR system (specifically [FHHJ17]) that are currently under development.At the time of writing the first version of our implementation, the "q-adics" interface or an interface to handle ramifed extensions of Q p did not exist.Thus, our article is written with Q p in mind to more closely reflect the actual software.Rapid progress is being made on the system, so we expect our implementation to evolve over time as well.

Matrix factorizations.
In this section, we very briefly review some matrix invariants and factorizations.Matrix factorizations relating to the topological structure of Q p are unsurprisingly related to the algebraic structure of the Z p -module spanned by the columns (or rows) of the matrix.We invite the reader to consider [Ked10, Chapter 4] for further details.
Proposition 2.6 (Iwasawa decomposition).Let k be either a finite extension of Q p , R, or C and let G ֒−→ GL n (k) be a linear algebraic group, let K be a maximal compact subgroup of G, and let B be a Borel subgroup.For any A ∈ G, there exists a Q ∈ K and an R ∈ B such that A = QR.
Note that a maximal compact subgroup of GL n (R) is the orthogonal group O n (R), and the subgroup of invertible real upper triangular matrices is a Borel subgroup.We have as a consequence: Corollary 2.7 (QR-factorization).Let A ∈ R n×m be a matrix.Then there exists a Q ∈ O n (R) and an upper triangular matrix R ∈ R n×m such that A = QR.
Of course, QR-factorization is really just a consequence of Gram-Schmidt.From the point of view of Proposition 2.6, we also have a p-adic analogue of QR-factorization.In GL n (Q p ), a maximal compact subgroup is GL n (Z p ), and the subgroup of invertible upper triangular matrices is a Borel subgroup.We decided to refer to the factorization above as a QR-factorization since many applications are similar to the real case.However, one important difference is that Q is generally not an element of O n (Q p ).As in the case over R, an elementary proof of Corollary 2.8 can be given by exhibiting an algorithm to compute a QR-decomposition.The algorithm is well-known, and in fact it is just the algorithm to compute a P LUdecomposition, though with each row-pivot chosen by taking the vector with the largest p-adic norm [Ked10, Chapter 4].We also note since Q ∈ GL n (Z p ), that R can be chosen to be the Hermite normal form of A. We summarize the discussion as: Remark 2.9.The p-adic QR-decomposition of A is computed using a P LU -decomposition, with pivots chosen by p-adic norm.
An essential property of the QR-factorization of a real matrix A is that Q has a well-conditioned inverse.The same is true p-adically.The columns of Q in the p-adic case are also "orthogonal" in the sense of [Sch84, Section 50].Unfortunately, this notion of orthogonality does not help to compute the inverse of Q like in the real case.As we might expect, QR-factorizations give us a way to compute the singular value decomposition of a matrix in the p-adic setting.
The singular value decomposition is the Smith normal form of the Z p -module generated by columns of A [Ked10, Chapter 4].Since the matrices U, V are well-conditioned, we can use a singular value decomposition to stabilize pnumerical computations (as in the real setting).

2.3.
The eigenvector problem: semantics.Let A ∈ M n (Z p ) be a matrix.In the setting of a finite p-adic precision computation, we may formulate two versions of the eigenvector problem: Problem 1 (Eigenvector problem, Version I).Let (λ, v) be an exact eigenpair for A over Z p .Given an approximation A to A such that Problem 2 (Eigenvector problem, Version II).Let (λ, v) be an exact eigenpair for A over Z p .Given an approximation A to A such that (ii) for any B = 0 + O(p N ), we have that M is maximal among such triples.
In Version II, the second criteria ensures that the problem is well-posed.
Lemma 2.11.Let A ∈ M n (Z p ), let A be an approximation such that A = A + O(p N ), and let (λ, v, M ) satisfy the first two criteria for a solution to Version II.Then M ≤ N . Proof.
The following example demonstrates the difference between Version I and Version II of the problem.
Example 2.12.Let p = 7, and consider the matrices Both matrices are rounded to the same result with precision N := 6.However, the respective eigenvectors are Version II exhibits pathological behavior whenever the characteristic polynomial of A is not square-free modulo p M .For instance, as X ranges over all Z n p matrices, we see that the eigenvectors of I + p M X could be anything.In other words, any solution to Version II of the problem for I + p M X has M ≤ 0. A similar difficulty arises in the eigenvector problem for real matrices.Because of this pathology, we restrict our attention to solving Version I.
Of course, we can give analogous formulations for computing a Jordan form of A or for computing a block Schur form for A. Generally, we will work with Version I of these problems as well.
2.4.Commentary on the classical algorithm.The naive algorithm to compute the eigenvalues of a matrix relies on the calculation and factorization of the characteristic polynomial.In the numerical setting over R it is desirable to avoid this step due to the instability of solving for the roots of a polynomial.Also p-adically, there is also a loss of precision in solving for the roots of a polynomial.Thus, in order to accurately compute the nullspace of A − λI the roots of the characteristic polynomial, and thus the characteristic polynomial itself, must be known to as much precision as possible.Presently, there are two algorithms to compute the characteristic polynomial at the maximum possible precision, both of which use more than O(n 3 ) arithmetic operations (asymptotically): (1) Use a division-free calculation.The asymptotic run-time is at least O(n 4 ).
(2) Use the algorithm of [CRV17].Compute the characteristic polynomial of a matrix with inflated precision using division, then truncate to the optimal precision computed from the precision lattice analysis.The runtime is O(ρn 3 ), where ρ represents the factor of performing arithmetic at the higher precision.Practically, the computation of the comatrix and the precision increase required for the precision analysis appears to impose a large cost [CRV17, Remark 5.3].
In the next subsection, we discuss how to avoid an expensive computation of the characteristic polynomial for most matrices in M n (Z p ).
Remark 2.13.To clarify the term "most", if A + O(p N ) is an n × n-matrix whose entries are chosen with the uniform probability distribution on [0, . . ., p N − 1], then the limit as n → ∞ of the probability that χ A is square-free is at least 1−p −5 1+p −3 [Ful02].

Algorithm 2.14 Eigvecs(A)
Input: A + O(p N ), an n × n-matrix over Z p .Output: The eigenpairs (λ i , v i ) of the matrix A.
Reset precision of V return Eigvecs(X i ) for i = 1, . . ., r 14: else 15: return ClassicalAlgorithm(A) Note in step 6, we may reset the precision since the input matrix on the previous step X was a multiple of p ν .If p is moderately large, then for most matrices the size of the Jordan blocks modulo p will be small (see Remark 2.13).The square matrix in step 13 is expected to be small, so we apply the costly method to compute the characteristic polynomial of this block.

Algorithm 2.15 PowerIterationDecomposition(A, χ A,p ) Input:
A + O(p N ), an n × n-matrix over Z p .χ A,p , the characteristic polynomial of A (mod p).
1: Compute eigenpairs (λ i , V i ) mod p 2: Set m to be the largest multiplicity of a linear factor of χ A,p .
Lift the eigenpair to an approximate eigenpair (λ i , V i ) in Z p 5: Set B = (A − λ i I).
To clarify, the lifting in step 4 is just the trivial operation of interpreting the elements {0, . . ., p − 1} as elements of Z p .2.6.Proof of correctness.We give a quick proof of correctness of the power iteration algorithm.
Lemma 2.16.Let λ be an approximate eigenvalue of A (mod p), let m be the multiplicity of λ as a root of χ A,p , and let M be the nullspace of (A − λI) N m .Then M is a free Z/p N Z-module with trivial annihilator.Letting V be a matrix whose columns are the generators of M , there exists an X ∈ M N (Z p ) with such that Proof.First, note that (A − λI) (mod p) has 0 as an eigenvalue with multiplicity m.Ordering the singular values and eigenvalues by size, we have that (A − λI) satisfies Thus, by the Hodge-Newton decomposition [Ked10, Theorem 4.3.11],there is a matrix U ∈ GL n (Z p ) such that B := U −1 (A − λI)U is block upper triangular, with the lower right block B 22 accounting for the m small eigenvalues and the upper left block B 11 accounting for the unit eigenvalues.Now, we see that B m 22 = 0 (mod p), so But as U ∈ GL n (Z p ), we see that the Smith normal form of (A − λI) N m is I 0 0 0 .Thus, the kernel of (A − λI) N m is free as a Z/p N Z-module with trivial annihilator in Z/p N Z.Next, note that we have so AM ⊆ M .Thus A restricted to M defines an endomorphism of M , represented by a square matrix.
As an immediate corollary, we have: Corollary 2.17.Let (V 1 , X 1 ), . . ., (V r , X r ) be the blocks returned from power iteration.With the horizontal join, we have We hope that step 15 can be improved, as we believe it is helpful in solving polynomial systems whose solutions have large valuations.Thus, we raise the question: Question 2.18.Given a topologically nilpotent matrix B ∈ M n (Z p ), is there an efficient iterative strategy to determine the eigenvectors of B + O(p N )? 2.7.The Schur form algorithm.Solving for the nullspace at the end of the iteration is the dominant cost in the power iteration algorithm.In the classical numerical setting, we can avoid this drawback by using the QR-iteration algorithm.We also have a p-adic version of the QR-algorithm.With the p-adic QR step replacing the usual QR factorization, our algorithm is simply the original LR-algorithm [Rut58] that inspired the QRalgorithm for real matrices.The LR-algorithm fell out of favor some time in the 1960's since it is numerically less stable than the QR-algorithm, but p-adically this situation is reversed!
We now give the algorithm.

the characteristic polynomial of A (mod p).
Output: A (block) triangular form T for A, and matrix V such that AV = V T + O(p N ).
1: Set λ 1 , . . ., λ ℓ to be the roots of χ A,p in F p , lifted to Z p .2: Set m 1 , . . ., m ℓ to be the multiplicities of the roots of χ A,p .3: Compute B, V such that AV = V B and B is in Hessenberg form.4: for i = 1, . . ., ℓ do 5:

Algorithm 2.20 BlockSchurForm(A)
Input: A + O(p N ), an n × n-matrix over Z p .Output: A block upper triangular matrix T with ℓ + 1 distinct blocks and matrix V such that AV = V T + O(p N ).Here, ℓ is the number of linear factors of χ A,p = χ A (mod p).
Reset precision of V Update B and V 19: return A, V We briefly comment on some aspects of this algorithm analogous to the archimedean setting (the proofs for the various statements are analogous as well).Note that each p-adic QR-iteration preserves the Hessenberg form, so the cost of applying the QR-step is O(n 2 ).From the discussion of the LR-algorithm in [Wil65, Sections 5-9], we see that the subdiagonal entry (i, i + 1) converges to 0 at a rate of O , where m is the largest multiplicity of an eigenvalue of B (mod p).Additionally, if some subdiagonal entries are identified to be zero during the iteration, the standard deflation strategies can be used to accelerate the algorithm.
Remark 2.21.In our version of the iteration we use a constant shift factor, and so the iteration converges in at most mN steps [Wil65].We believe that a p-adic analogue of the Rayleigh quotient strategy to refine the eigenvalue approximation will lead to a drastic speed-up in convergence.
Remark 2.22.It is possible to use the QR-iteration to compute the valuations of the eigenvalues iteratively.After applying a single round of the QR-iteration to a Hessenberg matrix A with no shift (i.e, λ i = 0), the valuations of the eigenvalues can be read from the diagonal blocks in the resulting matrix.If B ii is one of the diagonal blocks with no subdiagonal entry equal to 0, then all of the eigenvalues of B ii must have the same valuation, as otherwise some of the subdiagonal entries of B ii would have converged to 0. In fact, the valuations of the eigenvalues of A can be computed this way even if some of the eigenvalues are not defined over Q p .
2.8.Complexity analysis.We denote by M (n) the cost of matrix multiplication.
Proposition 2.23.Let N be the precision and let m be the largest size of a Jordan block of A (mod p).Then the number of Q p -arithmetic operations performed by Algorithm 2.14 is as most Proof.We tabulate the costs in a O(ℓn 3 + M (n) log 2 mN ).
The ρm 3 term corresponds to blocks which cannot be decomposed using power iteration, i.e, where the classical algorithm is invoked.Note that the cost of computing χ A (mod p) is negligible, since over a finite field the characteristic polynomial can be computed in O(M (n)(log n)(log log n)) finite field operations [Sto01].
If the linear factors occur with small multiplicities mod p, our algorithm significantly outperforms the naive strategy.If all linear factors mod p occur with multiplicity 1 the classical algorithm is never called, so the asymptotic run-time simplifies to O(n 3 + ℓM (n) log 2 (N )).The worst-case scenario of the algorithm is if χ A (mod p) ≡ (x − α) n , in which case the default algorithm is called immediately; the only wasted time in this case is the inexpensive computation of χ A (mod p).
We now study the complexity of the block Schur form algorithm.
Proposition 2.24.Let N be the precision, let m be the largest size of a Jordan block of A (mod p), and let ℓ be the number of linear factors of χ A (mod p).Then the block Schur form can be computed in O(n 3 + ℓmN n 2 + ρm 3 ) Q p -arithmetic operations.
Proof.We tabulate the costs in a O(n 3 + ℓmN n 2 ) 2.9.p-adic Householder reflections.We close the section with a quick observation about Householder reflections.
Lemma 2.25.Let p be an odd prime.Let x ∈ Z n p be a vector with exactly one coordinate with minimal valuation r.We further assume this coordinate is not x 1 .Let e 1 be the first standard basis vector.Then: (a) x T x is a non-zero square in Z p .(b) Let α := √ x T x ∈ Z p be one of the square roots and let v ′ := x − αe 1 .Choose (the unique) v ∈ Z n p such that v ′ = p r v. Then v T v p = 1, and the Householder transformation Proof.(a) Assume x i is the coordinate with minimal valuation.Write x = p r y.We have that x T x = p 2r (y T y) = p 2r (y 2 i + O(p)).As y i is a unit, we see by Hensel's lemma that this is a square.(b) Assume again x i is the coordinate with minimal valuation.Note that α 2 = n j=1 x 2 j , so as x i is the unique coordinate with minimal valuation, we have α = x i + O(p r+1 ).Now By our assumption on x, we have x i − x 1 = x i + O(p r+1 ).Thus, v T v is a unit and H is well-defined.That it is a reflection follows from the usual calculation, which is Finally, it is clear that all of the entries of H are integral.(c) The result is clear from direct calculation.
We now offer a Householder version of the Hessenberg algorithm.Our version of the Hessenberg algorithm is more a theoretical curiosity; for our practical implementations, we use [CRV17, Algorithm 1] as it appears to be more efficient.Swap the last (i.e,n-th) row with row i to ensure a nj has the minimal valuation of the sub-diagonal elements 5: Add multiples of the last (n-th row) to the i-th row, for j < i < n, to ensure a nj is the unique subdiagonal element with minimal valuation 6: Apply the corresponding column operations to preserve similarity.None of these operations affect columns 1, . . ., j

7:
Apply the Householder reflection to set the j-th subcolumn to a multiple of e 1 .Theorem 3.5.Let k be a field, let R := k[x 1 , . . ., x n ], and let I be an ideal such that dim k R/I is finite.Let p 1 , . . ., p ℓ be the geometric points of the affine scheme Spec(R/I) with multiplicities m 1 , . . ., m ℓ respectively.Let [f ] : R/I → R/I be the endomorphism of R/I induced by multiplication by f ∈ R. Then the eigenvalues of [f ] and their multiplicities are (f (p 1 ), m 1 ), . . ., (f (p ℓ ), m ℓ ).The i-th invariant subspace of [f ] is {g ∈ R/I : g(p j ) = 0 for j = i}.
Proof.Let k al be an algebraic closure for k.We may assume that R/I ⊗ k k al is local by the Chinese Remainder Theorem.Note R/I ⊗ k k al is Artinian, so for m the maximal ideal we have m(m r ) = m r for some r > 0. By Nakayama's lemma we have m r = 0. Thus, for any f ∈ m, [f ] is nilpotent and so has 0 as an eigenvalue with multiplicity dim k al (R/I ⊗ k k al ).Since R/I ⊗ k k al = k al ⊕ m as a k al -vector space, we now need only check the result for constant functions.For f = 1, the result is clear.
Note that if k is not algebraically closed, the eigenvalues and invariant subspaces may only be defined over a finite extension of k; these correspond to points whose coordinates lie in a finite extension of k.Critically, the eigenvalues of [x i ] are the i-th coordinates of the geometric points.We now give the p-adic solver algorithm: Algorithm 3.6 SolveQpSystem(f 1 , . . ., f m ) Input: Polynomials f 1 , . . ., f m ∈ Q p [x 1 , . . ., x n ] defining a 0-dimensional polynomial system.Output: A set of approximate solutions to the distinct Q p -solutions.The differences between the Q p and R versions of the algorithm are precisely the adaptations in the linear algebra, i.e, steps 2, 4, and 6.Note that if k is any finite precision field for which these steps can be implemented, the general recipe of [vBMT18, Algorithm 1] allows us to compute the approximate solutions.Examples include k = R, C, Q p , F q ((t)).
Finally, instead of computing a truncated normal form in steps 1 and 2, we may simply use a Gröbner basis instead.For systems where a Gröbner basis is already known or where the polynomials in the basis have small coefficients, the numerous optimizations of Gröbner basis algorithms are likely preferable from a practical standpoint.The major advantage of our implementation is the fact that the coefficient size is capped at p N .

ACKNOWLEDGEMENTS
I would like to thank Yue Ren, David Roe, and Simon Telen for helpful discussions about the article.I would also like to thank Marta Panizzut, Emre Sertöz, and Bernd Sturmfels for their helpful comments.
Corollary 2.8 (p-adic QR-factorization).Let A ∈ Z n×m p be a matrix.Then there exists a Q ∈ GL n (Z p ) and an upper triangular matrix R ∈ Z n×m p such that A = QR.

1:
Construct the resultant map [Res ij ] with D sufficiently large 2: Compute [π ij ] = (ker[Res ij ] T ) T 3: Extract the submatrix M of [π ij ] whose columns correspond to monomials m such that deg m < D. 4: Compute a square subblock of M that is pnumerically well conditioned, using the p-adic QR-factorization.Let b be the basis corresponding to this sub-block.5: Construct the multiplication by x i matrices [x 1 ] b , . . ., [x n ] b from the columns of [π ij ] b .6: Compute the eigenvectors of a random linear combination via iteration.7: return The eigenvalues for [x 1 ] b , . . ., [x n ] b indexed by the invariant subspaces.By [x i ] b , we mean the matrix corresponding to the multiplication-by-x i operator in b-coordinates.The choice of B := Span(b) in step 4 defines a subspace of V D of the same dimension as R/I, hence a section s : R/I → V D and thus a truncated normal form.Since [x i ]B ⊆ V D , we can easily read off the action of the operator [x i ] : B → B from the columns of [π ij ].For explicit details see [vBMT18, Theorem 4 and following text]. table. table.