Abdoulaye Maiga

Damien Robert

Towards computing canonical lifts of ordinary elliptic curves in medium characteristic

Keywords: Canonical lift of Elliptic curves, Isogeny computation, Point counting. 2010 Mathematics Subject Classification. Primary canonical lift, point counting. We

 developed several algorithms to compute the canonical lift of an ordinary elliptic curve E over F p n with j-invariant not in F p 2 . When p is constant, the best variant has complexity Õ(nm) Bit operations to lift E to p-adic precision m. As an application, lifting E to precision m = O(n) allows to recover its cardinality in time Õ(n 2 ). However, taking p into account the complexity is Õ(p 2 nm), so Satoh's algorithm can only be applied to small p.

We propose in this paper two variants of these algorithms, which do not rely on the modular polynomial, for computing the canonical lift of an ordinary curve. Our new method yield a complexity of Õ(pnm) to lift at precision m, and even Õ( √ pnm) when we are provided a rational point of p-torsion on the curve. This allows to extend Saoth's point counting algorithm to larger p.

Introduction

Let E/F q be an elliptic curve over a finite field. Schoof's method [START_REF] Schoof | Elliptic curves over finite fields and the computation of square roots mod p[END_REF] gives a polynomial time algorithm to count the number of point of E. The complexity was later improved by Atkin and Elkies to give the SEA algorithm [START_REF] Elkies | Explicit isogenies[END_REF][START_REF] Elkies | Elliptic and modular curves over finite fields and related computational issues[END_REF][START_REF] Bostan | Fast algorithms for computing isogenies between elliptic curves[END_REF]. The algorithm can be seen as an incarnation of ℓ-adic étale cohomology: if χ(t) is the characteristic polynomial of the Frobenius π q , χ(t) mod ℓ is computed modulo several primes ℓ by looking at the action of π q on (a subgroup of) the ℓ-torsion E [ℓ]. The CRT algorithm allows to reconstruct χ(t) once we have enough precision (as bounded by the Hasse-Weil bound). One can compute χ mod ℓ in Õ ((ℓ + log q)ℓ log q), hence reconstruct χ in Õ log 4 q) .

In 2000, a second class of algorithms was introduced by Satoh [Sat00], using the Lubin-Serre-Tate Theorem. Let q = p n , let Z q denotes the ring of Witt vectors of F q , and Q q = Frac(Z q ) the unique unramified extension of Q p of degree n. Then [START_REF] Lubin | Elliptic curves and formal groups[END_REF] establishes the existence of a unique (up to isomorphisms) elliptic curve E ↑ over Z q for every ordinary elliptic curves E/F q such that the modulo p reduction of E ↑ is E and End(E ↑ ) ∼ = End(E) as a ring. The curve E ↑ is called the canonical lift of E. Then the trace of the Frobenius morphism is deduced using crystalline cohomology. After improvements by Harley, Satoh's algorithm can compute the canonical lift to precision m in quasi-linear time Õp (nm). Here the notation Õp means that we assume that p is a constant. We can then recover the inversible eigenvalue of the Frobenius at precision m in the same time. By Hasse's bound, it suffices to work at precision m = O(n) to recover the full eigenvalue, so Satoh's algorithm gives a point counting algorithm of quasi quadratic complexity Õp (n 2 ).

We are interested in the dependency of p of the algorithm. We will now assume that p > 2 for simplicity For an ordinary elliptic curve E/F q , Satoh's algorithm and its improvements [Ver03, Gau04] proceeds in four steps:

(1) Compute the canonical lift E ↑ /Z q at p-adic precision m by solving the equation Φ p (j(E ↑ ), Σ(j(E ↑ ))) = 0 via a Newton lift. Here j is the j-invariant, Φ p the modular polynomial classifying p-isogenies, and Σ the (small) Frobenius on Z q . (2) Lift the kernel E [p] et of the Verschiebung to E ↑ via a Newton lifts. The kernel of the Verschiebung modulo p is defined by the x-coordinates of its points: H p (x) = P ∈E[p]et\0 E (x -x(P )), and its lift Hp is the unique étale lift dividing the p-division polynomial Ψ p (E ↑ ).

(3) Compute the isogeny E ↑ → E ↑ / Hp using Vélu's formula, and an isomorphism u between E ↑ / H and Σ(E ↑ ). Since Vélu's isogenies are normalised, applying Σ to this isomorphism u gives (up to a sign) the action λ 0 of π on the tangent spaces dx/y and Σ(dx/y) of E ↑ and E ↑σ . (4) Compute the norm λ = N Qq/Fq (λ 0 ). This recovers the inversible eigenvalue of the big

Frobenius π q at precision m, up to a sign. The correct sign is chosen using Hasse's invariant. The trace is then given by t = λ + q/λ, and if m ⩾ (n + 5)/2, the value of t at p-adic precision m is enough to recover t in Z. Then χ π (x) = x 2 -tx + q.

The modular polynomial Φ p (X, Y ) is of total degree p + 1 and its logarithm height is h(Φ p ) ⩽ 6p log p + 18p (see [START_REF] Bröker | An explicit height bound for the classical modular polynomial[END_REF]). Thus its total size is of Õ(p 3 ), and there are quasi-linear algorithms to compute it [START_REF] Enge | Computing modular polynomials in quasi-linear time[END_REF].

Step 1 is done via Newton iterations, the dominating step is evaluating Φ p at precision m in Z q , for a cost of Õ(p 2 m log q) = Õ(p 2 mn).

Step 2 is also done via Newton iterations, the dominating step is evaluating the division polynomial Ψ p (X), which is of degree (p 2 -1)/2 at precision m, for a total cost of Õ(p 2 m log q) = Õ(p 2 mn). Step 3 is dominated by Vélu's formula and costs Õ(pm log q) = Õ(pmn). In Step 4 the norm is done via a resultant, and also costs Õ(pm log q) = Õ(pmn). Since m = O(n), the final complexity of Satoh's algorithm is thus Õ(p 3 + p 2 m log q) = Õ(p 3 + p 2 n 2 ). By constrast, the SEA algorithm (in particular the version of [START_REF] Lercier | On elkies subgroups of ℓ-torsion points in elliptic curves defined over a finite field[END_REF] which works in all characteristic) has a complexity of Õ(n 4 ), so Satoh's algorithm has better complexity for small p and large n. We note that the complexity of Õ(p 3 ) comes from the computation of Φ p (x, y). This polynomial only depends on p, not on the elliptic curve, so this part may be seen as a precomputation, and the real complexity of Satoh's algorithm is Õ(p 2 n 2 ). Alternatively one could use the techniques of [Rob21, § 5.3.8] to evaluate Φ p directly.

In 2002, given an affine equation f (x, y) = 0 of E, Kedlaya proposed in [START_REF] Kedlaya | Counting points on hyperelliptic curves using Monsky-Washnitzer cohomology[END_REF] to use the Monsky-Washnitzer cohomology associated to A † = Q q ⟨⟨x, y⟩⟩/ f (x, y). The difference between these two p-adic methods is the unicity of the canonical lift in Satoh's method in contrast to Kedlaya's method where the lift is arbitrary. Kedlaya's approach [START_REF] Kiran | Counting points on hyperelliptic curves using monsky-washnitzer cohomology[END_REF] thus computes a non-specific lift with linear complexity in p and then reconstructs χ with complexity in time (and space) of O p (n 3+ϵ ). Havey in [START_REF] Harvey | Kedlaya's algorithm in larger characteristic[END_REF] improved the dependency on p of Kedlaya's algorithm. More precisely he shows that Kedlaya's original algorithm can compute the Frobenius to padic precision m with a complexity of Õ(pn 2 m), and Harvey improves the dependency on p to Õ( √ pn 5/2 m + n 4 m log p) (at the cost of a worse dependency on m).

It is such natural to ask whether there exists an algorithm that has the Õp (nm) quasi-linear complexity of Satoh's algorithm with respect to n and the precision m but improves the Õ(p 2 ) dependency on p (which is even Õ(p 3 ) if we take into account the precomputation of the modular polynomial when we don't use the direct evaluation strategy of [START_REF] Robert | Efficient algorithms for abelian varieties and their moduli spaces[END_REF]) to Harvey's Õ( √ p).

Isogeny based key exchange protocols rekindled the interest of the second author on computing canonical lifts to high precision m. (We stress that so far we are not aware of applications other than point counting, which only require a precision m = O(n).) He proposed in [START_REF] Robert | Efficient algorithms for abelian varieties and their moduli spaces[END_REF]Chapter 6] a new approach of Satoh's method which works by only using the modular polynomial Φ p to both lift the curve and the isogeny. This allows to dispense with the computation of the division polynomial Ψ p , but does not change the asymptotic because of the evaluation of the modular polynomial, so the algorithm is still in Õ(p 3 + p 2 nm log p) (although with better constants).

He proposed another method bypassing the need for the modular polynomial in [Rob21, Remark 6.6.2], assuming a point of p-torsion is given on E. More recently, he realized that by working modulo H p directly allows to define a "formal" point of p-torsion to which to apply [Rob21, Remark 6.6.2]. This allows to bypass the need to find a point of p-torsion. These two strategies (working with a rational point of E et [p] or with its formal point defined by H p ) were implemented in Pari/GP by the first author, who also carefully tracked the loss of precision. The implementation showed that the resulting algorithms not only improve the theoretical complexity but are also practical. This is the subject of the present work.

Indeed, we can keep the Õp (nm) complexity of Satoh's algorithm while improving the dependency on p.

Theorem 1.1. Let E/F q be an ordinary elliptic curve, with j(E) ̸ ∈ F p 2 . Then one can compute the canonical lift E ↑ /Z q and the trace of the Frobenius to p-adic precision m in time Õ(mnp).

In particular, for point counting where we need m = O(n), the complexity to compute χ π is Õ(pn 2 ).

The main idea behind Theorem 1.1, is that when doing a Newton lift to lift the root of a polynomial F (X), it is not necessary to be given F , one only needs to be able to evaluate it. We can thus circumvent computing the modular polynomial Φ p and the division polynomial Ψ p in Satoh's algorithm by directly evaluating isogenies (ie solving the equation j(E ν ) = j(E Σ) where E ν is computed via an isogeny) and the multiplication by [p] map.

Although we do not reach Harvey's Õ( √ p) complexity, in some cases a variant of our method achieve such a complexity.

Theorem 1.2. Let E/F q be an ordinary elliptic curve with j(E) ̸ ∈ F p 2 , and assume that we are given a rational étale point of p-torsion P . Let I(d, m, Z q ) be the cost of evaluating at precision m an isogeny of degree d on an elliptic curve E ′ over Z q given a generator P (defined over Z q ) of its kernel. Here by evaluating the isogeny, we only mean computing the equations of E ′ / < P > at precision m.

Then one can compute the canonical lift Ẽ/Z q and the trace of the Frobenius to p-adic precision m in time Õ(mn log p + I(p, m, Z q )).

Remark 1.3. We can also work on the Kummer line E/ ± 1, that is given only the x-coordinate x P of our point P , in which case I(d, m, Z q ) should be the cost of evaluating the induced isogeny E/ ± 1 → (E/ < P >)/ ± 1. Using Vélu's formula, we have I(p, m, Z q ) = Õ(pm log q) = Õ(pmn). A recent improvement of Vélu's formula [START_REF] Bernstein | Faster computation of isogenies of large prime degree[END_REF] improves this complexity to Õ( √ pm log q) = Õ( √ pmn).

