An elliptic curve is given by a Weierstrass model

y² + a₁ xy + a₃ y = x³ + a₂ x² + a₄ x + a₆,

whose discriminant is nonzero. One can also define an elliptic curve with a

y² = x³ + a₄ x + a₆.

Affine points on E are represented as two-component vectors [x,y]; the point at infinity, i.e. the identity element of the group law, is represented by the one-component vector [0].

Given a vector of coefficients [a₁,a₂,a₃,a₄,a₆] or [a₄,a₆], the function ellinit initializes and returns an ell structure. An additional optional argument allows to specify the base field in case it cannot be inferred from the curve coefficients. This structure contains data needed by elliptic curve related functions, and is generally passed as a first argument. Expensive data are skipped on initialization: they will be dynamically computed when (and if) needed, and then inserted in the structure. The precise layout of the ell structure is left undefined and should never be used directly. The following member functions are available, depending on the underlying domain.

All domains.

* a1, a2, a3, a4, a6: coefficients of the elliptic curve.

* b2, b4, b6, b8: b-invariants of the curve; in characteristic ! = 2, for Y = 2y + a₁x + a₃, the curve equation becomes Y² = 4 x³ + b₂ x² + 2b₄ x + b₆ = : g(x).

* c4, c6: c-invariants of the curve; in characteristic ! = 2,3, for X = x + b₂/12 and Y = 2y + a₁x + a₃, the curve equation becomes Y² = 4 X³ - (c₄/12) X - (c₆/216).

* disc: discriminant of the curve. This is only required to be nonzero, not necessarily a unit.

* j: j-invariant of the curve.

These are used as follows:

  ? E = ellinit([0,0,0, a4,a6]);
  ? E.b4
  %2 = 2*a4
  ? E.disc
  %3 = -64*a4^3 - 432*a6^2

Curves over ℂ.

This in particular includes curves defined over ℚ. All member functions in this section return data, as it is currently stored in the structure, if present; and otherwise compute it to the default accuracy, that was fixed at the time of ellinit (via a t_REAL D domain argument, or realprecision by default). The function ellperiods allows to recompute (and cache) the following data to current realprecision.

* area: volume of the complex lattice defining E.

* roots is a vector whose three components contain the complex roots of the right hand side g(x) of the attached b-model Y² = g(x). If the roots are all real, they are ordered by decreasing value. If only one is real, it is the first component.

* omega: [ω₁,ω₂], periods forming a basis of the complex lattice defining E. The first component ω₁ is the (positive) real period, in other words the integral of the Néron differential dx/(2y+a₁x+a₃) over the connected component of the identity component of E(ℝ). The second component ω₂ is a complex period, such that τ = (ω₁)/(ω₂) belongs to Poincaré's half-plane (positive imaginary part); not necessarily to the standard fundamental domain. It is normalized so that Im(ω₂) < 0 and either Re(ω₂) = 0, when E.disc > 0 (E(ℝ) has two connected components), or Re(ω₂) = ω₁/2

* eta is a row vector containing the quasi-periods η₁ and η₂ such that η_i = 2ζ(ω_i/2), where ζ is the Weierstrass zeta function attached to the period lattice; see ellzeta. In particular, the Legendre relation holds: η₂ω₁ - η₁ω₂ = 2π i.

Warning. As for the orientation of the basis of the period lattice, beware that many sources use the inverse convention where ω₂/ω₁ has positive imaginary part and our ω₂ is the negative of theirs. Our convention τ = ω₁/ω₂ ensures that the action of PSL₂ is the natural one: [a,b;c,d].τ = (aτ+b)/(cτ+d) = (a ω₁ + bω₂)/(cω₁ + dω₂), instead of a twisted one. (Our τ is -1/τ in the above inverse convention.)

Curves over ℚ_p.

We advise to input a model defined over ℚ for such curves. In any case, if you input an approximate model with t_PADIC coefficients, it will be replaced by a lift to ℚ (an exact model "close" to the one that was input) and all quantities will then be computed in terms of this lifted model.

For the time being only curves with multiplicative reduction (split or nonsplit), i.e. v_p(j) < 0, are supported by nontrivial functions. In this case the curve is analytically isomorphic to ℚ_p^*/q^ℤ := E_q(ℚ_p), for some p-adic integer q (the Tate period). In particular, we have j(q) = j(E).

* p is the residual characteristic

* roots is a vector with a single component, equal to the p-adic root e₁ of the right hand side g(x) of the attached b-model Y² = g(x). The point (e₁,0) corresponds to -1 ∈ ℚ_p^*/q^ℤ under the Tate parametrization.

* tate returns [u²,u,q,[a,b],Ei,L] in the notation of Henniart-Mestre (CRAS t. 308, p. 391–395, 1989): q is as above, u ∈ ℚ_p(sqrt{-c₆}) is such that φ^* dx/(2y + a₁x + a₃) = u dt/t, where φ: E_q → E is an isomorphism (well defined up to sign) and dt/t is the canonical invariant differential on the Tate curve; u² ∈ ℚ_p does not depend on φ. (Technicality: if u ∉ ℚ_p, it is stored as a quadratic t_POLMOD.) The parameters [a,b] satisfy 4u² b.agm(sqrt{a/b},1)² = 1 as in Theorem 2 (loc. cit.). Ei describes the sequence of 2-isogenous curves (with kernel generated by [0,0]) E_i: y² = x(x+A_i)(x+A_i-B_i) converging quadratically towards the singular curve E _oo . Finally, L is Mazur-Tate-Teitelbaum's ℒ-invariant, equal to log_p q / v_p(q).

Curves over 𝔽_q.

* p is the characteristic of 𝔽_q.

* no is #E(𝔽_q).

* cyc gives the cycle structure of E(𝔽_q).

* gen returns the generators of E(𝔽_q).

* group returns [no,cyc,gen], i.e. E(𝔽_q) as an abelian group structure.

Curves over ℚ.

All functions should return a correct result, whether the model is minimal or not, but it is a good idea to stick to minimal models whenever gcd(c₄,c₆) is easy to factor (minor speed-up). The construction

    E = ellminimalmodel(E0, &v)

replaces the original model E₀ by a minimal model E, and the variable change v allows to go between the two models:

    ellchangepoint(P0, v)
    ellchangepointinv(P, v)

respectively map the point P₀ on E₀ to its image on E, and the point P on E to its pre-image on E₀.

