COUNTING POINTS ON VARIETIES OVER FINITE FIELDS OF SMALL CHARACTERISTIC

ALAN G.B. LAUDER AND DAQING WAN

Abstract. We present a deterministic polynomial time algorithm for computing the zeta function of an arbitrary variety of fixed dimension over a finite field of small characteristic. One consequence of this result is an efficient method for computing the order of the group of rational points on the Jacobian of a smooth geometrically connected projective curve over a finite field of small characteristic.

1. Introduction

The purpose of this paper is to give an elementary and self-contained proof that one may efficiently compute zeta functions of arbitrary varieties of fixed dimension over finite fields of suitably small characteristic. This is achieved via the p-adic methods developed by Dwork in his proof of the rationality of the zeta function of a variety over a finite field [4, 5]. Dwork’s theorem shows that it is in principle possible to compute the zeta function. Our main contribution is to show how Dwork’s trace formula, Bombieri’s degree bound [3] and a semi-linear reduction argument yield an efficient algorithm for doing so. That p-adic methods may be used to efficiently compute zeta functions for small characteristic was first suggested in [20, 22], where Wan gives a simpler algorithm for counting the number of solutions to an equation over a finite field modulo small powers of the characteristic.

We now give more details of our results. For q = p^{a}, where p is a prime number and a a positive integer, let 𝔽_{q} denote a finite field with q elements. Let 𝔽̄_{q} denote an algebraic closure of 𝔽_{q}, and 𝔽_{q^{k}} the subfield of 𝔽̄_{q} of order q^{k}. Denote by 𝔽_{q}[X_{1}, …, X_{n}] the ring of all polynomials in n variables over 𝔽_{q}.

Date: October 28, 2002.

Alan Lauder gratefully acknowledges the support of the EPSRC (Grant GR/N35366/01) and St John’s College, Oxford, and thanks Richard Brent. Daqing Wan is partially supported by the NSF and the NSFC. The authors are pleased to thank Colin McDiarmid and Bernd Sturmfels for answering some questions on convex geometry. Many excellent suggestions were made for improving the paper by the anonymous referee, and these were incorporated by the authors during the revision. They are especially grateful for this help.

Mathematics Subject Classification 2000: 11Y16, 11T99, 14Q15 Key words and phrases:

variety, finite field, zeta function, algorithm.

For a polynomial f ∈ 𝔽_{q}[X_{1}, …, X_{n}], we denote by N_{k} the number of solutions to the equation f = 0 with coordinates in 𝔽_{q^{k}}. The zeta function of the variety defined by f is the formal power series in T with non-negative integer coefficients

Z(f/𝔽_{q})(T) = exp (∑_{k=1}^{∞} \frac{N_{k}T^{k}}{k} ) .

Dwork’s theorem asserts that Z(f/𝔽_{q}) is a rational function r(T)/s(T) with integer coefficients. From this it follows that knowledge of explicit bounds deg(r) ≤ D_{1} and deg(s) ≤ D_{2}, and of the values N_{k} for k = 1, 2, …, D_{1}+D_{2} is enough to efficiently determine Z(f/𝔽_{q}), see [22]. The Bombieri degree bound tells us that deg(r) + deg(s) ≤ (4d + 9)^{n+1}, see [3, 22], and so in particular we may take D_{1} = D_{2} = (4d + 9)^{n+1}. Each number N_{k} can be computed in a naive fashion by straightforward counting, using q^{nk} evaluations of the polynomial f. Thus one may compute the zeta function of a variety, but the naive method described requires a number of steps that is exponential in the parameters d^{n} and log q, where d is the total degree of f. Note that the dense input size of f is O((d + 1)^{n} log q), and the size of the zeta function is polynomial in O((d + 1)^{n} log q), by the Bombieri degree bound. We prove the following theorem. Theorem 1. There exist an explicit deterministic algorithm and an explicit polynomial P such that for any f ∈ 𝔽_{q}[X_{1}, …, X_{n}] of total degree d, where q = p^{a} and p is prime, the algorithm computes the zeta function Z(f/𝔽_{q})(T) of f in a number of bit operations which is bounded by P (p^{n}d^{n^{2}}a^{n}).

In particular, this computes the zeta function of a polynomial in a fixed number of variables over a finite field of “small characteristic” in deterministic polynomial time. Note that our result makes no assumption of non-singularity on the variety defined by the polynomial. Also, we shall explicitly describe all the algorithms in this paper, rather than just prove their existence, and we assume that the finite field 𝔽_{q} itself is presented as input via an irreducible polynomial of degree a over the prime field 𝔽_{p}, as explained in Section 3.

With regard to the exponents in the algorithm, we state these precisely in Theorem 37. For now we observe that if p, d and n are fixed then the time required to compute the number of points in 𝔽_{q} is O(a^{3n+7}), with space complexity O(a^{2n+4}). Here we are ignoring logarithmic factors (see also Proposition 36). Note that all our complexity estimates are made using standard methods for multiplication in various rings, and can be modestly reduced with faster methods.

We also present refinements to this result based upon the ideas of Adolphson and Sperber [2], and indeed from the outset will follow their approach, as it involves little extra complication. This refinement takes into account the terms which actually occur in the polynomial f rather than working solely with the total degree. We shall need more definitions: the support of a polynomial f is the set of exponents r = (r_{1}, …, r_{n}) of non-zero terms a_{r}X_{1}^{r_{1}} … X_{n}^{r_{n}} which occur in f, thought of as points in ℝ^{n}. The Newton polytope of f is defined to be the convex hull in ℝ^{n} of the support of f. Our refined version of Theorem 1 essentially replaces the parameter d^{n} with the normalised volume of the Newton polytope of f (see Proposition 35 and Section 6.3.3).

Zeta functions may also be defined for a finite collection of polynomials, and we next describe how our results may be extended to this case. An affine variety V over 𝔽_{q} is the set of common zeros in (𝔽̄_{q})^{n} of a set of polynomials f_{1}, f_{2}, …, f_{r}. An analogous zeta function Z(V/𝔽_{q}) may then be defined in terms of the number of solutions in each finite extension field of 𝔽_{q}. These numbers may be computed using an inclusion-exclusion argument involving the polynomials ∏_{i∈S} f_{i} where S is a subset of {1, 2, …, r}, see [22]. Theorem 1 then easily yields the following.

Corollary 2. There exist an explicit deterministic algorithm and an explicit polynomial Q with the following property. Let f_{1}, …, f_{r} ∈ 𝔽_{q}[X_{1}, …, X_{n}] have total degrees d_{1}, …, d_{r} respectively, where q = p^{a} and p is prime. Denote by V the affine variety defined by the common vanishing of these polynomials, and define d = ∑_{i=1}^{r} d_{i}. The algorithm computes the zeta function Z(V/𝔽_{q})(T) in a number of bit operations which is bounded by Q(p^{n}d^{n^{2}}a^{n}2^{r}).

Thus one has an efficient algorithm for computing the zeta function of an arbitrary affine variety over 𝔽_{q} assuming the characteristic, dimension and the number of defining polynomials are fixed. More generally still, an arbitrary variety is defined through patching together suitable affine varieties. Zeta functions for such general varieties may be defined. These zeta functions may be computed using the above ideas provided explicit data are given on how to construct them from affine patches. As an example, the zeta functions of arbitrary projective varieties or toric varieties may be computed in this way.

Conceptually our algorithm is rather straightforward. The zeta function can be expressed in terms of the “characteristic power series” (Fredholm determinant) of a certain “lifting of Frobenius” which acts on an infinite dimensional p-adic Banach space constructed by Dwork. Under modular reduction, we obtain an operator acting on a finite dimensional vector space. Unfortunately, this operator cannot be computed efficiently directly from its definition; however, it can be expressed as a product of certain semi-linear operators each of which can be computed efficiently if the characteristic p is small. This last step is thus of crucial importance in deriving an efficient algorithm. Note that the same idea is used in [20, 22] in a simpler situation. In more concrete language, the algorithm requires one to construct a certain “semi-linear” finite matrix and compute the “linear” matrix which is the product of the Galois conjugates of the semi-linear matrix. The number of points is then read off from the trace of the final linear matrix. The zeta function can be computed from the characteristic polynomial of the final linear matrix. The semi-linear matrix itself is defined over a certain finite “p-adic lifting” of the original finite field.

In the literature algorithms have already been developed for computing zeta functions of curves and abelian varieties [1, 6, 14, 17, 18]. These utilise the theory originally developed by Weil for abelian varieties, whereas we use Dwork’s more general and simpler p-adic theory. For example, in the case of a smooth geometrically irreducible projective plane curve of degree d over a field of size q = p^{a} these algorithms have a time complexity which grows as (log q)^{C_{d}}, where the constant C_{d} grows exponentially in the degree d. Given an absolutely irreducible bivariate polynomial f of degree d over 𝔽_{q}, the zeta function of the unique smooth projective curve birational to the affine curve defined by f may be computed in time polynomial in d, p and a using our approach. Thus our more general method is far better in terms of the degree d if the characteristic p is small, since the running time has polynomial growth in d (but much worse if p is large). The zeta function immediately gives the order of the group of rational points on the Jacobian, see [22].

Corollary 3. There exist an explicit deterministic algorithm and an explicit polynomial R with the following property. Let V be a geometrically irreducible affine curve defined by the vanishing of polynomials f_{1}, …, f_{r} ∈ 𝔽_{q}[X_{1}, …, X_{n}] of total degrees d_{1}, …, d_{r} respectively, where q = p^{a} and p is prime. Denote by Ṽ the unique smooth projective curve birational to the affine curve V, and let d = ∑_{i=1}^{r} d_{i}. The algorithm computes the order of the group of rational points on the Jacobian of Ṽ in a number of bit operations bounded by R(p^{n}d^{n^{2}}a^{n}2^{r}).

Thus one may compute the order of the group of rational points on the Jacobian of a smooth geometrically irreducible projective curve over a finite field of small characteristic in deterministic polynomial time, provided the number of variables and the number of defining equations are fixed. (Notice that the case r = 1 and n = 2 corresponds to that of being given a possibly singular plane model of the curve.) In particular, this answers a question posed in [15], which Poonen attributes to Katz and Sarnak. We note that recently a similar result for special classes of plane curves was independently obtained in [8, 10], using Monsky-Washnitzer’s method. Also, a different padic approach for elliptic curves has been developed in [16].

This paper is written primarily for theoretical computational interest, in obtaining a deterministic polynomial time algorithm for computing the zeta function in full generality if p is small. It can certainly be improved in many ways for practical computations. What we have done in this paper is to work on the easier but more flexible “chain level”. A general improvement in the smooth case, is to work on the cohomology level. Then there are several related p-adic cohomology theories available, each leading to a somewhat different version of the algorithm. Again, as indicated in [20, 22], these p-adic methods are expected to be practical only for small p. (The authors have recently applied Dwork’s cohomology theory to derive a practical algorithm in a special case [12].)