In general, the étale points of p-torsion will live in an extension of degree e ⩽ p -1 (which we can compute using Hasse's formula), and to find one of them require computing a root of the division polynomial Ψ p modulo p, which can be done in time Õ(p 2 log q + p log 2 q) = Õ(np 2 + pn 2 ). We explain how to improve this complexity to Õ(p log 2 q) in Section 4.2 and get: Corollary 1.4. Let e be the degree of the extension of F q where the étale points of p-torsion lives. Then one can compute the canonical lift Ẽ/Z q and the trace of the Frobenius to p-adic precision m in time Õ(p log 2 q + √ pme log q) = Õ(pn 2 + √ pmen). (If χ π is already known the complexity becomes Õ(log 2 (q e ) + √ pme log q) = Õ(e 2 n 2 + √ pmen).)

In the worst case, e = O(p) so the complexity is Õ(p log 2 p + p 3/2 m log q) = Õ(pn 2 + p 3/2 mn) is not better than Theorem 1.1. In the best cases, when e = O(log p); for instance if the trace t = 1 mod p (which implies e = 1); the complexity is Õ(p log 2 q + p 1/2 m log q) = Õ(pn 2 + p 1/2 mn). In general, to compute Ẽ at high p-adic precision, we improve on the complexity of Theorem 1.1 whenever e = O( √ p).

We organize this paper as follow. In Section 2, we recall the Serre-Tate theorem and Satoh's algorithm. We present our new approach to Newton lifts in Section 3. As a first application we explain how to lift the p-torsion in Section 4, then we give our canonical lift algorithm in Section 5.

1.1. Notation and Convention. In the following p is prime and q = p n with n ⩾ 1. We denote by Q q the unramified extension of the field of p-adic numbers Q p and by Z q is the valuation ring of Q q ; it is also the ring W (F q ) of the Witt vectors over F q . The extension Q q /Q p has a cyclic Galois group of order n, generated by an element Σ that reduces to the (small) Frobenius automorphism σ on the residue field F q . The large Frobenius (and its lift) will be denoted by σ q and Σ q respectively, and sometime we will denote σ by σ p to emphasize we work with the small Frobenius. As a convenience we let σ = σ -1 = σ n-1 , Σ = Σ n-1 , the "Verschiebung" Galois elements.

Explicitly

Q q = Q p [X]/M (X) and also Z q = Z p [X]/M (X)
with M is monic irreducible polynomial of degree n over Z p [X] with irreducible reduction modulo p. The complexity of an elementary operation require Õ(m log q) = Õ(mn) with Kronecker-Schönhage method at precision m. By p-adic precision m, we mean that we are working modulo p m Z q . Furthermore, fast modular composition [START_REF] Kiran | Fast polynomial factorization and modular composition[END_REF] allows to efficiently evaluate Σ and Σ in Õ(nm); it also allows to evaluate σ in Õ(log q) = Õ(n) rather than the slower Õ(n log q) = Õ(n 2 ) we get iterating the Frobenius n -1 times. It is also convenient to take for M the Teichmuller lift of an irreducible polynomial M (X) of degree n over F p , this allows for a fast computation of Σ without invoking modular composition.

We recall the Frobenius σ p induces an isogeny π p : E → E σ , P → P σ , and σ q induces an endomorphism π q . The Verschiebung π : E → E σ is the dual of π p : E σ → E (we warn that it is not given on points by P → P σ !). Both the Frobenius and Verschiebung lift uniquely to the canonical lifts, we denote them by π and π. In this article, E ↑ we always denote the canonical lift of E, while Ẽ will denote a candidate lift (which may or may not be canonical).

Background

Let E/k be an elliptic curve, and Ψ ℓ its polynomial of ℓ-torsion (or ℓ-division polynomial) associated with the equation of the curve. A point P = (x, y) on E is a point of ℓ-torsion if and only if its coordinates constitute a solution of Ψ ℓ . An isogeny ϕ is a non trivial morphism between elliptic curves which is also a group morphism. The multiplication morphism is identified with Z then Z ⊂ End(E). Furthermore when the base field k is F q we have: Z[π q ] ⊂ End(E) where π q is the Frobenius endomorphism. In the case where E is ordinary: χ(X) = X 2 -tX + q is the characteristic polynomial of π q where t is the trace of π q and verifies the relation |t|⩽ 2

√ q called Hasse's bound. Therefore, if we set

D πq = t 2 -4q < 0 then : #E(k) = q + 1 -t and Z[π q ] ⊂ End(E) ⊂ O K where K = Q[ D πq ].
2.1. Vélu's Algorithm. According to the inputs, the algorithms for calculating isogenies can be classified into two large groups. The first ones initiated by Vélu [START_REF] Vélu | Isogénies entre courbes elliptiques[END_REF] takes an elliptic curve E and a subgroup K of E then outputs an explicit form of isogeny ϕ : E -→ E/K and an equation of E/K. Then for every P ∈ E:

x ϕ(P ) = x P + Q∈K\{O} (x P +Q -x Q ) and y ϕ(P ) = y P + Q∈K\{O} (y P +Q -y Q ) .
Considering the improvements made by D. Kohel [START_REF] Kohel | Endomorphism rings of elliptic curves over finite fields[END_REF] we arrive at the same results when K = ker ϕ is represented by a polynomial h.

Example 2.1. When char(k) > 3 and an elliptic curve E over k is given by E :

y 2 = f (x) = x 3 + a 4 x + a 6 .
Let h be the polynomial defines the kernel K of a separable normalised isogeny ϕ of degree ℓ with domain E Set:

Q(x) = gcd (f (x), h(x)) D(x) = h(x) 2 /Q(x) = x ℓ-1 -d 1 x ℓ-2 + d 2 x ℓ-3 -d 3 x ℓ-4 + • • •
Then for every point P (x, y) in E we have:

ϕ(x, y) = (α(x), yα(x))
where

α(x) = ℓx -d 1 x -(3x 2 + a 4 ). D ′ (x) D(x) -2f (x). D ′ (x) D(x) ′ .
And E/K is given by the equation:

y 2 = x 3 + (a 4 -5v)x + (a 6 -7w) where v = a 4 (ℓ -1) + 3(d 2 1 -2d 2 ) and w = 3a 4 d 1 + 2a 6 (ℓ -1) + 5(d 3 1 -3d 1 d 2 + 3d 3 ).
On the other hand, the modular polynomial Φ p encodes directly the j-invariants of isgogeneous elliptic curves.

Lubin Serre Tate theory.

Theorem 2.2. (Lubin-Serre-Tate) Consider E an ordinary elliptic curve over F q , then there exist a unique elliptic curve up to isomorphism E ↑ over Z q such that.

• E is the reduction of E ↑ modulo p , • End(E ↑ ) ∼ = End(E), E ↑ is called the canonical lift of E,
and is also uniquely characterised by the fact that the Frobenius π q lift to E ↑ , or that π p lift to an isogeny E ↑ → E ↑Σ , ie by the equation

Φ p (j(E ↑ ), j(E ↑ ) Σ )) = 0. E ↑ E ↑Σ E E σ π π
We refer to [START_REF] Lubin | Elliptic curves and formal groups[END_REF] for the statements (without proofs) and [START_REF] Messing | The crystals associated to barsotti-tate groups, The crystals associated to Barsotti-Tate groups: with applications to abelian schemes[END_REF] for proofs.

Remark 2.3. If E ↑ /Z q is an elliptic curve, then it is the Néron model of its generic fiber E ↑ η . Furthermore, by the property of Néron models, E ↑ (Z q ) = E ↑ η (Q q ). Hence it is harmless to consider the curve over Q q .

2.3. Modular Equations. We present Satoh's (as improved by Harley) method to compute the canonical lift Ẽ using the modular relation Φ p (j(E ↑ ), j(E ↑ ) Σ )) = 0.

We know that the modular polynomial satisfy the Kronecker's relation:

Φ p (X, Y ) ≡ (X p -Y )(X -Y p ) mod p
Let j / ∈ F p 2 , the following statement is an immediate consequence of Kronecker's relation:

     ∂Φ p ∂X (j, j σ ) ≡ j p -j p ≡ 0 mod p ∂Φ p ∂Y (j, j σ ) ≡ j p 2 -j ̸ ≡ 0 mod p
Thus we deduce that the Frobenius σ has multiplicity 1 and the Verschiebung σ has multiplicity p. In fact we have over Q q two points (P, Q) on E ↑ [p] that reduce to P and 0 respectively. The p kernels ⟨P + kQ⟩ with 0 ⩽ k < p reduce to P (ie the kernel of the Vershiebung) and the last ⟨Q⟩ reduces to ⟨O⟩ (ie the kernel of the Frobenius). So we have p-isogenies on E ↑ which reduces to the Verschiebung π modulo p. A more detailed analysis show that they reduce to different isogenies modulo p 2 , hence:

Lemma 2.4. Let Ẽ/Z q be any lift of E/F q where j(E)

̸ ∈ F p 2 . Then ∂Φ p ∂X (j( Ẽ), j( Ẽ) Σ ) is of valuation 1. Proof. By [Nak93, Proposition 2], since j( Ẽ) ̸ = 0, 1728, Φ p (j( Ẽ), X) = (X -j( Ẽ) Σ )(X - j( Ẽ) Σ)G(X ) where G(X + j( Ẽ) Σ) is an Eisenstein polynomial. Since j( Ẽ) ̸ ∈ F p 2 , j( Ẽ) Σ ̸ = j( Ẽ)
Σ and the result follows.

□

This provides an algorithm to compute the lifted j-invariants of the p-isogenous curves Ẽ and ẼΣ .

We want to solve in Z q the equation Φ p ( j, jΣ ) = 0 knowing j modulo p. Suppose that we can compute efficiently the Frobenius Σ of Q q and j ∈ Z q is an approximation of j at precision k i.e j = j + p k e for some error e ∈ Z q that we want to find. Using the modular equation and Taylor expansion of Φ p we have:

0 = Φ p (j + p k e, j Σ + p k e Σ ) 0 = Φ p (j, j Σ ) + p k e ∂Φ p ∂X (j, j Σ ) + p k e Σ ∂Φ p ∂Y (j, j Σ ) + p 2k (...)
Dividing by p k , we get

u + e ∂Φ p ∂X (j, j Σ ) + e Σ ∂Φ p ∂Y (j, j Σ ) ≡ 0 mod p k .
If j / ∈ F p 2 , the Kronecker inequality implies that ∂Φ p ∂Y (j, j Σ ) ≡ 0 mod p and ∂Φ p ∂Y (j, j Σ ) ̸ ≡ 0 mod p. Then to have the error e we must solve over Z q the following equation:

e Σ + Ae + B = 0.
with A ≡ 0 mod p and B ̸ = 0 mod p called "Artin-Schreier equation"in [START_REF] Gaudry | Algorithmes de comptage de points d'une courbe définie sur un corps fini[END_REF]. Set e = x + p k α with α ∈ Z q , the error α can be determine using algorithm 2.3 (a general case of Harley's algorithm). In odd characteristic the short Weierstrass equations have two parameters that we denote A and B. Using the relation between the j-invariant and those parameters, given the j-invariant j( Ẽ), one can lift the equation of the elliptic curves defined over F q to Z q . Take an arbitrary lift of one parameter, then the equation between the lifted j-invariant and the second parameter provide a simple Newton algorithm to lift it. Furthermore in characteristic ⩾ 5 Skjernaa [START_REF] Skjernaa | Satoh's algorithm in characteristic 2[END_REF] has suggested to simply take A = 3λ and B = 2λ with

λ = j(E) 1728 -j(E)
This method is faster than the first. It needs only one inversion in Z q form the lifted j-invariant.

The division polynomial.

If E/k is an elliptic curve with a short Weierstrass equation, and P = (x, y), then ℓ.P = ( ξ ℓ (x)

ψ 2 ℓ (x) , ω ℓ (x,y) ψ 3 ℓ (x,y) )
where ξ ℓ and ω ℓ are expressible in terms of the ψ ℓ-2 , ψ ℓ-1 , ψ ℓ , ψ ℓ+1 , ψ ℓ+2 , and the ψ ℓ satisfy a recurrence relation expression ψ 2ℓ and ψ 2ℓ+1 in term of he ψ ℓ-2 , ψ ℓ-1 , ψ ℓ , ψ ℓ+1 , ψ ℓ+2 . In practice, the recurrence formula simply come from computing ℓ.P formally via the double and add algorithm. In particular, when ℓ is odd, the roots of ψ ℓ (x) are exactly the elements x(P ) for P ∈ E [ℓ].

In this article we will use a slightly different version of the division polynomial: we let Ψ ℓ (x) = ψ ℓ (x) when ℓ is odd, and Ψ ℓ (x) = ψ ℓ (x)/2y when ℓ is even. This reformulation is such that Ψ ℓ is always in k[x] whether ℓ is even or odd. It is easy to adapt the recurrence formula to compute the Ψ ℓ directly.

In the following, we will need to compute Ψ p (x) and Ψ ′ p (x) for an elliptic curve Ẽ/Z q (at precision m) modulo a polynomial H of degree d. In practice d will be equal to 1 when we want to evaluate Ψ p on a point x P (so H = (x -x P )), or d will be equal to (p -1)/2 when we want to evaluate Ψ p modulo Hp a candidate lift of H p .

We remark that we can evaluate Ψ p modulo H simply by evaluating the recurrence relation modulo H. Also from the recurrence relation on Ψ ℓ , we get a recurrence relation on Ψ ′ ℓ , so we can also evaluate it modulo H. We obtain Lemma 2.5. Given an elliptic curve Ẽ/Z q and a monic polynomial H(x) of degree d, we can evaluate Ψ Ẽ,p and Ψ ′ Ẽ,p modulo H at precision m in time Õ(dm log q log p) = Õ(dmn).

2.6. Lifting the Verschiebung. Since the Frobenius π q is inseparable, we lift the Verschiebung π p over Z q by lifting its kernel.

