Hensel Lifting
\gdef\p#1{\left({#1}\right)} \gdef\setn#1{\mathcal{#1}} \gdef\vec#1{\mathbf{#1}} \gdef\mat#1{\mathrm{#1}} \gdef\F{\mathbb{F}} \gdef\Z{\mathbb{Z}} \gdef\vd{π }
There are many variants of Hensel's lemma and the one derived here is for quadratic lifting of multivariate roots. This is useful solving systems of polynomial equations in modular and Galois rings. One starts by finding a solution in the residue field where other methods may be available (e.g. exhaustive search). Once found, these residue field solutions are 'lifted' to progressively larger rings until a solution in the target ring is found. See Mur05, Gro24 and Con24 for alternative lemmas and further elaboration of the multivariate case.
Algebraic Preliminaries
A ring is a set of numbers with operations of addition and multiplication (in this post assumed finite, unital and commutative). There is also a notion of division, but unlike fields this does not necessarily work for all non-zero numbers. An ubiquitous example of a ring is \Z_{2^{64}}, the ring of 64-bit unsigned integers that all computers use. In this ring, multiplicative inverses exist only for odd numbers.
An ideal for a ring R is a subset I that is closed under addition and multiplication, including multiplication by any element of R. If the values of I are all the R-multiples of a single element p, we write I = (p) and called it a principal ideal. The trivial ideals are (0) = \{0\} and (1) = R. The ring \Z_{2^{64}} has the ideals
(1), (2), (2^2), (2^3), β¦, (2^{63}), (0).
We can also define the product of two ideals as I β J = \{i β j \mid i \in I, j \in J\}. The product of two ideals is itself an ideal. We say an ideal I squares to zero if I β I = (0). This just means that for each a, b \in I we have a β b = 0. In \Z_{2^{64}} the ideals (2^k) for k β₯ 32 square to zero.
A quotient ring R/I is the set of all cosets r + I = \{r + i \mid i \in I\} with addition and multiplication defined as (r + I) + (s + I) = (r + s) + I and (r + I) β (s + I) = (r β s) + I. The quotient ring is a ring itself, and captures the idea of 'modulo I' arithmetic. In fact, \Z_{2^{64}}/(2^k) = \Z_{2^k} is the ring of integers modulo 2^k.
An interesting theorem is that the quotient ring R/I is a field if and only if I is maximal, meaning there is no ideal J such that I β J β R. In the ring \Z_{2^{64}} the ideal (2) is maximal, producing the field \Z_{2^{64}}/(2) = \F_2.
Ideal Taylor Expansion
Hensel lifting is analogous to Newton's method, and it starts with an analogue of Taylor expansion:
Lemma 1. Given ring R with ideal I that squares to zero and polynomial f \in R[X] then for x \in R and Ξ΄ \in I, we have
\begin{equation} f(x + Ξ΄) = f(x) + f'(x) β Ξ΄. \end{equation}
Proof: Assume f(x) = x^m and consider the binomial theorem for f(x + Ξ΄):
(x + Ξ΄)^m = \sum_{i \in [0,m]} \binom{m}{i} β x^{m-i} β Ξ΄^i
the terms i > 2 cancel due to Ξ΄^2 = 0, leaving only
x^m + m β x^{m-1} β Ξ΄ = f(x) + f'(x) β Ξ΄.
This proofs the case for f(x) = x^m. By linearity it thus holds for all polynomials f \in R[x].
Curiously this result is exact. Where in real numbers truncating the Taylor expansion results in an error with some bound, in the present setting the result is exact. An explanation for this can be found in the theory of p-adic numbers, where closeness is defined in terms of how small the ideal is that contains the difference. This can be made rigorous using a counterintuitive but entirely valid alternative norm, but such a treatment is out of scope here.
The lemma generalizes to the multi-variate case as one would expect:
Lemma 2. Given ring R with ideal I that squares to zero and polynomial map f: R^n \rightarrow R^m then for \vec x \in R^n and \vd \in I^n, we have
\begin{equation} f(\vec x + \vd) = f(\vec x) + \mat J_f (\vec x) β \vd \end{equation}
where the Jacobian \mat J_f: R^n \rightarrow R^{m Γ n} is given by
\begin{equation} \mat J_{i j} (\vec x) = \frac{\partial f_i}{\partial x_j} (\vec x). \end{equation}
Proof: Consider coordinates f_i and apply the multi-binomial theorem analogous to lemma 1.
Lifting Roots
Lifting is the process of taking a solution from a quotient ring R/I to R itself:
Definition 3. (Lift) Given ring R with ideal I, polynomial map f: R^n β R^m, and \vec x \in R^n such that
f(\vec x) = \vec 0 \pmod I
then lifts of \vec x are \vec x' \in R^n such that
\begin{aligned} \vec x' &= \vec x \pmod I \\ f(\vec x') &= \vec 0. \end{aligned}
Lemma 4. (Hensel Lifting) Suppose the ideal I squares to zero, then lifts of \vec x are given by \vec x' = \vec x + \vd with \vd \in I^n solutions of
\begin{equation} -f(\vec x) = \mat J _f (\vec x) β \vd. \end{equation}
In particular if \mat J _f (\vec x) is invertible there is a unique solution given by
\begin{equation} \vec x' = \vec x - \mat J_f^{-1} (\vec x) β f(\vec x). \end{equation}
Proof: Follows directly from lemma 2 and f(\vec x + \vd) = \vec 0.
As familiar from linear algebra, \mat J_f (\vec x) is invertible if and only if its determinant is invertible. It follows that if the determinant is invertible \operatorname{mod} I, it is also invertible in R. So if we have a chain of rings R/(p^k), and \mat J_f (\vec x) is invertible in the first, it is invertible in all rings. In this case the solution can be uniquely lifted to all the rings.
Implementation and Optimization
As an optimization, we often do not need to solve equation (4) in R but can instead pick a smaller ring:
Lemma 5. Given principal ideal I=(p) that squares to zero and ideal J such that I β J = (0), equation (4) can be solved in R/J by setting \vd = \vd' β p where \vd' is the solution of
\begin{equation} -\frac{f(\vec x)}{p} = \mat J _f (\vec x) β \vd' \pmod J. \end{equation}
Proof: Equation 4 is linear with both sides in I. As an R-module the ideal I is isomorhpic to R/J by the map x \mapsto x / p and the result follows.
Typically we use Hensel lifting over a ring R with ideals (p^s) for s \in [0,k] where (p) is maximal and (p^k) = (0). We first find an initial solution in the residue field R/(p) and then lift it to R/(p^2) using I = J = (p). Then we can either continue lifting to R/(p^4), R/(p^16), β¦, R/(p^k) = R using I_i = J_i = \p{p^{2^i}}. Or alternatively, we can fix J = R/(p) and lift to R/(p^3), R/(p^4), β¦ with I_i = (p^i). The former is known as quadratic lifting and the latter linear lifting. The quadratic method grows faster, but requires arithmetic in progressively larger rings. The linear method requires more iteration, but all arithmetic is in the residue field. Which one is faster depends on the specifics, but typically the quadratic method is preferred.
Other Functions
Note that the machinery of Hensel lifting works for any function f for which an f' can be defined such that lemma 1 holds. In particular it holds for rational functions with the familiar derivative:
Lemma 6. Given f(x) = p(x)/q(x) such that lemma 1 holds for p and q, then lemma 1 holds for f with derivative
\begin{equation} f'(x) = \frac{p'(x) β q(x) - p(x) β q'(x)}{(q(x))^2}. \end{equation}
Proof: Expand definitions,
f(x + Ξ΄) = \frac{p(x + Ξ΄)}{q(x + Ξ΄)} = \frac{p(x) + p'(x) β Ξ΄}{q(x) + q'(x) β Ξ΄}
drop the argument (x) for brevity,
\frac{p + p' β Ξ΄}{q + q' β Ξ΄} = \frac{(p + p' β Ξ΄)β (q - q' β Ξ΄)}{(q + q' β Ξ΄) β (q - q' β Ξ΄)} = \frac{pβ q + (p'β q - pβ q')β Ξ΄}{q^2}
and thus
f(x + Ξ΄) = \frac{p}{q} + \frac{p'β q - pβ q'}{q^2} β Ξ΄ = f(x) + f'(x) β Ξ΄.
With this we can construct an efficient method for computing inverses in R by lifting the root a^{-1} of f(x) = x^{-1} - a. This results in lifting equation x' = x β (2 - a β x), which is the basis of the fastest know inversion methods for \Z_{2^{64}} (see Rey22 for an overview). Implemented in Rust with further optimizations this looks like:
fn ring_inv(a: u64) -> u64 {
let mut x = a.wrapping_mul(3) ^ 2; // Initial guess correct to 5 bits
let mut y = 1.wrapping_sub(a.wrapping_mul(x));
for _ in 0..3 { // Successively lift to 10, 20 and 40 bits
x = x.wrapping_mul(y.wrappining_add(1));
y = y.wrapping_mul(y);
}
x.wrapping_mul(y.wrappining_add(1)) // Final lift to 80 bits (only 64 computed)
}