2. Additive character sums over finite fields

The most natural objects of study in Dwork’s theory are certain additive character sums over finite fields. In this section we introduce these sums, and explain their connection to varieties.

An additive character Ψ is a mapping from 𝔽_{q^{k}} to the group of units of some commutative ring S with identity 1 such that Ψ(x + y) = Ψ(x)Ψ(y) for x, y ∈ 𝔽_{q^{k}}. We say that it is non-trivial if Ψ(x) ≠ 1 for some x ∈ 𝔽_{q^{k}}. Let {Ψ_{k}}_{k≥1} be any family of mappings with each Ψ_{k} a non-trivial additive character from 𝔽_{q^{k}} to some extension ring of the integers whose image is a group of order p with elements summing to zero. We assume that the family {Ψ_{k}} forms a tower of characters in the sense that for each k ≥ 2,

Ψ_{k} = Ψ_{1} ∘ Tr_{𝔽_{q^{k}} /𝔽_{q}}.

In our application, the ring S will be taken to be a certain p-adic ring.

For the remainder of the paper we will use multi-index notation. Specifically, we let X^{u} represent the monomial X_{0}^{u_{0}}X_{1}^{u_{1}} … X_{n}^{u_{n}} for an integer vector u = (u_{0}, u_{1}, …, u_{n}); let x be an (n + 1)-tuple (x_{0}, x_{1}, …, x_{n}) of field elements; and X the list of indeterminates X_{0}, X_{1}, …, X_{n}. Observe here that we have introduced an extra indeterminate X_{0}. Even though the polynomial f whose zeta function we wish to compute is in the n variables X_{1}, …, X_{n}, in Dwork’s theory the extra indeterminate arises naturally, as we are about to see.

Lemma 4. Let f ∈ 𝔽_{q}[X_{1}, …, X_{n}] and N_{k}^{∗} denote the number of solutions to the equation f = 0 in the affine torus (𝔽_{q^{k}}^{∗})^{n}. Then

∑_{x∈(𝔽_{q^{k}}^{∗})^{n+1}} Ψ_{k}(x_{0}f(x_{1}, …, x_{n})) = q^{k}N_{k}^{∗} − (q^{k} − 1)^{n}

where x = (x_{0}, x_{1}, …, x_{n}).

Proof. We have that for any u ∈ 𝔽_{q^{k}}