We set E n-i = E σ i and π i is the isogeny between E i+1 and E i defined by (x, y) → (x σ , y σ ).

Then the Verschiebung π q decomposes as follow:

π q = π n-1 π n-2 • • • π 0 .
ker( π) is a subgroup of order p of E[p] defined by the monic separable factor H p of the p-division Ψ p given by :

H p (x) = P ∈ker π\{O} (x -x(P ))
Let Hp be the lift of H p over Z q , then Hp is a monic facteur of degree (p -1)/2 of Ψ p on Ẽ and Hp (x) = H p (x) mod p is square free. Furthermore Ψ p (x) ≡ H p (x) p mod p i.e modulo p, the factors H p (x) and Ψ p (x)/H p (x) are not coprime modulo p. T.Satoh introduced in [Sat00, § 2] a variant of Hensel's lift that compute Hp over Z q . Let p be an odd prime, and suppose that we have a polynomial G in Z q [X] and h ∈ F q [X] a monic factor of the reduction of G modulo p. We assume that h(x) is separable and relatively prime with p -t G ′ (x) where t = ord p (G ′ (x)). Let u ∈ N be such that G(x) ≡ q(x)h(x) mod p u+t . Then the polynomial :

H(x) = h(x) + G(x) G ′ (x) h ′ (x) mod h(x)
is a lift of h(x) at precision p 2u and G(x) ≡ Q(x)H(x) mod p v where v = 2u + min(t, u) (see [START_REF] Satoh | The canonical lift of an ordinary elliptic curve over a finite field and its point counting[END_REF]). This property provides an algorithm constructing a lift h with O((deg h + deg G) 2 ) arithmetic operations over Z q at precision O(n). Satoh then applies this construction to lift H p , by [Sat00, Lemma 3.7], in this case t = 1. An alternative method when we are provided an étale point P of p-torsion is to lift the equation (p ′ + 1).P = p ′ .P where p = 2p ′ + 1 as in [MR20, Proposition A.7.], or to work with only the x-coordinate to simply use the standard Newton method to lift Ψ p (x p ) = 0. This is faster than the euclidean extended GCD used in Satoh's formula above, we will revisit this in Section 4.

Application to point counting. When we have E ↑ at sufficient precision m (given by

Hasse-Weil bounds), one can evaluate the action of the Verschiebung on the differential form dx y as detailed by Satoh's diagram.

E ↑ E ↑ Σ E ↑ν = E ↑ / K Σ ν u
Here the isogeny ν is computed by Vélu's algorithm from the lift Hp of the kernel of the Verschiebung.

Since the isogeny ν is normalized, the action of the isogeny π on the differential form of E ↑ is given by the isomorphism ±u on E ↑ν ; let us denote it by λ 1 . Concretely, we have π = ±u • ν, and if u(x, y) = (u 2 x, u 3 y), λ 1 = ±u.

On the other hand, when we consider the q th -power Frobenius morphism decomposition:

E ↑ -→ E ↑Σ -→ • • • -→ E ↑Σ n-1
The action on the differential forms along the cycle will be given by the successive conjugates of λ 1 . Finally, by composition, the action of the dual endomorphism πq of πq on the main differential form of Ẽ is given by the product of all these conjugates, i.e. by the norm of λ 1 . On the other hand the norm of N Qq/Qp (λ 1 ) is simply given as the resultant of λ

1 modulo M (X) in Q p [X]
. This method due to Harley can be asymptotically done in quasi-linear time in the precision m using a fast GCD algorithm [CFA + 06]. A slower alternative is to use the formula Input Coefficients (A4, A6) of E an elliptic curve of F q with q = p n , n ∈ N. Output The Trace of Frobenius endomorphism of E.

Using algorithm 5.1, compute E ↑ at precision m = (n + 5)/2 ; Compute the action λ 1 of an isomorphism u : E ↑ν → E ↑ Σ ; Compute λ 2 = N Qq/Qp λ 2 1 ; Compute λ the correct square root from λ 2 and t = λ + q/λ mod q such that |t| < 2 √ q ;

Return χ(X) = X 2 -t • X + q .

Algorithm 2.2 Computing the characteristic polynomial of ordinary elliptic curve E N Qq/Qp (c) = exp Tr Qq/Qp (log c) using a specific implementation to compute it in time O(m 3/2 n) (available in [START_REF] Vercauteren | Pari/gp, version2.12.1[END_REF]). Since we only have λ 1 up to a sign, taking its norm λ and then computing the trace t = λ + q/λ only give t up to a sign. One can use Hasse's invariant to get the correct sign, see Section 4.2. Let χ(X) = X 2 -t • X + q be characteristic polynomial of the Frobenius of E, Hasse-Weil bound states |t| < 2 √ q. On the other hand we have #E(F q ) = χ(1). Then we deduce the following result:

Theorem 2.6. Let E/F q be an ordinary elliptic curve. Given the canonical lift E ↑ /Z q and the lift H ↑ p of the kernel of the Verschiebung to precision m, one can compute the trace of the Frobenius to p-adic precision m in time Õ(mnp).

In particular, for point counting where we need m = O(n), the complexity to compute χ π once we have E ↑ and H ↑ p to precision m is Õ(pn 2 ).

In the rest of this paper, we will explain how we can compute E ↑ and H ↑ p to precision m in time Õ(mnp) rather than in time Õ(mnp 2 ). By Theorem 2.6, this will show that we have a point counting algorithm in time Õ(n 2 p).

We also remark that we can bypass the computation of E ↑ Σ (since Σ is typically more expensive to compute than Σ) by applying the above method to E ↑Σ , the canonical lift of E σ instead.

Revisiting Newton's method

Let F (X) be a multivariate polynomial system defined over Z q , and suppose that we have a solution x modulo p (in other words, at precision 1) of the equation F (x) = 0 (modulo p). Assume furthermore that dF (x) is inversible modulo p. Then there is a unique lift x of x in Z q such that F (x) = 0 and x = x modulo p. Newton's method show that x can be approximated by the sequence (1)

x 0 = x, x 2k = x k -dF (x k ) -1 F (x k ).
A standard computation shows that x k approximates x to precision m = 2 k and that F (x k ) = 0 modulo p m .

Our trivial, but key remark which is at the core of this article, is that to use Newton's method we do not need to know F , we only need to be able to evaluate F at some precision m. Indeed from Equation (1) it is clear that we only need to be able to evaluate F and dF . But we can recover dF from evaluations of F at suitable points.

We illustrate this when F (X) is univariate. Then modulo p 2m , F (x + p m y) = F (x) + F ′ (x)p m y, hence F ′ (x) = (F (x + p m ) -F (x))/p m modulo p m . We can thus recover F ′ (x) modulo p m from two evaluations of F at precision 2m. The Newton process can thus be done as follow: given the solution x m at precision m, we evaluate F (x m ) and F (x m + p m ) at precision 2m. Then

x 2m = x m - F (xm) (F (xm+p m )-F (xm))/p m .
More generally, when F has N -variable, we can recover the Jacobian dF (x) at precision m in N + 1 evaluations of F at precision 2m.

We have proved:

Lemma 3.1. Given a multivariate polynomial system F (X) in N variables and N equations, and a solution x 0 modulo p of the equation F (x) = 0 modulo p such that dF (x 0 ) is inversible modulo p. Let C(m, Z q ) be the cost of evaluating F at a point x at precision m and M (m, Z q ) be the cost of doing the standard arithmetic operations in Z q at precision m, and assume that both C(m) and M (m) are superlinear.

Then one can compute the unique lift x of x 0 such that

F (x) = 0 to precision m in time O(N • C(2m, Z q ) + N • M (2m, Z q )).
Remark 3.2. We note that if we have an approximation x 0 of x to precision m, then for our method (and the convergence), we only need that F is analytic at x 0 on the ball of center x 0 and radius ||p m ||.

More generally, Newton's algorithm will converge whenever we have a x 0 modulo p e+1 such that f (x 0 ) = 0 modulo p 2e+1 and p e dF (x 0 ) is inversible. Iterating the Newton process then gives x modulo p e+2 k such that f (x) = 0 modulo p 2e+2 k .

When this is not the case, we need to push the Taylor expansion of F further:

F (x + ep k ) = p k dF (x). t e i + p 2k e i • d 2 F (x, x) • t e i + O(p 3k ).
Let J(x) = dF (x) be the Jacobian, and H(x) = d 2 F (x, x) be the Hessian matrix, we explain how to evaluate them to precision m. We assume here for simplicity that N = 2 and p > 2.

Set e 1 = (1, 0), e 2 = (0, 1) and e 5 = (1, 1), set

x 1 = x + e 1 p m , x 2 = x + e 2 p m , x 3 = x -e 1 p m , x 4 =
x -e 2 p m and x 5 = x + e 5 p m , and evaluate F (x i ) modulo p 3m . We have modulo p 2m :

J X (x) = F (x 1 ) -F (x 3 ) 2p m and J Y (x) = F (x 2 ) -F (P x )] 2p m
and modulo p m :

H X (x) = F (x 1 ) -J X (x)p m p 2m , H Y (x) = F (x) -J Y (x)p m p 2m H XY (x) = F (x 5 ) -F (x) -J X (x)p m -J Y (x)p m -H X (x)p 2m -H Y (x)p 2m p 2m .
More generally, in N variables we may compute d k F (x) at precision m by O(k + N k+1 ) evaluations of F at precision km when p is large enough.

Lifting the étale points of p-torsion

Let Ẽ/Z q be a (non necessarily canonic here) lift of an ordinary elliptic curve E/F q . In this section, we explain how to compute the polynomial H p which parametrizes E[p] ét (this is also the kernel of the Verschiebung) and how to lift it to Ẽ using Section 3. We also explain how to find an étale point P ∈ E[p] and how to lift it.

We first note that when Ẽ is an arbitrary lift of E, there is an obstruction to lifting an étale point of p-torsion P : in general Ẽ[p] may not have points living in an unramified extension of Z q , in particular even if P is rational, a lift of P to Ẽ[p](Z q ) will not exist. This obstruction vanish if and only if j( Ẽ) = j(E ↑ ) mod p 2 : Proposition 4.1. Let Ẽ/Z q be an arbitrary lift of E/F q , and let E ↑ /Z q be the canonical lift of E. Let Z ur q be the maximal unramified extension of Z q . The following are equivalent:

(1) Ẽ[p](Z ur q ) ̸ = {0 E }; (2) Ẽ[p](Z ur q ) is a lift of E[p] et ; (3) Z q ( Ẽ[p]) is tamely ramified; (4) j( Ẽ) = j(E ↑ ) mod p 2 .
If these conditions are satisfied, and F q e is the smallest extension of F q where the points of

E[p] et are defined, then Ẽ[p] = Ẽ0 [p] ⊕ Ẽ[p](Z q e )
where Ẽ0 is the relative connected component of E (ie the kernel of the reduction map Ẽ → E), and the points of Ẽ0 [p] live in the tamely ramified extension of Z q e given by adjoining a p-root of unity ζ p . Furthermore, if P ∈ Ẽ[p](Z ur q ), then

Ψ ′ P (P ) is of valuation 1.
Proof. We have a connected étale exact sequence [START_REF] Tate | Finite flat group schemes, Modular forms and Fermat's last theorem[END_REF]:

0 → Ẽ0 → Ẽ → Ẽet → 0.
This exact sequence commutes with specialisation, so since Spec Z q is connected, Ẽ0 is exactly the kernel of the projection map Ẽ → In particular, since Ẽet [p] is étale and Z q is complete hence Henselian, it is the unique étale lift of E[p] et , and Ẽ0 [p] is a lift of E[p] loc which is of multiplicative type (since it is the Cartier dual of E[p] et ), hence its points live in a ramified extension of Z q . So if P ∈ Ẽ[p](Z ur q ) ̸ = 0 Ẽ , the subgroup generated by P induces a a splitting Ẽ

[p] = Ẽ0 [p] ⊕ Ẽet [p], in particular P ̸ ∈ Ẽ0 [p]
. This proves the equivalence of (1) and (2). The rest of the equivalences are from [Sat00, Theorem 3.1]. Furthermore we have Ẽ[p](Z ur q ) = Ẽet [p](Z ur q ) = Ẽet [p](Z q e ), and since Ẽ0 [p] is the Cartier dual of Ẽet its points live in Q q (ζ p ).

Finally, Ψ ′ p (P ) = 1 by Satoh's lemma [Sat00, Lemma 3.7]. □ 4.1. Computing the kernel of the Verschiebung. To apply Section 2.6, we first need to compute the kernel H p of the Verschiebung (or a rational point in this kernel). We have Ψ p = H p p , so an easy method is to compute Ψ p using the recursive formula for division polynomials to get H p . But Ψ p is of degree p 2 , so this will cost Õ(p 2 n) operations.