A few routines — namely ellgenerators, ellidentify, ellsearch, forell — require the optional package elldata (John Cremona's database) to be installed. In that case, the function ellinit will allow alternative inputs, e.g. ellinit("11a1"). Functions using this package need to load chunks of a large database in memory and require at least 2MB stack to avoid stack overflows.

* gen returns the generators of E(ℚ), if known (from John Cremona's database)

Curves over number fields.

* nf return the nf structure attached to the number field over which E is defined.

* bnf return the bnf structure attached to the number field over which E is defined or raise an error (if only an nf is available).

* omega, eta, area: vectors of complex periods, quasi-periods and lattice areas attached to the complex embeddings of E, in the same order as E.nf.roots.

Let E be a curve defined over ℚ_p given by a p-integral model; if the curve has good reduction at p, we may define its reduction ~{E} over the finite field 𝔽_p:

  ? E = ellinit([-3,1], O(5^10));  \\   E/ℚ5
  ? Et = ellinit(E, 5)
  ? ellcard(Et)  \\   ~{E}/𝔽5 has 7 points
  %3 = 7
  ? ellinit(E, 7)
   ***   at top-level: ellinit(E,7)
   ***                 ^ —  —  —  — 
   *** ellinit: inconsistent moduli in ellinit: 5 != 7

Likewise, if a curve is defined over a number field K and 𝔭 is a maximal ideal with finite residue field 𝔽_q, we define the reduction ~{E}/𝔽_q provided E has good reduction at 𝔭. E/ℚ is an important special case:

  ? E = ellinit([-3,1]);
  ? factor(E.disc)
  %2 =
  [2 4]
  
  [3 4]
  ? Et = ellinit(E, 5);
  ? ellcard(Et) \\   ~{E} / 𝔽5 has 7 points
  %4 = 7
  ? ellinit(E, 3)  \\ bad reduction at 3
  %5 = []

General number fields are similar:

  ? K = nfinit(x^2+1); E = ellinit([x,x+1], K);
  ? idealfactor(K, E.disc)  \\ three primes of bad reduction
  %2 =
  [  [2, [1, 1]~, 2, 1, [1, -1; 1, 1]] 10]
  
  [ [5, [-2, 1]~, 1, 1, [2, -1; 1, 2]]  2]
  
  [[5, [2, 1]~, 1, 1, [-2, -1; 1, -2]]  2]
  ? P = idealprimedec(K, 3); \\ a prime of good reduction
  ? idealnorm(K, P)
  %4 = 9
  ? Et = ellinit(E, P);
  ? ellcard(Et)  \\   ~{E} / 𝔽9 has 4 points
  %6 = 4

If the model is not locally minimal at 𝔭, the above will fail: elllocalred and ellchangecurve allow to reduce to that case.

Some functions such as ellap, ellcard, ellgroup and ellissupersingular even implicitly replace the given equation by a local minimal model and consider the group of nonsingular points ~{E}^ns so they make sense even when the curve has bad reduction.

If E is an elliptic curve over ℚ, returns a basis of the set of everywhere locally soluble 2-covers of the curve E. For each cover a pair [R,P] is returned where y²-R(x) is a quartic curve and P is a point on E(k), where k = ℚ(x)[y] / (y²-R(x)). E can also be given as the output of ellrankinit(E), or as a pair [e, f], where e is an elliptic curve given by ellrankinit and f is a quadratic twist of e. We then look for points on f.

  ? E = ellinit([-25,4]);
  ? C = ell2cover(E); #C
  %2 = 2
  ? [R,P] = C[1]; R
  %3 = 64*x^4+480*x^2-128*x+100
  ? P[1]
  %4 = -320/y^2*x^4 + 256/y^2*x^3 + 800/y^2*x^2 - 320/y^2*x - 436/y^2
  ? ellisoncurve(E, Mod(P, y^2-R))
  %5 = 1
  ? H = hyperellratpoints(R,10)
  %6 = [[0,10], [0,-10], [1/5,242/25], [1/5,-242/25], [2/5,282/25],
        [2/5,-282/25]]
  ? A = substvec(P,[x,y],H[1])
  %7 = [-109/25, 686/125]

The library syntax is GEN ell2cover(GEN E, long prec).

Returns the value at s = 1 of the derivative of order r of the L-function of the elliptic curve E/ℚ.

  ? E = ellinit("11a1"); \\ order of vanishing is 0
  ? ellL1(E)
  %2 = 0.2538418608559106843377589233
  ? E = ellinit("389a1");  \\ order of vanishing is 2
  ? ellL1(E)
  %4 = -5.384067311837218089235032414 E-29
  ? ellL1(E, 1)
  %5 = 0
  ? ellL1(E, 2)
  %6 = 1.518633000576853540460385214

The main use of this function, after computing at low accuracy the order of vanishing using ellanalyticrank, is to compute the leading term at high accuracy to check (or use) the Birch and Swinnerton-Dyer conjecture:

  ? \p18
    realprecision = 18 significant digits
  ? E = ellinit("5077a1"); ellanalyticrank(E)
  time = 8 ms.
  %1 = [3, 10.3910994007158041]
  ? \p200
    realprecision = 202 significant digits (200 digits displayed)
  ? ellL1(E, 3)
  time = 104 ms.
  %3 = 10.3910994007158041387518505103609170697263563756570092797 [...]

Analogous and more general functionalities for E defined over general number fields are available through lfun.

The library syntax is GEN ellL1(GEN E, long r, long bitprec).

Sum of the points z1 and z2 on the elliptic curve corresponding to E.

The library syntax is GEN elladd(GEN E, GEN z1, GEN z2).

Computes the coefficient a_n of the L-function of the elliptic curve E/ℚ, i.e. coefficients of a newform of weight 2 by the modularity theorem (Taniyama-Shimura-Weil conjecture). E must be an ell structure over ℚ as output by ellinit. E must be given by an integral model, not necessarily minimal, although a minimal model will make the function faster.

  ? E = ellinit([1,-1,0,4,3]);
  ? ellak(E, 10)
  %2 = -3
  ? e = ellchangecurve(E, [1/5,0,0,0]); \\ made not minimal at 5
  ? ellak(e, 10) \\ wasteful but works
  %3 = -3
  ? E = ellminimalmodel(e); \\ now minimal
  ? ellak(E, 5)
  %5 = -3

If the model is not minimal at a number of bad primes, then the function will be slower on those n divisible by the bad primes. The speed should be comparable for other n:

  ? for(i=1,10^6, ellak(E,5))
  time = 699 ms.
  ? for(i=1,10^6, ellak(e,5)) \\ 5 is bad, markedly slower
  time = 1,079 ms.
  
  ? for(i=1,10^5,ellak(E,5*i))
  time = 1,477 ms.
  ? for(i=1,10^5,ellak(e,5*i)) \\ still slower but not so much on average
  time = 1,569 ms.

The library syntax is GEN akell(GEN E, GEN n).

Computes the vector of the first n Fourier coefficients a_k corresponding to the elliptic curve E defined over a number field. If E is defined over ℚ, the curve may be given by an arbitrary model, not necessarily minimal, although a minimal model will make the function faster. Over a more general number field, the model must be locally minimal at all primes above 2 and 3.

The library syntax is GEN ellan(GEN E, long n). Also available is GEN ellanQ_zv(GEN e, long n), which returns a t_VECSMALL instead of a t_VEC, saving on memory.

Returns the order of vanishing at s = 1 of the L-function of the elliptic curve E/ℚ and the value of the first nonzero derivative. To determine this order, it is assumed that any value less than eps is zero. If eps is omitted, 2^-b/2 is used, where b is the current bit precision.

  ? E = ellinit("11a1"); \\ rank 0
  ? ellanalyticrank(E)
  %2 = [0, 0.2538418608559106843377589233]
  ? E = ellinit("37a1"); \\ rank 1
  ? ellanalyticrank(E)
  %4 = [1, 0.3059997738340523018204836835]
  ? E = ellinit("389a1"); \\ rank 2
  ? ellanalyticrank(E)
  %6 = [2, 1.518633000576853540460385214]
  ? E = ellinit("5077a1"); \\ rank 3
  ? ellanalyticrank(E)
  %8 = [3, 10.39109940071580413875185035]

Analogous and more general functionalities for E defined over general number fields are available through lfun and lfunorderzero.

The library syntax is GEN ellanalyticrank(GEN E, GEN eps = NULL, long bitprec).

Let E be an ell structure as output by ellinit, attached to an elliptic curve E/K. If the field K = 𝔽_q is finite, return the trace of Frobenius t, defined by the equation #E(𝔽_q) = q+1 - t.

For other fields of definition and p defining a finite residue field 𝔽_q, return the trace of Frobenius for the reduction of E: the argument p is best left omitted if K = ℚ_ℓ (else we must have p = ℓ) and must be a prime number (K = ℚ) or prime ideal (K a general number field) with residue field 𝔽_q otherwise. The equation need not be minimal or even integral at p; of course, a minimal model will be more efficient.

For a number field K, the trace of Frobenius is the a_p coefficient in the Euler product defining the curve L-series, whence the function name: L(E/K,s) = ∏_{bad p} (1-a_p (Np)^-s)^-1 ∏_{good p} (1-a_p (Np)^-s + (Np)^1-2s)^-1.

When the characteristic of the finite field is large, the availability of the seadata package will speed up the computation.

  ? E = ellinit([0,1]);  \\ y^2 = x^3 + 0.x + 1, defined over Q
  ? ellap(E, 7) \\ 7 necessary here
  %2 = -4       \\ #E(F7) = 7+1-(-4) = 12
  ? ellcard(E, 7)
  %3 = 12       \\ OK
  
  ? E = ellinit([0,1], 11);  \\ defined over F_11
  ? ellap(E)       \\ no need to repeat 11
  %4 = 0
  ? ellap(E, 11)   \\ ... but it also works
  %5 = 0
  ? ellgroup(E, 13) \\ ouch, inconsistent input!
     ***   at top-level: ellap(E,13)
     ***                 ^ —  —  — --
     *** ellap: inconsistent moduli in Rg_to_Fp:
       11
       13
  ? a = ffgen(ffinit(11,3), 'a); \\ defines Fq := F11^3
  ? E = ellinit([a+1,a]);  \\ y^2 = x^3 + (a+1)x + a, defined over Fq
  ? ellap(E)
  %8 = -3

If the curve is defined over a more general number field than ℚ, the maximal ideal p must be explicitly given in idealprimedec format. There is no assumption of local minimality at p.

  ? K = nfinit(a^2+1); E = ellinit([1+a,0,1,0,0], K);
  ? fa = idealfactor(K, E.disc)
  %2 =
  [  [5, [-2, 1]~, 1, 1, [2, -1; 1, 2]] 1]
  
  [[13, [5, 1]~, 1, 1, [-5, -1; 1, -5]] 2]
  ? ellap(E, fa[1,1])
  %3 = -1 \\ nonsplit multiplicative reduction
  ? ellap(E, fa[2,1])
  %4 = 1  \\ split multiplicative reduction
  ? P17 = idealprimedec(K,17)[1];
  ? ellap(E, P17)
  %6 = 6  \\ good reduction
  ? E2 = ellchangecurve(E, [17,0,0,0]);
  ? ellap(E2, P17)
  %8 = 6  \\ same, starting from a nonminimal model
  
  ? P3 = idealprimedec(K,3)[1];
  ? ellap(E, P3)  \\ OK: E is minimal at P3
  %10 = -2
  ? E3 = ellchangecurve(E, [3,0,0,0]);
  ? ellap(E3, P3) \\ not integral at P3
   ***   at top-level: ellap(E3,P3)
   ***                 ^ —  —  —  — 
   *** ellap: impossible inverse in Rg_to_ff: Mod(0, 3).

Algorithms used. If E/𝔽_q has CM by a principal imaginary quadratic order we use a fast explicit formula (involving essentially Kronecker symbols and Cornacchia's algorithm), in O(log q)² bit operations. Otherwise, we use Shanks-Mestre's baby-step/giant-step method, which runs in time Õ(q^1/4) using Õ(q^1/4) storage, hence becomes unreasonable when q has about 30 digits. Above this range, the SEA algorithm becomes available, heuristically in Õ(log q)⁴, and primes of the order of 200 digits become feasible. In small characteristic we use Mestre's (p = 2), Kohel's (p = 3,5,7,13), Satoh-Harley (all in Õ(p² n²)) or Kedlaya's (in Õ(p n³)) algorithms.

The library syntax is GEN ellap(GEN E, GEN p = NULL).

Deprecated alias for ellheight(E,P,Q).

The library syntax is GEN bilhell(GEN E, GEN z1, GEN z2, long prec).

E being an elliptic curve over a number field, returns a real number c such that the Birch and Swinnerton-Dyer conjecture predicts that L_E^(r)(1)/r! = c R S, where r is the rank, R the regulator and S the cardinal of the Tate-Shafarevich group.

  ? e = ellinit([0,-1,1,-10,-20]); \\ rank 0
  ? ellbsd(e)
  %2 = 0.25384186085591068433775892335090946105
  ? lfun(e,1)
  %3 = 0.25384186085591068433775892335090946104
  ? e = ellinit([0,0,1,-1,0]); \\ rank 1
  ? P = ellheegner(e);
  ? ellbsd(e)*ellheight(e,P)
  %6 = 0.30599977383405230182048368332167647445
  ? lfun(e,1,1)
  %7 = 0.30599977383405230182048368332167647445
  ? e = ellinit([1+a,0,1,0,0],nfinit(a^2+1)); \\ rank 0
  ? ellbsd(e)
  %9 = 0.42521832235345764503001271536611593310
  ? lfun(e,1)
  %10 = 0.42521832235345764503001271536611593309

The library syntax is GEN ellbsd(GEN E, long prec).

Let E be an ell structure as output by ellinit, attached to an elliptic curve E/K. If K = 𝔽_q is finite, return the order of the group E(𝔽_q).

  ? E = ellinit([-3,1], 5); ellcard(E)
  %1 = 7
  ? t = ffgen(3^5,'t); E = ellinit([t,t^2+1]); ellcard(E)
  %2 = 217

For other fields of definition and p defining a finite residue field 𝔽_q, return the order of the reduction of E: the argument p is best left omitted if K = ℚ_ℓ (else we must have p = ℓ) and must be a prime number (K = ℚ) or prime ideal (K a general number field) with residue field 𝔽_q otherwise. The equation need not be minimal or even integral at p; of course, a minimal model will be more efficient. The function considers the group of nonsingular points of the reduction of a minimal model of the curve at p, so also makes sense when the curve has bad reduction.

  ? E = ellinit([-3,1]);
  ? factor(E.disc)
  %2 =
  [2 4]
  
  [3 4]
  ? ellcard(E, 5)  \\ as above !
  %3 = 7
  ? ellcard(E, 2) \\ additive reduction
  %4 = 2

When the characteristic of the finite field is large, the availability of the seadata package will speed the computation. See also ellap for the list of implemented algorithms.

The library syntax is GEN ellcard(GEN E, GEN p = NULL). Also available is GEN ellcard(GEN E, GEN p) where p is not NULL.

Changes the data for the elliptic curve E by changing the coordinates using the vector v = [u,r,s,t], i.e. if x' and y' are the new coordinates, then x = u²x'+r, y = u³y'+su²x'+t. E must be an ell structure as output by ellinit. The special case v = 1 is also used instead of [1,0,0,0] to denote the trivial coordinate change.

The library syntax is GEN ellchangecurve(GEN E, GEN v).

Changes the coordinates of the point or vector of points x using the vector v = [u,r,s,t], i.e. if x' and y' are the new coordinates, then x = u²x'+r, y = u³y'+su²x'+t (see also ellchangecurve).

  ? E0 = ellinit([1,1]); P0 = [0,1]; v = [1,2,3,4];
  ? E = ellchangecurve(E0, v);
  ? P = ellchangepoint(P0,v)
  %3 = [-2, 3]
  ? ellisoncurve(E, P)
  %4 = 1
  ? ellchangepointinv(P,v)
  %5 = [0, 1]

The library syntax is GEN ellchangepoint(GEN x, GEN v). The reciprocal function GEN ellchangepointinv(GEN x, GEN ch) inverts the coordinate change.

Changes the coordinates of the point or vector of points x using the inverse of the isomorphism attached to v = [u,r,s,t], i.e. if x' and y' are the old coordinates, then x = u²x'+r, y = u³y'+su²x'+t (inverse of ellchangepoint).

  ? E0 = ellinit([1,1]); P0 = [0,1]; v = [1,2,3,4];
  ? E = ellchangecurve(E0, v);
  ? P = ellchangepoint(P0,v)
  %3 = [-2, 3]
  ? ellisoncurve(E, P)
  %4 = 1
  ? ellchangepointinv(P,v)
  %5 = [0, 1]  \\ we get back P0

The library syntax is GEN ellchangepointinv(GEN x, GEN v).

Converts an elliptic curve name, as found in the elldata database, from a string to a triplet [conductor, isogeny class, index]. It will also convert a triplet back to a curve name. Examples:

  ? ellconvertname("123b1")
  %1 = [123, 1, 1]
  ? ellconvertname(%)
  %2 = "123b1"

The library syntax is GEN ellconvertname(GEN name).

n-division polynomial f_n for the curve E in the variable v. In standard notation, for any affine point P = (X,Y) on the curve and any integer n ≥ 0, we have [n]P = (φ_n(P)ψ_n(P) : ω_n(P) : ψ_n(P)³) for some polynomials φ_n,ω_n,ψ_n in ℤ[a₁,a₂,a₃,a₄,a₆][X,Y]. We have f_n(X) = ψ_n(X) for n odd, and f_n(X) = ψ_n(X,Y) (2Y + a₁X+a₃) for n even. We have f₀ = 0, f₁ = 1, f₂ = 4X³ + b₂X² + 2b₄ X + b₆, f₃ = 3 X⁴ + b₂ X³ + 3b₄ X² + 3 b₆ X + b8, f₄ = f₂(2X⁶ + b₂ X⁵ + 5b₄ X⁴ + 10 b₆ X³ + 10 b₈ X² + (b₂b₈-b₄b₆)X + (b₈b₄ - b₆²)), ... When n is odd, the roots of f_n are the X-coordinates of the affine points in the n-torsion subgroup E[n]; when n is even, the roots of f_n are the X-coordinates of the affine points in E[n] \ E[2] when n > 2, resp. in E[2] when n = 2. For n < 0, we define f_n := - f_-n.

The library syntax is GEN elldivpol(GEN E, long n, long v = -1) where v is a variable number.

k being an even positive integer, computes the numerical value of the Eisenstein series of weight k at the lattice w, as given by ellperiods, namely

(2i π/ω₂)^k (1 + 2/ζ(1-k) ∑_{n ≥ 1} n^k-1qⁿ / (1-qⁿ)),

where q = exp(2iπ τ) and τ := ω₁/ω₂ belongs to the complex upper half-plane. It is also possible to directly input w = [ω₁,ω₂], or an elliptic curve E as given by ellinit.

  ? w = ellperiods([1,I]);
  ? elleisnum(w, 4)
  %2 = 2268.8726415508062275167367584190557607
  ? elleisnum(w, 6)
  %3 = -3.977978632282564763 E-33
  ? E = ellinit([1, 0]);
  ? elleisnum(E, 4)
  %5 = -48.000000000000000000000000000000000000

When flag is nonzero and k = 4 or 6, returns the elliptic invariants g₂ or g₃, such that y² = 4x³ - g₂ x - g₃ is a Weierstrass equation for E.

  ? g2 = elleisnum(E, 4, 1)
  %6 = -4.0000000000000000000000000000000000000
  ? g3 = elleisnum(E, 6, 1)  \\ ~ 0
  %7 = 0.E-114 - 3.909948178422242682 E-57*I

The library syntax is GEN elleisnum(GEN w, long k, long flag, long prec).

Returns the quasi-periods [η₁,η₂] attached to the lattice basis w = [ω₁, ω₂]. Alternatively, w can be an elliptic curve E as output by ellinit, in which case, the quasi periods attached to the period lattice basis E.omega (namely, E.eta) are returned.

  ? elleta([1, I])
  %1 = [3.141592653589793238462643383, 9.424777960769379715387930149*I]

The library syntax is GEN elleta(GEN w, long prec).

Let ω := dx / (2y+a₁x+a₃) be the invariant differential form attached to the model E of some elliptic curve (ellinit form), and η := x(t)ω. Return n terms (seriesprecision by default) of f(t),g(t) two power series in the formal parameter t = -x/y such that ω = f(t) dt, η = g(t) dt: f(t) = 1+a₁ t + (a₁² + a₂) t² +..., g(t) = t^-2 +...

   ? E = ellinit([-1,1/4]); [f,g] = ellformaldifferential(E,7,'t);
   ? f
   %2 = 1 - 2*t^4 + 3/4*t^6 + O(t^7)
   ? g
   %3 = t^-2 - t^2 + 1/2*t^4 + O(t^5)

The library syntax is GEN ellformaldifferential(GEN E, long precdl, long n = -1) where n is a variable number.

The elliptic formal exponential Exp attached to E is the isomorphism from the formal additive law to the formal group of E. It is normalized so as to be the inverse of the elliptic logarithm (see ellformallog): Exp o L = Id. Return n terms of this power series:

  ? E=ellinit([-1,1/4]); Exp = ellformalexp(E,10,'z)
  %1 = z + 2/5*z^5 - 3/28*z^7 + 2/15*z^9 + O(z^11)
  ? L = ellformallog(E,10,'t);
  ? subst(Exp,z,L)
  %3 = t + O(t^11)

The library syntax is GEN ellformalexp(GEN E, long precdl, long n = -1) where n is a variable number.

The formal elliptic logarithm is a series L in t K[[t]] such that d L = ω = dx / (2y + a₁x + a₃), the canonical invariant differential attached to the model E. It gives an isomorphism from the formal group of E to the additive formal group.

  ? E = ellinit([-1,1/4]); L = ellformallog(E, 9, 't)
  %1 = t - 2/5*t^5 + 3/28*t^7 + 2/3*t^9 + O(t^10)
  ? [f,g] = ellformaldifferential(E,8,'t);
  ? L' - f
  %3 = O(t^8)

The library syntax is GEN ellformallog(GEN E, long precdl, long n = -1) where n is a variable number.

If E is an elliptic curve, return the coordinates x(t), y(t) in the formal group of the elliptic curve E in the formal parameter t = -x/y at oo : x = t^-2 -a₁ t^-1 - a₂ - a₃ t +... y = - t^-3 -a₁ t^-2 - a₂t^-1 -a₃ +... Return n terms (seriesprecision by default) of these two power series, whose coefficients are in ℤ[a₁,a₂,a₃,a₄,a₆].

  ? E = ellinit([0,0,1,-1,0]); [x,y] = ellformalpoint(E,8,'t);
  ? x
  %2 = t^-2 - t + t^2 - t^4 + 2*t^5 + O(t^6)
  ? y
  %3 = -t^-3 + 1 - t + t^3 - 2*t^4 + O(t^5)
  ? E = ellinit([0,1/2]); ellformalpoint(E,7)
  %4 = [x^-2 - 1/2*x^4 + O(x^5), -x^-3 + 1/2*x^3 + O(x^4)]

The library syntax is GEN ellformalpoint(GEN E, long precdl, long n = -1) where n is a variable number.

Return the formal power series w attached to the elliptic curve E, in the variable t: w(t) = t³(1 + a₁ t + (a₂ + a₁²) t² +...+ O(tⁿ)), which is the formal expansion of -1/y in the formal parameter t := -x/y at oo (take n = seriesprecision if n is omitted). The coefficients of w belong to ℤ[a₁,a₂,a₃,a₄,a₆].

  ? E=ellinit([3,2,-4,-2,5]); ellformalw(E, 5, 't)
  %1 = t^3 + 3*t^4 + 11*t^5 + 35*t^6 + 101*t^7 + O(t^8)

The library syntax is GEN ellformalw(GEN E, long precdl, long n = -1) where n is a variable number.

Given a genus 1 plane curve, defined by the affine equation f(x,y) = 0, return the coefficients [a₁,a₂,a₃,a₄,a₆] of a Weierstrass equation for its Jacobian. This allows to recover a Weierstrass model for an elliptic curve given by a general plane cubic or by a binary quartic or biquadratic model. The function implements the f ⟼ f^* formulae of Artin, Tate and Villegas (Advances in Math. 198 (2005), pp. 366–382).

In the example below, the function is used to convert between twisted Edwards coordinates and Weierstrass coordinates.

  ? e = ellfromeqn(a*x^2+y^2 - (1+d*x^2*y^2))
  %1 = [0, -a - d, 0, -4*d*a, 4*d*a^2 + 4*d^2*a]
  ? E = ellinit(ellfromeqn(y^2-x^2 - 1 +(121665/121666*x^2*y^2)),2^255-19);
  ? isprime(ellcard(E) / 8)
  %3 = 1

The elliptic curve attached to the sum of two cubes is given by

  ? ellfromeqn(x^3+y^3 - a)
  %1 = [0, 0, -9*a, 0, -27*a^2]

Congruent number problem. Let n be an integer, if a²+b² = c² and a b = 2 n, then by substituting b by 2 n/a in the first equation, we get ((a²+(2 n/a)²)-c²) a² = 0. We set x = a, y = a c.

  ? En = ellfromeqn((x^2 + (2*n/x)^2 - (y/x)^2)*x^2)
  %1 = [0, 0, 0, -16*n^2, 0]

For example 23 is congruent since the curve has a point of infinite order, namely:

  ? ellheegner( ellinit(subst(En, n, 23)) )
  %2 = [168100/289, 68053440/4913]

The library syntax is GEN ellfromeqn(GEN P).

Returns the coefficients [a₁,a₂,a₃,a₄,a₆] of a fixed elliptic curve with j-invariant j. The given model is arbitrary; for instance, over the rationals, it is in general not minimal nor even integral.

  ? v = ellfromj(1/2)
  %1 = [0, 0, 0, 10365/4, 11937025/4]
  ? E = ellminimalmodel(ellinit(v)); E[1..5]
  %2 = [0, 0, 0, 41460, 190992400]
  ? F = ellminimalmodel(elltwist(E, 24)); F[1..5]
  %3 = [1, 0, 0, 72, 13822]
  ? [E.disc, F.disc]
  %4 = [-15763098924417024000, -82484842750]

For rational j, the following program returns the integral curve of minimal discriminant and given j invariant:

  ellfromjminimal(j)=
  { my(E = ellinit(ellfromj(j)));
    my(D = ellminimaltwist(E));
  
    ellminimalmodel(elltwist(E,D));
  }
  ? e = ellfromjminimal(1/2); e.disc
  %1 = -82484842750

Using flag = 1 in ellminimaltwist would instead return the curve of minimal conductor. For instance, if j = 1728, this would return a different curve (of conductor 32 instead of 64).

The library syntax is GEN ellfromj(GEN j).

If E is an elliptic curve over the rationals, return a ℤ-basis of the free part of the Mordell-Weil group attached to E. This relies on the elldata database being installed and referencing the curve, and so is only available for curves over ℤ of small conductors. If E is an elliptic curve over a finite field 𝔽_q as output by ellinit, return a minimal set of generators for the group E(𝔽_q).

Caution. When the group is not cyclic, of shape ℤ/d₁ℤ x ℤ/d₂ℤ with d₂ | d₁, the points [P,Q] returned by ellgenerators need not have order d₁ and d₂: it is true that P has order d₁, but we only know that Q is a generator of E(𝔽_q)/ <P> and that the Weil pairing w(P,Q) has order d₂, see ??ellgroup. If you need generators [P,R] with R of order d₂, find x such that R = Q-[x]P has order d₂ by solving the discrete logarithm problem [d₂]Q = [x]([d₂]P) in a cyclic group of order d₁/d₂. This will be very expensive if d₁/d₂ has a large prime factor.

The library syntax is GEN ellgenerators(GEN E).

Let E be an ell structure as output by ellinit attached to an elliptic curve defined over a number field. This function calculates the arithmetic conductor and the global Tamagawa number c. The result [N,v,c,F,L] is slightly different if E is defined over ℚ (domain D = 1 in ellinit) or over a number field (domain D is a number field structure, including nfinit(x) representing ℚ !):

* N is the arithmetic conductor of the curve,

* v is an obsolete field, left in place for backward compatibility. If E is defined over ℚ, v gives the coordinate change for E to the standard minimal integral model (ellminimalmodel provides it in a cheaper way); if E is defined over another number field, v gives a coordinate change to an integral model (ellintegralmodel provides it in a cheaper way).

* c is the product of the local Tamagawa numbers c_p, a quantity which enters in the Birch and Swinnerton-Dyer conjecture,

* F is the factorization of N,

* L is a vector, whose i-th entry contains the local data at the i-th prime ideal divisor of N, i.e. L[i] = elllocalred(E,F[i,1]). If E is defined over ℚ, the local coordinate change has been deleted and replaced by a 0; if E is defined over another number field the local coordinate change to a local minimal model is given relative to the integral model afforded by v (so either start from an integral model so that v be trivial, or apply v first).

The library syntax is GEN ellglobalred(GEN E).

Let E be an ell structure as output by ellinit, attached to an elliptic curve E/K. We first describe the function when the field K = 𝔽_q is finite, it computes the structure of the finite abelian group E(𝔽_q):

* if flag = 0, returns the structure [] (trivial group) or [d₁] (nontrivial cyclic group) or [d₁,d₂] (noncyclic group) of E(𝔽_q) ~ ℤ/d₁ℤ x ℤ/d₂ℤ, with d₂ | d₁.

* if flag = 1, returns a triple [h,cyc,gen], where h is the curve cardinality, cyc gives the group structure as a product of cyclic groups (as per flag = 0). More precisely, if d₂ > 1, the output is [d₁d₂, [d₁,d₂], [P,Q]] where P is of order d₁ and [P,Q] generates the curve. Caution. It is not guaranteed that Q has order d₂, which in the worst case requires an expensive discrete log computation. Only that ellweilpairing(E, P, Q, d₁) has order d₂.

For other fields of definition and p defining a finite residue field 𝔽_q, returns the structure of the reduction of E: the argument p is best left omitted if K = ℚ_ℓ (else we must have p = ℓ) and must be a prime number (K = ℚ) or prime ideal (K a general number field) with residue field 𝔽_q otherwise. The curve is allowed to have bad reduction at p and in this case we consider the (cyclic) group of nonsingular points for the reduction of a minimal model at p.

If flag = 0, the equation need not be minimal or even integral at p; of course, a minimal model will be more efficient.

If flag = 1, the requested generators depend on the model, which must then be minimal at p, otherwise an exception is thrown. Use ellintegralmodel and/or ellocalred first to reduce to this case.

  ? E = ellinit([0,1]);  \\ y^2 = x^3 + 0.x + 1, defined over Q
  ? ellgroup(E, 7)
  %2 = [6, 2] \\ Z/6 x Z/2, noncyclic
  ? E = ellinit([0,1] * Mod(1,11));  \\ defined over F_11
  ? ellgroup(E)   \\ no need to repeat 11
  %4 = [12]
  ? ellgroup(E, 11)   \\ ... but it also works
  %5 = [12]
  ? ellgroup(E, 13) \\ ouch, inconsistent input!
     ***   at top-level: ellgroup(E,13)
     ***                 ^ —  —  —  — --
     *** ellgroup: inconsistent moduli in Rg_to_Fp:
       11
       13
  ? ellgroup(E, 7, 1)
  %6 = [12, [6, 2], [[Mod(2, 7), Mod(4, 7)], [Mod(4, 7), Mod(4, 7)]]]

Let us now consider curves of bad reduction, in this case we return the structure of the (cyclic) group of nonsingular points, satisfying #E_ns(𝔽_p) = p - a_p:

  ? E = ellinit([0,5]);
  ? ellgroup(E, 5, 1)
  %2 = [5, [5], [[Mod(4, 5), Mod(2, 5)]]]
  ? ellap(E, 5)
  %3 = 0 \\ additive reduction at 5
  ? E = ellinit([0,-1,0,35,0]);
  ? ellgroup(E, 5, 1)
  %5 = [4, [4], [[Mod(2, 5), Mod(2, 5)]]]
  ? ellap(E, 5)
  %6 = 1 \\ split multiplicative reduction at 5
  ? ellgroup(E, 7, 1)
  %7 = [8, [8], [[Mod(3, 7), Mod(5, 7)]]]
  ? ellap(E, 7)
  %8 = -1 \\ nonsplit multiplicative reduction at 7

The library syntax is GEN ellgroup0(GEN E, GEN p = NULL, long flag). Also available is GEN ellgroup(GEN E, GEN p), corresponding to flag = 0.

Let E be an elliptic curve over the rationals, assumed to be of (analytic) rank 1. This returns a nontorsion rational point on the curve, whose canonical height is equal to the product of the elliptic regulator by the analytic Sha.

This uses the Heegner point method, described in Cohen GTM 239; the complexity is proportional to the product of the square root of the conductor and the height of the point (thus, it is preferable to apply it to strong Weil curves).

  ? E = ellinit([-157^2,0]);
  ? u = ellheegner(E); print(u[1], "\n", u[2])
  69648970982596494254458225/166136231668185267540804
  538962435089604615078004307258785218335/67716816556077455999228495435742408
  ? ellheegner(ellinit([0,1]))         \\ E has rank 0 !
   ***   at top-level: ellheegner(E=ellinit
   ***                 ^ —  —  —  —  —  — --
   *** ellheegner: The curve has even analytic rank.

The library syntax is GEN ellheegner(GEN E).

Let E be an elliptic curve defined over K = ℚ or a number field, as output by ellinit; it needs not be given by a minimal model although the computation will be faster if it is.

* Without arguments P,Q, returns the Faltings height of the curve E using Deligne normalization. For a rational curve, the normalization is such that the function returns -(1/2)*log(ellminimalmodel(E).area).

* If the argument P ∈ E(K) is present, returns the global Néron-Tate height h(P) of the point, using the normalization in Cremona's Algorithms for modular elliptic curves.

* If the argument Q ∈ E(K) is also present, computes the value of the bilinear form (h(P+Q)-h(P-Q)) / 4.

The library syntax is GEN ellheight0(GEN E, GEN P = NULL, GEN Q = NULL, long prec). Also available is GEN ellheight(GEN E, GEN P, long prec) (Q omitted).

x being a vector of points, this function outputs the Gram matrix of x with respect to the Néron-Tate height, in other words, the (i,j) component of the matrix is equal to ellheight(E,x[i],x[j]). The rank of this matrix, at least in some approximate sense, gives the rank of the set of points, and if x is a basis of the Mordell-Weil group of E, its determinant is equal to the regulator of E. Note our height normalization follows Cremona's Algorithms for modular elliptic curves: this matrix should be divided by 2 to be in accordance with, e.g., Silverman's normalizations.

The library syntax is GEN ellheightmatrix(GEN E, GEN x, long prec).

Look up the elliptic curve E, defined by an arbitrary model over ℚ, in the elldata database. Return [[N, M, G], C] where N is the curve name in Cremona's elliptic curve database, M is the minimal model, G is a ℤ-basis of the free part of the Mordell-Weil group E(ℚ) and C is the change of coordinates from E to M, suitable for ellchangecurve.

The library syntax is GEN ellidentify(GEN E).

Initialize an ell structure, attached to the elliptic curve E. E is either

* a 5-component vector [a₁,a₂,a₃,a₄,a₆] defining the elliptic curve with Weierstrass equation Y² + a₁ XY + a₃ Y = X³ + a₂ X² + a₄ X + a₆,

* a 2-component vector [a₄,a₆] defining the elliptic curve with short Weierstrass equation Y² = X³ + a₄ X + a₆,

* a single-component vector [j] giving the j-invariant for the curve, with the same coefficients as given by ellfromj.

* a character string in Cremona's notation, e.g. "11a1", in which case the curve is retrieved from the elldata database if available.

The optional argument D describes the domain over which the curve is defined:

* the t_INT 1 (default): the field of rational numbers ℚ.

* a t_INT p, where p is a prime number: the prime finite field 𝔽_p.

* an t_INTMOD Mod(a, p), where p is a prime number: the prime finite field 𝔽_p.

* a t_FFELT, as returned by ffgen: the corresponding finite field 𝔽_q.

* a t_PADIC, O(pⁿ): the field ℚ_p, where p-adic quantities will be computed to a relative accuracy of n digits. We advise to input a model defined over ℚ for such curves. In any case, if you input an approximate model with t_PADIC coefficients, it will be replaced by a lift to ℚ (an exact model "close" to the one that was input) and all quantities will then be computed in terms of this lifted model, at the given accuracy.

* a t_REAL x: the field ℂ of complex numbers, where floating point quantities are by default computed to a relative accuracy of precision(x). If no such argument is given, the value of realprecision at the time ellinit is called will be used.

* a number field K, given by a nf or bnf structure; a bnf is required for ellminimalmodel.

* a prime ideal 𝔭, given by a prid structure; valid if x is a curve defined over a number field K and the equation is integral and minimal at 𝔭.

This argument D is indicative: the curve coefficients are checked for compatibility, possibly changing D; for instance if D = 1 and an t_INTMOD is found. If inconsistencies are detected, an error is raised:

  ? ellinit([1 + O(5), 1], O(7));
   ***   at top-level: ellinit([1+O(5),1],O
   ***                 ^ —  —  —  —  —  — --
   *** ellinit: inconsistent moduli in ellinit: 7 != 5

If the curve coefficients are too general to fit any of the above domain categories, only basic operations, such as point addition, will be supported later.

If the curve (seen over the domain D) is singular, fail and return an empty vector [].

  ? E = ellinit([0,0,0,0,1]); \\ y^2 = x^3 + 1, over Q
  ? E = ellinit([0,1]);       \\ the same curve, short form
  ? E = ellinit("36a1");      \\ sill the same curve, Cremona's notations
  ? E = ellinit([0]);         \\ a curve of j-invariant 0
  ? E = ellinit([0,1], 2)     \\ over F2: singular curve
  %4 = []
  ? E = ellinit(['a4,'a6] * Mod(1,5));  \\ over F5[a4,a6], basic support !

Note that the given curve of j-invariant 0 happens to be 36a1 but a priori any model for an arbitrary twist could have been returned. See ellfromj.

The result of ellinit is an ell structure. It contains at least the following information in its components:

a₁,a₂,a₃,a₄,a₆,b₂,b₄,b₆,b₈,c₄,c₆, Δ,j.

All are accessible via member functions. In particular, the discriminant is E.disc, and the j-invariant is E.j.

  ? E = ellinit([a4, a6]);
  ? E.disc
  %2 = -64*a4^3 - 432*a6^2
  ? E.j
  %3 = -6912*a4^3/(-4*a4^3 - 27*a6^2)

Further components contain domain-specific data, which are in general dynamic: only computed when needed, and then cached in the structure.

  ? E = ellinit([2,3], 10^60+7);  \\ E over Fp, p large
  ? ellap(E)
  time = 4,440 ms.
  %2 = -1376268269510579884904540406082
  ? ellcard(E);  \\ now instantaneous !
  time = 0 ms.
  ? ellgenerators(E);
  time = 5,965 ms.
  ? ellgenerators(E); \\ second time instantaneous
  time = 0 ms.

See the description of member functions related to elliptic curves at the beginning of this section.

The library syntax is GEN ellinit(GEN x, GEN D = NULL, long prec).

Let E be an ell structure over a number field K or ℚ_p. This function returns an integral model. If v is present, sets v = [u,0,0,0] to the corresponding change of variable: the return value is identical to that of ellchangecurve(E, v).

  ? e = ellinit([1/17,1/42]);
  ? e = ellintegralmodel(e,&v);
  ? e[1..5]
  %3 = [0, 0, 0, 15287762448, 3154568630095008]
  ? v
  %4 = [1/714, 0, 0, 0]

The library syntax is GEN ellintegralmodel(GEN E, GEN *v = NULL).

Let E an elliptic curve over a number field. Return 0 if E is not CM, otherwise return the discriminant of its endomorphism ring.

  ? E = ellinit([0,0,-5,-750,7900]);
  ? D = elliscm(E)
  %2 = -27
  ? w = quadgen(D, 'w);
  ? P = ellheegner(E)
  %4 = [10,40]
  ? Q = ellmul(E,P,w)
  %5 = [110/7-5/49*w,85/49-225/343*w]

An example over a number field:

  ? nf=nfinit(a^2-5);
  ? E = ellinit([261526980*a-584793000,-3440201839360*a+7692525148000],nf);
  ? elliscm(E)
  %3 = -20
  ? ellisomat(E)[2]
  %4 = [1,2,5,10;2,1,10,5;5,10,1,2;10,5,2,1]

The library syntax is long elliscm(GEN E).

Given E/K a number field and P in E(K) return 1 if P = [n]R for some R in E(K) and set Q to one such R; and return 0 otherwise.

  ? K = nfinit(polcyclo(11,t));
  ? E = ellinit([0,-1,1,0,0], K);
  ? P = [0,0];
  ? ellorder(E,P)
  %4 = 5
  ? ellisdivisible(E,P,5, &Q)
  %5 = 1
  ? lift(Q)
  %6 = [-t^7-t^6-t^5-t^4+1, -t^9-2*t^8-2*t^7-3*t^6-3*t^5-2*t^4-2*t^3-t^2-1]
  ? ellorder(E, Q)
  %7 = 25

We use a fast multimodular algorithm over ℚ whose complexity is essentially independent of n (polynomial in log n). Over number fields, we compute roots of division polynomials and the algebraic complexity of the underlying algorithm is in O(p⁴), where p is the largest prime divisor of n. The integer n ≥ 0 may be given as ellxn(E,n), if many points need to be tested; this provides a modest speedup over number fields but is likely to slow down the algorithm over ℚ.

The library syntax is long ellisdivisible(GEN E, GEN P, GEN n, GEN *Q = NULL).

Return 0 if the elliptic curves E and F defined over the same number field are not isomorphic, otherwise return [u,r,s,t] suitable for ellchangecurve, mapping E to F.

  ? E = ellinit([1,2]);
  ? ellisisom(E, ellinit([1,3]))
  %2 = 0
  ? F = ellchangecurve(E, [-1,1,3,2]);
  ? ellisisom(E,F)
  %4 = [1, 1, -3, -2]

  ? nf = nfinit(a^3-2); E = ellinit([a^2+1,2*a-5], nf);
  ? F = ellchangecurve(E,Mod([a, a+1, a^2, a^2+a-3], nf.pol));
  ? v = ellisisom(E,F)
  %3 = [Mod(-a, a^3 - 2), Mod(a + 1, a^3 - 2), Mod(-a^2, a^3 - 2),
        Mod(-a^2 - a + 3, a^3 - 2)]
  ? ellchangecurve(E,v) == F
  %4 = 1

The library syntax is GEN ellisisom(GEN E, GEN F).

Given an elliptic curve E, a finite subgroup G of E is given either as a generating point P (for a cyclic G) or as a polynomial whose roots vanish on the x-coordinates of the nonzero elements of G (general case and more efficient if available). This function returns the [a₁,a₂,a₃,a₄,a₆] invariants of the quotient elliptic curve E/G and (if only_image is zero (the default)) a vector of rational functions [f, g, h] such that the isogeny E → E/G is given by (x,y) ⟼ (f(x)/h(x)², g(x,y)/h(x)³).

  ? E = ellinit([0,1]);
  ? elltors(E)
  %2 = [6, [6], [[2, 3]]]
  ? ellisogeny(E, [2,3], 1)  \\ Weierstrass model for E/<P>
  %3 = [0, 0, 0, -135, -594]
  ? ellisogeny(E,[-1,0])
  %4 = [[0,0,0,-15,22], [x^3+2*x^2+4*x+3, y*x^3+3*y*x^2-2*y, x+1]]

The library syntax is GEN ellisogeny(GEN E, GEN G, long only_image, long x = -1, long y = -1) where x, y are variable numbers.

Given an isogeny of elliptic curves f:E' → E (being the result of a call to ellisogeny), apply f to g:

* if g is a point P in the domain of f, return the image f(P);

* if g:E" → E' is a compatible isogeny, return the composite isogeny f o g: E" → E.

  ? one = ffgen(101, 't)^0;
  ? E = ellinit([6, 53, 85, 32, 34] * one);
  ? P = [84, 71] * one;
  ? ellorder(E, P)
  %4 = 5
  ? [F, f] = ellisogeny(E, P);  \\ f: E->F = E/<P>
  ? ellisogenyapply(f, P)
  %6 = [0]
  ? F = ellinit(F);
  ? Q = [89, 44] * one;
  ? ellorder(F, Q)
  %9 = 2
  ? [G, g] = ellisogeny(F, Q); \\  g: F->G = F/<Q>
  ? gof = ellisogenyapply(g, f); \\ gof: E -> G

The library syntax is GEN ellisogenyapply(GEN f, GEN g).

Given an elliptic curve E defined over a number field K, computes representatives of the set of isomorphism classes of elliptic curves defined over K and K-isogenous to E, assuming it is finite (see below). For any such curve E_i, let f_i: E → E_i be a rational isogeny of minimal degree and let g_i: E_i → E be the dual isogeny; and let M be the matrix such that M_i,j is the minimal degree for an isogeny E_i → E_j.

The function returns a vector [L,M] where L is a list of triples [E_i, f_i, g_i] (flag = 0), or simply the list of E_i (flag = 1, which saves time). The curves E_i are given in [a₄,a₆] form and the first curve E₁ is isomorphic to E by f₁.

The set of isomorphism classes is finite except when E has CM over a quadratic order contained in K. In that case the function only returns the discriminant of the quadratic order.

If p is set, it must be a prime number; in this which case only isogenies of degree a power of p are considered.

Over a number field, the possible isogeny degrees are determined by Billerey's algorithm.

  ? E = ellinit("14a1");
  ? [L,M] = ellisomat(E);
  ? LE = apply(x->x[1], L)  \\ list of curves
  %3 = [[215/48,-5291/864],[-675/16,6831/32],[-8185/48,-742643/864],
       [-1705/48,-57707/864],[-13635/16,306207/32],[-131065/48,-47449331/864]]
  ? L[2][2]  \\ isogeny f2
  %4 = [x^3+3/4*x^2+19/2*x-311/12,
        1/2*x^4+(y+1)*x^3+(y-4)*x^2+(-9*y+23)*x+(55*y+55/2),x+1/3]
  ? L[2][3]  \\ dual isogeny g2
  %5 = [1/9*x^3-1/4*x^2-141/16*x+5613/64,
        -1/18*x^4+(1/27*y-1/3)*x^3+(-1/12*y+87/16)*x^2+(49/16*y-48)*x
        +(-3601/64*y+16947/512),x-3/4]
  ? apply(E->ellidentify(ellinit(E))[1][1], LE)
  %6 = ["14a1","14a4","14a3","14a2","14a6","14a5"]
  ? M
  %7 =
  [1  3  3 2  6  6]
  
  [3  1  9 6  2 18]
  
  [3  9  1 6 18  2]
  
  [2  6  6 1  3  3]
  
  [6  2 18 3  1  9]
  
  [6 18  2 3  9  1]

The library syntax is GEN ellisomat(GEN E, long p, long flag).

Gives 1 (i.e. true) if the point z is on the elliptic curve E, 0 otherwise. If E or z have imprecise coefficients, an attempt is made to take this into account, i.e. an imprecise equality is checked, not a precise one. It is allowed for z to be a vector of points in which case a vector (of the same type) is returned.

The library syntax is GEN ellisoncurve(GEN E, GEN z). Also available is int oncurve(GEN E, GEN z) which does not accept vectors of points.

Given an elliptic curve E defined over ℚ or a set of ℚ-isogenous curves as given by ellisomat, return a pair [L,M] where

* L lists the minimal models of the isomorphism classes of elliptic curves ℚ-isogenous to E (or in the set of isogenous curves),

* M is the adjacency matrix of the prime degree isogenies tree: there is an edge from E_i to E_j if there is an isogeny E_i → E_j of prime degree such that the Néron differential forms are preserved.

  ? E = ellinit("14a1");
  ? [L,M] = ellisotree(E);
  ? M
  %3 =
  [0 0 3 2 0 0]
  
  [3 0 0 0 2 0]
  
  [0 0 0 0 0 2]
  
  [0 0 0 0 0 3]
  
  [0 0 0 3 0 0]
  
  [0 0 0 0 0 0]
  ? [L2,M2] = ellisotree(ellisomat(E,2,1));
  %4 =
  [0 2]
  
  [0 0]
  ? [L3,M3] = ellisotree(ellisomat(E,3,1));
  ? M3
  %6 =
  [0 0 3]
  
  [3 0 0]
  
  [0 0 0]

Compare with the result of ellisomat.

  ? [L,M]=ellisomat(E,,1);
  ? M
  %7 =
  [1  3  3 2  6  6]
  
  [3  1  9 6  2 18]
  
  [3  9  1 6 18  2]
  
  [2  6  6 1  3  3]
  
  [6  2 18 3  1  9]
  
  [6 18  2 3  9  1]

The library syntax is GEN ellisotree(GEN E).

Return 1 if the elliptic curve E defined over a number field, ℚ_p or a finite field is supersingular at p, and 0 otherwise. If the curve is defined over ℚ or a number field, p must be explicitly given, and must be a prime number, resp. a maximal ideal; we return 1 if and only if E has supersingular good reduction at p.

Alternatively, E can be given by its j-invariant in a finite field. In this case p must be omitted.

  ?  g = ffprimroot(ffgen(7^5))
  %1 = 4*x^4+5*x^3+6*x^2+5*x+6
  ?  [g^n | n <- [1 .. 7^5 - 1], ellissupersingular(g^n)]
  %2 = [6]
  ?  j = ellsupersingularj(2^31-1)
  %3 = 1618591527*w+1497042960
  ?  ellissupersingular(j)
  %4 = 1
  
  ?  K = nfinit(y^3-2); P = idealprimedec(K, 2)[1];
  ?  E = ellinit([y,1], K);
  ?  ellissupersingular(E, P)
  %7 = 1
  ?  Q = idealprimedec(K,5)[1];
  ?  ellissupersingular(E, Q)
  %9 = 0

The library syntax is GEN ellissupersingular(GEN E, GEN p = NULL). Also available is int elljissupersingular(GEN j) where j is a j-invariant of a curve over a finite field.

Elliptic j-invariant. x must be a complex number with positive imaginary part, or convertible into a power series or a p-adic number with positive valuation.

The library syntax is GEN jell(GEN x, long prec).

Calculates the Kodaira type of the local fiber of the elliptic curve E at p. E must be an ell structure as output by ellinit, over ℚ_ℓ (p better left omitted, else equal to ℓ) over ℚ (p a rational prime) or a number field K (p a maximal ideal given by a prid structure). The result is a 4-component vector [f,kod,v,c]. Here f is the exponent of p in the arithmetic conductor of E, and kod is the Kodaira type which is coded as follows:

1 means good reduction (type I₀), 2, 3 and 4 mean types II, III and IV respectively, 4+ν with ν > 0 means type I_ν; finally the opposite values -1, -2, etc. refer to the starred types I₀^*, II^*, etc. The third component v is itself a vector [u,r,s,t] giving the coordinate changes done during the local reduction; u = 1 if and only if the given equation was already minimal at p. Finally, the last component c is the local Tamagawa number c_p.

The library syntax is GEN elllocalred(GEN E, GEN p = NULL).

Given two points P and G on the elliptic curve E/𝔽_q, returns the discrete logarithm of P in base G, i.e. the smallest nonnegative integer n such that P = [n]G. See znlog for the limitations of the underlying discrete log algorithms. If present, o represents the order of G, see Section se:DLfun; the preferred format for this parameter is [N, factor(N)], where N is the order of G.

If no o is given, assume that G generates the curve. The function also assumes that P is a multiple of G.

  ? a = ffgen(ffinit(2,8),'a);
  ? E = ellinit([a,1,0,0,1]);  \\ over F2^8
  ? x = a^3; y = ellordinate(E,x)[1];
  ? P = [x,y]; G = ellmul(E, P, 113);
  ? ord = [242, factor(242)]; \\ P generates a group of order 242. Initialize.
  ? ellorder(E, G, ord)
  %4 = 242
  ? e = elllog(E, P, G, ord)
  %5 = 15
  ? ellmul(E,G,e) == P
  %6 = 1

The library syntax is GEN elllog(GEN E, GEN P, GEN G, GEN o = NULL).

This function is deprecated, use lfun(E,s) instead.

E being an elliptic curve, given by an arbitrary model over ℚ as output by ellinit, this function computes the value of the L-series of E at the (complex) point s. This function uses an O(N^1/2) algorithm, where N is the conductor.

The optional parameter A fixes a cutoff point for the integral and is best left omitted; the result must be independent of A, up to realprecision, so this allows to check the function's accuracy.

The library syntax is GEN elllseries(GEN E, GEN s, GEN A = NULL, long prec).

Let E be an elliptic curve over Q given by ellinit or a rational isogeny class given by ellisomat. Return the Manin constant of the curve, see ellweilcurve. The algorithm is slow but unconditional. The function also accepts the output of ellisomat and returns the list of Manin constants for all the isogeny class.

  ? E = ellinit("11a3");
  ? ellmaninconstant(E)
  %2 = 5
  ? L=ellisomat(E,,1);
  ? ellmaninconstant(L)
  %4 = [5,1,1]

The library syntax is GEN ellmaninconstant(GEN E).

E being an elliptic curve defined over a number field output by ellinit, return the minimal discriminant ideal of E.

The library syntax is GEN ellminimaldisc(GEN E).

Let E be an ell structure over a number field K. This function determines whether E admits a global minimal integral model. If so, it returns it and sets v = [u,r,s,t] to the corresponding change of variable: the return value is identical to that of ellchangecurve(E, v).

Else return the (nonprincipal) Weierstrass class of E, i.e. the class of ∏ 𝔭^(v_𝔭{Δ - δ_𝔭) / 12} where Δ = E.disc is the model's discriminant and 𝔭^δ_𝔭 is the local minimal discriminant. This function requires either that E be defined over the rational field ℚ (with domain D = 1 in ellinit), in which case a global minimal model always exists, or over a number field given by a bnf structure. The Weierstrass class is given in bnfisprincipal format, i.e. in terms of the K.gen generators.

The resulting model has integral coefficients and is everywhere minimal, the coefficients a₁ and a₃ are reduced modulo 2 (in terms of the fixed integral basis K.zk) and a₂ is reduced modulo 3. Over ℚ, we further require that a₁ and a₃ be 0 or 1, that a₂ be 0 or ± 1 and that u > 0 in the change of variable: both the model and the change of variable v are then unique.

  ? e = ellinit([6,6,12,55,233]);  \\ over Q
  ? E = ellminimalmodel(e, &v);
  ? E[1..5]
  %3 = [0, 0, 0, 1, 1]
  ? v
  %4 = [2, -5, -3, 9]

  ? K = bnfinit(a^2-65);  \\ over a nonprincipal number field
  ? K.cyc
  %2 = [2]
  ? u = Mod(8+a, K.pol);
  ? E = ellinit([1,40*u+1,0,25*u^2,0], K);
  ? ellminimalmodel(E) \\ no global minimal model exists over ZK
  %6 = [1]~

The library syntax is GEN ellminimalmodel(GEN E, GEN *v = NULL).

Let E be an elliptic curve defined over ℚ, return a discriminant D such that the twist of E by D is minimal among all possible quadratic twists, i.e. if flag = 0, its minimal model has minimal discriminant, or if flag = 1, it has minimal conductor.

In the example below, we find a curve with j-invariant 3 and minimal conductor.

  ? E = ellminimalmodel(ellinit(ellfromj(3)));
  ? ellglobalred(E)[1]
  %2 = 357075
  ? D = ellminimaltwist(E,1)
  %3 = -15
  ? E2 = ellminimalmodel(elltwist(E,D));
  ? ellglobalred(E2)[1]
  %5 = 14283

In the example below, flag = 0 and flag = 1 give different results.

  ? E = ellinit([1,0]);
  ? D0 = ellminimaltwist(E,0)
  %7 = 1
  ? D1 = ellminimaltwist(E,1)
  %8 = 8
  ? E0 = ellminimalmodel(elltwist(E,D0));
  ? [E0.disc, ellglobalred(E0)[1]]
  %10 = [-64, 64]
  ? E1 = ellminimalmodel(elltwist(E,D1));
  ? [E1.disc, ellglobalred(E1)[1]]
  %12 = [-4096, 32]

The library syntax is GEN ellminimaltwist0(GEN E, long flag). Also available are GEN ellminimaltwist(E) for flag = 0, and GEN ellminimaltwistcond(E) for flag = 1.

e being an elliptic curve defined over ℚ output by ellinit, compute the modular degree of e divided by the square of the Manin constant c. It is conjectured that c = 1 for the strong Weil curve in the isogeny class (optimal quotient of J₀(N)) and this can be proven using ellweilcurve when the conductor N is moderate.

  ? E = ellinit("11a1"); \\ from Cremona table: strong Weil curve and c = 1
  ? [v,smith] = ellweilcurve(E); smith \\ proof of the above
  %2 = [[1, 1], [5, 1], [1, 1/5]]
  ? ellmoddegree(E)
  %3 = 1
  ? [ellidentify(e)[1][1] | e<-v]
  %4 = ["11a1", "11a2", "11a3"]
  ? ellmoddegree(ellinit("11a2"))
  %5 = 5
  ? ellmoddegree(ellinit("11a3"))
  %6 = 1/5

The modular degree of 11a1 is 1 (because ellweilcurve or Cremona's table prove that the Manin constant is 1 for this curve); the output of ellweilcurve also proves that the Manin constants of 11a2 and 11a3 are 1 and 5 respectively, so the actual modular degree of both 11a2 and 11a3 is 5.

The library syntax is GEN ellmoddegree(GEN e).

Given a prime N < 500, return a vector [P,t] where P(x,y) is a modular equation of level N, i.e. a bivariate polynomial with integer coefficients; t indicates the type of this equation: either canonical (t = 0) or Atkin (t = 1). This function requires the seadata package and its only use is to give access to the package contents. See polmodular for a more general and more flexible function.

Let j be the j-invariant function. The polynomial P satisfies the functional equation, P(f,j) = P(f | W_N, j | W_N) = 0 for some modular function f = f_N (hand-picked for each fixed N to minimize its size, see below), where W_N(τ) = -1 / (N τ) is the Atkin-Lehner involution. These two equations allow to compute the values of the classical modular polynomial Φ_N, such that Φ_N(j(τ), j(Nτ)) = 0, while being much smaller than the latter. More precisely, we have j(W_N(τ)) = j(N τ); the function f is invariant under Γ₀(N) and also satisfies

* for Atkin type: f | W_N = f;

* for canonical type: let s = 12/gcd(12,N-1), then f | W_N = N^s / f. In this case, f has a simple definition: f(τ) = N^s (η(N τ) / η(τ) )^{2 s}, where η is Dedekind's eta function.

The following GP function returns values of the classical modular polynomial by eliminating f_N(τ) in the above functional equation, for N ≤ 31 or N ∈ {41,47,59,71}.

  classicaleqn(N, X='X, Y='Y)=
  {
    my([P,t] = ellmodulareqn(N), Q, d);
    if (poldegree(P,'y) > 2, error("level unavailable in classicaleqn"));
    if (t == 0, \\ Canonical
      my(s = 12/gcd(12,N-1));
      Q = 'x^(N+1) * substvec(P,['x,'y],[N^s/'x,Y]);
      d = N^(s*(2*N+1)) * (-1)^(N+1);
    , \\ Atkin
      Q = subst(P,'y,Y);
      d = (X-Y)^(N+1));
    polresultant(subst(P,'y,X), Q) / d;
  }

The library syntax is GEN ellmodulareqn(long N, long x = -1, long y = -1) where x, y are variable numbers.

Computes [n]z, where z is a point on the elliptic curve E. The exponent n is in ℤ, or may be a complex quadratic integer if the curve E has complex multiplication by n (if not, an error message is issued).

  ? Ei = ellinit([1,0]); z = [0,0];
  ? ellmul(Ei, z, 10)
  %2 = [0]     \\ unsurprising: z has order 2
  ? ellmul(Ei, z, I)
  %3 = [0, 0]  \\ Ei has complex multiplication by Z[i]
  ? ellmul(Ei, z, quadgen(-4))
  %4 = [0, 0]  \\ an alternative syntax for the same query
  ? Ej  = ellinit([0,1]); z = [-1,0];
  ? ellmul(Ej, z, I)
    ***   at top-level: ellmul(Ej,z,I)
    ***                 ^ —  —  —  — --
    *** ellmul: not a complex multiplication in ellmul.
  ? ellmul(Ej, z, 1+quadgen(-3))
  %6 = [1 - w, 0]

The simple-minded algorithm for the CM case assumes that we are in characteristic 0, and that the quadratic order to which n belongs has small discriminant.

The library syntax is GEN ellmul(GEN E, GEN z, GEN n).

Opposite of the point z on elliptic curve E.

The library syntax is GEN ellneg(GEN E, GEN z).

Given an elliptic curve E/ℚ (more precisely, a model defined over ℚ of a curve) and a rational point P ∈ E(ℚ), returns the pair [R,n], where n is the least positive integer such that R := [n]P has good reduction at every prime. More precisely, its image in a minimal model is everywhere nonsingular.

  ? e = ellinit("57a1"); P = [2,-2];
  ? ellnonsingularmultiple(e, P)
  %2 = [[1, -1], 2]
  ? e = ellinit("396b2"); P = [35, -198];
  ? [R,n] = ellnonsingularmultiple(e, P);
  ? n
  %5 = 12

The library syntax is GEN ellnonsingularmultiple(GEN E, GEN P).

Gives the order of the point z on the elliptic curve E, defined over a finite field or a number field. Return (the impossible value) zero if the point has infinite order.

  ? E = ellinit([-157^2,0]);  \\ the "157-is-congruent" curve
  ? P = [2,2]; ellorder(E, P)
  %2 = 2
  ? P = ellheegner(E); ellorder(E, P) \\ infinite order
  %3 = 0
  ? K = nfinit(polcyclo(11,t)); E=ellinit("11a3", K); T = elltors(E);
  ? ellorder(E, T.gen[1])
  %5 = 25
  ? E = ellinit(ellfromj(ffgen(5^10)));
  ? ellcard(E)
  %7 = 9762580
  ? P = random(E); ellorder(E, P)
  %8 = 4881290
  ? p = 2^160+7; E = ellinit([1,2], p);
  ? N = ellcard(E)
  %9 = 1461501637330902918203686560289225285992592471152
  ? o = [N, factor(N)];
  ? for(i=1,100, ellorder(E,random(E)))
  time = 260 ms.

The parameter o, is now mostly useless, and kept for backward compatibility. If present, it represents a nonzero multiple of the order of z, see Section se:DLfun; the preferred format for this parameter is [ord, factor(ord)], where ord is the cardinality of the curve. It is no longer needed since PARI is now able to compute it over large finite fields (was restricted to small prime fields at the time this feature was introduced), and caches the result in E so that it is computed and factored only once. Modifying the last example, we see that including this extra parameter provides no improvement:

  ? o = [N, factor(N)];
  ? for(i=1,100, ellorder(E,random(E),o))
  time = 260 ms.

The library syntax is GEN ellorder(GEN E, GEN z, GEN o = NULL). The obsolete form GEN orderell(GEN e, GEN z) should no longer be used.

Gives a 0, 1 or 2-component vector containing the y-coordinates of the points of the curve E having x as x-coordinate.

The library syntax is GEN ellordinate(GEN E, GEN x, long prec).

Returns the value (or r-th derivative) on a character χ^s of ℤ_p^* of the p-adic L-function of the elliptic curve E/ℚ, twisted by D, given modulo pⁿ.

Characters. The set of continuous characters of Gal(ℚ(μ_{p ^oo} )/ ℚ) is identified to ℤ_p^* via the cyclotomic character χ with values in ℚ_p^*. Denote by τ:ℤ_p^* → ℤ_p^* the Teichmüller character, with values in the (p-1)-th roots of 1 for p ! = 2, and {-1,1} for p = 2; finally, let <χ>= χ τ^-1, with values in 1 + 2pℤ_p. In GP, the continuous character of Gal(ℚ(μ_{p ^oo} )/ ℚ) given by <χ>^s₁ τ^s₂ is represented by the pair of integers s = (s₁,s₂), with s₁ ∈ ℤ_p and s₂ mod p-1 for p > 2, (resp. mod 2 for p = 2); s may be also an integer, representing (s,s) or χ^s.

The p-adic L function. The p-adic L function L_p is defined on the set of continuous characters of Gal(ℚ(μ_{p ^oo} )/ ℚ), as ∫_{ℤ_{p}^*} χ^s d μ for a certain p-adic distribution μ on ℤ_p^*. The derivative is given by L_p^(r)(E, χ^s) = ∫_{ℤ_{p}^*} log_p^r(a) χ^s(a) dμ(a). More precisely:

* When E has good supersingular reduction, L_p takes its values in D := H¹_dR(E/ℚ) ⨂ _ℚ ℚ_p and satisfies (1-p^-1 F)^-2 L_p(E, χ⁰) = (L(E,1) / Ω).ω where F is the Frobenius, L(E,1) is the value of the complex L function at 1, ω is the Néron differential and Ω the attached period on E(ℝ). Here, χ⁰ represents the trivial character.

The function returns the components of L_p^(r)(E,χ^s) in the basis (ω, F ω).

* When E has ordinary good reduction, this method only defines the projection of L_p(E,χ^s) on the α-eigenspace, where α is the unit eigenvalue for F. This is what the function returns. We have (1- α^-1)^-2 L_p,α(E,χ⁰) = L(E,1) / Ω.

Two supersingular examples:

  ? cxL(e) = bestappr( ellL1(e) / e.omega[1] );
  
  ? e = ellinit("17a1"); p=3; \\ supersingular, a3 = 0
  ? L = ellpadicL(e,p,4);
  ? F = [0,-p;1,ellap(e,p)]; \\ Frobenius matrix in the basis (omega,F(omega))
  ? (1-p^(-1)*F)^-2 * L / cxL(e)
  %5 = [1 + O(3^5), O(3^5)]~ \\ [1,0]~
  
  ? e = ellinit("116a1"); p=3; \\ supersingular, a3 != 0~
  ? L = ellpadicL(e,p,4);
  ? F = [0,-p; 1,ellap(e,p)];
  ? (1-p^(-1)*F)^-2*L~ / cxL(e)
  %9 = [1 + O(3^4), O(3^5)]~

Good ordinary reduction:

  ? e = ellinit("17a1"); p=5; ap = ellap(e,p)
  %1 = -2 \\ ordinary
  ? L = ellpadicL(e,p,4)
  %2 = 4 + 3*5 + 4*5^2 + 2*5^3 + O(5^4)
  ? al = padicappr(x^2 - ap*x + p, ap + O(p^7))[1];
  ? (1-al^(-1))^(-2) * L / cxL(e)
  %4 = 1 + O(5^4)

Twist and Teichmüller:

  ? e = ellinit("17a1"); p=5; \\ ordinary
  \\ 2nd derivative at tau^1, twist by -7
  ? ellpadicL(e, p, 4, [0,1], 2, -7)
  %2 = 2*5^2 + 5^3 + O(5^4)

We give an example of non split multiplicative reduction (see ellpadicbsd for more examples).

  ? e=ellinit("15a1"); p=3; n=5;
  ? L = ellpadicL(e,p,n)
  %2 = 2 + 3 + 3^2 + 3^3 + 3^4 + O(3^5)
  ? (1 - ellap(e,p))^(-1) * L / cxL(e)
  %3 = 1 + O(3^5)

This function is a special case of mspadicL and it also appears as the first term of mspadicseries:

  ? e = ellinit("17a1"); p=5;
  ? L = ellpadicL(e,p,4)
  %2 = 4 + 3*5 + 4*5^2 + 2*5^3 + O(5^4)
  ? [M,phi] = msfromell(e, 1);
  ? Mp = mspadicinit(M, p, 4);
  ? mu = mspadicmoments(Mp, phi);
  ? mspadicL(mu)
  %6 = 4 + 3*5 + 4*5^2 + 2*5^3 + 2*5^4 + 5^5 + O(5^6)
  ? mspadicseries(mu)
  %7 = (4 + 3*5 + 4*5^2 + 2*5^3 + 2*5^4 + 5^5 + O(5^6))
        + (3 + 3*5 + 5^2 + 5^3 + O(5^4))*x
        + (2 + 3*5 + 5^2 + O(5^3))*x^2
        + (3 + 4*5 + 4*5^2 + O(5^3))*x^3
        + (3 + 2*5 + O(5^2))*x^4 + O(x^5)

These are more cumbersome than ellpadicL but allow to compute at different characters, or successive derivatives, or to twist by a quadratic character essentially for the cost of a single call to ellpadicL due to precomputations.

The library syntax is GEN ellpadicL(GEN E, GEN p, long n, GEN s = NULL, long r, GEN D = NULL).

Given an elliptic curve E over ℚ, its quadratic twist E_D and a prime number p, this function is a p-adic analog of the complex functions ellanalyticrank and ellbsd. It calls ellpadicL with initial accuracy pⁿ and may increase it internally; it returns a vector [r, L_p] where

* L_p is a p-adic number (resp. a pair of p-adic numbers if E has good supersingular reduction) defined modulo p^N, conjecturally equal to R_p S, where R_p is the p-adic regulator as given by ellpadicregulator (in the basis (ω, F ω)) and S is the cardinal of the Tate-Shafarevich group for the quadratic twist E_D.

* r is an upper bound for the analytic rank of the p-adic L-function attached to E_D: we know for sure that the i-th derivative of L_p(E_D,.) at χ⁰ is O(p^N) for all i < r and that its r-th derivative is nonzero; it is expected that the true analytic rank is equal to the rank of the Mordell-Weil group E_D(ℚ), plus 1 if the reduction of E_D at p is split multiplicative; if r = 0, then both the analytic rank and the Mordell-Weil rank are unconditionnally 0.

Recall that the p-adic BSD conjecture (Mazur, Tate, Teitelbaum, Bernardi, Perrin-Riou) predicts an explicit link between R_p S and (1-p^-1 F)^-2.L_p^(r)(E_D, χ⁰) / r! where r is the analytic rank of the p-adic L-function attached to E_D and F is the Frobenius on H¹_dR; see ellpadicL for definitions.

  ? E = ellinit("11a1"); p = 7; n = 5; \\ good ordinary
  ? ellpadicbsd(E, 7, 5) \\ rank 0,
  %2 = [0, 1 + O(7^5)]
  
  ? E = ellinit("91a1"); p = 7; n = 5; \\ non split multiplicative
  ? [r,Lp] = ellpadicbsd(E, p, n)
  %5 = [1, 2*7 + 6*7^2 + 3*7^3 + 7^4 + O(7^5)]
  ? R = ellpadicregulator(E, p, n, E.gen)
  %6 = 2*7 + 6*7^2 + 3*7^3 + 7^4 + 5*7^5 + O(7^6)
  ? sha = Lp/R
  %7 = 1 + O(7^4)
  
  ? E = ellinit("91b1"); p = 7; n = 5; \\ split multiplicative
  ? [r,Lp] = ellpadicbsd(E, p, n)
  %9 = [2, 2*7 + 7^2 + 5*7^3 + O(7^4)]
  ? ellpadicregulator(E, p, n, E.gen)
  %10 = 2*7 + 7^2 + 5*7^3 + 6*7^4 + 2*7^5 + O(7^6)
  ? [rC, LC] = ellanalyticrank(E);
  ? [r, rC]
  %12 = [2, 1]  \\ r = rC+1 because of split multiplicative reduction
  
  ? E = ellinit("53a1"); p = 5; n = 5; \\ supersingular
  ? [r, Lp] = ellpadicbsd(E, p, n);
  ? r
  %15 = 1
  ? Lp
  %16 = [3*5 + 2*5^2 + 2*5^5 + O(5^6), \
         5 + 3*5^2 + 4*5^3 + 2*5^4 + 5^5 + O(5^6)]
  ? R = ellpadicregulator(E, p, n, E.gen)
  %17 = [3*5 + 2*5^2 + 2*5^5 + O(5^6), 5 + 3*5^2 + 4*5^3 + 2*5^4 + O(5^5)]
  \\ expect Lp = R*#Sha, hence (conjecturally) #Sha = 1
  
  ? E = ellinit("84a1"); p = 11; n = 6; D = -443;
  ? [r,Lp] = ellpadicbsd(E, 11, 6, D) \\ Mordell-Weil rank 0, no regulator
  %19 = [0, 3 + 2*11 + O(11^6)]
  ? lift(Lp)  \\ expected cardinal for Sha is 5^2
  %20 = 25
  ? ellpadicbsd(E, 3, 12, D)  \\ at 3
  %21 = [1, 1 + 2*3 + 2*3^2 + O(3^8)]
  ? ellpadicbsd(E, 7, 8, D)   \\ and at 7
  %22 = [0, 4 + 3*7 + O(7^8)]

The library syntax is GEN ellpadicbsd(GEN E, GEN p, long n, GEN D = NULL).

If p > 2 is a prime and E is an elliptic curve on ℚ with good reduction at p, return the matrix of the Frobenius endomorphism ϕ on the crystalline module D_p(E) = ℚ_p ⨂ H¹_dR(E/ℚ) with respect to the basis of the given model (ω, η = x ω), where ω = dx/(2 y+a₁ x+a₃) is the invariant differential. The characteristic polynomial of ϕ is x² - a_p x + p. The matrix is computed to absolute p-adic precision pⁿ.

  ? E = ellinit([1,-1,1,0,0]);
  ? F = ellpadicfrobenius(E,5,3);
  ? lift(F)
  %3 =
  [120 29]
  
  [ 55  5]
  ? charpoly(F)
  %4 = x^2 + O(5^3)*x + (5 + O(5^3))
  ? ellap(E, 5)
  %5 = 0

The library syntax is GEN ellpadicfrobenius(GEN E, long p, long n).

Cyclotomic p-adic height of the rational point P on the elliptic curve E (defined over ℚ), given to n p-adic digits. If the argument Q is present, computes the value of the bilinear form (h(P+Q)-h(P-Q)) / 4.

Let D := H¹_dR(E) ⨂ _ℚ ℚ_p be the ℚ_p vector space spanned by ω (invariant differential dx/(2y+a₁x+a₃) related to the given model) and η = x ω. Then the cyclotomic p-adic height h_E associates to P ∈ E(ℚ) an element f ω + g η in D. This routine returns the vector [f, g] to n p-adic digits. If P ∈ E(ℚ) is in the kernel of reduction mod p and if its reduction at all finite places is non singular, then g = -(log_E P)², where log_E is the logarithm for the formal group of E at p.

If furthermore the model is of the form Y² = X³ + a X + b and P = (x,y), then f = log_p(denominator(x)) - 2 log_p(σ(P)) where σ(P) is given by ellsigma(E,P).

Recall (Advanced topics in the arithmetic of elliptic curves, Theorem 3.2) that the local height function over the complex numbers is of the form λ(z) = -log (|E.disc|) / 6 + Re(z η(z)) - 2 log( σ(z)). (N.B. our normalization for local and global heights is twice that of Silverman's).

   ? E = ellinit([1,-1,1,0,0]); P = [0,0];
   ? ellpadicheight(E,5,3, P)
   %2 = [3*5 + 5^2 + 2*5^3 + O(5^4), 5^2 + 4*5^4 + O(5^5)]
   ? E = ellinit("11a1"); P = [5,5]; \\ torsion point
   ? ellpadicheight(E,19,6, P)
   %4 = [0, 0]
   ? E = ellinit([0,0,1,-4,2]); P = [-2,1];
   ? ellpadicheight(E,3,3, P)
   %6 = [2*3^2 + 2*3^3 + 3^4 + O(3^5), 2*3^2 + 3^4 + O(3^5)]
   ? ellpadicheight(E,3,5, P, elladd(E,P,P))
   %7 = [3^2 + 2*3^3 + O(3^7), 3^2 + 3^3 + 2*3^4 + 3^5 + O(3^7)]

* When E has good ordinary reduction at p or non split multiplicative reduction, the "canonical" p-adic height is given by

  s2 = ellpadics2(E,p,n);
  ellpadicheight(E, p, n, P) * [1,-s2]~

Since s₂ does not depend on P, it is preferable to compute it only once:

  ? E = ellinit("5077a1"); p = 5; n = 7;  \\ rank 3
  ? s2 = ellpadics2(E,p,n);
  ? M = ellpadicheightmatrix(E,p, n, E.gen) * [1,-s2]~;
  ? matdet(M)   \\ p-adic regulator on the points in E.gen
  %4 = 5 + 5^2 + 4*5^3 + 2*5^4 + 2*5^5 + 2*5^6 + O(5^7)

* When E has split multiplicative reduction at p (Tate curve), the "canonical" p-adic height is given by

  Ep = ellinit(E[1..5], O(p^(n))); \\ E seen as a Tate curve over Qp
  [u2,u,q] = Ep.tate;
  ellpadicheight(E, p, n, P) * [1,-s2 + 1/log(q)/u2]]~

where s₂ is as above. For example,

  ? E = ellinit("91b1"); P =[-1, 3]; p = 7; n = 5;
  ? Ep = ellinit(E[1..5], O(p^(n)));
  ? s2 = ellpadics2(E,p,n);
  ? [u2,u,q] = Ep.tate;
  ? H = ellpadicheight(E,p, n, P) * [1,-s2 + 1/log(q)/u2]~
  %5 = 2*7 + 7^2 + 5*7^3 + 6*7^4 + 2*7^5 + O(7^6)

These normalizations are chosen so that p-adic BSD conjectures are easy to state, see ellpadicbsd.

The library syntax is GEN ellpadicheight0(GEN E, GEN p, long n, GEN P, GEN Q = NULL).

Q being a vector of points, this function returns the "Gram matrix" [F,G] of the cyclotomic p-adic height h_E with respect to the basis (ω, η) of D = H¹_dR(E) ⨂ _ℚ ℚ_p given to n p-adic digits. In other words, if ellpadicheight(E,p,n, Q[i],Q[j]) = [f,g], corresponding to f ω + g η in D, then F[i,j] = f and G[i,j] = g.

  ? E = ellinit([0,0,1,-7,6]); Q = [[-2,3],[-1,3]]; p = 5; n = 5;
  ? [F,G] = ellpadicheightmatrix(E,p,n,Q);
  ? lift(F)  \\ p-adic entries, integral approximation for readability
  %3 =
  [2364 3100]
  
  [3100 3119]
  
  ? G
  %4 =
  [25225 46975]
  
  [46975 61850]
  
  ? [F,G] * [1,-ellpadics2(E,p,n)]~
  %5 =
  [4 + 2*5 + 4*5^2 + 3*5^3 + O(5^5)           4*5^2 + 4*5^3 + 5^4 + O(5^5)]
  
  [    4*5^2 + 4*5^3 + 5^4 + O(5^5) 4 + 3*5 + 4*5^2 + 4*5^3 + 5^4 + O(5^5)]
  

The library syntax is GEN ellpadicheightmatrix(GEN E, GEN p, long n, GEN Q).

Let p be a prime number and let E/ℚ be a rational elliptic curve with good or bad multiplicative reduction at p. Return the Iwasawa invariants λ and μ for the p-adic L function L_p(E), twisted by (D/.) and the i-th power of the Teichmüller character τ, see ellpadicL for details about L_p(E).

Let χ be the cyclotomic character and choose γ in Gal(ℚ_p(μ_{p ^oo} )/ℚ_p) such that χ(γ) = 1+2p. Let ^{L(i), D} ∈ ℚ_p[[X]] ⨂ D_cris such that ( < χ > ^s τⁱ) (^{L(i), D}(γ-1)) = L_p(E, < χ > ^sτⁱ (D/.)).

* When E has good ordinary or bad multiplicative reduction at p. By Weierstrass's preparation theorem the series ^{L(i), D} can be written p^μ (X^λ + p G(X)) up to a p-adic unit, where G(X) ∈ ℤ_p[X]. The function returns [λ,μ].

* When E has good supersingular reduction, we define a sequence of polynomials P_n in ℚ_p[X] of degree < pⁿ (and bounded denominators), such that ^{L(i), D} = P_n ϕⁿ⁺¹ω_E - ξ_n P_n-1ϕⁿ⁺²ω_E mod ((1+X)^{p^{n}}-1) ℚ_p[X] ⨂ D_cris, where ξ_n = polcyclo(pⁿ, 1+X). Let λ_n,μ_n be the invariants of P_n. We find that

* μ_n is nonnegative and decreasing for n of given parity hence μ_2n tends to a limit μ⁺ and μ_2n+1 tends to a limit μ⁻ (both conjecturally 0).

* there exists integers λ⁺, λ⁻ in ℤ (denoted with a ~ in the reference below) such that lim_{n → oo} λ_2n + 1/(p+1) = λ⁺ and lim_{n → oo} λ_2n+1 + p/(p+1) = λ⁻. The function returns [[λ⁺, λ⁻], [μ⁺,μ⁻]].

Reference: B. Perrin-Riou, Arithmétique des courbes elliptiques à réduction supersinguli\`ere en p, Experimental Mathematics, 12, 2003, pp. 155-186.

The library syntax is GEN ellpadiclambdamu(GEN E, long p, long D, long i).

Given E defined over K = ℚ or ℚ_p and P = [x,y] on E(K) in the kernel of reduction mod p, let t(P) = -x/y be the formal group parameter; this function returns L(t) to relative p-adic precision pⁿ, where L denotes the formal logarithm (mapping the formal group of E to the additive formal group) attached to the canonical invariant differential: dL = dx/(2y + a₁x + a₃).

  ? E = ellinit([0,0,1,-4,2]); P = [-2,1];
  ? ellpadiclog(E,2,10,P)
  %2 = 2 + 2^3 + 2^8 + 2^9 + 2^10 + O(2^11)
  ? E = ellinit([17,42]);
  ? p=3; Ep = ellinit(E,p); \\ E mod p
  ? P=[114,1218]; ellorder(Ep,P) \\ the order of P on (E mod p) is 2
  %5 = 2
  ? Q = ellmul(E,P,2) \\ we need a point of the form 2*P
  %6 = [200257/7056, 90637343/592704]
  ? ellpadiclog(E,3,10,Q)
  %7 = 3 + 2*3^2 + 3^3 + 3^4 + 3^5 + 3^6 + 2*3^8 + 3^9 + 2*3^10 + O(3^11)

The library syntax is GEN ellpadiclog(GEN E, GEN p, long n, GEN P).

Let E/ℚ be an elliptic curve. Return the determinant of the Gram matrix of the vector of points S = (S₁,..., S_r) with respect to the "canonical" cyclotomic p-adic height on E, given to n (p-adic) digits.

When E has ordinary reduction at p, this is the expected Gram deteterminant in ℚ_p.

In the case of supersingular reduction of E at p, the definition requires care: the regulator R is an element of D := H¹_dR(E) ⨂ _ℚ ℚ_p, which is a two-dimensional ℚ_p-vector space spanned by ω and η = x ω (which are defined over ℚ) or equivalently but now over ℚ_p by ω and Fω where F is the Frobenius endomorphism on D as defined in ellpadicfrobenius. On D we define the cyclotomic height h_E = f ω + g η (see ellpadicheight) and a canonical alternating bilinear form [.,.]_D such that [ω, η]_D = 1.

For any ν ∈ D, we can define a height h_ν := [ h_E, ν ]_D from E(ℚ) to ℚ_p and <.,.>_ν the attached bilinear form. In particular, if h_E = f ω + gη, then h_η = [ h_E, η ]_D = f and h_ω = [ h_E, ω ]_D = - g hence h_E = h_η ω - h_ω η. Then, R is the unique element of D such that [ω,ν]_D^r-1 [R, ν]_D = det(<S_i, S_j >_ν) for all ν ∈ D not in ℚ_p ω. The ellpadicregulator function returns R in the basis (ω, Fω), which was chosen so that p-adic BSD conjectures are easy to state, see ellpadicbsd.

Note that by definition [R, η]_D = det(<S_i, S_j >_η) and [R, ω+η]_D = det(<S_i, S_j >_ω+η).

The library syntax is GEN ellpadicregulator(GEN E, GEN p, long n, GEN S).

If p > 2 is a prime and E/ℚ is an elliptic curve with ordinary good reduction at p, returns the slope of the unit eigenvector of ellpadicfrobenius(E,p,n), i.e., the action of Frobenius ϕ on the crystalline module D_p(E) = ℚ_p ⨂ H¹_dR(E/ℚ) in the basis of the given model (ω, η = x ω), where ω is the invariant differential dx/(2 y+a₁ x+a₃). In other words, η + s₂ω is an eigenvector for the unit eigenvalue of ϕ.

  ? e=ellinit([17,42]);
  ? ellpadics2(e,13,4)
  %2 = 10 + 2*13 + 6*13^3 + O(13^4)

This slope is the unique c ∈ 3^-1ℤ_p such that the odd solution σ(t) = t + O(t²) of - d((1)/(σ) (d σ)/(ω)) = (x(t) + c) ω is in tℤ_p[[t]].

It is equal to b₂/12 - E₂/12 where E₂ is the value of the Katz p-adic Eisenstein series of weight 2 on (E,ω). This is used to construct a canonical p-adic height when E has good ordinary reduction at p as follows

  s2 = ellpadics2(E,p,n);
  h(E,p,n, P, s2) = ellpadicheight(E, [p,[1,-s2]],n, P);

Since s₂ does not depend on the point P, we compute it only once.

The library syntax is GEN ellpadics2(GEN E, GEN p, long n).

Let w describe a complex period lattice (w = [w₁,w₂] or an ellinit structure). Returns normalized periods [W₁,W₂] generating the same lattice such that τ := W₁/W₂ has positive imaginary part and lies in the standard fundamental domain for SL₂(ℤ).

If flag = 1, the function returns [[W₁,W₂], [η₁,η₂]], where η₁ and η₂ are the quasi-periods attached to [W₁,W₂], satisfying η₂ W₁ - η₁ W₂ = 2 i π.

The output of this function is meant to be used as the first argument given to ellwp, ellzeta, ellsigma or elleisnum. Quasi-periods are needed by ellzeta and ellsigma only.

  ? L = ellperiods([1,I],1);
  ? [w1,w2] = L[1]; [e1,e2] = L[2];
  ? e2*w1 - e1*w2
  %3 = 6.2831853071795864769252867665590057684*I
  ? ellzeta(L, 1/2 + 2*I)
  %4 = 1.5707963... - 6.283185307...*I
  ? ellzeta([1,I], 1/2 + 2*I) \\ same but less efficient
  %4 = 1.5707963... - 6.283185307...*I

The library syntax is GEN ellperiods(GEN w, long flag, long prec).

If E/ℂ ~ ℂ/Λ is a complex elliptic curve (Λ = E.omega), computes a complex number z, well-defined modulo the lattice Λ, corresponding to the point P; i.e. such that P = [℘_Λ(z),℘'_Λ(z)] satisfies the equation y² = 4x³ - g₂ x - g₃, where g₂, g₃ are the elliptic invariants.

If E is defined over ℝ and P ∈ E(ℝ), we have more precisely, 0 ≤ Re(t) < w1 and 0 ≤ Im(t) < Im(w2), where (w1,w2) are the real and complex periods of E.

  ? E = ellinit([0,1]); P = [2,3];
  ? z = ellpointtoz(E, P)
  %2 = 3.5054552633136356529375476976257353387
  ? ellwp(E, z)
  %3 = 2.0000000000000000000000000000000000000
  ? ellztopoint(E, z) - P
  %4 = [2.548947057811923643 E-57, 7.646841173435770930 E-57]
  ? ellpointtoz(E, [0]) \\ the point at infinity
  %5 = 0

If E is defined over a general number field, the function returns the values corresponding to the various complex embeddings of the curve and of the point, in the same order as E.nf.roots:

  ? E=ellinit([-22032-15552*x,0], nfinit(x^2-2));
  ? P=[-72*x-108,0];
  ? ellisoncurve(E,P)
  %3 = 1
  ? ellpointtoz(E,P)
  %4 = [-0.52751724240790530394437835702346995884*I,
        -0.090507650025885335533571758708283389896*I]
  ? E.nf.roots
  %5 = [-1.4142135623730950488016887242096980786, \\ x-> -sqrt(2)
         1.4142135623730950488016887242096980786] \\ x->  sqrt(2)

If E/ℚ_p has multiplicative reduction, then E/ℚ_p is analytically isomorphic to ℚ_p^*/q^ℤ (Tate curve) for some p-adic integer q. The behavior is then as follows:

* If the reduction is split (E.tate[2] is a t_PADIC), we have an isomorphism φ: E(ℚ_p) ~ ℚ_p^*/q^ℤ and the function returns φ(P) ∈ ℚ_p.

* If the reduction is not split (E.tate[2] is a t_POLMOD), we only have an isomorphism φ: E(K) ~ K^*/q^ℤ over the unramified quadratic extension K/ℚ_p. In this case, the output φ(P) ∈ K is a t_POLMOD; the function is not fully implemented in this case and may fail with a "u not in ℚ_p" exception:

  ? E = ellinit([0,-1,1,0,0], O(11^5)); P = [0,0];
  ? [u2,u,q] = E.tate; type(u) \\ split multiplicative reduction
  %2 = "t_PADIC"
  ? ellmul(E, P, 5)  \\ P has order 5
  %3 = [0]
  ? z = ellpointtoz(E, [0,0])
  %4 = 3 + 11^2 + 2*11^3 + 3*11^4 + 6*11^5 + 10*11^6 + 8*11^7 + O(11^8)
  ? z^5
  %5 = 1 + O(11^9)
  ? E = ellinit(ellfromj(1/4), O(2^6)); x=1/2; y=ellordinate(E,x)[1];
  ? z = ellpointtoz(E,[x,y]); \\ t_POLMOD of t_POL with t_PADIC coeffs
  ? liftint(z) \\ lift all p-adics
  %8 = Mod(8*u + 7, u^2 + 437)
  ? x=33/4; y=ellordinate(E,x)[1]; z = ellpointtoz(E,[x,y])
  ***   at top-level: ...;y=ellordinate(E,x)[1];z=ellpointtoz(E,[x,y])
  ***                                             ^ —  —  —  —  —  — --
  *** ellpointtoz: sorry, ellpointtoz when u not in Qp is not yet implemented.

The library syntax is GEN zell(GEN E, GEN P, long prec).

Deprecated alias for ellmul.

The library syntax is GEN ellmul(GEN E, GEN z, GEN n).

If E is an elliptic curve over ℚ, attempts to compute the Mordell-Weil group attached to the curve. The output is [r₁,r₂,s,L], where r₁ ≤ rank(E) ≤ r₂, s gives informations on the Tate-Shafarevic group (see below), and L is a list of independent, non-torsion rational points on the curve. E can also be given as the output of ellrankinit(E).

If points is provided, it must be a vector of rational points on the curve, which are not computed again.

The parameter effort is a measure of the time employed to find rational points before giving up. If effort is not 0, the search is randomized, so rerunning the function might yield different or even a different number of rational points. Values up to 10 or so are reasonable but the parameter can be increased futher, with running times increasing roughly like the cube of the effort value.

  ? E = ellinit([-127^2,0]);
  ? ellrank(E)
  %2 = [1, 1, 0, []] \\ rank is 1 but no point has been found.
  ? ellrank(E,4) \\ with more effort we find a point.
  %3 = [1, 1, 0, [[38902300445163190028032/305111826865145547009,
       680061120400889506109527474197680/5329525731816164537079693913473]]]

In addition to the previous calls, the first argument E can be a pair [e,f], where e is an elliptic curve given by ellrankinit and f is a quadratic twist of e. We then look for points on f. Note that the ellrankinit initialization is independent of f, so this can speed up computations significantly!

Technical explanation. The algorithm, which computes the 2-descent and the 2-part of the Cassels pairings has an intrinsic limitation: r₁ = r₂ never holds when the Tate-Shafarevic group G has 4-torsion. Thus, in this case we cannot determine the rank precisely. The algorithm computes unconditionally three quantities:

* the rank C of the 2-Selmer group.

* the rank T of the 2-torsion subgroup.

* the (even) rank s of G[2]/2G[4]; then r₂ is defined by r₂ = C - T - s.

The following quantities are also relevant:

* the rank R of the free part of E(ℚ); it always holds that r₁ ≤ R ≤ r₂.

* the rank S of G[2] (conjecturally even); it always holds that s ≤ S and that C = T + R + S. Then r₂ = C - T - s ≥ R.

When the conductor of E is small, the BSD conjecture can be used to (conditionally) find the true rank:

  ? E=ellinit([-113^2,0]);
  ? ellrootno(E) \\ rank is even (parity conjecture)
  %2 = 1
  ? ellrank(E)
  %3 = [0, 2, 0, []] \\ rank is either 0 or 2, $2$-rank of $G$ is
  ? ellrank(E, 3) \\ try harder
  %4 = [0, 2, 0, []] \\ no luck
  ? [r,L] = ellanalyticrank(E) \\ assume BSD
  %5 = [0, 3.9465...]
  ? L / ellbsd(E) \\ analytic rank is 0, compute Sha
  %6 = 16.0000000000000000000000000000000000000

We find that the rank is 0 and the cardinal of the Tate-Shafarevich group is 16 (assuming BSD!). Moreover, since s = 0, it is isomorphic to (ℤ/4ℤ)².

When the rank is 1 and the conductor is small, ellheegner can be used to find a non-torsion point:

   ? E = ellinit([-157^2,0]);
   ? ellrank(E)
   %2 = [1, 1, 0, []] \\ rank is 1, no point found
   ? ellrank(E, 5) \\ Try harder
   time = 1,094 ms.
   %3 = [1, 1, 0, []] \\ No luck
   ? ellheegner(E) \\ use analytic method
   time = 492 ms.
   %4 = [69648970982596494254458225/166136231668185267540804, ...]

In this last example, an effort about 10 would also (with probability about 80%) find a random point, not necessarily the Heegner point, in about 5 seconds.

The library syntax is GEN ellrank(GEN E, long effort, GEN points = NULL, long prec).

If E is an elliptic curve over ℚ, initialize data to speed up further calls to ellrank.

  ? E = ellinit([0,2429469980725060,0,275130703388172136833647756388,0]);
  ? rk = ellrankinit(E);
  ? [r, R, s, P] = ellrank(rk)
  %3 = [12, 14, 0, [...]]
  ? [r, R, s, P] = ellrank(rk, 1, P) \\ more effort, using known points
  %4 = [14, 14, 0, [...]] \\ this time all points are found

The library syntax is GEN ellrankinit(GEN E, long prec).

E being an integral model of elliptic curve , return a vector containing the affine rational points on the curve of naive height less than h. If flag = 1, stop as soon as a point is found; return either an empty vector or a vector containing a single point. See hyperellratpoints for how h can be specified.

  ? E=ellinit([-25,1]);
  ? ellratpoints(E,10)
  %2 = [[-5,1],[-5,-1],[-3,7],[-3,-7],[-1,5],[-1,-5],
        [0,1],[0,-1],[5,1],[5,-1],[7,13],[7,-13]]
  ? ellratpoints(E,10,1)
  %3 = [[-5,1]]

The library syntax is GEN ellratpoints(GEN E, GEN h, long flag).

E being an ell structure over ℚ as output by ellinit, this function computes the local root number of its L-series at the place p (at the infinite place if p = 0). If p is omitted, return the global root number and in this case the curve can also be defined over a number field.

Note that the global root number is the sign of the functional equation and conjecturally is the parity of the rank of the Mordell-Weil group. The equation for E needs not be minimal at p, but if the model is already minimal the function will run faster.

The library syntax is long ellrootno(GEN E, GEN p = NULL).

Let E be an elliptic curve over ℚ and and V be a set of independent non-torsion rational points on E of infinite order that generate a subgroup G of E(ℚ) of finite index. Return a new set W of the same length that generate a subgroup H of E(ℚ) containing G and such that [E(ℚ):H] is not divisible by any prime number less than B. The running time is roughly quadratic in B.

  ? E = ellinit([0,0, 1, -7, 6]);
  ? [r,R,s,V] = ellrank(E)
  %2 = [3, 3, 0, [[-1,3], [-3,0], [11,35]]]
  ? matdet(ellheightmatrix(E, V))
  %3 = 3.7542920288254557283540759015628405708
  ? W = ellsaturation(E, V, 2) \\ index is now odd
  time = 1 ms.
  %4 = [[-1, 3], [-3, 0], [11, 35]]
  ? W = ellsaturation(E, W, 10) \\ index not divisible by p <= 10
  %5 = [[1, -1], [2, -1], [0, -3]]
  time = 2 ms.
  ? W = ellsaturation(E, V, 100) \\ looks OK now
  time = 171 ms.
  %6 = [[1, -1], [2, -1], [0, -3]]
  ? matdet(ellheightmatrix(E,V))
  %7 = 0.41714355875838396981711954461809339675
  ? lfun(E,1,3)/3! / ellbsd(E) \\ conductor is small, check assuming BSD
  %8 = 0.41714355875838396981711954461809339675

The library syntax is GEN ellsaturation(GEN E, GEN V, long B, long prec).

Let E be an ell structure as output by ellinit, defined over a finite field 𝔽_q. This low-level function computes the order of the group E(𝔽_q) using the SEA algorithm; compared to the high-level function ellcard, which includes SEA among its choice of algorithms, the tors argument allows to speed up a search for curves having almost prime order and whose quadratic twist may also have almost prime order. When tors is set to a nonzero value, the function returns 0 as soon as it detects that the order has a small prime factor not dividing tors; SEA considers modular polynomials of increasing prime degree ℓ and we return 0 as soon as we hit an ℓ (coprime to tors) dividing #E(𝔽_q):

  ? ellsea(ellinit([1,1], 2^56+3477), 1)
  %1 = 72057594135613381
  ? forprime(p=2^128,oo, q = ellcard(ellinit([1,1],p)); if(isprime(q),break))
  time = 6,571 ms.
  ? forprime(p=2^128,oo, q = ellsea(ellinit([1,1],p),1);if(isprime(q),break))
  time = 522 ms.

In particular, set tors to 1 if you want a curve with prime order, to 2 if you want to allow a cofactor which is a power of two (e.g. for Edwards's curves), etc. The early exit on bad curves yields a massive speedup compared to running the cardinal algorithm to completion.

When tors is negative, similar checks are performed for the quadratic twist of the curve.

The following function returns a curve of prime order over 𝔽_p.

  cryptocurve(p) =
  {
    while(1,
      my(E, N, j = Mod(random(p), p));
      E = ellinit(ellfromj(j));
      N = ellsea(E, 1); if (!N, continue);
      if (isprime(N), return(E));
      \\ try the quadratic twist for free
      if (isprime(2*p+2 - N), return(elltwist(E)));
    );
  }
  ? p = randomprime([2^255, 2^256]);
  ? E = cryptocurve(p); \\ insist on prime order
  %2 = 47,447ms

The same example without early abort (using ellcard(E) instead of ellsea(E, 1)) runs for about 5 minutes before finding a suitable curve.

The availability of the seadata package will speed up the computation, and is strongly recommended. The generic function ellcard should be preferred when you only want to compute the cardinal of a given curve without caring about it having almost prime order:

* If the characteristic is too small (p ≤ 7) or the field cardinality is tiny (q ≤ 523) the generic algorithm ellcard is used instead and the tors argument is ignored. (The reason for this is that SEA is not implemented for p ≤ 7 and that if q ≤ 523 it is likely to run into an infinite loop.)

* If the field cardinality is smaller than about 2⁵⁰, the generic algorithm will be faster.

* Contrary to ellcard, ellsea does not store the computed cardinality in E.

The library syntax is GEN ellsea(GEN E, long tors).

This function finds all curves in the elldata database satisfying the constraint defined by the argument N:

* if N is a character string, it selects a given curve, e.g. "11a1", or curves in the given isogeny class, e.g. "11a", or curves with given conductor, e.g. "11";

* if N is a vector of integers, it encodes the same constraints as the character string above, according to the ellconvertname correspondance, e.g. [11,0,1] for "11a1", [11,0] for "11a" and [11] for "11";

* if N is an integer, curves with conductor N are selected.

If N codes a full curve name, for instance "11a1" or [11,0,1], the output format is [N, [a₁,a₂,a₃,a₄,a₆], G] where [a₁,a₂,a₃,a₄,a₆] are the coefficients of the Weierstrass equation of the curve and G is a ℤ-basis of the free part of the Mordell-Weil group attached to the curve.

  ? ellsearch("11a3")
  %1 = ["11a3", [0, -1, 1, 0, 0], []]
  ? ellsearch([11,0,3])
  %2 = ["11a3", [0, -1, 1, 0, 0], []]

If N is not a full curve name, then the output is a vector of all matching curves in the above format:

  ? ellsearch("11a")
  %1 = [["11a1", [0, -1, 1, -10, -20], []],
        ["11a2", [0, -1, 1, -7820, -263580], []],
        ["11a3", [0, -1, 1, 0, 0], []]]
  ? ellsearch("11b")
  %2 = []

The library syntax is GEN ellsearch(GEN N). Also available is GEN ellsearchcurve(GEN N) that only accepts complete curve names (as t_STR).

Computes the value at z of the Weierstrass σ function attached to the lattice L as given by ellperiods(,1): including quasi-periods is useful, otherwise there are recomputed from scratch for each new z. σ(z, L) = z ∏_{ω ∈ L^*} (1 - (z)/(ω)) e^{(z)/(ω) + (z^{2})/(2ω²)}. It is also possible to directly input L = [ω₁,ω₂], or an elliptic curve E as given by ellinit (L = E.omega).

  ? w = ellperiods([1,I], 1);
  ? ellsigma(w, 1/2)
  %2 = 0.47494937998792065033250463632798296855
  ? E = ellinit([1,0]);
  ? ellsigma(E) \\ at 'x, implicitly at default seriesprecision
  %4 = x + 1/60*x^5 - 1/10080*x^9 - 23/259459200*x^13 + O(x^17)

If flag = 1, computes an arbitrary determination of log(σ(z)).

The library syntax is GEN ellsigma(GEN L, GEN z = NULL, long flag, long prec).

Difference of the points z1 and z2 on the elliptic curve corresponding to E.

The library syntax is GEN ellsub(GEN E, GEN z1, GEN z2).

Return a random supersingular j-invariant defined over 𝔽_p² as a t_FFELT in the variable w, if p is a prime number, or over the field of definition of p if p is a t_FFELT. The field must be of even degree. The random distribution is close to uniform except when 0 or 1728 are supersingular j-invariants, in which case they are less likely to be returned. This bias becomes negligible as p grows.

  ?  j = ellsupersingularj(1009)
  %1 = 12*w+295
  ?  ellissupersingular(j)
  %2 = 1
  ?  a = ffgen([1009,2],'a);
  ?  j = ellsupersingularj(a)
  %4 = 867*a+721
  ?  ellissupersingular(j)
  %5 = 1
  ?  E = ellinit([j]);
  ?  F = elltwist(E);
  ?  ellissupersingular(F)
  %8 = 1
  ?  ellap(E)
  %9 = 2018
  ?  ellap(F)
  %10 = -2018

The library syntax is GEN ellsupersingularj(GEN p).

The object E being an elliptic curve over a number field, returns the global Tamagawa number of the curve (including the factor at infinite places).

  ? e = ellinit([1, -1, 1, -3002, 63929]); \\ curve "90c6" from elldata
  ? elltamagawa(e)
  %2 = 288
  ? [elllocalred(e,p)[4] | p<-[2,3,5]]
  %3 = [6, 4, 6]
  ? vecprod(%)  \\ since e.disc > 0 the factor at infinity is 2
  %4 = 144
  ? ellglobalred(e)[4] \\ product without the factor at infinity
  %5 = 144

The library syntax is GEN elltamagawa(GEN E).

Computes the modular parametrization of the elliptic curve E/ℚ, where E is an ell structure as output by ellinit. This returns a two-component vector [u,v] of power series, given to n significant terms (seriesprecision by default), characterized by the following two properties. First the point (u,v) satisfies the equation of the elliptic curve. Second, let N be the conductor of E and Φ: X₀(N) → E be a modular parametrization; the pullback by Φ of the Néron differential du/(2v+a₁u+a₃) is equal to 2iπ f(z)dz, a holomorphic differential form. The variable used in the power series for u and v is x, which is implicitly understood to be equal to exp(2iπ z).

The algorithm assumes that E is a strong Weil curve and that the Manin constant is equal to 1: in fact, f(x) = ∑_{n > 0} ellak(E, n) xⁿ.

The library syntax is GEN elltaniyama(GEN E, long precdl).

Let E be an elliptic curve defined over a finite field k and m ≥ 1 be an integer. This function computes the (nonreduced) Tate pairing of the points P and Q on E, where P is an m-torsion point. More precisely, let f_m,P denote a Miller function with divisor m[P] - m[O_E]; the algorithm returns f_m,P(Q) ∈ k^*/(k^*)^m.

The library syntax is GEN elltatepairing(GEN E, GEN P, GEN Q, GEN m).

If E is an elliptic curve defined over a number field or a finite field, outputs the torsion subgroup of E as a 3-component vector [t,v1,v2], where t is the order of the torsion group, v1 gives the structure of the torsion group as a product of cyclic groups (sorted by decreasing order), and v2 gives generators for these cyclic groups. E must be an ell structure as output by ellinit.

  ?  E = ellinit([-1,0]);
  ?  elltors(E)
  %1 = [4, [2, 2], [[0, 0], [1, 0]]]

Here, the torsion subgroup is isomorphic to ℤ/2ℤ x ℤ/2ℤ, with generators [0,0] and [1,0].

The library syntax is GEN elltors(GEN E).

Let E be an elliptic curve over a base field and a point P defined over an extension field using t_POLMOD constructs. Returns the sum of the Galois conjugates of P. The field over which P is defined must be specified, even in the (silly) case of a trivial extension:

  ? E = ellinit([1,15]);  \\ y^2 = x^3 + x + 15, over Q
  ? P = Mod([a/8-1, 1/32*a^2-11/32*a-19/4], a^3-135*a-408);
  ? ellisoncurve(E,P) \\ P defined over a cubic extension
  %3 = 1
  ? elltrace(E,P)
  %4 = [2,-5]

  ? E = ellinit([-13^2, 0]);
  ? P = Mod([13,0], a^2-2); \\ defined over Q, seen over a quadratic extension
  ? elltrace(E,P) == ellmul(E,P,2)
  %3 = 1
  ? elltrace(E,[13,0]) \\ number field of definition of the point unspecified!
    ***   at top-level: elltrace(E,[13,0])
    ***                 ^ —  —  —  —  —  — 
    *** elltrace: incorrect type in elltrace (t_INT).
  ? elltrace(E,Mod([13,0],a)) \\ trivial extension
  %5 = [Mod(13, a), Mod(0, a)]
  ? P = Mod([-10*x^3+10*x-13, -16*x^3+16*x-34], x^4-x^3+2*x-1);
  ? ellisoncurve(E,P)
  %7 = 1
  ? Q = elltrace(E,P)
  %8 = [11432100241 / 375584400, 1105240264347961 / 7278825672000]
  ? ellisoncurve(E,Q)
  %9 = 1

  ? E = ellinit([2,3], 19); \\ over F_19
  ? T = a^5+a^4+15*a^3+16*a^2+3*a+1; \\ irreducible
  ? P = Mod([11*a^3+11*a^2+a+12,15*a^4+9*a^3+18*a^2+18*a+6], T);
  ? ellisoncurve(E, P)
  %4 = 1
  ? Q = elltrace(E, P)
  %5 = [Mod(1,19), Mod(14,19)]
  ? ellisoncurve(E, Q)
  %6 = 1

The library syntax is GEN elltrace(GEN E, GEN P).

Returns an ell structure (as given by ellinit) for the twist of the elliptic curve E by the quadratic extension of the coefficient ring defined by P (when P is a polynomial) or quadpoly(P) when P is an integer. If E is defined over a finite field, then P can be omitted, in which case a random model of the unique nontrivial twist is returned. If E is defined over a number field, the model should be replaced by a minimal model (if one exists).

The elliptic curve E can be given in some of the formats allowed by ellinit: an ell structure, a 5-component vector [a₁,a₂,a₃,a₄,a₆] or a 2-component vector [a₄,a₆].

Twist by discriminant -3:

  ? elltwist([0,a2,0,a4,a6], -3)[1..5]
  %1 = [0, -3*a2, 0, 9*a4, -27*a6]
  ? elltwist([a4,a6], -3)[1..5]
  %2 = [0, 0, 0, 9*a4, -27*a6]

Twist by the Artin-Schreier extension given by x²+x+T in characteristic 2:

  ? lift(elltwist([a1,a2,a3,a4,a6]*Mod(1,2), x^2+x+T)[1..5])
  %1 = [a1, a2+a1^2*T, a3, a4, a6+a3^2*T]

Twist of an elliptic curve defined over a finite field:

  ? E = elltwist([1,7]*Mod(1,19)); lift([E.a4, E.a6])
  %1 = [11, 12]

The library syntax is GEN elltwist(GEN E, GEN P = NULL).

If E' is an elliptic curve over ℚ, let L_E' be the sub-ℤ-module of Hom_{Γ_{0}(N)}(Δ₀,ℚ) attached to E' (It is given by x[3] if [M,x] = msfromell(E').)

On the other hand, if N is the conductor of E and f is the modular form for Γ₀(N) attached to E, let L_f be the lattice of the f-component of Hom_{Γ_{0}(N)}(Δ₀,ℚ) given by the elements φ such that φ({0,γ^-1 0}) ∈ ℤ for all γ ∈ Γ₀(N) (see mslattice).

Let E' run through the isomorphism classes of elliptic curves isogenous to E as given by ellisomat (and in the same order). This function returns a pair [vE,vS] where vE contains minimal models for the E' and vS contains the list of Smith invariants for the lattices L_E' in L_f. The function also accepts the output of ellisomat, i.e. the isogeny class. If the optional argument ms is present, it contains the output of msfromell(vE, 0), i.e. the new modular symbol space M of level N and a vector of triples [x⁺,x⁻, L] attached to each curve E'.

In particular, the strong Weil curve amongst the curves isogenous to E is the one whose Smith invariants are [c,c], where c is the Manin constant, conjecturally equal to 1.

  ? E = ellinit("11a3");
  ? [vE, vS] = ellweilcurve(E);
  ? [n] = [ i | i<-[1..#vS], vS[i]==[1,1] ]  \\ lattice with invariant [1,1]
  %3 = [2]
  ? ellidentify(vE[n]) \\ ... corresponds to strong Weil curve
  %4 = [["11a1", [0, -1, 1, -10, -20], []], [1, 0, 0, 0]]
  
  ? [vE, vS] = ellweilcurve(E, &ms); \\ vE,vS are as above
  ? [M, vx] = ms; msdim(M) \\ ... but ms contains more information
  %6 = 3
  ? #vx
  %7 = 3
  ? vx[1]
  %8 = [[1/25, -1/10, -1/10]~, [0, 1/2, -1/2]~, [1/25,0; -3/5,1; 2/5,-1]]
  ? forell(E, 11,11, print(msfromell(ellinit(E[1]), 1)[2]))
  [1/5, -1/2, -1/2]~
  [1, -5/2, -5/2]~
  [1/25, -1/10, -1/10]~

The last example prints the modular symbols x⁺ in M⁺ attached to the curves 11a1, 11a2 and 11a3.

The library syntax is GEN ellweilcurve(GEN E, GEN *ms = NULL).

Let E be an elliptic curve defined over a finite field and m ≥ 1 be an integer. This function computes the Weil pairing of the two m-torsion points P and Q on E, which is an alternating bilinear map. More precisely, let f_m,R denote a Miller function with divisor m[R] - m[O_E]; the algorithm returns the m-th root of unity ϵ(P,Q)^m.f_m,P(Q) / f_m,Q(P), where f(R) is the extended evaluation of f at the divisor [R] - [O_E] and ϵ(P,Q) ∈ {±1} is given by Weil reciprocity: ϵ(P,Q) = 1 if and only if P, Q, O_E are not pairwise distinct.

The library syntax is GEN ellweilpairing(GEN E, GEN P, GEN Q, GEN m).

Computes the value at z of the Weierstrass ℘ function attached to the lattice w as given by ellperiods. It is also possible to directly input w = [ω₁,ω₂], or an elliptic curve E as given by ellinit (w = E.omega).

  ? w = ellperiods([1,I]);
  ? ellwp(w, 1/2)
  %2 = 6.8751858180203728274900957798105571978
  ? E = ellinit([1,1]);
  ? ellwp(E, 1/2)
  %4 = 3.9413112427016474646048282462709151389

One can also compute the series expansion around z = 0:

  ? E = ellinit([1,0]);
  ? ellwp(E)              \\ 'x implicitly at default seriesprecision
  %5 = x^-2 - 1/5*x^2 + 1/75*x^6 - 2/4875*x^10 + O(x^14)
  ? ellwp(E, x + O(x^12)) \\ explicit precision
  %6 = x^-2 - 1/5*x^2 + 1/75*x^6 + O(x^9)

Optional flag means 0 (default): compute only ℘(z), 1: compute [℘(z),℘'(z)].

For instance, the Dickson elliptic functions sm and sn can be implemented as follows

   smcm(z) =
   { my(a, b, E = ellinit([0,-1/(4*27)])); \\ ell. invariants (g2,g3)=(0,1/27)
     [a,b] = ellwp(E, z, 1);
     [6*a / (1-3*b), (3*b+1)/(3*b-1)];
   }
   ? [s,c] = smcm(0.5);
   ? s
   %2 = 0.4898258757782682170733218609
   ? c
   %3 = 0.9591820206453842491187464098
   ? s^3+c^3
   %4 = 1.000000000000000000000000000
   ? smcm('x + O('x^11))
   %5 = [x - 1/6*x^4 + 2/63*x^7 - 13/2268*x^10 + O(x^11),
         1 - 1/3*x^3 + 1/18*x^6 - 23/2268*x^9 + O(x^10)]

The library syntax is GEN ellwp0(GEN w, GEN z = NULL, long flag, long prec). For flag = 0, we also have GEN ellwp(GEN w, GEN z, long prec), and GEN ellwpseries(GEN E, long v, long precdl) for the power series in variable v.

For any affine point P = (t,u) on the curve E, we have [n]P = (φ_n(P)ψ_n(P) : ω_n(P) : ψ_n(P)³) for some φ_n,ω_n,ψ_n in ℤ[a₁,a₂,a₃,a₄,a₆][t,u] modulo the curve equation. This function returns a pair [A,B] of polynomials in ℤ[a₁,a₂,a₃,a₄,a₆][v] such that [A(t),B(t)] = [φ_n(P),ψ_n(P)²] in the function field of E, whose quotient give the abscissa of [n]P. If P is an n-torsion point, then B(t) = 0.

  ? E = ellinit([17,42]); [t,u] = [114,1218];
  ? T = ellxn(E, 2, 'X)
  %2 = [X^4 - 34*X^2 - 336*X + 289, 4*X^3 + 68*X + 168]
  ? [a,b] = subst(T,'X,t);
  %3 = [168416137, 5934096]
  ? a / b == ellmul(E, [t,u], 2)[1]
  %4 = 1

The library syntax is GEN ellxn(GEN E, long n, long v = -1) where v is a variable number.

Computes the value at z of the Weierstrass ζ function attached to the lattice w as given by ellperiods(,1): including quasi-periods is useful, otherwise there are recomputed from scratch for each new z. ζ(z, L) = (1)/(z) + z²∑_{ω ∈ L^*} (1)/(ω²(z-ω)). It is also possible to directly input w = [ω₁,ω₂], or an elliptic curve E as given by ellinit (w = E.omega). The quasi-periods of ζ, such that ζ(z + aω₁ + bω₂) = ζ(z) + aη₁ + bη₂ for integers a and b are obtained as η_i = 2ζ(ω_i/2). Or using directly elleta.

  ? w = ellperiods([1,I],1);
  ? ellzeta(w, 1/2)
  %2 = 1.5707963267948966192313216916397514421
  ? E = ellinit([1,0]);
  ? ellzeta(E, E.omega[1]/2)
  %4 = 0.84721308479397908660649912348219163647

One can also compute the series expansion around z = 0 (the quasi-periods are useless in this case):

  ? E = ellinit([0,1]);
  ? ellzeta(E) \\ at 'x, implicitly at default seriesprecision
  %4 = x^-1 + 1/35*x^5 - 1/7007*x^11 + O(x^15)
  ? ellzeta(E, x + O(x^20)) \\ explicit precision
  %5 = x^-1 + 1/35*x^5 - 1/7007*x^11 + 1/1440257*x^17 + O(x^18)

The library syntax is GEN ellzeta(GEN w, GEN z = NULL, long prec).

E being an ell as output by ellinit, computes the coordinates [x,y] on the curve E corresponding to the complex or p-adic parameter z. Hence this is the inverse function of ellpointtoz.

* If E is defined over a p-adic field and has multiplicative reduction, then z is understood as an element on the Tate curve Q_p^* / q^ℤ.

  ? E = ellinit([0,-1,1,0,0], O(11^5));
  ? [u2,u,q] = E.tate; type(u)
  %2 = "t_PADIC" \\ split multiplicative reduction
  ? z = ellpointtoz(E, [0,0])
  %3 = 3 + 11^2 + 2*11^3 + 3*11^4 + 6*11^5 + 10*11^6 + 8*11^7 + O(11^8)
  ? ellztopoint(E,z)
  %4 = [O(11^9), O(11^9)]
  
  ? E = ellinit(ellfromj(1/4), O(2^6)); x=1/2; y=ellordinate(E,x)[1];
  ? z = ellpointtoz(E,[x,y]); \\ nonsplit: t_POLMOD with t_PADIC coefficients
  ? P = ellztopoint(E, z);
  ? P[1] \\ y coordinate is analogous, more complicated
  %8 = Mod(O(2^4)*x + (2^-1 + O(2^5)), x^2 + (1 + 2^2 + 2^4 + 2^5 + O(2^7)))

* If E is defined over the complex numbers (for instance over ℚ), z is understood as a complex number in ℂ/Λ_E. If the short Weierstrass equation is y² = 4x³ - g₂x - g₃, then [x,y] represents the Weierstrass ℘-function and its derivative. For a general Weierstrass equation we have x = ℘(z) - b₂/12, y = ℘'(z)/2 - (a₁ x + a₃)/2. If z is in the lattice defining E over ℂ, the result is the point at infinity [0].

  ? E = ellinit([0,1]); P = [2,3];
  ? z = ellpointtoz(E, P)
  %2 = 3.5054552633136356529375476976257353387
  ? ellwp(E, z)
  %3 = 2.0000000000000000000000000000000000000
  ? ellztopoint(E, z) - P
  %4 = [2.548947057811923643 E-57, 7.646841173435770930 E-57]
  ? ellztopoint(E, 0)
  %5 = [0] \\ point at infinity

The library syntax is GEN pointell(GEN E, GEN z, long prec).

Let PQ be a polynomial P, resp. a vector [P,Q] of polynomials defined over a field F of characteristic ! = 2. Returns the Igusa invariants [J₂,J₄,J₆,J₈,J₁₀] of the hyperelliptic curve C/F, defined by the equation y² = P(x), resp. y² + Q(x)*y = P(x). If k is given, only return the invariant of degree k (k must be even between 2 and 10).

  ? genus2igusa(x^5+3*x^2-4)
  %1 = [0, 9600, 20736, -23040000, 177926144]
  ? genus2igusa([x^6+x^5-x^4+3*x^3+x^2-2*x+1,x^3-x^2+x-1])
  %2 = [-788, 1958, 341220, -68178781, -662731520]
  ? genus2igusa([x^6+x^5-x^4+3*x^3+x^2-2*x+1,x^3-x^2+x-1],4)
  %3 = 1958
  ? genus2igusa(x^5+3*Mod(a,a^2-3)*x^2-4) \\  {over ℚ(sqrt{3})}
  %4 = [Mod(0, a^2 - 3), Mod(9600*a, a^2 - 3), Mod(186624, a^2 - 3),
        Mod(-69120000, a^2 - 3), Mod(-241864704*a + 204800000, a^2 - 3)]
  ? a = ffgen(3^4,'a); \\  {over 𝔽3^4 = 𝔽3[a]}
  ? genus2igusa(x^6+a*x^5-a*x^4+2*x^3+a*x+a+1)
  %6 = [2*a^2, a^3 + a^2 + a + 1, a^2 + a + 2, 2*a^3 + 2*a^2 + a + 1,
        2*a^2 + 2]
  ? a = ffgen(2^4,'a); \\  {𝔽2^4 = 𝔽2[a]}
  ? genus2igusa(x^6+a*x^5+a*x^4+a*x+a+1) \\ doesn't work in characteristic 2
    ***   at top-level: genus2igusa(x^6+a*x^5+a*x^4+a*x+a+1)
    ***                 ^ —  —  —  —  —  —  —  —  —  —  —  — 
    *** genus2igusa: impossible inverse in FF_mul2n: 2.

The library syntax is GEN genus2igusa(GEN PQ, long k).

Let PQ be a polynomial P, resp. a vector [P,Q] of polynomials, with rational coefficients. Determines the reduction at p > 2 of the (proper, smooth) genus 2 curve C/ℚ, defined by the hyperelliptic equation y² = P(x), resp. y² + Q(x)*y = P(x). (The special fiber X_p of the minimal regular model X of C over ℤ.)

If p is omitted, determines the reduction type for all (odd) prime divisors of the discriminant.

This function was rewritten from an implementation of Liu's algorithm by Cohen and Liu (1994), genus2reduction-0.3, see https://www.math.u-bordeaux.fr/~liu/G2R/.

CAVEAT. The function interface may change: for the time being, it returns [N,FaN, [P_m, Q_m], V] where N is either the local conductor at p or the global conductor, FaN is its factorization, y² +Q_m y = P_m defines a minimal model over ℤ and V describes the reduction type at the various considered p. Unfortunately, the program is not complete for p = 2, and we may return the odd part of the conductor only: this is the case if the factorization includes the (impossible) term 2^-1; if the factorization contains another power of 2, then this is the exact local conductor at 2 and N is the global conductor.

  ? default(debuglevel, 1);
  ? genus2red(x^6 + 3*x^3 + 63, 3)
  (potential) stable reduction: [1, []]
  reduction at p: [III{9}] page 184, [3, 3], f = 10
  %1 = [59049, Mat([3, 10]), x^6 + 3*x^3 + 63, [3, [1, []],
         ["[III{9}] page 184", [3, 3]]]]
  ? [N, FaN, T, V] = genus2red(x^3-x^2-1, x^2-x);  \\ X1(13), global reduction
  p = 13
  (potential) stable reduction: [5, [Mod(0, 13), Mod(0, 13)]]
  reduction at p: [I{0}-II-0] page 159, [], f = 2
  ? N
  %3 = 169
  ? FaN
  %4 = Mat([13, 2])   \\ in particular, good reduction at 2 !
  ? T
  %5 = x^6 + 58*x^5 + 1401*x^4 + 18038*x^3 + 130546*x^2 + 503516*x + 808561
  ? V
  %6 = [[13, [5, [Mod(0, 13), Mod(0, 13)]], ["[I{0}-II-0] page 159", []]]]

We now first describe the format of the vector V = V_p in the case where p was specified (local reduction at p): it is a triple [p, stable, red]. The component stable = [type, vecj] contains information about the stable reduction after a field extension; depending on types, the stable reduction is

* 1: smooth (i.e. the curve has potentially good reduction). The Jacobian J(C) has potentially good reduction.

* 2: an elliptic curve E with an ordinary double point; vecj contains j mod p, the modular invariant of E. The (potential) semi-abelian reduction of J(C) is the extension of an elliptic curve (with modular invariant j mod p) by a torus.

* 3: a projective line with two ordinary double points. The Jacobian J(C) has potentially multiplicative reduction.

* 4: the union of two projective lines crossing transversally at three points. The Jacobian J(C) has potentially multiplicative reduction.

* 5: the union of two elliptic curves E₁ and E₂ intersecting transversally at one point; vecj contains their modular invariants j₁ and j₂, which may live in a quadratic extension of 𝔽_p and need not be distinct. The Jacobian J(C) has potentially good reduction, isomorphic to the product of the reductions of E₁ and E₂.

* 6: the union of an elliptic curve E and a projective line which has an ordinary double point, and these two components intersect transversally at one point; vecj contains j mod p, the modular invariant of E. The (potential) semi-abelian reduction of J(C) is the extension of an elliptic curve (with modular invariant j mod p) by a torus.

* 7: as in type 6, but the two components are both singular. The Jacobian J(C) has potentially multiplicative reduction.

The component red = [NUtype, neron] contains two data concerning the reduction at p without any ramified field extension.

The NUtype is a t_STR describing the reduction at p of C, following Namikawa-Ueno, The complete classification of fibers in pencils of curves of genus two, Manuscripta Math., vol. 9, (1973), pages 143-186. The reduction symbol is followed by the corresponding page number or page range in this article.

The second datum neron is the group of connected components (over an algebraic closure of 𝔽_p) of the Néron model of J(C), given as a finite abelian group (vector of elementary divisors).

If p = 2, the red component may be omitted altogether (and replaced by [], in the case where the program could not compute it. When p was not specified, V is the vector of all V_p, for all considered p.

Notes about Namikawa-Ueno types.

* A lower index is denoted between braces: for instance, [I{2}-II-5] means [I2-II-5].

* If K and K' are Kodaira symbols for singular fibers of elliptic curves, then [K-K'-m] and [K'-K-m] are the same.

We define a total ordering on Kodaira symbol by fixing I < I* < II < II*,.... If the reduction type is the same, we order by the number of components, e.g. I₂ < I₄, etc. Then we normalize our output so that K ≤ K'.

* [K-K'--1] is [K-K'-α] in the notation of Namikawa-Ueno.

* The figure [2I0-m] in Namikawa-Ueno, page 159, must be denoted by [2I0-(m+1)].

The library syntax is GEN genus2red(GEN PQ, GEN p = NULL).

C being a nonsingular hyperelliptic model of a curve, apply the change of coordinate given by m = [e, [a,b;c,d], H].

If (x,y) is a point on the new model, the corresponding point (X,Y) on C is given by

X = (a*x + b) / (c*x + d), Y = e (y + H(x)) / (c*x + d)^g+1.

C can be given either by a squarefree polynomial P such that C: y² = P(x) or by a vector [P,Q] such that C: y² + Q(x) y = P(x) and Q²+4 P is squarefree.

The library syntax is GEN hyperellchangecurve(GEN C, GEN m).

X being a nonsingular hyperelliptic curve defined over a finite field, return the characteristic polynomial of the Frobenius automorphism. X can be given either by a squarefree polynomial P such that X: y² = P(x) or by a vector [P,Q] such that X: y² + Q(x) y = P(x) and Q²+4 P is squarefree.

The library syntax is GEN hyperellcharpoly(GEN X).

X being a nonsingular hyperelliptic model of a curve, defined over a field of characteristic distinct from 2, returns its discriminant. X can be given either by a squarefree polynomial P such that X has equation y² = P(x) or by a vector [P,Q] such that X has equation y² + Q(x) y = P(x) and Q²+4 P is squarefree.

  ? hyperelldisc([x^3,1])
  %1 = -27
  ? hyperelldisc(x^5+1)
  %2 = 800000

The library syntax is GEN hyperelldisc(GEN X).

X being a nonsingular hyperelliptic model of a curve, test whether the point p is on the curve.

X can be given either by a squarefree polynomial P such that X: y² = P(x) or by a vector [P,Q] such that X: y² + Q(x) y = P(x) and Q²+4 P is squarefree.

  ? W = [2*x^6+3*x^5+x^4+x^3-x,x^3+1]; p = [px, py] = [1/3,-14/27];
  ? hyperellisoncurve(W, p)
  %2 = 1
  ? [Px,Qx]=subst(W,x,px); py^2+py*Qx == Px
  %3 = 1

The library syntax is GEN hyperellisoncurve(GEN X, GEN p).

C being a nonsingular integral hyperelliptic model of a curve, return the minimal discriminant of an integral model of C. If pr is given, it must be a list of primes and the discriminant is then only garanteed minimal at the elements of pr. C can be given either by a squarefree polynomial P such that C: y² = P(x) or by a vector [P,Q] such that C: y² + Q(x) y = P(x) and Q²+4 P is squarefree.

  ? W = [x^6+216*x^3+324,0];
  ? D = hyperelldisc(W)
  %2 = 1828422898924853919744000
  ? M = hyperellminimaldisc(W)
  %4 = 29530050606000

The library syntax is GEN hyperellminimaldisc(GEN C, GEN pr = NULL).

C being a nonsingular integral hyperelliptic model of a curve, return an integral model of C with minimal discriminant. If pr is given, it must be a list of primes and the model is then only garanteed minimal at the elements of pr. If present, m is set to the mapping from the original model to the new one: a three-component vector [e,[a,b;c,d],H] such that if (x,y) is a point on W, the corresponding point on C is given by x_C = (a*x+b)/(c*x+d) y_C = (e*y+H(x))/(c*x+d)^g+1 where g is the genus. C can be given either by a squarefree polynomial P such that C: y² = P(x) or by a vector [P,Q] such that C: y² + Q(x) y = P(x) and Q²+4 P is squarefree.

  ? W = [x^6+216*x^3+324,0];
  ? D = hyperelldisc(W)
  %2 = 1828422898924853919744000
  ? Wn = hyperellminimalmodel(W,&M)
  %3 = [2*x^6+18*x^3+1,x^3];
  ? M
  %4 = [18, [3, 0; 0, 1], 9*x^3]
  ? hyperelldisc(Wn)
  %5 = 29530050606000
  ? hyperellchangecurve(W, M)
  %6 = [2*x^6+18*x^3+1,x^3]

The library syntax is GEN hyperellminimalmodel(GEN C, GEN *m = NULL, GEN pr = NULL).

Gives a 0, 1 or 2-component vector containing the y-coordinates of the points of the curve H having x as x-coordinate.

The library syntax is GEN hyperellordinate(GEN H, GEN x).

Let X be the curve defined by y² = Q(x), where Q is a polynomial of degree d over ℚ and q ≥ d is a prime such that X has good reduction at q. Return the matrix of the Frobenius endomorphism ϕ on the crystalline module D_p(X) = ℚ_p ⨂ H¹_dR(X/ℚ) with respect to the basis of the given model (ω, x ω,...,x^g-1 ω), where ω = dx/(2 y) is the invariant differential, where g is the genus of X (either d = 2 g+1 or d = 2 g+2). The characteristic polynomial of ϕ is the numerator of the zeta-function of the reduction of the curve X modulo q. The matrix is computed to absolute q-adic precision qⁿ.

Alternatively, q may be of the form [T,p] where p is a prime, T is a polynomial with integral coefficients whose projection to 𝔽_p[t] is irreducible, X is defined over K = ℚ[t]/(T) and has good reduction to the finite field 𝔽_q = 𝔽_p[t]/(T). The matrix of ϕ on D_q(X) = ℚ_q ⨂ H¹_dR(X/K) is computed to absolute p-adic precision pⁿ.

  ? M=hyperellpadicfrobenius(x^5+'a*x+1,['a^2+1,3],10);
  ? liftall(M)
  [48107*a + 38874  9222*a + 54290  41941*a + 8931 39672*a + 28651]
  
  [ 21458*a + 4763  3652*a + 22205 31111*a + 42559 39834*a + 40207]
  
  [ 13329*a + 4140 45270*a + 25803  1377*a + 32931 55980*a + 21267]
  
  [15086*a + 26714  33424*a + 4898 41830*a + 48013  5913*a + 24088]
  ? centerlift(simplify(liftpol(charpoly(M))))
  %8 = x^4+4*x^2+81
  ? hyperellcharpoly((x^5+Mod(a,a^2+1)*x+1)*Mod(1,3))
  %9 = x^4+4*x^2+81

The library syntax is GEN hyperellpadicfrobenius0(GEN Q, GEN q, long n). The functions GEN hyperellpadicfrobenius(GEN H, ulong p, long n) and GEN nfhyperellpadicfrobenius(GEN H, GEN T, ulong p, long n) are also available.

X being a nonsingular hyperelliptic curve given by an rational model, return a vector containing the affine rational points on the curve of naive height less than h. If flag = 1, stop as soon as a point is found; return either an empty vector or a vector containing a single point.

X is given either by a squarefree polynomial P such that X: y² = P(x) or by a vector [P,Q] such that X: y²+Q(x) y = P(x) and Q²+4 P is squarefree.

The parameter h can be

* an integer H: find the points [n/d,y] whose abscissas x = n/d have naive height ( = max(|n|, d)) less than H;

* a vector [N,D] with D ≤ N: find the points [n/d,y] with |n| ≤ N, d ≤ D.

* a vector [N,[D₁,D₂]] with D₁ < D₂ ≤ N find the points [n/d,y] with |n| ≤ N and D₁ ≤ d ≤ D₂.

The library syntax is GEN hyperellratpoints(GEN X, GEN h, long flag).

Let C be a nonsingular integral hyperelliptic model of a curve of positive genus g > 0. Return an integral model of C with the same discriminant but small coefficients, using Cremona-Stoll reduction.

The optional argument m is set to the mapping from the original model to the new one, given by a three-component vector [1,[a,b;c,d],H] such that a*d-b*c = 1 and if (x,y) is a point on W, the corresponding point (X,Y) on C is given by

X = (a*x + b) / (c*x + d), Y = (y + H(x)) / (c*x + d)^g+1.

C can be given either by a squarefree polynomial P such that C: y² = P(x) or by a vector [P,Q] such that C: y² + Q(x) y = P(x) and Q²+4 P is squarefree.

  ? P = 1001*x^4 + 3704*x^3 + 5136*x^2 + 3163*x + 730;
  ? hyperellred(P, &m)
  %2 = [x^3 + 1, 0]
  ? hyperellchangecurve(P, m)
  %3 = [x^3 + 1, 0]

The library syntax is GEN hyperellred(GEN C, GEN *m = NULL).

Also available is GEN ZX_hyperellred(GEN P, GEN *M) where C: y² = P(x) and *M is set to [a,b;c,d]