∑_{x_{0}∈𝔽_{q^{k}}} Ψ_{k}(x_{0}u) = \{ \begin{matrix}0 & if u ∈ 𝔽_{q^{k}}^{∗} \\ q^{k} & if u = 0. \end{matrix}

This is a standard result from the theory of additive character sums [13, page 168]. Thus ∑ Ψ_{k}(x_{0}f(x_{1}, …, x_{n})), where the sum is taken over points x ∈ 𝔽_{q^{k}} × (𝔽_{q^{k}}^{∗})^{n}, equals q^{k}N_{k}^{∗}. Removing the contribution of (q^{k} − 1)^{n} from the terms with x_{0} = 0 in this sum gives the required result. □

Our next step will be to find an alternative formula for the lefthand side of the equation in Lemma 4. This is achieved in Proposition 11, which leads us eventually to Dwork’s Trace Formula (Theorem 26).

3. p-Adic Theory

3.1. p-adic rings. We first introduce notation for the p-adic rings we shall need, before explaining how to construct them and compute in them. Let ℚ_{p} be the field of p-adic rationals, and ℤ_{p} the ring of p-adic integers (see [11]). Denote by Ω the completion of an algebraic closure of ℚ_{p}. Select π ∈ Ω with π^{p−1} = −p and define R_{1} = ℤ_{p}[π], a totally ramified extension of ℤ_{p} of degree p − 1. By binomial expansion and Hensel’s lifting lemma, one sees that the equation (1 + πt)^{p} = 1 has exactly p distinct solutions t in R_{1}. In particular, R_{1} contains all p-th roots of unity. The motivation behind the introduction of R_{1} is that to define an additive character of order p, we need a small p-adic ring which contains a primitive pth root of unity.

Let R_{0} denote the ring of integers of the unique unramified extension of ℚ_{p} in Ω of degree a, where q = p^{a}. Finally, let R be the compositum ring of R_{1} and R_{0}. We have the diagram of ring extensions

Ω | R

⧸  ⧹
R_{0}  R_{1}
⧹  ⧸
ℤ_{p}

The residue class ring of R_{0} is 𝔽_{q}, and we shall “lift” the coefficients of the polynomial f ∈ 𝔽_{q}[X_{1}, …, X_{n}] to this characteristic zero ring, using the Teichmüller lifting. Precisely, the Teichmüller lift ω(x) of a non-zero element x ∈ 𝔽̄_{q} is defined as the unique root of unity in the maximal unramified extension of ℚ_{p} which is congruent to x modulo p and has order coprime to p. We define ω(0) = 0. Thus the compositum R contains both the lifting of the coefficients of f and the image of an additive character we shall construct.

3.2. Algorithmic aspects.

3.2.1. Construction and lifting Frobenius. We assume that 𝔽_{q} is presented as the quotient 𝔽_{p}[y]/(h) where h(y) is a monic, irreducible polynomial of degree a over the prime field 𝔽_{p}. For any positive integer N we describe how the quotient ring R/(p^{N}) may be constructed: Lift the polynomial h to an integer polynomial ĥ whose coefficients lie in the open interval (−p/2, (p + 1)/2). Take the unramified extension R_{0} to be ℤ_{p}[μ] = ℤ_{p}[y]/(ĥ). Elements in R_{0}/(p^{N}) can now be represented as linear combinations over ℤ_{p}/(p^{N}) of the basis elements 1, μ, …, μ^{a−1}. (Recall that ℤ_{p}/(p^{N}) can be identified with ℤ/(p^{N}).) The extension R/(p^{N}) = R_{0}[π]/(p^{N}) is easily constructed by adjoining an element π and specifying the relation π^{p−1} = −p.

We define a lifting of the Frobenius automorphism on 𝔽_{q} to an automorphism of R which is the identity on R_{1}. Define the map τ : R → R by setting τ(μ) to be the unique root of ĥ which is congruent to μ^{p} modulo p. Define τ(π) = π, and extend to the whole of R by insisting τ is an automorphism. We extend τ to act on R[[X]] coefficient-wise, fixing monomials. Here R[[X]] is the ring of all power series in the indeterminates X = X_{0}, …, X_{n} with coefficients from R.

3.2.2. Complexity of arithmetic. In this section we bound the complexity of the basic arithmetic operations in the ring R/(p^{N}), along with that of computing the map τ and Teichmüller lifts. Note that these estimates are all the simplest possible, and can be improved using more advanced methods. The reader may wish to skip the proof of the next lemma, and refer back when required in Section 6.3.2. Lemma 5. Elements in R/(p^{N}) can be represented using O(paN log p) bits. Addition and subtraction can be performed in O(paN log p) bit operations, and multiplication and inversion of units in O((paN log p)^{2}) bit operations. The Teichmüller lifting to R/(p^{N}) of a finite field element can be computed in O((a log p)^{3}N^{2}) bit operations. For any 1 ≤ i ≤ a − 1, the map τ^{i} on R/(p^{N}) may be evaluated using O(p(aN log p)^{2}) bit operations. (For the powers of τ we require a total of O(a^{4}N^{2}(log p)^{3}) bits of precomputation.)

Proof. Elements in R/(p^{N}) can be written as

(1)  ∑_{i=0}^{a−1} ∑_{j=0}^{p−2} c_{ij}π^{j}μ^{i}

where the coefficients c_{ij} belong to the ring ℤ_{p}/(p^{N}) of size p^{N}. The bit size of such an expression is a(p − 1) log(p^{N}) = O(paN log p). Addition and subtraction are straightforward, just involving the addition of integers and reduction modulo p^{N}. Likewise, multiplication of two expansions of the form (1) is straightforward, using the reduction relations π^{p−1} = −p and ĥ(μ) = 0.

For Teichmüller lifting and inversion we shall use Newton iteration with quadratic convergence. Specifically, we shall use Newton lifting with respect to the prime π in R, and define l = (p − 1)N so that π^{l} = (−p)^{N}. Given a polynomial ϕ(Y) ∈ R[Y] and an element g_{0} ∈ R such that ϕ(g_{0}) ≡ 0 mod π and ϕ′(g_{0}) is invertible (modulo π) with inverse s_{0} modulo π, this algorithm computes an element g ∈ R/(π^{l}) such that ϕ(g) ≡ 0 mod π^{l} and g ≡ g_{0} mod π (compare with [7, Algorithm 9.22]). For i ≥ 1, assuming that g_{i−1} and s_{i−1} have been found, we define g_{i} = g_{i−1} − ϕ(g_{i−1})s_{i−1} mod π^{2^{i}} and s_{i} = 2s_{i−1} −ϕ′(g_{i})s_{i−1}^{2} mod π^{2^{i}}. As in [7, Theorem 9.23] one checks that g_{i} ≡ g_{0} mod π, ϕ(g_{i}) ≡ 0 mod π^{2^{i}} and s_{i} ≡ ϕ′(g_{i})^{−1} mod π^{2^{i}} at the ith step.

Thus after ⌈log_{2}(l)⌉ steps we shall have found the required approximate root g = g_{⌈log_{2}(l)⌉}. The ith step involves three additions/subtractions and three multiplications in the ring R/(π^{2^{i}}), which requires O((2^{i}a log p)^{2}) bit operations, along with evaluation of the polynomials ϕ and ϕ′ modulo π^{2^{i}}. Let c(ϕ, i) be the complexity of these latter operations. Thus the total complexity is O((paN log p)^{2} + c(ϕ)) bit operations where c(ϕ) = ∑_{i=1}^{⌈log_{2}(l)⌉} c(ϕ, i). Note that to compute approximate roots in the ring R_{0} rather than R, one can lift using the prime p rather than π and obtain a complexity of O((aN log p)^{2} + c(ϕ)).

Suppose now we are given a unit x ∈ R/(p^{N}). We compute a Newton lifting starting from an element g_{0} ∈ R/(π) such that xg_{0} ≡ 1 mod π; that is, we use the equation ϕ(Y) = xY − 1. Now R/(π) is just the finite field 𝔽_{q} and so an inverse g_{0} of x modulo π can be computed in O((a log p)^{2}) bit operations [7, Corollary 4.6]. (Here ϕ′(Y) = x and s_{0} = g_{0}, and so in fact we only need to iterate the formula for s_{i}.) In this case c(ϕ, i) is just one multiplication and a subtraction modulo π^{2^{i}}. By the above paragraph the total complexity is then O((paN log p)^{2}) for inversion. For Teichmüller lifts we use the same approach, only with the polynomial ϕ(Y) = Y^{q−1} − 1 and lifting in R_{0} via the element p. Here using a fast exponentiation routine we find that c(ϕ, i) involves O(log q) multiplications and a subtraction in R_{0}/(p^{2^{i}}). Thus c(ϕ) = O((aN log p)^{2} log q) which gives the Teichmüller lifting estimate.

We precompute τ on the basis elements 1, μ, …, μ^{a−1}. To do this, recall that τ(μ) ∈ R_{0} is defined as the unique root of ĥ which is congruent modulo p to μ^{p}. This may be approximated modulo p^{N} by Newton lifting in R_{0} with respect to p using the polynomial ϕ(Y) = ĥ(Y) and the initial value g_{0} = μ^{p}. (Finding μ^{p} takes O(a^{2}(log p)^{3}) bit operations and this is absorbed in the stated precomputation estimate.) Here c(ϕ, i) is O(a) additions and multiplications in R_{0}/(p^{2^{i}}), using Horner’s method for polynomial evaluation [7, Page 93]. Thus the Newton lifting estimate gives a complexity of O(a(aN log p)^{2}) bit operations. Using τ(μ^{i}) = (τ(μ))^{i} one can now find the image of all basis elements in a further O(a(aN log p)^{2}) bit operations. One stores this information as a matrix for τ acting on R_{0}/(p^{N}) as an ℤ_{p}/(p^{N})-module with basis the powers of μ. The map τ can now be computed on any element in R_{0}/(p^{N}) in O(a^{2}(N log p)^{2}) bit operations using linear algebra. By taking powers of the matrix, matrices for the maps τ^{i} for 1 ≤ i ≤ a−1 can also be found in O(a^{3}a(N log p)^{2}) bit operations. Thus the total precomputation is bounded by O(a^{4}N^{2}(log p)^{3}), and each evaluation of τ^{i} on R_{0}/(p^{N}) takes O((aN log p)^{2}) bit operations. Finally, to evaluate τ on R/(p^{N}) one writes elements of R on the R_{0}-basis 1, π, …, π^{p−2} and applies τ component-wise. □

3.3. p-adic valuations and convergence of power series. Denote by ord the additive valuation on Ω normalised so that ord (p) = 1. Thus ord (π) = 1/(p − 1). Define a p-adic norm |.|_{p} on Ω by |x|_{p} = p^{−ord (x)}. The set of all x ∈ Ω with |x|_{p} ≤ 1 (equivalently ord (x) ≥ 0) is called the closed unit disk. Given any formal power series ∑_{r} A_{r}X^{r} where A_{r} ∈ Ω, we say it converges at a point x = (x_{0}, …, x_{n}) ∈ Ω^{n+1} if the sequence of partial sums ∑_{r, |r|<e} A_{r}x^{r} tends to a limit under the p-adic norm. (Here |r| = ∑_{i=0}^{n} r_{i}.) This sequence of partial sums will converge if and only if the summands A_{r}x^{r} tend to zero p-adically (that is, are divisible in the ring of integers of Ω by increasingly large powers of p), as |r| goes to infinity. In particular, if there is a real number c > 0 such that ord (A_{r}) ≥ c|r| for all r, then the series will certainly converge for all points x which are Teichmüller liftings of points over 𝔽̄_{q}. Note that throughout the paper we shall use the additive valuation ord rather than the p-adic norm |.|_{p} itself.

4. Analytic representation of characters

4.1. Dwork’s splitting functions. We now need to find a suitable p-adic expression for a non-trivial additive character from 𝔽_{q} to R_{1}. In the case of complex characters, this is done via the exponential function. However, the radius of convergence of the exponential function in Ω is not large enough; in particular, it does not converge on the Teichmüller lifting of all the points in 𝔽_{q}. Instead we use the power series constructed by Dwork using the exponential function. (The reader may find the discussion of the related Artin-Hasse function on [11, Pages 92-93] helpful.)

Let μ be the Möbius function. Taking the logarithmic derivative, one checks that the exponential function has the following product expansion

exp(z) = ∑_{k=0}^{∞} \frac{z^{k}}{k!} = ∏_{k=1}^{∞}(1 − z^{k})^{−μ(k)/k}.

This can be rewritten as

exp(z) = ∏_{(k,p)=1}^{∞}(1 − z^{k})^{−μ(k)/k}(1 − z^{kp})^{μ(k)/(kp)}.

It follows that

(2)  exp (z + \frac{z^{p}}{p} ) = ∏_{(k,p)=1}^{∞}(1 − z^{k})^{−μ(k)/k}(1 − z^{kp^{2}})^{μ(k)/(kp^{2})}.

Replacing z by πz in the above relation and noting that π^{p} = −pπ, we define a power series in z by

θ(z) = exp(πz − πz^{p}).

Writing θ(z) = ∑_{r=0}^{∞} λ_{r}z^{r} we see that λ_{r} = π^{r}/r! for r < p, and we shall shortly show that all λ_{r} lie in R_{1}. From (2) we get the product expansion

(3)  θ(z) = ∏_{(k,p)=1}^{∞}(1 − π^{k}z^{k})^{−μ(k)/k}(1 − π^{kp^{2}}z^{kp^{2}})^{μ(k)/(kp^{2})}.

By the binomial expansion, the first factor

(1 − π^{k}z^{k})^{−μ(k)/k} = ∑_{j=0}^{∞}(−1)^{j}\binom{−μ(k)/k}{j}π^{jk}z^{jk} = ∑_{j=0}^{∞} b_{j}(k)z^{kj}

is a power series in z^{k}. Now for (k, p) = 1, we have [11, Page 82]

−μ(k)/k ∈ ℤ_{p}, \binom{−μ(k)/k}{j} ∈ ℤ_{p}.

Thus for j > 0 the coefficient of z^{jk} satisfies

ord (b_{j}(k)) ≥ \frac{jk}{p − 1} > \frac{p − 1}{p^{2}}jk.

Similarly, for the second factor, we write

(1 − π^{kp^{2}}z^{kp^{2}})^{μ(k)/(kp^{2})} = ∑_{j=0}^{∞}(−1)^{j}\binom{μ(k)/kp^{2}}{j}(πz)^{jkp^{2}} = ∑_{j=0}^{∞} c_{j}(kp^{2})z^{jkp^{2}}.

Now for j > 0 we have ord (j!) < \frac{j}{p−1}, see [11, Page 79]. Thus for j > 0 and k positive and coprime to p, the coefficient of z^{jkp^{2}} satisfies

ord (c_{j}(kp^{2})) > \frac{jkp^{2}}{p − 1} − 2j − \frac{j}{p − 1} ≥ \frac{p − 1}{p^{2}}jkp^{2}.

Putting the above two inequalities together, we conclude that for r > 0 the coefficients λ_{r} of θ(z) satisfy

(4)  ord (λ_{r}) > \frac{(p − 1)r}{p^{2}} , λ_{r} ∈ R_{1}.

This shows that the power series θ(z) is convergent in the disk |z|_{p} < 1 + ϵ for some ϵ > 0. In particular, θ(z) converges on the closed unit disk, and Definition 6 makes sense.

In the proof of the next lemma we shall use the fact that ord (λ_{r}) ≥ 2/(p − 1) for r ≥ 2, and so

(5)  θ(1) ≡ 1 + (πz) mod (πz)^{2}.

This can be seen as follows: The Artin-Hasse exponential function [11, Page 93]

E(z) := ∏_{(k,p)=1}^{∞}(1 − z^{k})^{−μ(k)/k} = exp (z + \frac{z^{p}}{p} + \frac{z^{p^{2}}}{p^{2}} + …)

has coefficients in ℤ_{p}, because each factor in the product expansion does. Since E(πz) ≡ θ(z) mod z^{p^{2}} we see

ord (λ_{r}) ≥ \frac{r}{p − 1}

for 0 ≤ r < p^{2}. This estimate combined with (4) gives (5).

Definition 6. [Dwork’s splitting function] Let

Φ_{k}(z) = ∏_{i=0}^{ak−1} θ(z^{p^{i}}) ∈ R_{1}[[z]],

and

Ψ_{k} = Φ_{k} ∘ ω : 𝔽_{q^{k}} → R_{1},

where ω is the Teichmüller map. (Recall that q = p^{a}.)

(That Ψ_{k} has image in R_{1} can be seen as follows: Let R_{1}[ω(𝔽_{q^{k}})] be R_{1} adjoined the image of ω on 𝔽_{q^{k}}. Then R_{1}[ω(𝔽_{q^{k}})] is an unramified extension of R_{1} of degree k. The Galois group of the corresponding quotient field extension is generated by τ. The map τ acts on a Teichmüller point ω(x) as τ(ω(x)) = ω(x)^{p}. Hence it fixes the element Ψ_{k}(x) for x ∈ 𝔽_{q^{k}}, and so Ψ_{k}(x) ∈ R_{1}.) Lemma 7. The maps Ψ_{k} form a tower of non-trivial additive characters from the fields 𝔽_{q^{k}} to the ring R_{1}.

Proof. (This is the case “s = 1” on [5, Pages 55-57].) We first show that θ(1) is a primitive pth root of unity. By (5) we see θ(1) ≠ 1. As a formal power series in z,

θ(z)^{p} = exp(pπz) exp(−pπz^{p}).

Now, θ(z), exp(pπz) and exp(−pπz^{p}) are all convergent in |z|_{p} < 1 + ε for some ε > 0. We can thus substitute z = 1 and find that θ(1)^{p} = exp(pπ) exp(−pπ) = 1. Thus θ(1) is a primitive p-th root of unity in R_{1}.

Next, for γ ∈ R_{0} with γ^{p^{ak}} = γ, we claim that

∏_{i=0}^{ak−1} θ(γ^{p^{i}}) = θ(1)^{γ+γ^{p}+···+γ^{p^{ak−1}}} .

Using (5) it is clear that both sides are congruent to

1 + π(γ + γ^{p} + ⋯ + γ^{p^{ak−1}})

modulo π^{2}. To prove the claim, it remains to prove that both sides are p-th roots of unity. The right side is a p-th root of unity since θ(1) is a p-th root of unity. The p-th power of the left side is

∏_{i=0}^{ak−1} θ(γ^{p^{i}})^{p} = exp(pπ ∑_{i=0}^{ak−1}(γ^{p^{i}} − γ^{p^{i+1}})) = exp(pπγ) exp(−pπγ^{p^{ak}}) = 1.

Thus, the left side is also a p-th root of unity. The claim is proved. Note that the individual factor θ(γ) is not necessarily a p-th root of unity. We conclude that

Ψ_{k}(x) = θ(1)^{Tr_{k}(x)}

for any x ∈ 𝔽_{q^{k}}, where Tr_{k} is the trace function from 𝔽_{q^{k}} to 𝔽_{p}, and the exponent is thought of as an integer. □

Note 8 The infinite sum

exp(π(z − z^{p})) = ∑_{k=0}^{∞} \frac{(π(z − z^{p}))^{k}}{k!}

is convergent for |z|_{p} < 1, but not necessarily convergent for |z|_{p} = 1. If |z|_{p} = 1, then it is possible that |z − z^{p}|_{p} = 1 and for such z the above infinite sum does not converge. Since the infinite sum does not converge everywhere on the disk |z|_{p} ≤ 1, one cannot simply substitute z = 1 into the above infinite sum and get the contradiction that θ(1) = 1. There is no contradiction here!

We now define another power series related to our original polynomial f whose relevance will become apparent in Proposition 11. Definition 9. Let f be the polynomial whose zeta function we wish to compute, and write

X_{0}f = ∑_{j∈J} ā_{j}X^{j},

where J is the support of X_{0}f. Let a_{j} be the Teichmüller lifting of ā_{j}. Let F be the formal power series in the indeterminates X with coefficients in R given by

F (X) = ∏_{j∈J} θ(a_{j}X^{j}).

Let F^{(a)}(X) be the formal power series in the indeterminates X with coefficients in R given by

F^{(a)}(X) = ∏_{j∈J} ∏_{s=0}^{a−1} θ((a_{j}X^{j})^{p^{s}}).

The relation between F (X) and F^{(a)}(X) is clear. Lemma 10. Let the power series F^{(a)} and F be as in Definition 9. Then

F^{(a)}(X) = ∏_{i=0}^{a−1} τ^{i}(F (X^{p^{i}}))

where the map τ acts coefficient-wise on the power series F.

The power series F^{(a)}(X) relates to rational point counting in the following way. Proposition 11. Let f ∈ 𝔽_{q}[X_{1}, …, X_{n}] and F^{(a)} ∈ R[[X]] be as in Definition 9. Denoting by N_{k}^{∗} the number of solutions to the equation f = 0 in the affine torus (𝔽_{q^{k}}^{∗})^{n}, we have

q^{k}N_{k}^{∗} − (q^{k} − 1)^{n} = ∑_{x^{q^{k}−1}=1} F^{(a)}(x)F^{(a)}(x^{q}) … F^{(a)}(x^{q^{k−1}})

where the summation is over the Teichmüller lifting of points on the torus (𝔽_{q^{k}}^{∗})^{n+1}.

Proof. For any point x̄ in (𝔽_{q^{k}})^{n+1} with Teichmüller lifting x we have that

Ψ_{k}(x̄_{0}f(x̄_{1}, …, x̄_{n})) = Ψ_{k}(∑_{j∈J} ā_{j}x̄^{j})
= ∏_{j∈J} Ψ_{k}(ā_{j}x̄^{j})
= ∏_{j∈J} Φ_{k}(a_{j}x^{j})
= ∏_{j∈J}{∏_{i=0}^{ak−1} θ((a_{j}x^{j})^{p^{i}})}
= ∏_{j∈J}{∏_{i=0}^{k−1} ∏_{s=0}^{a−1} θ((a_{j}x^{j})^{q^{i}p^{s}})}
= ∏_{i=0}^{k−1}{∏_{j∈J}(∏_{s=0}^{a−1} θ((a_{j}^{q^{i}}x^{jq^{i}})^{p^{s}}))}
= ∏_{i=0}^{k−1}{∏_{j∈J}(∏_{s=0}^{a−1} θ((a_{j}(x^{q^{i}})^{j})^{p^{s}}))}
= F^{(a)}(x)F^{(a)}(x^{q}) … F^{(a)}(x^{q^{k−1}}),

where F^{(a)}(X) is given by

F^{(a)}(X) = ∏_{j∈J} ∏_{s=0}^{a−1} θ((a_{j}X^{j})^{p^{s}}).

(We pause to justify the steps above: the first four equalities follow straight from definitions and from the homomorphic property of Ψ_{k}; the fifth and sixth by rearrangement; and the seventh since a_{j} satisfies a_{j}^{q} = a_{j}.)

Thus we have

∑_{x̄∈(𝔽_{q^{k}}^{∗})^{n+1}} Ψ_{k}(x̄_{0}f(x̄_{1}, …, x̄_{n})) = ∑_{x^{q^{k}−1}=1} F^{(a)}(x)F^{(a)}(x^{q}) … F^{(a)}(x^{q^{k−1}})

where the latter sum is over the Teichmüller lifting in Ω^{n+1} of points in (𝔽_{q^{k}}^{∗})^{n+1}. Combining this with Lemma 4 gives us the result. □

4.2. Decay rates and weight functions. We now describe the decay rates of the coefficients of the power series F^{(a)} and F. Specifically, we obtain lower bounds for the p-adic order of the coefficients of the power series F expressed in terms of a certain weight function on integer vectors.

Write F = ∑_{r} F_{r}X^{r} where the sum is over non-negative integer vectors in ℤ_{≥0}^{n+1}. Let A be the (n + 1) × |J| matrix whose columns are j = (j_{0}, j_{1}, …, j_{n}) ∈ J. Then from Definition 9 one sees

(6)  F_{r} = ∑_{u}(∏_{j∈J} λ_{u_{j}}a_{j}^{u_{j}})

where the outer sum is over all |J|-tuples u = (u_{j}) of non-negative integers such that

Au = r, (7)

thinking of u and r as column vectors. Since j_{0} = 1 for all (j_{0}, …, j_{n}) ∈ J, the first row of the matrix A is the vector (1, 1, …, 1). The first equation in the above linear system is then

∑_{j∈J} u_{j} = r_{0}. (8)

Now F_{r} is zero if (7) has no solutions. Otherwise, since ord (λ_{u_{j}}a_{j}^{u_{j}}) = ord (λ_{u_{j}}) from (4),(6) and (8) we get

ord (F_{r}) ≥ inf_{u} \{∑_{j∈J} \frac{(p − 1)u_{j}}{p^{2}} \} = \frac{p − 1}{p^{2}}r_{0}, (9)

where the inf is over all non-negative integer vector solutions u of (7). We now define a weight function w such that (9) gives estimates on ord (F_{r}) in terms of this weight function.

Let δ_{1} ⊂ ℝ^{n} denote the convex hull of the support of f (the set of exponents of non-zero terms). Let δ_{2} ⊂ ℝ^{n} be the convex hull of the origin and the n points (d, 0, …, 0), (0, d, …, 0), …, (0, …, 0, d) where d is the total degree of f. We call δ_{1} the Newton polytope of f; the polytope δ_{2} is just a simplex containing δ_{1}. Definition 12. Let δ be any convex polytope with integer vertices such that δ_{1} ⊆ δ ⊆ δ_{2}. Denote by Δ the convex polytope in ℝ^{n+1} obtained by embedding δ in ℝ^{n+1} via the map x ↦ (1, x) for x ∈ ℝ^{n}, and taking the convex hull with the origin. Denote by C(Δ) the cone generated in ℝ^{n+1} as the positive hull of Δ. So C(Δ) is the union of all rays emanating from the origin and passing through Δ.

Ultimately, in Sections 6.3.3 and 6.4, we shall only be interested in the simplest choice of polytope δ = δ_{2} (although the choice δ = δ_{1} leads to the most refined algorithm). Letting Δ_{1} denote the polytope in ℝ^{n+1} obtained by choosing δ = δ_{1} we see that C(Δ_{1}) is the cone generated by the exponents of non-zero terms in X_{0}f. Equation (7) has no non-negative integer (or even real) solutions when r does not lie in C(Δ_{1}). Thus for any choice of δ(⊇ δ_{1}) and corresponding Δ(⊇ Δ_{1}), all exponents of F lie in the cone C(Δ).

Definition 13. Define a weight function w from ℝ^{n+1} to ℝ ∪ {∞} in the following way: For r = (r_{0}, r_{1}, …, r_{n}) ∈ ℝ^{n+1} define

w(r) = \{ \begin{matrix}r_{0} & if r ∈ C(Δ), \\ ∞ & otherwise. \end{matrix}

In particular w(r) is a non-negative integer for any r ∈ C(Δ) ∩ ℤ^{n+1}.

Note 14 Choosing δ = δ_{1} corresponds to working with the weight function of Adolphson and Sperber [2], and taking δ = δ_{2} to Dwork’s original weight function [4]. Note that Dwork’s weight function can equivalently be defined as w(r) = r_{0} if r_{1} + ⋯ + r_{n} ≤ r_{0}d and ∞ otherwise, where d is the total degree of f.

The weight function has a simple geometric interpretation: For a real number c define cΔ = {cx | x ∈ Δ}. The next lemma is straightforward. Lemma 15. For any point r ∈ C(Δ) we have that w(r) is the smallest non-negative number c such that r ∈ cΔ. If r ∉ C(Δ) then w(r) = ∞.

By (9) and the sentence preceding Definition 13, we have the following result.

Lemma 16. We have the inequality

ord (F_{r}) ≥ w(r)\frac{p − 1}{p^{2}}.

We also note the following simple property whose proof is straightforward.

Lemma 17. Let r, r′ ∈ ℝ^{n+1} and k a non-negative integer. Then w(kr) = kw(r) and w(r + r′) ≤ w(r) + w(r′). In particular when w(r′) ≠ ∞,

w(kr − r′) ≥ kw(r) − w(r′).

We shall work in certain subrings of R[[X]] defined in terms of the weight function.

Definition 18. Define L_{Δ} to be the subring of R[[X]] given by

L_{Δ} = { ∑_{r∈C(Δ)∩ℤ^{n+1}} A_{r}X^{r} | A_{r} ∈ R}.

Thus, L_{Δ} is just the ring of all power series over R whose terms have exponents in the cone C(Δ). Certainly F ∈ L_{Δ} and from Lemma 10 we see easily that F^{(a)} ∈ L_{Δ}.

Lemma 19. The power series F and F^{(a)} both belong to L_{Δ}.

This concludes all results in this section which shall be essential to the proof of our modular version of Dwork’s Trace Formula. We conclude with a definition and some comments which we will refer to in the analysis of the running time of our algorithm.

Definition 20. For any positive real number b, we define a set of power series by

L_{Δ}(b) = {∑_{r} A_{r}X^{r} ∈ L_{Δ} | ord A_{r} ≥ bw(r)}.

The set L_{Δ}(b) is easily seen to be a subring of L_{Δ}. Elements in L_{Δ}(b) for large b can be thought of as having fast decaying coefficients. Such rings will reduce to rings of small dimension modulo small powers of p.

We have that

(10)  F ∈ L_{Δ} (\frac{p − 1}{p^{2}})
(11)  F^{(a)} ∈ L_{Δ} (\frac{p − 1}{qp}) .

The first is immediate from the inequality in Lemma 16 and the second follows since τ^{i}(F (X^{p^{i}})) ∈ L_{Δ}(\frac{p−1}{p^{i}p^{2}} ) for each 0 ≤ i ≤ a − 1. Thus for a > 1 the coefficients of F^{(a)} decay more slowly than those of F itself.

5. Dwork’s trace formula

5.1. Lifting Frobenius. We now introduce Dwork’s “left inverse of Frobenius” mapping ψ_{p} on the ring R[[X]]. Definition 21. Let ψ_{p} be defined on the monomials in R[[X]] by

ψ_{p}(X^{r}) = \{ \begin{matrix}X^{r/p} & if p|r \\ 0 & otherwise \end{matrix}

and extend ψ_{p} by τ^{−1}-linearity to all of R[[X]]. That is

ψ_{p}(∑_{r} A_{r}X^{r}) = ∑_{r} τ^{−1}(A_{r})ψ_{p}(X^{r}) = ∑_{r} τ^{−1}(A_{pr})X^{r}.

Here p|r means that p divides all of the entries in the integer vector r. This map is a left inverse of the “Frobenius” map on the ring R[[X]] which takes a power series ∑_{r} A_{r}X^{r} to ∑_{r} τ(A_{r})X^{pr}. Definition 22. Let α_{a} be the map from R[[X]] to itself defined as

α_{a} = ψ_{p}^{a} ∘ F^{(a)}.

Precisely, this is the map which is the composition of multiplication by the power series F^{(a)} followed by the mapping ψ_{p}^{a} on the ring R[[X]]. (Notice that ψ_{p}^{a} just acts as

ψ_{p}^{a}(∑_{r} A_{r}X^{r}) = ∑_{r} A_{qr}X^{r}

since τ^{−a}(A_{r}) = A_{r}.) Let the map α from R[[X]] to itself be defined as

α = ψ_{p} ∘ F.

Thus α is multiplication by F followed by the mapping ψ_{p}.

We have the following result relating these two maps, which shall be of crucial importance in our derivation of an efficient algorithm. Lemma 23. With α_{a} and α as in Definition 22 we have

α_{a} = α^{a}

where the second exponent is a power under composition.

Proof. Firstly let H ∈ R[[X]] and denote by ψ_{p} ∘ H(X^{p}) the map composed of multiplication by H(X^{p}) followed by ψ_{p}. We claim that

(12)  ψ_{p} ∘ H(X^{p}) = τ^{−1}(H(X)) ∘ ψ_{p}

To see this write H = ∑_{r} H_{r}X^{pr}. Then we have that

ψ_{p} ∘ H(X^{p}) = ∑_{r} τ^{−1}(H_{r})(ψ_{p} ∘ X^{pr}).

Here the infinite series is interpreted as a mapping. Now ψ_{p} ∘ X^{pr} = X^{r} ∘ ψ_{p} as these two maps are τ^{−1}-linear and agree on monomials. Hence we have ψ_{p} ∘ H(X^{p}) = ∑_{r} τ^{−1}(H_{r})(X^{r} ∘ ψ_{p}) = τ^{−1}(H(X)) ∘ ψ_{p}.

Next we claim that for any b ≥ 1 and power series H we have

ψ_{p}^{b} ∘ ∏_{i=0}^{b−1} τ^{i}(H(X^{p^{i}})) = (ψ_{p} ∘ H(X))^{b}.

We prove this by induction, the result trivially holding if b = 1. For b > 1, by b − 1 applications of (12) we get

ψ_{p}^{b} ∘ ∏_{i=0}^{b−1} τ^{i}(H(X^{p^{i}})) = (ψ_{p} ∘ H(X)) ∘ (ψ_{p}^{b−1} ∘ ∏_{i=0}^{b−2} τ^{i}(H(X^{p^{i}})).

The second claim then follows by induction. Putting H = F and b = a we get the required result. □

The map α_{a} is linear and continuous, in the sense that

α_{a}(∑_{r} A_{r}X^{r}) = ∑_{r} A_{r}α_{a}(X^{r})

for any element ∑_{r} A_{r}X^{r} ∈ R[[X]]. The map α is τ^{−1}-linear and continuous, in the sense that

α(∑_{r} A_{r}X^{r}) = ∑_{r} τ^{−1}(A_{r})α(X^{r}).

From Lemma 19 the next lemma follows easily. Lemma 24. The subring L_{Δ} is stable under both α and α_{a}.

Both maps when restricted to the subring L_{Δ} are determined by their action on the monomials X^{r} (with w(r) < ∞), which we now consider.

5.2. Matrix representations of mappings. Recall that C(Δ) is the cone in ℝ^{n+1} from Definition 12. The set

Γ_{Δ} = {X^{u} | u ∈ C(Δ) ∩ ℤ_{≥0}^{n+1}}

written as a row vector, is a formal basis for the space L_{Δ}. Precisely, this means that any power series in L_{Δ} may be written in exactly one way as an infinite sum ∑_{u} A_{u}X^{u} with X^{u} in the above set. Notice that this is different from the usual notion of a basis in linear algebra, since we allow infinite combinations of basis elements. It is also different from the notion of an “orthonormal basis” in the literature, where one requires that the coefficient A_{u} goes to zero as w(u) goes to ∞ (see [21] for a more detailed discussion of these notions).

By Lemma 24, both α and α_{a} send the ring L_{Δ} to itself. We define certain matrices associated to the maps α and α_{a} restricted to L_{Δ} with regard to the formal row basis Γ_{Δ} of monomials. Definition 25. Let the infinite matrices M and M_{a} have columns describing the images of the monomials X^{v} ∈ L_{Δ} under the maps α and α_{a} with respect to our formal row basis Γ_{Δ}:

α(Γ_{Δ}) = Γ_{Δ}M, α_{a}(Γ_{Δ}) = Γ_{Δ}M_{a}.

Specifically, the (u, v)th entries of M and M_{a} for u, v ∈ C(Δ) are m_{uv} and m_{uv}^{(a)} respectively where

\begin{aligned}m_{uv} & = & τ^{−1}(F_{pu−v}) \\ m_{uv}^{(a)} & = & F_{qu−v}^{(a)}. \end{aligned}

Here F^{(a)} = ∑_{r} F_{r}^{(a)}X^{r} and as before F = ∑_{r} F_{r}X^{r}, and we take the coefficients of exponents r with negative entries to be zero.

(Note that as yet we have not ordered the basis; however, we shall choose a convenient ordering in the proof of Theorem 28.) We have by Lemmas 16 and 17

ord (F_{pu−v}) ≥ \frac{p − 1}{p^{2}}w(pu − v)
(13)  ≥ \frac{p − 1}{p} (w(u) − \frac{1}{p}w(v))

and certainly ord (m_{uv}) = ord (τ^{−1}(F_{pu−v})) = ord (F_{pu−v}).

Notice that the matrix powers M_{a}^{k} and M^{k} are defined for every positive integer k since the entries in M_{a}^{k}, say, are just finite sums of the entries in M_{a} and M_{a}^{k−1}. This follows for M_{a}, say, since all entries m_{uv}^{(a)} in M_{a} are zero when the vector qu − v contain negative entries. We define the trace of an infinite matrix to be the sum of its diagonal entries, when this sum converges, and ∞ when the sum does not converge. We shall see shortly that the trace of the infinite matrix M_{a}^{k} is finite. We write this as Tr(M_{a}^{k}).

We can now prove the following result of Dwork.

Theorem 26 (Dwork’s Trace Formula). Let f ∈ 𝔽_{q}[X_{1}, …, X_{n}] and let α_{a} be the mapping on the ring R[[X]] given as α_{a} = ψ_{p}^{a} ∘ F^{(a)}, where F^{(a)} and ψ_{p} are described in Definitions 9 and 21. Let M_{a} denote the infinite matrix representing the map α_{a} restricted to the subring L_{Δ} as described in Definition 25. For k ≥ 1, denote by N_{k}^{∗} the number of solutions to the equation f = 0 in the torus (𝔽_{q^{k}}^{∗})^{n}. Then

(q^{k} − 1)^{n+1}Tr(M_{a}^{k}) = q^{k}N_{k}^{∗} − (q^{k} − 1)^{n}.

Proof. (The following is a hybrid of the matrix proof given in [19] and Dwork’s original argument [4].)

By Proposition 11 we have that

q^{k}N_{k}^{∗} − (q^{k} − 1)^{n} = ∑_{x^{q^{k}−1}=1} F^{(a)}(x)F^{(a)}(x^{q}) … F^{(a)}(x^{q^{k−1}})

where the sum is over all (n + 1)-tuples of (q^{k} − 1)st roots of unity in Ω (namely the Teichmüller lifting in Ω^{n+1} of points on the torus (𝔽_{q^{k}}^{∗})^{n+1}).

We first consider the case k = 1. Since F^{(a)} ∈ L_{Δ}, we can write F^{(a)}(X) = ∑_{r} F_{r}^{(a)}X^{r} where the sum is over all lattice vectors r which belong to C(Δ).

Then the latter sum is

∑_{x^{q−1}=1} F^{(a)}(x) = ∑_{x^{q−1}=1}(∑_{r} F_{r}^{(a)}x^{r})
= ∑_{r} F_{r}^{(a)}(∑_{x^{q−1}=1} x^{r})
= (q − 1)^{n+1} ∑_{r,(q−1)|r} F_{r}^{(a)}
= (q − 1)^{n+1} ∑_{s} F_{(q−1)s}^{(a)}
= (q − 1)^{n+1}Tr(M_{a}).

Here by (q − 1)|r we mean that the integer q − 1 divides every entry in the vector r. Also, we use the fact that for any integer r_{i} [11, page 120]

∑_{x_{i}^{q−1}=1} x_{i}^{r_{i}} = \{ \begin{matrix}q − 1 & if (q − 1)|r_{i} \\ 0 & otherwise. \end{matrix}

For k > 1, define M_{ak} to be the matrix for the map ψ_{p}^{ak} ∘ ∏_{i=0}^{k−1} F^{(a)}(X^{q^{i}}) with respect to the formal basis Γ_{Δ}. Then by an analogous argument to that in the case k = 1 we see

∑_{x^{q^{k}−1}=1} ∏_{i=0}^{k−1} F^{(a)}(x^{q^{i}}) = (q^{k} − 1)^{n+1}Tr(M_{ak}).

Using (12) in the proof of Lemma 23 one sees that

ψ_{p}^{ak} ∘ ∏_{i=0}^{k−1} F^{(a)}(X^{q^{k}}) = (ψ_{p}^{a} ∘ F^{(a)}(X))^{k}.

Since both maps are linear it follows that M_{ak} = M_{a}^{k}, and the theorem is proved. □

To compute the trace of the matrix in Dwork’s formula we shall use the following matrix identity derived from Lemma 23. Lemma 27. Let M_{a} and M denote the matrices for the maps α_{a} and α described above. Then

M_{a} = ∏_{i=0}^{a−1} τ^{−i}(M) = MM^{τ^{−1}} ⋯ M^{τ^{−(a−1)}}

where the map τ^{−i} acts entry wise on the matrix M.

Proof. First suppose that N_{1} and N_{2} are matrices representing maps β_{1} and β_{2} on some subspace of R[[X]]. We assume that β_{1} is τ^{−1}-linear and β_{2} is τ^{−j}-linear for some j ∈ ℤ. Then it is not difficult to prove that the matrix for the τ^{−(j+1)}-linear map β_{1} ∘ β_{2} is just N_{1}τ^{−1}(N_{2}).

By Lemma 23 we have α_{a} = α^{a}. We claim that for any positive integer b the matrix for the map α^{b} is ∏_{i=0}^{b−1} τ^{−i}(M). The result is trivially true if b = 1. For b > 1 we have α^{b} = α ∘ α^{b−1}. By induction the matrix for α^{b−1} is ∏_{i=0}^{b−2} τ^{−i}(M). By τ^{−1}-linearity of α it follows from the observations in the preceding paragraph that the matrix for α^{b} is Mτ^{−1}(∏_{i=0}^{b−2} τ^{−i}(M)) = ∏_{i=0}^{b−1} τ^{−i}(M). The required result now follows by taking b = a. □

We note in passing that M_{a} in Lemma 27 also equals

τ^{a−1}(S)τ^{a−2}(S) … S,

where S = τ(M) is the matrix with (u, v)th entry simply F_{pu−v}. Although this is slightly more desirable from a practical point of view we shall not use this expression for M_{a}.

5.3. Modular reduction of the trace formula. We now examine the reduction of the Dwork trace formula modulo a power p^{N} of p. Observe that both sides in the trace formula, and all the entries in the matrices M and M_{a}, are elements of R, and so this reduction is defined.

We first recall some notation: Let C(Δ) be the cone in ℝ^{n+1} from Definition 12, and L_{Δ} denote the ring of power series over R whose monomials have exponents lying in C(Δ) (Definition 18). Let M denote the matrix for the τ^{−1}-linear map α = ψ_{p} ∘ F with respect to the formal row basis Γ_{Δ} (Definition 25). Here ψ_{p} is the “left inverse of Frobenius” given in Definition 21 and F is the power series obtained from the polynomial f as in Definition 9. Theorem 28. Let N denote any positive integer and A_{N} the finite square matrix over the finite ring R/(p^{N}) obtained by reducing modulo p^{N} all those entries in M whose rows and columns are indexed by vectors u ∈ C(Δ) ∩ ℤ_{≥0}^{n+1} with w(u) < (p/(p − 1))^{2}N. Then

(q^{k} − 1)^{n+1}Tr((∏_{i=0}^{a−1} τ^{−i}(A_{N}))^{k}) = q^{k}N_{k}^{∗} − (q^{k} − 1)^{n} mod p^{N}

where N_{k}^{∗} is the number of solutions to the equation f = 0 in the affine torus (𝔽_{q^{k}}^{∗})^{n}. Moreover, the size of A_{N} is W = #(tΔ), where t = ⌈p^{2}N/(p−1)^{2}⌉− 1 and #(tΔ) is number of lattice points in a dilation by a factor t of the polytope Δ.

Proof. For any finite or infinite matrix L with coefficients in R, we define \overline{L} to be the matrix obtained by reducing all its entries modulo p^{N}. Thus \overline{L} has entries in R/(p^{N}).

The theorem will follow from the Dwork Trace Formula (Theorem 26) once we find a suitable expression for the reduction modulo p^{N} of the trace of the matrix M_{a}^{k}. This is equal to the trace of the matrix \overline{M_{a}^{k}} = (\overline{M_{a}})^{k}. By Lemma 27 this matrix can be computed as a matrix product from the matrix \overline{M}.

By inequality (13) every entry m_{uv} = τ^{−1}(F_{pu−v}) in M satisfies

ord (m_{uv}) ≥ \frac{p − 1}{p} (w(u) − \frac{w(v)}{p}) .

Define t to be the greatest integer less than (p/(p − 1))^{2}N. When w(u) ≥ w(v) and w(u) > t we have

(14)  \frac{p − 1}{p} (w(u) − \frac{w(v)}{p}) ≥ \frac{(p − 1)^{2}}{p^{2}} w(u) ≥ \frac{(p − 1)^{2}}{p^{2}} (t + 1) ≥ N.

Recall that we have not yet ordered the basis. Now choose any total ordering on the basis set such that for distinct lattice points u, v ∈ C(Δ), the monomial X^{v} comes before X^{u} if w(v) < w(u). Thus, by the inequalities in (14) and the choice of ordering of the basis, the matrix \overline{M} is of the form

(15)  ( \begin{matrix}A_{N} & B_{N} \\ 0 & C_{N} \end{matrix} )

where C_{N} is a strictly upper triangular infinite matrix and A_{N} is the finite square matrix indexed by lattice points u ∈ C(Δ) such that w(u) ≤ t. The size of A_{N} is the number of lattice points u with w(u) ≤ t. By Lemma 15 this is exactly the number of lattice points in the polytope tΔ.

By Lemma 27 and modular reduction, one has

(\overline{M_{a}})^{k} = (∏_{i=0}^{a−1} τ^{−i}(\overline{M}))^{k}.

By (15) we see that (\overline{M_{a}})^{k} is of the form

( \begin{matrix}(∏_{i=0}^{a−1} τ^{−i}(A_{N}))^{k} & B′_{N} \\ 0 & C′_{N} \end{matrix} )

where C′_{N} = (∏_{i=0}^{a−1} τ^{−i}(C_{N}))^{k} is strictly upper triangular.

Hence the trace of (\overline{M_{a}})^{k} equals the trace of the finite matrix

(16)  (∏_{i=0}^{a−1} τ^{−i}(A_{N}))^{k}.

The theorem now follows from (16) and Theorem 26. □

6. Algorithms

In this section we present an algorithm for counting points based upon Theorem 28, and complete the proofs of the results in the introduction. This is a relatively straightforward matter, although the precise complexity estimates require a little care.

6.1. Toric point counting algorithm. We first give the algorithm.

Algorithm 29 Toric point counting Input: Positive integers a, k, n, d and a prime p; a polynomial f; a polytope Δ. (Here f is a polynomial in n variables of total degree d with coefficients in the field 𝔽_{q} where q = p^{a}. The polytope Δ is as in Definition 12. We assume a model of 𝔽_{q} is given as in Section 3.2.1.) Output: The number of solutions N_{k}^{∗} to the equation f = 0 in the torus (𝔽_{q^{k}}^{∗})^{n}.

Step 0: Set N = (n + 1)ak, where q = p^{a}.

Step 1: Compute the polynomial F mod p^{N} in the ring (R/(p^{N}))[X] where F is the power series in Definition 9. Specifically, writing X_{0}f = ∑_{j∈J} ā_{j}X^{j} we have F = ∏_{j∈J} θ(a_{j}X^{j}) mod p^{N}. Here a_{j} is the Teichmüller lifting of ā_{j} and θ(z) = exp(π(z − z^{p})). (See Section 3 for a description of the ring R/(p^{N}) and the element π.)

Step 2: Construct the matrix A_{N} which occurs in the statement of Theorem 28. Specifically, the matrix A_{N} is indexed by pairs (u, v) where u and v are lattice points in the dilation by a factor t of the polytope Δ, and t = ⌈p^{2}N/(p − 1)^{2}⌉ − 1. The (u, v)th entry of A_{N} is τ^{−1} of the coefficient of X^{pu−v} in the polynomial F mod p^{N}. The action of τ is as described in Section 3.

Step 3: Compute the product (∏_{i=0}^{a−1} τ^{−i}(A_{N}))^{k}. Let T denote the trace of this product.

Step 4: Output

N_{k}^{∗} = q^{−k}[((q^{k} − 1)^{n+1}T + (q^{k} − 1)^{n})mod p^{N}],

where the square brackets denote the smallest non-negative residue modulo p^{N}.

In the algorithm we assume that the polynomial is presented as input explicitly via its list of non-zero terms. The manner of presentation of Δ only affects the time required to find all lattice points in the dilated polytope tΔ. For concreteness, let us say it is presented via its list of vertices, although any other reasonable presentation would suffice.

6.2. Proof of correctness of the algorithm. We know by Theorem 28 that

q^{k}N_{k}^{∗} − (q^{k} − 1)^{n} = (q^{k} − 1)^{n+1}Tr((∏_{i=0}^{a−1} τ^{−i}(A_{N}))^{k}) mod p^{N}

in the ring R/(p^{N}). Thus,

q^{k}N_{k}^{∗} = (q^{k} − 1)^{n} + (q^{k} − 1)^{n+1}T mod p^{N}.

The lefthand side is a non-negative integer. The first term on the righthand side is an integer. The second term on the righthand side is the reduction modulo p^{N} of the trace of M_{a}^{k}, which is known to be an integer by Dwork’s Trace Formula (Theorem 26), and thus it is also an integer. Since N_{k}^{∗} ≤ (q^{k} − 1)^{n} it follows that the lefthand side is smaller than q^{k(n+1)}. Hence in the case N ≥ ak(n + 1) we must have

q^{k}N_{k}^{∗} = [((q^{k} − 1)^{n+1}T + (q^{k} − 1)^{n})mod p^{N}].

The proof is complete.

6.3. Complexity analysis. We shall use Big-Oh and Soft-Oh notation in our analysis of the complexity of the above algorithm. If C_{1} and C_{2} are real functions we write C_{1} = O(C_{2}) if |C_{1}| ≤ c(|C_{2}| + 1) for some positive constant c. We write C_{1} = Õ(C_{2}) if C_{1} = O(C_{2} log(|C_{2}| + 1)^{c′}) for some constant c′. Thus in the latter notation one ignores logarithmic factors.

6.3.1. Ring operations. We shall first of all count the number of operations in the ring R/(p^{N}) required in Steps 1, 2, 3, ignoring for the time being any other auxiliary computations. More precisely, because of the complexity bounds in Lemma 5 it is convenient for our analysis to define a “ring operation” to be either arithmetic or the evaluation of the map τ^{i}, for 1 ≤ i ≤ a−1, in the ring R/(p^{N}) (excluding precomputation). In Section 6.3.2 we shall add back in the small contribution from computing Teichmüller liftings and also the precomputation required for the maps τ^{i}. Similarly, in Section 6.3.3 we shall account for the remaining operations in Steps 1, 2, 3, arising mainly from computations with the exponents of polynomials (at this stage we will restrict the input polytope Δ to avoid complications from convex geometry). The contributions from Steps 0 and 4 are easily seen to be absorbed into the other estimates, and we shall not mention them again.

Our running time will be in terms of the parameters t, t̃, W, W̃. Here

\begin{matrix}t = ⌈(\frac{p}{p−1})^{2} N⌉ − 1, & W = #(tΔ), \\ t̃= ⌈\frac{p^{2}}{p−1}N⌉ − 1, & W̃ = #(t̃Δ), \end{matrix}

where N = ak(n+1) and the # operator counts lattice points in convex sets. The sizes of the sets W and W̃ are the number of lattice points in certain “truncated” cones. Since t̃ is about p times as large as t, the integer W̃ will be around p^{n+1} times as large as W, since we are working in n+1 dimensional space. The integer W is precisely the size of the matrix which occurs in Step 2. The integer W̃ will turn out to be the maximum number of terms in the polynomial we compute in Step 1. Define L_{Δ}((p − 1)/p^{2}) mod p^{N} to be the ring of polynomials obtained by reducing the coefficients of power series in L_{Δ}((p − 1)/p^{2}) ⊆ R[[X]] modulo p^{N}. Then W̃ is the number of monomials which occur in the finite ring L_{Δ}((p − 1)/p^{2}) mod p^{N}.

For Step 1 we have the following estimate. Lemma 30. Let F be the power series given in Definition 9 and N a positive integer. The polynomial F mod p^{N} may be computed in

O(|J|W̃^{2})

operations in the ring R/(p^{N}).

Proof. By (4) we have that θ(z) mod p^{N} is a polynomial of degree not greater than p^{2}N/(p − 1). Thus we can obtain θ(z) via the formula

θ(z) = exp(πz) exp(−πz^{p})

by computing the first O(pN) terms in the expansion for exp(πz), substituting z = −z^{p}, and one multiplication of polynomials. Note that exp(πz) has p-adic integral coefficients. Thus θ(z) can be found in time O((pN)^{2}) operations in the ring R/(p^{N}) using standard polynomial arithmetic.

Now by (10) we have that F ∈ L_{Δ}((p−1)/p^{2}). Thus F mod p^{N} ∈ L_{Δ}((p− 1)/p^{2}) mod p^{N}. One may then compute F mod p^{N} directly from Definition 9 in |J| − 1 multiplications of polynomials of the form θ(a_{j}X^{j}) mod p^{N}. Each such polynomial lies in the ring L_{Δ}((p − 1)/p^{2}) mod p^{N}, because ord (a_{j}) = 0, w(j) = 1 and the coefficients of θ decay at a suitable rate. Hence all computations required in computing F mod p^{N} involve polynomials in this ring. Such polynomials have at most W̃ terms. Exactly |J| − 1 multiplications are required. Thus the complexity is O(|J|W̃^{2}) ring operations. Noting that pN = O(W̃) we have the result. □

Note 31 It is crucial here that we only need to compute F mod p^{N} and not F^{(a)} mod p^{N}, as one might attempt to do using a more naive approach. The latter polynomial has very high degree (O(q)) because of the slow decay rate of the coefficients of F^{(a)}.

With regard to Step 2, given that the polynomial F mod p^{N} has already been computed the only task required is to identify those pairs of points (u, v) such that u, v ∈ tΔ, compute the integer point pu − v, and copy τ^{−1} of the term F_{pu−v} mod p^{N} from F mod p^{N} into the correct position in the matrix. Thus no arithmetic operations in the ring are required here, except W^{2} computations of τ^{−1} = τ^{a−1}. These arithmetic operations can safely be ignored since in Lemma 30 we have already counted O(|J|W̃^{2}) ring operations. Computation of the appropriate indices (u, v) does require one to find all lattice points in certain polytopes and we return to that in Section 6.3.3.

Finally, for Step 3 we have the following estimate. Lemma 32. With the notation as in the statement of Theorem 28, the product

(∏_{i=0}^{a−1} τ^{−i}(A_{N}))^{k}

can be computed given the matrix A_{N} in

O(W^{3} log (ak))

operations in the ring R/(p^{N}).

Proof. A fast square-and-multiply style algorithm may be used to compute the power α^{a} in O(log a) “matrix ring operations”. Specifically, working with the matrix representations one may compute α^{r+s} from α^{r} and α^{s} using

∏_{i=0}^{r+s−1} τ^{−i}(A_{N}) = ∏_{i=0}^{r−1} τ^{−i}(A_{N}) (τ^{−r}(∏_{i=0}^{s−1} τ^{−i}(A_{N}))) .

Now the case r = s = 2^{c} for some c gives us the “square” step (computing α^{2^{c+1}} from α^{2^{c}}) and the case r = 2^{c} and s arbitrary the “multiply” step. These two operations may be combined to give a fast exponentiation method in a straightforward way. The time required to compute a matrix for α^{a} from one for α is thus O(log a) “matrix ring operations”. By matrix ring operations we mean multiplication of matrices of size W over R/(p^{N}), and also computing τ^{−i}(B) for some 1 ≤ i ≤ a − 1 and matrix B of this form. The former requires O(W^{3}) operations in the ring R/(p^{N}) using standard algorithms. Since τ^{−i} = τ^{a−i} the latter may be computed in O(W^{2}) ring operations (that is, applications of a power of τ). Thus the time for computing a matrix for α^{a} from one for α is O(W^{3} log a).

Having obtained a matrix for α^{a} one may then compute a matrix for α^{ak} using the standard square-and-multiply algorithm. This requires O(log k) matrix multiplications, that is O(W^{3} log k) operations in R/(p^{N}). Thus the total time required is as claimed. □

Note that the exponent 3 for multiplication of matrices can be improved to around 2.4 using faster methods [7, page 330].

Gathering these results we find the following.

Lemma 33. The running time of Algorithm 29 is

O(W̃^{2}|J| + W^{3} log (ak))

ring operations. Here |J| is the number of non-zero terms in f, and W̃ and W are defined as at the start of Section 6.3.1, with N = ak(n + 1).

6.3.2. Bit complexity arising from ring operations. Using Lemma 5 one may now calculate the number of bit operations in the algorithm which arise from operations in the ring R/(p^{N}). In this section we shall also count the small contribution from computing the Teichmüller lifting of the coefficients of f, and also the precomputation required for powers of τ, which it was convenient to ignore in Section 6.3.1,

Definition 34. Let the polytope Δ from Definition 12 have dimension ñ ≤ n + 1. Let V (Δ) be the ñ-dimensional volume of Δ. Denote by v = ñ!V (Δ) the “normalised” volume of Δ.

Since f is non-zero we have ñ ≥ 1 and certainly v > 0.

Proposition 35. The running time of Algorithm 29 is

Õ(a^{3n+7}k^{3n+5}n^{3n+5}v^{3}p^{2n+4})

bit operations, plus the contribution from operations outside of the ring R/(p^{N}). The space complexity in bits is

Õ(a^{2n+4}k^{2n+3}n^{2n+3}v^{2}p)

plus the contribution from operations outside of R/(p^{N}).

Proof. To compute the complexity first observe that t = ⌈p^{2}N/(p−1)^{2}⌉−1 ≤ 4N. Thus t, N = O(akn) in the algorithm. Also t̃= O(pN) = O(pt). Now

W = #(tΔ)
≤ ñ!V (tΔ) + ñ
= ñ!tñV (Δ) + ñ
= O(vt^{n+1}).

Here we use the Blichfeldt bound #(P) ≤ m!V (P)+m for any m-dimensional polytope P from [9, page 144], the fact V (tΔ) = tñV (Δ) since Δ is ñdimensional, and also that ñ ≤ n + 1. Similarly W̃ = #(t̃Δ) = O(vt̃^{n+1}) = O(vt^{n+1}p^{n+1}).

Thus from Lemma 33 the number of ring operations is

O((vt^{n+1}p^{n+1})^{2}(v + n) + (t^{n+1}v)^{3} log (ak))

since |J| ≤ #(Δ) − 1 ≤ v + n. Thus the bit complexity which arises from ring operations is by Lemma 5

O(((vt^{n+1}p^{n+1})^{2}(v + n) + (t^{n+1}v)^{3} log (ak))(paN log p)^{2}).

Tidying up and ignoring logarithmic factors we get

Õ(v^{3}t^{2n+2}p^{2n+4}a^{2}N^{2} + v^{3}t^{3n+3}p^{2}a^{2}N^{2}).

The second term is dominant in all factors except p. For simplicity we take the estimate of

Õ(v^{3}t^{3n+3}p^{2n+4}a^{2}N^{2}).

Now putting t, N = O(ank) we get

Õ(v^{3}a^{3n+7}k^{3n+5}n^{3n+5}p^{2n+4}).

This is the total bit complexity which arises from “ring operations”, as defined at the start of Section 6.3.1. There remains the contribution from computing the Teichmüller liftings of the coefficients of f in Step 1, and also the precomputation required for the map τ. By Lemma 5, this is easily seen to be absorbed in the above estimate.

With regard to the space complexity, this is dominated by the space required to store the matrix A_{N}, which is O(W^{2}) ring elements. Putting W = O(v(ank)^{n+1}) and using Lemma 5 gives us the result. □

One may replace n in the exponents in Proposition 35 by ñ − 1; however, this only gives an improvement when the Newton polytope is not fulldimensional.

6.3.3. Bit complexity arising from auxiliary operations. It remains to bound the complexity which arises from operations outside of the ring R/(p^{N}) in Steps 1 and 2. In Step 1 manipulation of exponents of polynomials will add an extra term O(log(pNd)) to the running time, which can safely be ignored.

In Step 2 one is required to find all lattice points which lie in tΔ, for t = O(ank). The complexity of this step will depend upon the input polytope Δ. For a “general” Δ one requires methods from computational convex geometry which are not in the spirit of the present exposition. Thus our total bit complexity estimate for Algorithm 29 will just be as in Proposition 35 “plus the contribution from finding all lattice points in tΔ”.

At this stage for simplicity we shall restrict to the choice of δ = δ_{2} in Definition 12. Thus we take Δ = Δ_{2} as the convex hull in ℝ^{n+1} of the origin and the n + 1 points

(1, 0, …, 0), (1, d, 0, …, 0), …, (1, 0, …, 0, d).

For this case the required set of lattice points is

{(r_{0}, r_{1}, …, r_{n}) | r_{1} + ⋯ + r_{n} ≤ dr_{0} ≤ dt}

and so no computations are required here. Also, now v equals d^{n} and directly from Proposition 35 we get the following result. Proposition 36. Let Algorithm 29′ be exactly as Algorithm 29 only with input restricted to the choice of polytope Δ = Δ_{2} described in the preceding paragraph. The total running time of Algorithm 29′ is

Õ(a^{3n+7}k^{3n+5}n^{3n+5}d^{3n}p^{2n+4})

bit operations. The total space complexity in bits is

Õ(a^{2n+4}k^{2n+3}n^{2n+3}d^{2n}p).

We shall use this restricted version of Algorithm 29 in the proofs of our main results in the next section.

6.4. Proofs of the results in the Introduction. To compute the number of points on the affine variety defined by a polynomial f one simply uses the torus decomposition of 𝔽_{q^{k}}^{n}. Specifically, for any subset S ⊆ {1, 2, …, n} let G_{k}^{S} denote the set of points

{(x_{1}, …, x_{n}) | x_{i} ∈ 𝔽_{q^{k}}, x_{i} = 0 ⇔ x ∈ S}.

Denote by f^{S} the polynomial obtained from f by setting to zero all indeterminates X_{i} which occur in f for i ∈ S. Denote by N_{k}^{S} the number of solutions of f^{S} = 0 in the torus G_{k}^{S} of dimension n−|S|. Then N_{k} = ∑_{S} N_{k}^{S} where the sum is over all subsets of {1, 2, …, n}. Each number N_{k}^{S} can be computed using Algorithm 29′. (If some f^{S} is identically zero or has degree 0 then N_{k}^{S} = (q^{k} − 1)^{n−|S|} or 0, respectively, and Algorithm 29′ is not required!) Thus by 2^{n} applications of this algorithm we obtain N_{k} as desired.

Now to obtain the whole zeta function Z(f/𝔽_{q}) it suffices to count N_{k} for all k = 1, …, deg(r) + deg(s), where

Z(f/𝔽_{q})(T) = \frac{r(T)}{s(T)}

with r and s coprime polynomials in 1 + Tℤ[T]. More precisely, it is enough to know upper bounds deg(r) ≤ D_{1} and deg(s) ≤ D_{2}, and compute N_{k} for k = 1, …, D_{1} + D_{2}. Then use the linear algebra method described prior to [22, Corollary 2.8], which we now supplement with further details.

Let u(T) = 1 + ∑_{i=1}^{D_{1}} u_{i}T^{i} and v(T) = 1 + ∑_{i=1}^{D_{2}} v_{i}T^{i} have indeterminate coefficients. Write Z(f/𝔽_{q})(T) = 1 + z_{1}T + z_{2}T^{2} + . . .. This power series has non-negative integer coefficients and it can easily be computed modulo T^{D_{1}+D_{2}+1} given N_{k} for k = 1, …, D_{1} + D_{2}. The equation v(T)Z(f/𝔽_{q}) ≡ u(T) mod T^{D_{1}+D_{2}+1} defines a linear system Ax = y where A is a known square D_{1} + D_{2} integer matrix and y a known integer column vector. The entries in A and y are just coefficients from the power series Z(f/𝔽_{q}) mod T^{D_{1}+D_{2}+1}. Let b be a bound on their bit length. The unknown entries in x are the coefficients of u and v. By [22] the set of all solutions to this system consists of precisely those vectors x derived by specialising the coefficients of u(T) and v(T) to equal those of d(T)r(T) and d(T)s(T), respectively, for some d(T) ∈ 1 + Tℤ[T] with degree at most min(D_{1} − deg(r), D_{2} − deg(s)). In particular, the system has a unique solution (i.e. det(A) ≠ 0) if and only if either deg(r) = D_{1} or deg(s) = D_{2} (or both). The determinant det(A) can be computed using the small primes method in [7, Algorithm 5.10] in a number of bit operations bounded by Õ((D_{1} + D_{2})^{4}b^{2}) [7, Theorem 5.12]. Let us assume now that det(A) ≠ 0 and so the system has a unique solution, namely the unknown vector containing the integer coefficients of r and s. Let B be a bound on the bit length of these coefficients. Find the unique solution to the linear system modulo enough small primes which do not divide det(A), and recover this integer solution using the Chinese remainder theorem. Precisely, work modulo a collection of such primes whose product has bit length greater than B. This second step requires Õ((D_{1} + D_{2})^{4}B^{2}) bit operations using Gaussian elimination (this can be improved with a Padé Approximation algorithm [7, Section 5.9]). Values for b and B may be deduced from the bound N_{k} ≤ q^{nk}. Specifically, one may show from this that the absolute values of the reciprocal zeros of r and s are all ≤ q^{n}, and so we can take B = O((D_{1} + D_{2})n log q). Also, we can take b = O((D_{1} + D_{2})^{2}n log q). If in the above we find det(A) = 0 then we must have that deg(r) < D_{1} and deg(s) < D_{2}. In this case one must first reduce D_{2}, say, and compute determinants until the correct value D_{2} = deg(s) is found (then det(A) ≠ 0 and the above method works).

By the refinement of Bombieri’s degree bound [3] from [2, Equation (1.13)], the “total degree” deg(r) + deg(s) is bounded by 2^{n+1}6^{n+1}(n + 1)!V (Δ_{1}), where Δ_{1} is the polytope in ℝ^{n+1} derived from the Newton polytope of f (see the paragraph following Definition 12). Certainly (n + 1)!V (Δ_{1}) ≤ d^{n}. Hence we may take D_{1}, D_{2} = 2^{4n+4}d^{n} and so D_{1} + D_{2} = 2^{4n+5}d^{n}. Theorem 37. Let f be a polynomial in n variables of total degree d > 0 over 𝔽_{q}, where q = p^{a}. The full zeta function Z(f/𝔽_{q}) can be computed deterministically in

Õ(2^{13n^{2}}a^{3n+7}d^{3n^{2}+9n}p^{2n+4})

bit operations. (Here we use soft-Oh notation which ignores logarithmic factors, as defined at the start of Section 6.3.)

Proof. From Proposition 36 and the torus decomposition method, the bit complexity of computing N_{k} for k = 1, …, 2^{4n+5}d^{n} is

Õ( ∑_{k=1}^{2^{4n+5}d^{n}}{a^{3n+7}k^{3n+5}n^{3n+5}d^{3n}p^{2n+4}}2^{n}).

(The contribution from recovering Z(f/𝔽_{q}) from the N_{k} is absorbed in this estimate.) Tidying up the factor in n we get the claimed result. □

Since we may assume that d > 1 we have 2^{n^{2}} = O(d^{n^{2}}), and Theorem 1 now follows.

The proof of Corollary 2 was explained in the introduction, and we finish with some comments on Corollary 3. By Weil’s theorem, the zeta function of the smooth projective curve Ṽ from Corollary 3 is of the form

Z(Ṽ)(T) = \frac{P (T)}{(1 − T)(1 − qT)}

for some polynomial P (T) whose reciprocal roots have complex absolute value q^{1/2}. Since the (possibly singular) affine curve V and the smooth projective curve Ṽ differ in only finitely many closed points, we deduce that the zeta function of V is of the form

Z(V)(T) = \frac{P (T)Q(T)}{1 − qT},

where Q(T) is a rational function whose zeros and poles are roots of unity. This zeta function, and in particular the rational function P (T)Q(T), may be computed within the time bound in Corollary 3 by Corollary 2. In terms of the pure weight decomposition [22], the polynomial P (T) (respectively, Q(T)) is exactly the pure weight 1 (respectively, weight 0) part of the product P (T)Q(T), and can be recovered quickly from P (T)Q(T) via the LLL polynomial factorization algorithm. In our current special case, one can proceed directly without using the LLL-factorization algorithm. By repeatedly removing the common factor of the numerator of P (T)Q(T) with T^{s} − 1 for ϕ(s) (Euler totient function) not greater than the total degree of P (T)Q(T), the desired polynomial P (T) can be recovered. The order of the group of rational points on Ṽ is simply P (1), see [22]. Thus for fixed dimension and finite field one can compute the order of the group of rational points on the Jacobian of a smooth projective curve in time polynomial in the degree d.

References

[1] L. Adleman and M.D. Huang, Counting rational points on curves and abelian varieties over finite fields, in “Algorithmic Number Theory II” (H. Cohen), Lecture Notes in Computer Science 1122, Springer-Verlag, 1996, 1-16.

[2] A. Adolphson and S. Sperber, Newton polyhedra and the degree of the L-functions associated to an exponential sum, Invent. Math., 88 (1987), 555-569.

[3] E. Bombieri, On exponential sums in finite fields, II, Invent. Math. 47, (1978), 29-39.

[4] B. Dwork, On the rationality of the zeta function of an algebraic variety, Amer. J. Math., 82, (1960), 631-648.

[5] B. Dwork, On the zeta function of a hypersurface, Pub. IHES No. 12, 1962.

[6] N. Elkies, Elliptic and modular curves over finite fields and related computational issues, in “Computational perspectives in number theory: Proceedings of a conference in honour of A.O.L. Atkin” , (D.A. Buell and J.T. Teitelbaum), American Mathematical Society International Press 7, 1998, 21-76.

[7] J. von zur Gathen and J. Gerhard, Modern Computer Algebra, Cambridge University Press, 1999.

[8] P. Gaudry and N. Gürel, An extension of Kedlaya’s algorithm for counting points on superelliptic curves, in “Advances in Cryptology - ASIACRYPT 2001” (C. Boyd), Lecture Notes in Computer Science 2248, Springer, (2001), 480-494.

[9] J.E. Goodman and J. O’Rourke (Eds), Handbook of Discrete and Computational Geometry, CRC Press, 1997.

[10] K.S. Kedlaya, Counting points on hyperelliptic curves using Monsky-Washnitzer cohomology, J. Ramanujan Math. Soc. 16 (2001), 323-338.

[11] N. Koblitz, p-Adic Numbers, p-Adic Analysis and Zeta Functions, 2nd Edition, Springer, 1984.

[12] A.G.B. Lauder and D. Wan, Computing zeta functions of Artin-Schreier curves over finite fields, London Math. Soc. JCM, Vol. 5, (2002), 34-55.

[13] R. Lidl and H. Niederreiter, Introduction to Finite Fields and Their Applications, Cambridge University Press, 1986.

[14] J. Pila, Frobenius maps of Abelian varieties and finding roots of unity in finite fields, Math. Comp 55 (1990), 745-763.

[15] B. Poonen, Computational aspects of curves of genus at least 2, in “Algorithmic Number Theory II” (H. Cohen), Lecture Notes in Computer Science 1122, Springer, 1996, 283-306.

[16] T. Satoh, The canonical lift of an ordinary elliptic curve over a finite field and its point counting, J. Ramanujan Math. Soc. 15, (2000), 247-270.

[17] R. Schoof, Elliptic curves over finite fields and the computation of square roots mod p, Math. Comp. 44 (1985), no. 170, 483-494.

[18] R. Schoof, Counting points on elliptic curves over finite fields, J. Théorie des Nombres de Bordeaux 7, 1998, 219-254.

[19] D. Wan, Meromorphic continuation of L-functions of p-adic representations, Ann. Math., 143 (1996), 469-498.

[20] D. Wan, Computing zeta functions over finite fields, Contemporary Mathematics, 225 (1999), 131-141.

[21] D. Wan, Rank one case of Dwork’s conjecture, J. Amer. Math. Soc., 13(2000), 853-908.

[22] D. Wan, Algorithmic theory of zeta functions over finite fields, preprint.

Computing Laboratory, Oxford University, Oxford OX1 3QD, UK E-mail

address: alan.lauder@comlab.ox.ac.uk

Department of Mathematics, University of California, Irvine, CA 92697,

USA E-mail address: dwan@math.uci.edu