Let π be the Verschiebung. By definition [p] = ππ = ππ, so we have π(π(P )) = [p].P . In particular we can efficiently evaluate the Verschiebung on the point π(P ). We can thus recover the Verschiebung by interpolation, from which we get the kernel.

More precisely we only need to work with x-coordinates. We can then sample p-random points x p ∈ E σ (F q )/ ± 1, and compute the values p.x p in x-coordinates only. Let R(x) be the rational fraction of degree O(p) interpolating the points (π(x p ), p.x p ). Then the kernel H p of the Verschiebung is simply the denominator of R.

In summary:

Lemma 4.2. Let E/F q be an ordinary elliptic curve. The kernel H p (x) of the (small) Verschiebung can be computed in time Õ(p log q) = Õ(pn).

Finding an étale point of torsion.

If we furthermore need the x-coordinate of an étale point P of p-torsion, we need to find a root of H p . First we need to compute the degree e of the extension where the étale points of p-torsion live. Assume that we know λ, the inversible eigenvalue of the Frobenius modulo p. Then σ(P ) = λ.P , so e is the order of λ.

There are two methods to find λ to precision 1. The first one is to use Hasse's formula. Using the recurrence formula to compute the Hasse invariant A q (see [Sil86, V.4.1]), this costs Õ(n 2 +np) operations: Õ(np) to compute A p , then Õ(log 2 q) = Õ(n 2 ) to compute A q .

The other approach evaluates the Verschiebung from its kernel H p using Vélu's formula, and look at the action on the differentials (ie we apply Satoh's algorithm at precision m = 1, so without lifting), as in Section 2.7. This costs Õ(np) operations, but this only recovers ±λ.

Indeed, we compute an isomorphism u : E/H p ≃ E σ , (x, y) → (u 2 x, u 3 y), so if ϕ : E → E/H p is given by Vélu's formula, the Verschiebung is equal to ±u • ϕ. To know the correct sign, we need to stop working in x-coordinate only and take a random point P ∈ E(F q ) and check whether [p]σ(P ) = u • ϕ(P ) or [p]σ(P ) = -u • ϕ(P ). Then replacing u by -u if necessary, we have that λ = N Fq/Fp (u) since ϕ is normalised. To take P we need to compute a square root, so this costs Õ(log 2 q) = Õ(n 2 ), and the total cost to recover e exactly (rather than potentially 2e) is Õ(n 2 + np) operations, like the computation of the Hasse invariant.

The factorisation of H p using an equal degree factorisation algorithm then costs Õ(p log 2 q) = Õ(pn 2 ). We note that without the knowledge of e, we would need to use a distinct degree factorisation algorithm insted, which would cost Õ(p 1.5 n + pn 2 ) by [START_REF] Kiran | Fast polynomial factorization and modular composition[END_REF]. In summary: Lemma 4.3. Let E/F q be an ordinary elliptic curve. The kernel H p (x) of the (small) Verschiebung can be computed and factorized in time Õ(pn 2 ).

There is a faster method when we already know N = #E(F q ). Compute e as above, and N e = #E(F q e ). Take a random point Q ∈ E(F q e ) and multiply by the cofactor: P = N e /p • Q. If P ̸ = 0 0 we have found a point of p-torsion. A random point Q can be taken by taking a random x Q and trying to find a square root of x 3 Q + ax Q + b. We can also work in x-coordinates only, this gains a square root. In any case, the total cost of this method is Õ(log 2 q e ) = Õ(e 2 n 2 ).

Lifting a point of p-torsion.

We now assume that we are given a lift Ẽ that satisfy the equivalent conditions of Proposition 4.1.

Given a point P of p-torsion on E, to lift it to Ẽ we apply Lemma 3.1 to the equation p. P = 0 E . To stay in affine coordinates, we can rewrite this equation as (p ′ + 1). P = -p ′ . P for p = 2p ′ + 1.

Evaluating this equation by a double and add algorithm takes O(log p) operations in Z q (at a given precision m), hence by Lemma 3.1 we can compute P to precision m in time Õ(nm).

Remark that the p-torsion P points is defined equivalently by systems of the form:

f (x, y) = 0 Ψ p (x) = 0 or f (x, y) = 0 g(x, y) = 0 where g(x, y) is one of the equation [p ′ + 1]P = -[p ′ ]P such that p = 2p ′ + 1.

Since p ̸ = 2, we have ∂f ∂y (P ) non null modulo p. The Jacobian of the system is given by:

  ∂f ∂x (x, y) ∂f ∂y (x, y) Ψ ′ p (x) 0  
Then using Satoh's lemma [Sat00, Lemma 3.7], we conclude that the valuation the determinant of the Jacobian of those system at P is 1: it has the form ⋆ ⋆ p 0 at P Thus in the Newton's lifting steps, we lost 1 precision on the coordinates of P . Also to boostrap at precision 1, it seems like we would need to compute the Hessian. Fortunately, the situation simplifies if we only try to lift the x-coordinate of P , the system then becomes Ψ p (x) = 0. We never compute Ψ p but evaluate it on x directly via the double and add formula for x-coordinates (in other words via the standard recurrence formula for the division polynomials) by Lemma 2.5. However, our lift is such that Ψ ′ p (x P ) is of valuation 1 and This subgroup is defined by a polynomial H G (x) (say when ℓ is odd) of degree (d -1)/2. The standard method would be to lift H G (X) as a factor of χ ℓ, Ẽ (X), which would cost Õ(ℓ 2 m log q) at precision m. Our method only computes χ ℓ modulo (potential) lifts of H G , hence only cost Õ(d log ℓm log q) = Õ(d log ℓmn) (where the log ℓ comes from the recurrence formula for Ψ ℓ may not be absorbed in the Õ notation here).

Computing the Canonical Lift without the use Modular Polynomials

In this section we will focus on the case where p is odd for simplicity.

Lemma 5.1. Let E an ordinary elliptic curve over F q , then E ↑ is the unique elliptic curve up to isomorphism over Z q such that.

• E is the reduction of

E ↑ modulo p , • Let K ⊂ E ↑ (Q un q ) be such that K reduces to E[p] et modulo p and ν : E ↑ → E ↑ /K. Then j(E ↑ν ) = j(E ↑ ) Σ.
Proof. Immediate by Section 2 and Theorem 2.2, □

We can then apply Lemma 3.1 to the equation of Lemma 5.1 to compute the j-invariant J ↑ of the canonical lift. We first note that Proposition 4.1 gives a convenient criteria to compute the canonical lift E ↑ to precision 2. Lemma 5.2. Let E/F q be an ordinary elliptic curve, P a point of p-torsion on E, and H p the kernel of the Verschiebung. Let Ẽ/Z q be a lift of E. Then j( Ẽ) = j(E ↑ ) mod p 2 if and only if Ψ Ẽ,p (P ) = 0 modulo p 2 , if and only if Ψ Ẽ,p = 0 modulo (p 2 , H p ).

Proof. We first note that the value of Ψ Ẽ,p (P ) does does not depend on the choice of lift P of P to precision 2 since Ψ ′ p (P ) = 0 mod p. The same hold for Ψ Ẽ,p modulo H p . By Section 4, is Ψ Ẽ,p (P ) = 0 modulo p 2 then Newton's method lifts P to a point of p-torsion on Ẽ, alternatively the existence of a point of p-torsion on Ẽ is given by [Sat00, Theorem 3.1].

The lemma is then a direct application of Proposition 4.1. □ So Ẽ mod p 2 correspond to the unique elliptic curve (up to isomorphism) Ẽ mod p 2 such that Ψ Ẽ,p (P ) = 0 Such we look for Ẽ : y 2 = x 3 + ãx + b mod p 2 such that Ψ Ẽ,p (P ) = 0. Taking ã an arbitrary lift of a, we look for b = b 0 + pb 1 , and we solve for b 1 by using the methods of Section 3. If we have H p instead, we do the same computation using the equation Ψ Ẽ,p %H p = 0 mod p 2 .

Assume that we have J at p-adic precision k ⩾ 2, we want to find it at precision 2k. We assume here that we are given H p , we explain how to adjust the algorithm when we are given a point of p-torsion P instead afterwards.

We let F (X) be the following process (at precision 2k): given x such that x = j(E) modulo p, we construct the elliptic curve E with j-invariant x, we let Hp be the lift of H p to E, and E ν the isogenous variety E/ Hp . Then x = j(E) is the lift we look for whenever F (x) = x Σ.

We can evaluate F (X) using Vélu's formula and Section 4, hence we can also evaluate F ′ (X) by Section 3. Lemma 5.3. Let J satisfy F (J) = J Σ at precision k ⩾ 2, and take an arbitrary lift at precision 2k. Let A = F (J) and B = F ′ (J). Then F (J + ep k ) = A + Bp k e, where B is of valuation -1.

Proof. By definition of F , if J ′ = J + p k e, we have Φ p (J ′ , F (J ′ )) = 0 modulo p 2k . Write F (J ′ ) = F (J) + p k e ′ , then Φ p (J, F (J)) + ∂Φ p /∂ x (J, F (J))p k e + ∂Φ p /∂ y (J, F (J))p k e ′ = 0 modulo p 2k , Input E an elliptic curve of F q with q = p n , n ∈ N. Output The canonical lift of E at precision m.

Compute H p over F q using Section 4.1 ; Compute E ↑ at precision 2 using the equation Ψ ↑ p mod H p = 0 at precision 2 ; Compute H ↑ p mod p 2+1 . k = 1 while k < ⌈(m -1)/2⌉ ; a. Compute at precision 2(k + 1) two lifts of the curve E ↑ mod p k+1 ; b. Compute Hp mod p 2k+1 on these two curves using Section 4.4. c. Solve the equation j( Ẽν ) = j( Ẽ) Σ to recover the curve E ↑ mod p 2k+1 ; d. k = 2k; Return E ↑ at precision m . Algorithm 5.1 Computing the canonical lift by lifting the étale p-torsion hence B = ∂Φ p /∂ x (J, F (J))/∂Φ p /∂ y (J, F (J)) is of valuation -1 by by Kronecker's formula (see Section 2.3 and Lemma 2.4). □

We look for a lift of the form J + p k e, and we want:

F (J + p k e) = A + B.e.p k = J Σ + e Σ.p k = A Σ + B Σp k e Σ.
Since B is of valuation -1, evaluating F (X) only make sense modulo p 2k-1 . Concretely this stems from the fact that given Ẽ at precision 2k, we can only compute Hp at precision 2k -1, so the corresponding isogeny at precision 2k -1. So we solve F (J + p k e) = A + (Bp)p k-1 e = J Σ + e Σp k mod p 2k-1 . By applying the Frobenius Σ, we get: where a = 0 modulo p. We then solve this equation using algorithm Section 2.3. This proves Theorem 1.1. The resulting algorithm is given in (Algorithm 5).

A Σ + B Σ .
Computing the canonical lift from lifting a point of p-torsion. Instead of lifting H p to compute the isogeny, we could also lift a point of p-torsion P directly and use an isogeny algorithm that takes a point of the kernel as input to compute the isogenous curve E ν . This second strategy gives the complexity stated by Theorem 1.2, the algorithms are summarized in algorithms 5.2 and 5.3.

To illustrate the flexibility of Section 3, rather than working with the j-invariant, we also illustrate a variant which works directly with the coefficients of Ẽ.

Suppose at precision k we have E : y 2 = x 3 + Ax + B and P (x p , y p ) ∈ E[p] -{O}.

Let (e, r) be the couple of error such that: x P -x P = e.p k and y P -y P = r.p k , then we can obtain a lift of the coefficients A and B by taking A and B + θ.p k . These three errors e, r and θ

  Input a, b ∈ Z q and the precision m. Output e such that e σ + ae + b ≡ 0 mod p N with a = 0 mod p If N = 1 Return e the unique root of e σ + b ≡ 0 mod p.x ← ArtinSchreier(a, b, N/2).Lift arbitrarilyx at precision p N . b ′ ← (x Σ + ax + b)/p N/2 . e ← ArtinSchreier(a, b ′ , N/2). Return x + p N/2 e. Algorithm 2.1 Artin-Schreier 2.4. Lift of the Weierstrass Equation.

  e Σ .p k = J + e.p k So dividing by (pB)p k-1 (recall that pB is inversible), we obtain an equation of the form: e Σ + a.e + b = 0 mod p k .

  Computing zeta functions of curves over finite fields, PhD thesis, Katholieke Universiteit Leuven (2003). [Vil18] Gilles Villard, On computing the resultant of generic bivariate polynomials, Proceedings of the 2018 acm international symposium on symbolic and algebraic computation, 2018, pp. 391-398. Cheikh Anta Diop University, Dakar, SENEGAL Équipe FAST, LIRIMA (Laboratoire International de Recherche en Informatique et Mathématiques Appliquées) Email address: abdoulaye.maiga@aims-senegal.org INRIA Bordeaux-Sud-Ouest, 200 avenue de la Vieille Tour, 33405 Talence Cedex FRANCE Email address: damien.robert@inria.fr URL: http://www.normalesup.org/~robert/ Institut de Mathématiques de Bordeaux, 351 cours de la liberation, 33405 Talence cedex FRANCE Équipe FAST, LIRIMA (Laboratoire International de Recherche en Informatique et Mathématiques Appliquées)

	[Ver03]	F. Vercauteren,

Acknowledgements. We thank Djiby SOW (FST-UCAD-Senegal) for his comments and suggestion on ealier version of this paper. We are supported by the FAST project and ANR CIAO. Experiments presented in this paper were carried out using PARI/GP [Dev19].

Ψ p (x P ) is only of valuation 2 (by Proposition 4.1), not 3 as needed to bootstrap the Newton method, see Remark 3.2. But since Ψ" p (x P ) = 0 modulo p (because Ψ p = H p p , the first iteration of the Newton method does still allow to go to precision 3, as remarked in [Rob21, Aside 6.2.2]. We can then apply Remark 3.2: we compute x P modulo p k such that Ψ p (x P ) = 0 modulo p k+1 . We can then lift y P by solving the square root via Newton's algorithm.

In summary:

Lemma 4.4. Let Ẽ/Z q be a lift of E that satisfy the conditions of Proposition 4.1, and let P be an étale point of p-torsion on E which lives in F q e . Then P can be lifted to a point of p-torsion P ∈ Ẽ[p](Z q e ) to precision m in time Õ(m log q log p) = Õ(mn).

4.4.

Lifting all the étale p-Torsion. An alternative is to lift directly the kernel of the Verschiebung H p . Suppose that we are given Hp at precision k, and we want to compute it at precision 2k. First we note that by employing the same strategy as in Section 4.3 but working over the algebra A p = F q [u]/H p (u) rather than over F q e , we can find a lift P = u + p.a 1 (u) + p 2 .a 2 (u) + • • • of the formal point of p-torsion P : x = u. Notably, P encodes simultaneously the lifts of all points of p-torsion: if P λ is the point of p-torsion with x-coordinate given by the root λ of H p , its lift Pλ is given by λ + p.u 1 (λ) + p 2 .u 2 (λ) + • • • .

Then Hp (x) is given by the resultant Res u (x -u, P ). But it is not clear if this resultant can be computed in quasi-linear time (the best generic algorithm in [Vil18] is not quasi-linear, but in our situation the roots of Hp are deformations of the roots of H p so there may be more efficient algorithms.) So rather than lifting the formal point P : x = u over A p , we simply lift Hp directly. We give two methods.

The first is to use Section 3 applied to the equation Ψ p mod Hp = 0 mod p 2k (where 2k is our target precision).

Indeed by Lemma 2.5, Ψ p can be evaluated modulo our candidate polynomial H * p via the recurrence formula in quasi-linear time (so we never need to compute it fully, only modulo a polynomial of degree p).

The Newton formula is as follow: take an arbitrary lift H * p , let a = Ψ p mod H * p and b = Ψ p mod ( H * p + p). Then the derivative of our Newton process is given by c = (b -a)/p k , and we solve the equation a + cp k Q = 0 mod ( H * p , p 2k ) (since the equation is valid at precision k, this equation does not depend on the choice of H * p ). The correct lift is then Hp = H * p + p k Q. The second one is to use the strategy of Section 2.6. Given Hp at precision k, take an arbitrary lift H * p to precision 2k + 1, then we have

We can use Lemma 2.5 to compute e in quasi-linear time.

In summary both methods give:

Proposition 4.5. Given an ordinary elliptic curve E/F q , and a lift (not necessarily canonical) Ẽ at precision m that satisfy the conditions of Proposition 4.1, the kernel H p of the Verschiebung can be lifted to precision m in quasi-linear time Õ(mp log q) = Õ(pmn).

Remark 4.6. As for Section 4.3, when Ẽ is given at precision m, Hp is only determined to precision m -1.

It is easy to see that we can extend the methods of this section to lift to Ẽ a subgroup G of degree d of E[ℓ], when p ∤ ℓ. (In this case there is no restriction on Ẽ since E[ℓ] is étale.)

Input Coefficients (a 4 , a 6 ), p-torsion P and integer m the precision. Output Coefficients (a ↑ 4 , a ↑ 6 ) at precision m. k = 1; Use Initialization Phase (Lemma 5.2) to compute J ↑ at precision 2; While (2 ⩽ k ⩽ (m + 1)/2) ; a. Use algorithm 5.3 to compute J ↑ at precision 2k -1; b. Compute the lift (a ↑ 4 , a ↑ 6 ) of coefficients at 2k -1 using J ↑ and the method Section 2.4 ;

Algorithm 5.2 Lifting the coefficients (a 4 , a 6 ) of an ordinary elliptic curve E.

constitute the unique root of system given by:

Then using the first equation we get: r = 1 2y P (x 3 P + A.x P + B -y 2 P )/p k + (3.e.x 2 P + A.e + θ) . And using the second equation, one can extract p.e from Ψ p x p + e.p k , A, B + θ.p k = 0 (since from Satoh's lemma we have Ψ ′ p (x p + e.p k ) has valuation 1). On other hand, let A v and B v be the coefficients of the p-isogenous curve given by the Vélu's formula from the input (A, B + θ.p k , P e,r ) at the precision 2k where P e,r = (x P + e.p k , y P + r.p k ). Since the curves E ν and E Σ must be isomorphic, the equality between the corresponding jinvariants gives the third equation:

After replacing in this equation r and p.e by theirs values, we obtain at precision 2k an equation of the form:

that can be solved using algorithm 2.3 to obtain θ, e and r. This variant is summarized in ?? 5.4. We note that the advantage of using all constraints at the same time, is that it automatically solve the initialisation problem (see Lemma 5.2).

Example 5.4. Let E be an ordinary elliptic curve over F 5 [T ]/M (T ) with M (T ) = T 10 +3T 6 +3T 5 +T 2 +2T +4 given by the equation y 2 = x 3 + a 4 x + a 6 where :

its j-inariant is

and P 0 = (x 0 , y 0 ) is its nonzero 5-torsion where:

The Teichmuller polynomial M of M at precision 13 is

Input Coefficients (a ↑ 4 , a ↑ 6 ) and p-torsion P at precision k. Output J ↑ = J + e.p k at precision 2k -1.

Choose nonzero r1 and r2 ; Set A6 r1 = a ↑ 6 + r1.p k and A6 r2 = a ↑ 6 + r2.p k ; Compute the Lift of P r1 and P r2 in both cases (a ↑ 4 , A6 r1 ) (a ↑ 4 , A6 r2 ) ; Compute the j-invariants Jv r1 and Jv r2 of the curves (a ↑ 4 , A6 r1 ) ν and (a ↑ 4 , A6 r2 ) ν at precision p 2k ; Set J r1 =Jinvariant(a ↑ 4 , A6 r1 ) and J r2 =Jinvariant(a ↑ 4 , A6 r2 ) Then R 1 = J r1 -J and R 2 = J r2 -J at precision p 2k ; Solve the system of equations


Joel E Cohen

Thierry E Huillet

TAYLOR'S LAW FOR SOME INFINITELY DIVISIBLE PROBABILITY DISTRIBUTIONS FROM POPULATION MODELS

Keywords: branching process, compound distribution, infinite divisibility, Markov process, Ornstein-Uhlenbeck process, Pólya-Aeppli distribution, power law, self-decomposability, Taylor

In a family of random variables, Taylor's law or Taylor's power law of fluctuation scaling is a variance function that gives the variance σ 2 > 0 of a random variable (rv) X with expectation µ > 0 as a power of µ: σ 2 = Aµ b for finite real A > 0, b that are the same for all rvs in the family. Equivalently, TL holds when log σ 2 = a + b log µ, a = log A, for all rvs in some set. Here we analyze the possible values of the TL exponent b in five families of infinitely divisible two-parameter distributions and show how the values of b depend on the parameters of these distributions. The five families are Tweedie-Bar-Lev-Enis, negative binomial, compound Poisson-geometric, compound geometric-Poisson (or Pólya-Aeppli), and gamma distributions. These families arise frequently in empirical data and population models, and they are limit laws of Markov processes that we exhibit in each case.

Introduction

For practical purposes like pest control, conservation, and evaluating the yields of alternative agricultural treatments, agronomists, ecologists, and agricultural statisticians invented a concept today called the "variance function." In a set of samples of some quantitative observations, the sample variance function describes how the sample variance varies from sample to sample as a function of the mean of each sample [START_REF] Bartlett | Some notes on insecticide tests in the laboratory and in the field[END_REF][START_REF] Beall | Methods of estimating the population of insects in a field[END_REF][START_REF] Beall | The transformation of data from entomological field experiments so that the analysis of variance becomes applicable[END_REF][START_REF] Bliss | Statistical problems in estimating populations of Japanese beetle larvae[END_REF].

Over several decades in the mid-twentieth century, multiple ecologists, apparently independently, observed that the log of the sample variance was approximated well by a linear function of the log of the sample mean, and that this model of the sample variance function was superior to several alternatives [START_REF] Bliss | Statistical problems in estimating populations of Japanese beetle larvae[END_REF][START_REF] Fracker | Measuring the local distribution of Ribes[END_REF][START_REF] Hayman | The transformation of counts of the cabbage aphid (Brevicoryne brassicae (L.))[END_REF][START_REF] Taylor | Aggregation, variance and the mean[END_REF]. The sample variance function that specifies the log sample variance as a linear function of the log sample mean has become known as Taylor's law (TL), after the last ecologist who discovered it [START_REF] Taylor | Aggregation, variance and the mean[END_REF]. The mathematically equivalent statement that the sample variance is approximately proportional to some power of the sample mean 1 has become known as Taylor's power law, or fluctuation scaling in physical applications. Thousands of empirical examples of TL have been found in many different fields of science and finance [START_REF] Eisler | Fluctuation scaling in complex systems: Taylor's law and beyond[END_REF][START_REF] Taylor | Taylor's Power Law: Order and Pattern in Nature[END_REF].

In theoretical studies of TL, a set of samples is modeled by a family (set) of random variables (rvs) {X(θ)} θ∈Θ indexed by some scalar-or vector-valued label θ ∈ Θ in an index set Θ with at least two elements. The sample mean of empirical studies is replaced by the population mean 0 < EX(θ) := µ(θ) ≤ ∞, and the sample variance of empirical studies is replaced by the population variance 0 < VarX(θ) := σ 2 (θ) ≤ ∞. Then, for finite real a, b which may or may not depend on θ, and for all θ ∈ Θ, TL asserts that log σ 2 (θ) = a + b log µ(θ) (log-linear form) or σ 2 (θ) = A[µ(θ)] b , A = exp(a) (power form) or [START_REF] Applebaum | Lévy Processes and Stochastic Calculus 2d ed[END_REF] b = log σ 2 (θ)a log µ(θ) , ∀θ ∈ Θ.

The connections between the TL parameters a, b and the distributions of X(θ) with moments µ(θ), σ 2 (θ) have long been of interest, and have sometimes been misunderstood. The best known example is the single-parameter family of Poisson(θ) distributions with expectation θ ∈ (0, ∞), which are widely used as a model of pure randomness in integer counts. Here σ 2 (θ) = µ(θ) = θ, so TL holds with a = 0, b = 1 for every θ ∈ (0, ∞). It does not follow that if TL holds with a = 0, b = 1, then any of the distributions is Poisson(θ) [START_REF] Cohen | Every variance function, including Taylor's power law of fluctuation scaling, can be produced by any location-scale family of distributions with positive mean and variance[END_REF]. If a = 0 or b = 1, then it is safe to conclude that at least one of the distributions is not a Poisson(θ) distribution.

Another well known single-parameter family of distributions is the negative exponential family Exp(θ), θ ∈ (0, ∞). Here σ 2 (θ) = [µ(θ)] 2 , so TL holds with a = 0, b = 2, regardless of θ. Again, it does not follow that if TL holds with a = 0, b = 2, then any of the distributions is Exp(θ).

Here we analyze the possible values of the TL exponent b in five families of infinitely divisible two-parameter distributions and show how the values of b depend on the parameters of these distributions. The shape of b (θ) appears to be specific to each family. A rv X is defined to be infinitely divisible if and only if, for every positive integer n, there exist n independent and identically distributed (iid) rvs

X n,1 , . . . , X n,n such that X n,1 + • • • + X n,n d = X, where d = means "has the same distribution as". If X is infinitely divisible and E(e itXn,1 ) = E(e itX ) 1/n is the characteristic function of (say) X n,1 , then for all t ∈ R exp -n 1 -E(e itXn,1 ) → E e itX as n → ∞,
showing that X is a weak limit of a compound Poisson sequence. More general iid sequences (of size k n → ∞) converging weakly to X can be found in Theorem 5.2 of [START_REF] Steutel | Infinite Divisibility of Probability Distributions on the Real Line[END_REF].

If X is infinitely divisible and satisfies Taylor's power law with exponent b, then each of its constitutive summands X n,m also satisfies Taylor's power law with the same exponent b. To see this, suppose every rv in a family

{X(θ)} θ∈Θ , Θ = ∅, satisfies (2) X(θ) d = X n,1 (θ) + • • • + X n,n (θ), n ≥ 1,
where X n,1 (θ), . . . , X n,n (θ) are iid as X n,1 (θ). Then X(θ) satisfies TL [START_REF] Applebaum | Lévy Processes and Stochastic Calculus 2d ed[END_REF] with exponent b if and only if X n,1 (θ) satisfies TL (1) with the same exponent b. Indeed, taking the expectation of both sides of (2) gives µ X (θ) = nµ Xn,1 (θ) and σ 2

X (θ) = nσ 2 Xn,1 (θ). If σ 2 X (θ) = A[µ X (θ)] b , then nσ 2 Xn,1 (θ) = An b (µ Xn,1 (θ)) b or σ 2 Xn,1 (θ) = An b-1 (µ Xn,1 (θ)) b , which is TL with the same exponent b. The coefficient A of X(θ) leads to the coefficient A n = An b-1 of X n,1 (θ).
The converse is obvious.

The five families to be analyzed here are Tweedie-Bar-Lev-Enis, negative binomial, compound Poisson-geometric, compound geometric-Poisson (or Pólya-Aeppli), and gamma distributions. (Kendal [27] gives another interesting family.) These families arise frequently in empirical data and population models. They all have the special property that it is possible to express the dependence of the mean and variance on their parameters in such a way that a = 0 or σ

2 (θ) = [µ(θ)] b , ∀θ ∈ Θ or (3) b = log σ 2 (θ) log µ(θ) , ∀θ ∈ Θ.
Non-zero a can arise from the rescaling described in section 7.

In addition to being infinitely divisible, some of these distributions are also selfdecomposable (SD). A distribution or rv X is defined to be SD if, for every c ∈ (0, 1), there exists an independent rv X c such that ( 4)

X(θ) d = c(θ)X(θ) + X c(θ) (θ).
If X takes values in N 0 := {0, 1, 2, . . .}, cX is to be interpreted as the c-Bernoulli thinning of X, which is defined as the sum of X iid Bernoulli rvs with success probability c.

If every rv in a family {X(θ)} θ∈Θ , Θ = ∅, satisfies TL (1) and ( 4), then the family {X c (θ) | θ ∈ Θ} also satisfies TL with the same value of b as the family {X(θ) | θ ∈ Θ}; and conversely. To see this, take the expectation of both sides of (4). Then µ(θ) = c(θ)µ(θ) + EX c(θ) (θ), hence µ(θ) = EX c(θ) (θ)/(1c(θ)). Taking the variance of both sides of (4) gives σ 2 (θ) = c(θ) 2 • σ 2 (θ) + VarX c(θ) (θ), hence, using TL for the second equality,

VarX c(θ) (θ) = σ 2 (θ) • (1 -c(θ) 2 ) = A[µ(θ)] b • (1 -c(θ) 2 ) = A[EX c(θ) (θ)/(1 -c(θ))] b • (1 -c(θ) 2 ) = {A(1 -c(θ) 2 )(1 -c(θ)) -b }[EX c(θ) (θ)] b .
This is TL with the same b, as claimed. The coefficient A of X(θ) leads to the coefficient

A(1 -c(θ) 2 )(1 -c(θ)) -b of X c(θ) (θ).
The proof of the converse is routine, given the above equalities.

Every SD rv is unimodal. If X is SD and not discrete-valued, then it has a density. If X is SD and continuous, it is a weak limit of a continuous-time Lévydriven Ornstein-Uhlenbeck (OU) process. If X is SD and N 0 -valued, it is a weak limit of a pure-death branching process with immigration. See the Appendix.

Thus we consider three kinds of infinitely divisible distributions with two parameters. First, "bare" infinitely divisible distributions arise as specific limit laws, as we will describe. Second, infinitely divisible distributions that are compound Poisson geometric are the limit laws of immigration processes with total disasters. Third, infinitely divisible distributions that are SD are the limit laws of a pure-death branching process with immigration (if discrete) or of a Lévy-driven Ornstein-Uhlenbeck process (if continuous). In the latter two cases, there is a balance between events of birth (or immigration) and events of death, resulting in an equilibrium distribution. For some of the families under study, depending on the parameter range, the rv can switch from "bare" infinitely divisible to SD, perhaps suggesting some kind of phase transition. There may be other examples of infinitely divisible distributions that fall into these three categories, and it is of interest to know whether they would satisfy Taylor's law.

Tweedie-Bar-Lev-Enis distributions

The definition of Tweedie-Bar-Lev-Enis distributions is elaborate [START_REF] Jørgensen | Exponential dispersion models[END_REF][START_REF] Jørgensen | The Theory of Dispersion Models[END_REF][START_REF] Bar-Lev | Reproducibility and natural exponential families with power variance functions[END_REF][START_REF] Bar-Lev | Characterizations of natural exponential families with power variance functions by zero regression properties[END_REF] and will not be attempted here. An accessible expository account, with historical background, is [2]. For brevity, we will refer to these distributions and associated random variables as the TweBLE family.

With parameter α ∈ [-∞, 0) ∪ {0} ∪ (0, 1) ∪ (1, 2) ∪ {2} and ( 5)

k (θ) = 1 -α α θ 1 -α α
for the values of parameter θ for which θ 1-α α is well-defined, the bilateral probability Laplace-Stieltjes transform (PLSt) of a TweBLE random variable (rv) [START_REF] Jørgensen | Exponential dispersion models[END_REF][START_REF] Jørgensen | The Theory of Dispersion Models[END_REF] is

Φ(λ) = e -[k(θ+λ)-k(θ)] .
The expression of k(θ) is extended to cover the boundary cases α = -∞ and α = 0, respectively, by 1e -θ and log θ. TweBLE rvs include Poisson (α = -∞), compound Poisson-gamma (α ∈ (-∞, 0)) (which is the sum of a Poisson-distributed number of iid gamma rvs all independent of the Poisson), negative exponential (α = 0), tempered one-sided stable (α ∈ (0, 1)), tempered two-sided extreme stable (α ∈ (1, 2)) and Gaussian (α = 2) rvs. Considering ( 5), the log-Laplace transform (LLt) is

L(λ) := -log Φ(λ) = 1 -α α θ + λ 1 -α α - θ 1 -α α .
Hence the mean and variance of a TweBLE rv are

µ = L ′ (0) = θ 1 -α α-1 , σ 2 = -L ′′ (0) = θ 1 -α α-2 .
All TweBLE distributions satisfy σ 2 = µ b , which is TL with a = 0. Taylor's law holds with exponent A TweBLE rv has a negative TL exponent b < 0 only when its support is the whole real line, unlike some of the following examples, where b < 0 can occur when the support is the nonnegative half line or the nonnegative integers.

(6) b = log σ 2 log µ = (2 -α) / (1 -α) .
No TweBLE distribution has b ∈ (0, 1) because there is no non-degenerate Twe-BLE PLSt with θ > 2. Hence a population model that obeys TL with b ∈ (0, 1) cannot have the distribution of, and cannot be explained by using, a TweBLE rv. Examples of TL with b ∈ (0, 1) include a multiplicative population process in a Markovian environment [14, p. 33] and an infinity of examples with arbitrary distributions with finite means and finite variances [START_REF] Cohen | Every variance function, including Taylor's power law of fluctuation scaling, can be produced by any location-scale family of distributions with positive mean and variance[END_REF]Example 3]. Therefore TweBLE distributions are not sufficient for modeling and understanding all instances of TL.

When α < 1, TweBLE rvs are infinitely divisible because Φ(λ) = exp -L(λ) where L ′ (λ) is completely monotone, obeying (-1) n L (n+1) (λ) ≥ 0 for all λ > 0. When α ∈ [1, 2], TweBLE rvs are tempered α-stable. In all cases, the Lévy jump measure can be found in [START_REF] Bar-Lev | Natural exponential families and selfdecomposability[END_REF]. Moreover, [5, Example 3] shows that they are SD as well except when α ∈ (-∞, 0). For TweBLE rvs that are SD, the PLSt is

(7) Φ (λ) = exp λ 0 log Φ 0 λ ′ λ ′ dλ ′
for some PLSt Φ 0 (λ) of an infinitely divisible rv [30, Theorem 2.9, Eq. 2.12]. TweBLE distributions as defined in this Section belongs to the class of exponential families [2]. All such TweBLE distributions satisfy σ 2 (θ) = [µ(θ)] b , which is TL with A = 1, a = 0, and b given by [START_REF] Barndorff-Nielsen | The multivariate supOU stochastic volatility model[END_REF]. Introducing the scaling ( 14) below, we switch to the class of exponential dispersion models, as defined by Jorgensen [START_REF] Jørgensen | The Theory of Dispersion Models[END_REF]. For this new class, σ 2 (θ) = A[µ(θ)] b with A = e a and a = 0, independent of b.

Negative binomial distributions

We shall show that the next two-parameter probability distribution families can be formulated in such a way that they also obey TL with the TweBLE family's special property that a = 0, and we shall express b in terms of their parameters. As for TweBLE models, the scaling [START_REF] Cohen | Stochastic population dynamics in a Markovian environment implies Taylor's power law of fluctuation scaling[END_REF] will yield the version of TL with a = 0, independent of b.

With the parameters α > 0, p ∈ (0, 1), q := 1p, a negative binomial rv

X d = N B(α, p) has probability mass function P [X = k] = k + α -1 k q α p k , k = 0, 1, 2, . . . .

It has PLSt and LLt

Φ (λ) = q 1 -pe -λ α , L (λ) = -log Φ (λ) = -α log q + α log 1 -pe -λ , L ′ (λ) = αpe -λ 1 -pe -λ ,
and mean and variance

µ = L ′ (0) = αp q , σ 2 = -L ′′ (0) = αp q 2 = µ + 1 α µ 2 . ( 8
)
Overdispersion holds, in the sense that, unlike the Poisson distribution, where the variance equals the mean, the negative binomial variance strictly exceeds the mean. The smaller α is, the larger the variance compared to the mean. In ecology, overdispersion in population counts is sometimes interpreted as a resulting from heterogeneous field conditions (environmental variation making some places more favorable for individuals of the population than other places) or from aggregation or clustering (individuals behaving in a way that brings them near other individuals) [START_REF] Bliss | Statistical problems in estimating populations of Japanese beetle larvae[END_REF][START_REF] Taylor | Aggregation, variance and the mean[END_REF].

Taylor's law holds with exponent

(9) b = log σ 2 log µ = 1 - log q log αp q . Then b < 1 ⇔ αp q < 1 ⇔ p < p c := 1/ (1 + α) , b < 0 ⇔ log q log αp q > 1 ⇔ q 2 + αq -α < 0 ⇔ q < q 0 := -α + √ α 2 + 4α 2 ⇔ p > p 0 := 2 + α - √ α 2 + 4α 2 > 0.
We always have [START_REF] Applebaum | Lévy Processes and Stochastic Calculus 2d ed[END_REF]2). In this case, the negative binomial distribution reduces to the geometric distribution, and p c = 1/2 and p 0 = 3 -√ 5 /2 = φ 2 , where φ = ( √ 5 -1)/2 is a golden ratio. The formula (9) identifies a manifold of (α, p) points on which the TL exponent b is constant. To maintain a constant b, (9) requires both α and p to vary. This result is surprising in light of the above negative binomial variance function (8), σ 2 = µ + α -1 µ 2 . For suppose α is constant. Then as µ → 0, µ 2 ≪ µ and asymptotically σ 2 is proportional to µ, which is asymptotically TL with a = 0, b = 1. However, as µ → ∞, µ 2 ≫ µ and asymptotically σ 2 is proportional to α -1 µ 2 , which is asymptotically TL with e a = α -1 , b = 2. Moreover, log σ 2 is a strictly convex increasing function of log µ [16, Supp.Mat., lines 29-42], instead of a linear function as in TL. If the expectation µ varies over only a small range in the data or the family of rvs, then log σ 2 can be closely approximated by (may be statistically indistinguishable from, in the case of empirical data) a locally linear function of log µ with a local slope that increases smoothly from 1 (for very small µ) to 2 (for very large µ). A wide range of means µ would be required to display the curvature in the relation of log σ 2 to log µ when the underlying distributions are negative binomial with constant α.

p 0 < p c = 1 1 + α .
if α > 1 (respectively α ≤ 1). Hence b ∈ (0, 1) if 0 < p < p 0 , b ≤ 0 if p 0 ≤ p < p c , b > b min if p > p c . The range (1, b min ) for b is excluded. If α = 1, then b min = 2 and b cannot fall in the interval
These observations illustrate that the form of the variance function of a twoparameter family of probability distributions depends strongly on how the parameters and their relation are constrained. We illustrate this point again in analyzing the gamma distribution below.

The negative binomial distribution is discrete and infinitely divisible (compound Poisson) because its probability generating function (pgf)

φ (z) := E z X is (10) φ (z) = q 1 -pz α = exp{-[-α log q] (1 -c (z))},
where c(z) = log(1pz)/ log q is the pgf of a logarithmic series distribution for its clusters' sizes.

The negative binomial distribution is also SD because its pgf may be written as (compare with [START_REF] Bartlett | Some notes on insecticide tests in the laboratory and in the field[END_REF] in the continuum)

φ (z) = exp -r 1 z 1 -h (z ′ ) 1 -z ′ dz ′
for some rate r > 0 and pgf h(z) such that h(0) = 0 [30, Theorem 4.11, Eq. 4.13]. Indeed, r = αp > 0 and h(z) = qz/(1pz) is a geometric pgf. See the Appendix. When α ∈ (0, 1), then π := 1q α ∈ (0, 1) and the pgf φ(z) in [START_REF] Bertoin | Subordinators, Lévy processes with no negative jumps, and branching processes[END_REF] equals

φ(z) = 1 -π 1 -πϕ(z) ,
where

ϕ(z) := 1 -(1 -pz) α
1q α satisfies ϕ(0) = 0. This ϕ(z) is a pgf because it is an absolutely monotone function: it obeys ϕ (n) (z) ≥ 0 for all 0 < z < 1 and has positive coefficients in its power series expansion. Indeed, the coefficient of

z n in ϕ(z) is [z n ]ϕ(z) = p n 1 -q α (-1) n-1 (α) n n! = p n 1 -q α α[α] n-1 n! , n ≥ 1 where (α) n := Γ (α + 1) /Γ (α + 1 -n) = α (α -1) ... (α -n + 1) is the falling fac- torial and [α] n := α (α + 1) ... (α + n -1) , n ≥ 1, is the rising factorial of α with [α] 0 := 1, α = 1 -α.
Because the negative binomial distribution is SD for all α > 0, it follows from our Appendix that the negative binomial distribution is the limit law of a pure-death branching process with immigration. But when α ∈ (0, 1), the negative binomial distribution is also compound-geometric, so it is also the limit law of a Markov chain with total disasters, as defined in the next section.

Compound Poisson-geometric distributions

With parameters α > 0, p ∈ (0, 1), q := 1p, a compound Poisson-geometric rv is the sum of N iid Poisson(α)-distributed rvs, where N + 1 is geometrically distributed with parameter p. The PLSt and LLt of a compound Poisson-geometric rv are

Φ(λ) = q 1 -pe -α(1-e -λ ) , L(λ) := -log Φ(λ) = -log q + log 1 -pe -α(1-e -λ ) .
The mean and variance are

µ = L ′ (0) = αp q , σ 2 = -L ′′ (0) = αp(α + q) q 2 = α + q q µ > µ.
Overdispersion holds. Taylor's law holds with exponent

(11) b = log σ 2 log µ = 1 + log 1 + α q log αp q .
The TL exponent b diverges when q = q c := α/(1 + α) or α = q/p. It vanishes when q = q 0 , defined as the solution (0, 1) of

f (q) := (1 + α) q 2 -qα (1 -α) -α 2 = 0.
Further, because f (q c ) = -α 2 /(1 + α) < 0, we have q 0 > q c . So the graph of b has two branches for a given fixed α:

-one is increasing and convex from b = 2 to b = ∞ as q varies from 0 + to q - c . -one is increasing from b = -∞ to b = 1 as q increases from q + c to 1, passing through 0 as q = q 0 > q c . Hence

b ∈ (2, ∞) if 0 < q < q c , b < 1 if q > q c .
The TL exponent b is excluded from [START_REF] Applebaum | Lévy Processes and Stochastic Calculus 2d ed[END_REF]2).

If α > 0 is small enough that e -α < (1α)/α (e.g., if α < 0.659) and if p ≥ p * := α/(1e -α α) ∈ (0, 1), then the compound Poisson-geometric rv is SD.

To prove this claim, we observe that, with ϕ(z) = e -α(1-z)

φ(z) = q 1 -pϕ(z) = e -(-log q)[1-c(ϕ(z))]
where c(z) is the pgf of a logarithmic series distribution. So φ(z) is compound Poisson with clusters' size pgf c(ϕ(z)).

The probability mass function P [X = k] = [z k ]φ(z) of a compound Poissongeometric rv X is in principle explicitly given by the Faa di Bruno formula for compositions of two pgfs [18, p. 146

]. It involves the ordinary Bell polynomials in the variables b

k = [z k ]ϕ(z).
Next, we try to write the pgf φ(z) of X in the form

φ(z) = q 1 -pϕ(z) = exp{-r 1 z 1 -h(z ′ ) 1 -z ′ dz ′ }
for some rate r > 0 and pgf h(z) obeying h(0) = 0. Forcing h(0) = 0 yields

r = pϕϕ ′ (0) 1 -pϕ(0) = pαe -α 1 -pe -α . Hence h(z) = 1 - 1 r (1 -z) pϕ ′ (z) 1 -pϕ(z) = 1 r r(1 -pϕ(z)) -p(1 -z)ϕ ′ (z) 1 -pϕ(z) .
Denoting the numerator by N (z), a sufficient condition for h to be a pgf is that

[z k ]N (z) ≥ 0 for all k ≥ 1,
where [z k ]N (z) denotes the coefficient of z k in the power series expansion of N (z). But, with b k = e -α α k /k!, we have

N (z) = p k≥1 z k [(k -r)b k -(k + 1)b k+1 ]. So h is a pgf if b k+1 b k = α k + 1 ≤ p(k -b 1 ) k + 1 for any k ≥ 1 or equivalently α ≤ p(1 -b 1 ) = p(1 -e -α α) or p ≥ p * := α/(1 -e -α α).
This completes the proof that, under the above condition, the compound Poissongeometric rv is SD.

A compound Poisson-geometric rv is the limiting rv of a Markov chain with total disasters. To see this, let (β n ) n≥1 be an iid sequence taking values in N 0 := {0, 1, 2, ...}. Let ϕ(z) := E(z β ) be the pgf of the βs. The Markov chain X n evolves by: X n+1 = X n + β n+1 with probability p, X n+1 = 0 with probability q := 1p.

This simple population growth model alternates periods of births of amplitude β n+1 with one or more total disasters where population size X n is instantaneously reset to 0. We assume without loss of generality that X 0 = 0 and that q does not depend on X n = x. We define p * X to equal the product of X times an independent Bernoulli rv B(p) with success parameter p. (This product is not to be confused with Bernoulli thinning.) Because X := X ∞ solves the distributional equation

X d = B(p)(X ′ + β) =: p * (X ′ + β)
, where (B(p), X ′ d = X, β) are mutually independent and B(p) is a Bernoulli rv with success parameter p, the above Markov chain {X n | n = 0, 1, 2, . . .} with total disasters is clearly ergodic. The limiting rv X := X ∞ exists and it has the compound geometric pgf φ(z) = q/(1pϕ(z))

(shifted to the left by one unit). In our example, β is Poisson-distributed with pgf ϕ(z) = e -α(1-z) .

Compound geometric-Poisson distributions

With parameters α > 0, p ∈ (0, 1), q := 1p, a compound geometric-Poisson (or Pólya-Aeppli) rv is the sum of a Poisson(α)-distributed number of geometric rvs with parameter p. The PLSt and LLt are

Φ(λ) = e -α 1-qe -λ 1-pe -λ = e -α 1-e -λ 1-pe -λ , L(λ) := -log Φ(λ) = α
1e -λ 1pe -λ . Hence the mean and variance are

µ = L ′ (0) = α q > 0, σ 2 = -L ′′ (0) = α(1 + p) q 2 = 1 + p q µ > µ.
Its probability mass function is given by

P [X = k] = e -α k l=1 k -1 l -1 α l l! p k-l q l , if k ≥ 1 = e -α if k = 0
Overdispersion holds. Taylor's law holds with exponent

(12) b = log σ 2 log µ = 1 + log 1+p q log α q .
If α ≥ 2, the graph of b has only one decreasing branch as q varies from 0 to 1, with b → 2 -as q → 0 + and b = 1 at q = 1. Hence b ∈ [START_REF] Applebaum | Lévy Processes and Stochastic Calculus 2d ed[END_REF]2].

If 1 < α < 2, the graph of b increases from b = 2 to a maximum less than 3 and then drops to 1 as q increases from 0 to 1.

If α = 1, then b = log(1 + p) = log(2q), so b falls from 2 to 0 as q increases from 0 to 1.

If 0 < α < 1, b diverges when q = q c := α, so b has two branches: -one is increasing and convex from b = 2 + to b = +∞ as q increases from 0 + to q - c .

-one is increasing and concave from b = -∞ to b = 1 as q increases from q + c to 1 -. This branch passes through b = 0 when q = q 0 := (

√ α 2 + 8α -α)/2 > q c . Hence the full range of b is covered: b ∈ (2, +∞) if 0 < q < q c , b ∈ (-∞, 2) if q > q c . The pgf φ(z) of a geometric-Poisson or Pólya-Aeppli rv is φ(z) = exp -α 1 -z 1 -pz = exp -r 1 z 1 -h (z ′ ) 1 -z ′ dz ′
for some rate r > 0 and pgf h(z) obeying h(0) = 0, only if p > 1/2. The necessary condition p > 1/2 arises because h(0) = 0 yields r = αq and

h(z) = z 1 -2p + p 2 z (1 -pz) 2 ,
which is absolutely monotone only if p > 1/2, and then h(z) satisfies [z n ]h(z) ≥ 0 for all n ≥ 1.

Gamma distributions

A gamma rv X with shape parameter α > 0 and scale parameter β > 0 has probability law

P (X ∈ dx) = x α-1 e -x/β Γ(α)β α dx, x ∈ (0, ∞). The PLSt and LLt are Φ(λ) = 1 1 + λβ α , L(λ) = -log Φ(λ) = α log(1 + λβ).
The mean and variance are

µ = L ′ (0) = αβ > 0, σ 2 = -L ′′ (0) = αβ 2 = 1 α µ 2 = βµ.
Overdispersion holds if and only if β > 1.

As a reviewer pointed out, a gamma rv can satisfy TL in multiple ways. For example, if we fix b = 2, then σ 2 /µ 2 = 1/α = A and a =log α. If we fix b = 1, then σ 2 /µ = β = A and a = log β. If we fix b = 3/2, then σ 2 /µ 3/2 = (β/α) 1/2 and a = (1/2) log(β/α). Here we fix a = 0, A = 1 so that TL holds with exponent

→ 1 + if α → 0 + , b → +∞ if α → α - c and b → -∞ if α → α + c , b → 1 -if α → ∞. If β > 1,
the graph of b is a symmetric image of the previous one with respect to the horizontal line b = 0, so with two decreasing branches.

If α = 1, giving an exponential distribution, then b = 2 is a fixed point. The gamma rv is SD because [30, p. 4, Theorem 2.9]

Φ(λ) = exp λ 0 log Φ 0 λ ′ λ ′ dλ ′
for some PLSt Φ 0 (λ) of an infinitely divisible rv, which is here found to be

Φ 0 (λ) = e αβλ/(1+λβ) = e -α(1-1/(1+λβ)) .
This Φ 0 (λ) is the PLSt of a compound Poisson(α) rv with Exp(β) distribution for the random size ∆ of its clusters. Consider the Ornstein-Uhlenbeck process

dX t = -X t dt + dL t , X 0 = 0,
driven by the Lévy-process L t (here a rate-α compound-Poisson exponential process) for which

Ee -λLt = Φ 0 (λ) t , t ≥ 0. With φ ∆ (λ) = Ee -λ∆ = (1 + λβ) -1 , Φ t (λ) := E(e -λXt | X 0 = 0) = exp -α t 0 (1 -φ ∆ (λe -s ))ds = 1 + λβ 1 + λβe -t -α → (1 + λβ) -α = E(e -λX∞ | X 0 = 0) = Φ(λ) as t → ∞.
The function Φ t (λ) is the PLSt of some rv X t with an atom at 0 with probability mass e -αt .

Scaling

In each example above, scaling the log-Laplace transform according to ( 14)

L(λ) → L 1 (λ) = 1 σ 2 1 L(σ 2 1 λ)
defines the law of a new rv. Indeed, the scaled function L 1 (λ) defined from the LLt L(λ) is itself the LLt of some rv if and only L(λ) is the LLt of some infinitelydivisible rv (which is true in our examples here). Under this scaling ( 14), the mean µ = L ′ 1 (0) remains invariant while the variance is rescaled from -L ′′ (0) = σ 2 to -L ′′ 1 (0) = σ 2 1 σ 2 . (For instance, a Γ(α, β) rv transforms to a Γ(α/A, βA) rv with A = σ 2 1 .) Then TL transforms according to

σ 2 = µ b → σ 2 = σ 2
1 µ b where σ 2 1 (previously named A) is the variance of the new scaled rv when its mean is 1. The log-linear version of TL now includes a non-zero constant term a = log σ 2 1 : log σ 2 = b log µ → log σ 2 = a + b log µ.

If X has support N 0 := {0, 1, 2, . . .}, then under the scaling [START_REF] Cohen | Stochastic population dynamics in a Markovian environment implies Taylor's power law of fluctuation scaling[END_REF] its support becomes e a N 0 := {0, e a , 2e a , . . .}. Scaling in the discrete case changes the original 'counting' support N 0 to e a N 0 , emphasizing that A = e a plays the role of a scale parameter of the underlying distribution. In this discrete dispersion model, A = e a could be interpreted as the individual mass of each individual constituting the population [START_REF] Jørgensen | Discrete dispersion models and their Tweedie asymptotics[END_REF].

If the support of X is the whole real line or the nonnegative half line, the above rescaling by e a does not change the support.

Conclusion and open questions

Empirical observations that sample means and sample variances are approximately consistent with TL cannot specify the underlying distribution, although they can reject some possibilities, as when the best estimates of a = 0 or b = 1 reject a Poisson distribution. As [START_REF] Cohen | Every variance function, including Taylor's power law of fluctuation scaling, can be produced by any location-scale family of distributions with positive mean and variance[END_REF] observed, for every µ > 0 and any rv Y with mean 0 and variance 1, the rv X = µ + µ b/2 Y obeys TL with exponent b and intercept a = 0. X can then be scaled following [START_REF] Cohen | Stochastic population dynamics in a Markovian environment implies Taylor's power law of fluctuation scaling[END_REF] to include the affine term a, resulting in a model obeying TL with any a and b.

For each of the infinitely divisible two-parameter models in this note, we showed that the TL exponent b depends on the parameters in very specific ways. We showed that the admissible range of b depends on the distribution, and in some cases on the parameters of the distribution, in ways that could help identify the underlying law from data approximately obeying TL, or at least could exclude some possible underlying laws. By contrast, for samples of increasing size from nonnegative distributions with regularly varying upper tails and tail index α ∈ (0, 1), TL holds asymptotically, but with different formulas for the exponent b, when the sample variance is replaced by the sample lower semivariance (b = 2, regardless of α ∈ (0, 1)), the sample lower local semivariance (b = 2, regardless of α ∈ (0, 1)), or the sample local upper semivariance (b = (2α 2 )/(1α)) [START_REF] Brown | Taylor's law of fluctuation scaling for semivariances and higher moments of heavytailed data[END_REF].

Can four appearances of ( 6), one here and three earlier, be coincidences? Or is some underlying process or mechanism common to all these different appearances?

Appendix

Here we briefly sketch that if X is SD, it is a weak limit of a pure-death branching process with immigration (if N 0 -valued), or a weak limit of a continuous-time Lévydriven Ornstein-Uhlenbeck process (if continuous). The analysis concerns the rvs X obeying TL with a = 0. From the scaling transform ( 14) introducing a = log σ 2 1 , the modifications for a = 0 could be readily obtained. In both discrete and continuous cases, a population that is randomly annihilated is randomly regenerated by the recurrent arrivals of random quantities of immigrants, yielding a stationary or invariant distribution of population size. 9.1. Discrete self-decomposable rvs and pure-death branching processes with immigration in continuous time. van Harn et al. [START_REF] Van Harn | Self-decomposable discrete distributions and branching processes[END_REF] construct a regenerative process in continuous time that produces discrete SD distributions in the long run. Consider a continuous-time homogeneous compound Poisson process P r (t), t ≥ 0, P r (0) = 0, having rate r > 0, with pgf

(15) E Pr (0)=0 (z Pr(t) ) = exp{-rt(1 -h(z))},
where h(z), with h(0) = 0, is the pgf of the sizes of the clones or immigrant clusters arriving at the jump times of P r (t). Let (16)

ϕ t (z) = 1 -e -t (1 -z)
be the pgf of a pure-death branching process started with one particle at t = 0. (More general subcritical branching processes could be considered.) This expression of ϕ t (z) is easily seen to solve . ϕ t (z) = f (ϕ t (z)) = 1ϕ t (z), ϕ 0 (z) = z, as is usual for a pure-death continuous-time Bellman-Harris branching process [START_REF] Harris | The theory of branching processes[END_REF] with affine branching mechanism f (z) = r d (1z) with fixed death rate r d = 1. The distribution function of the lifetime of the initial particle is 1e -t . Let X t with X 0 = 0 be a random process counting the current size of some population for which a random number of individuals (determined by h(z)) immigrate at the jump times of P r (t). Each newly arrived individual is independently and immediately subject to the pure death process above. We have [START_REF] Cohen | Heavy-tailed distributions, correlations, kurtosis, and Taylor's law of fluctuation scaling[END_REF] φ t (z) := E(z Xt ) = exp -r t 0

[1h(ϕ s (z))]ds , φ 0 (z) = 1, with φ t (0) = P (X t = 0) = exp{-r t 0 (1h(1e -s ))ds}, the probability that the population is extinct at t. As t → ∞,

φ t (z) → φ ∞ (z) = exp -r ∞ 0 [1 -h(1 -e -s (1 -z))]ds = exp -r 1 z 1 -h (u) 1 -u du . (18)
So X := X ∞ , as the limiting population size of this pure-death process with immigration, is a SD rv [START_REF] Van Harn | Self-decomposable discrete distributions and branching processes[END_REF]. Define the rv X c implicitly by requiring that X d = cX ′ +X c , where X ′ is an iid copy of X and 0 < c < 1. Then

φ Xc (z) = φ ∞ (z) φ ∞ (1 -c(1 -z)) = exp -r 1-c(1-z) z 1 -h(u)
1u du is a pgf. In such models typically, a decaying subcritical branching population is regenerated by a random number of incoming immigrants at random Poissonian times. 9.2. Continuous self-decomposable rvs and Ornstein-Uhlenbeck process in continuous-time. When X is continuous and SD, X is the limiting distribution of population size as t → ∞ of some Ornstein-Uhlenbeck process X t : dX t = -X t dt + dL t , X 0 = 0, driven by the Lévy process L t for which

Ee -λLt = Φ 0 (λ) t , t ≥ 0, where Φ 0 (λ) is the PLSt of an infinitely divisible rv appearing in the representation (7) of Φ(λ) = Ee -λX . See [START_REF] Jurek | Remarks on the selfdecomposability and new examples[END_REF]. We now show that for a TweBLE rv with α ∈ (-∞, 0), there is no L 0 (λ) = log Φ 0 (λ) such that L ′ 0 (λ) is completely monotone on (0, ∞). This result means that the TweBLE rv for α ∈ (-∞, 0) is not SD, just infinitely divisible. Indeed, with L(λ) =log Φ(λ), L 0 (λ) = λL ′ (λ) = (1α) 1-α λ(θ + λ) α-1 , L ′ 0 (λ) = (1α) 1-α (θ + λ) α-2 [θ + λα], with L ′ 0 (λ) > 0 only if λ > λ c = -θ/α > 0, so not in the full range λ ∈ (0, ∞). So Φ 0 (λ) is not completely monotone on (0, ∞), and is therefore not an infinitely divisible PLSt. This Poisson-gamma regime for which the limiting distribution is a Poisson sum P of iid gamma-distributed clusters of size ∆ was studied by [19, p. 17, section 3.3.2], who underline what they call its "impact inhomogeneity": E(∆) = C(α) • E(P ) -1/α for some constant C(α) > 0.

By contrast, if α ∈ (0, 1), then L 0 (λ) may be written as

L 0 (λ) = ∞ 0
(1e -λx )π(x)dx where π(x)dx = 1 Γ(1α)

x -(α+1) (α + θx)e -θx dx is a tempered Lévy measure integrating 1 ∧ x. The driving process L t of X t , with PLSt Ee -λLt = e -tL0(λ) , is a subordinator and X = X ∞ is a SD TweBLE rv obtained as the limiting distribution of the corresponding Ornstein-Uhlenbeck process. One could extend this construction to other SD subordinated Lévy families, such as those in [START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF][START_REF] Applebaum | Lévy Processes and Stochastic Calculus 2d ed[END_REF][START_REF] Bertoin | Subordinators, Lévy processes with no negative jumps, and branching processes[END_REF][START_REF] Barndorff-Nielsen | The multivariate supOU stochastic volatility model[END_REF][START_REF] Schoutens | Lévy Processes in Finance: Pricing Financial Derivatives[END_REF].

  The extensions of TweBLE rvs to α = -∞ and α = 0 must have b = 1 and b = 2, respectively, to be consistent, and these values of b follow directly from the moments of Poisson and exponential distributions, respectively. Although a TweBLE rv has two parameters α, θ, the TL exponent b depends only on α, unlike the two-parameter distributions in the following sections.The graph of b versus α when α ∈ [0, 2] diverges at α c := 1 and has two hyperbolic branches:-one increasing from α = -∞ to α = 1 -, with b ranging from 1 + (Poisson) to b = +∞, through b = 2 (gamma) when α = 0.-one increasing from α = 1 + to α = 2, with b ranging from b = -∞ to b = 0 (Gaussian when b = 0, α = 2).

  So for any fixed α, the graph of b against p has two branches: -one is concave and decreasing from b = 1 to b = -∞ as p varies from 0 + to p - c as a diverging point according as b → -∞ as p → p - c , passing through 0 as p = p 0 . -one is varying from b = +∞ as p → p + c to b = 2 as p → 1 -possibly passing through a minimum b min > 1 with b min < 2 (= 2) if and only

  If β = 1, then b = 1 constant (equidispersion, or variance equal to the mean). If β = 1, the graph of b versus α shows a singularity at α c := 1/β, with the full range of b covered. If β < 1, the graph of b versus α has two increasing branches. We have b

  The expression b = (2α)/(1α) for the TL exponent in[START_REF] Barndorff-Nielsen | The multivariate supOU stochastic volatility model[END_REF] arose in at least three different prior examples of TL. First, TL holds asymptotically (as sample size increases toward infinity) with exponent b = (2α)/(1α) for the sample variance and sample mean of samples from nonnegative stable distributions with tail index α ∈ (0, 1) [12, page 663, Proposition 2]. Such distributions have infinite mean. Second, more generally, TL holds asymptotically with exponent b = (2α)/(1α) for the sample variance and sample mean of samples of increasing size from nonnegative distributions with regularly varying upper tails and tail index α ∈ (0, 1) [17, page 6, their eq. (3.2)]. Again, such distributions have infinite mean. Third, TL holds asymptotically with exponent b = (2α)/(1α) for samples of increasing size from nonnegative distributions with regularly varying upper tails and tail index α ∈ (0, 1) when the sample variance is replaced by the sample upper semivariance [13, page 4, Theorem 2, their eq. 19].

Acknowledgments

We thank two reviewers for excellent constructive comments. T.H. acknowledges partial support from the "Chaire Modélisation mathématique et biodiversité" of Veolia-Ecole Polytechnique-MNHN-Fondation X, and support from the labex MME-DII Center of Excellence (Modèles mathématiques et économiques de la dynamique, de l'incertitude et des interactions, ANR-11-LABX-0023-01 project). This work was also funded by CY Initiative of Excellence (grant "Investissements d'Avenir "ANR-16-IDEX-0008), Project "EcoDep" PSI-AAP2020-0000000013. The ECODEP project is organized by Paul Doukhan. We are grateful for this opportunity to collaborate.