From d431329f0000bf36cdfd6af0f45b93ca99160bf3 Mon Sep 17 00:00:00 2001 From: Kevin Lui Date: Wed, 14 Sep 2016 15:35:29 -0700 Subject: [PATCH] Copied over files from Github repo. Made small docstring changes. --- src/sage/modular/abvar/abvar.py | 376 ++++++++++++++++- .../modular/abvar/abvar_ambient_jacobian.py | 30 +- src/sage/modular/abvar/lseries.py | 152 ++++++- src/sage/modular/abvar/torsion_subgroup.py | 389 ++++++++++++++---- 4 files changed, 849 insertions(+), 98 deletions(-) diff --git a/src/sage/modular/abvar/abvar.py b/src/sage/modular/abvar/abvar.py index 54c08161f37..0d7b68e606e 100644 --- a/src/sage/modular/abvar/abvar.py +++ b/src/sage/modular/abvar/abvar.py @@ -41,15 +41,25 @@ from sage.rings.all import ZZ, QQ, QQbar, Integer from sage.arith.all import LCM, divisors, prime_range, next_prime from sage.rings.ring import is_Ring +from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing +from sage.rings.polynomial.polynomial_element import Polynomial +from sage.rings.infinity import infinity +from sage.rings.fraction_field import FractionField from sage.modules.free_module import is_FreeModule from sage.modular.arithgroup.all import is_CongruenceSubgroup, is_Gamma0, is_Gamma1, is_GammaH from sage.modular.modsym.all import ModularSymbols from sage.modular.modsym.space import ModularSymbolsSpace +from sage.modular.modform.constructor import Newform from sage.matrix.all import matrix, block_diagonal_matrix, identity_matrix from sage.modules.all import vector from sage.groups.all import AbelianGroup from sage.databases.cremona import cremona_letter_code from sage.misc.all import prod +from sage.arith.misc import is_prime +from sage.databases.cremona import CremonaDatabase +from sage.schemes.elliptic_curves.constructor import EllipticCurve +from sage.sets.primes import Primes +from sage.functions.other import real from copy import copy @@ -174,6 +184,50 @@ def groups(self): """ return self.__groups + def is_J0(self): + """ + Returns whether of not self is of the form J0(N). + + OUTPUT: bool + + EXAMPLES:: + + sage: J0(23).is_J0() + True + sage: J1(11).is_J0() + False + sage: (J0(23) * J1(11)).is_J0() + False + sage: J0(37)[0].is_J0() + False + sage: (J0(23) * J0(21)).is_J0() + False + """ + return len(self.groups()) == 1 and is_Gamma0(self.groups()[0]) \ + and self.is_ambient() + + def is_J1(self): + """ + Returns whether of not self is of the form J1(N). + + OUTPUT: bool + + EXAMPLES:: + + sage: J1(23).is_J1() + True + sage: J0(23).is_J1() + False + sage: (J1(11) * J1(13)).is_J1() + False + sage: (J1(11) * J0(13)).is_J1() + False + sage: J1(23)[0].is_J1() + False + """ + return len(self.groups()) == 1 and is_Gamma1(self.groups()[0]) \ + and self.is_ambient() + ############################################################################# # lattice() *must* be defined by every derived class!!!! def lattice(self): @@ -446,6 +500,61 @@ def label(self): degen = str(self.degen_t()).replace(' ','') return '%s%s'%(self.newform_label(), degen) + def newform(self, names=None): + """ + Return the newform `f` such that this abelian variety is isogenous to + the newform abelian variety `A_f`. If this abelian variety is not + simple, raise a ValueError. + + INPUT: + + - ``names`` -- (default: None) If the newform has coefficients in a + number field, then a generator name must be specified. + + OUTPUT: A newform `f` so that self is isogenous to `A_f`. + + EXAMPLES:: + + sage: J0(11).newform() + q - 2*q^2 - q^3 + 2*q^4 + q^5 + O(q^6) + + sage: f = J0(23).newform(names='a') + sage: AbelianVariety(f) == J0(23) + True + + sage: J = J0(33) + sage: [s.newform('a') for s in J.decomposition()] + [q - 2*q^2 - q^3 + 2*q^4 + q^5 + O(q^6), + q - 2*q^2 - q^3 + 2*q^4 + q^5 + O(q^6), + q + q^2 - q^3 - q^4 - 2*q^5 + O(q^6)] + + The following fails since `J_0(33)` is not simple:: + + sage: J0(33).newform() + Traceback (most recent call last): + ... + ValueError: self must be simple + """ + return Newform(self.newform_label(), names=names) + + def newform_decomposition(self, names=None): + """ + Return the newforms of the simple subvarieties in the decomposition of + self as a product of simple subvarieties, up to isogeny. + + OUTPUT: + + - an array of newforms + + EXAMPLES:: + + sage: J = J1(11) * J0(23) + sage: J.newform_decomposition('a') + [q - 2*q^2 - q^3 + 2*q^4 + q^5 + O(q^6), + q + a0*q^2 + (-2*a0 - 1)*q^3 + (-a0 - 1)*q^4 + 2*a0*q^5 + O(q^6)] + """ + return [S.newform(names=names) for S in self.decomposition()] + def newform_label(self): """ Return the label [level][isogeny class][group] of the newform @@ -478,6 +587,78 @@ def newform_label(self): group = 'GH%s'%(str(G._generators_for_H()).replace(' ','')) return '%s%s%s'%(N, cremona_letter_code(self.isogeny_number()), group) + def elliptic_curve(self): + """ + Return an elliptic curve isogenous to self. If self is not dimension 1 + with rational base ring, raise a ValueError. + + The elliptic curve is found by looking it up in the CremonaDatabase. + The CremonaDatabase contains all curves up to some large conductor. If + a curve is not found in the CremonaDatabase, a RuntimeError will be + raised. In practice, only the most committed users will see this + RuntimeError. + + OUTPUT: an elliptic curve isogenous to self. + + EXAMPLES:: + + sage: J = J0(11) + sage: J.elliptic_curve() + Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field + + sage: J = J0(49) + sage: J.elliptic_curve() + Elliptic Curve defined by y^2 + x*y = x^3 - x^2 - 2*x - 1 over Rational Field + + sage: A = J0(37)[1] + sage: E = A.elliptic_curve() + sage: A.lseries()(1) + 0.725681061936153 + sage: E.lseries()(1) + 0.725681061936153 + + Elliptic curves are of dimension 1. :: + + sage: J = J0(23) + sage: J.elliptic_curve() + Traceback (most recent call last): + ... + ValueError: self must be of dimension 1 + + This is only implemented for curves over QQ. :: + + sage: J = J0(11).change_ring(CC) + sage: J.elliptic_curve() + Traceback (most recent call last): + ... + ValueError: base ring must be QQ + """ + + if self.dimension() > 1: + raise ValueError("self must be of dimension 1") + if self.base_ring() != QQ: + raise ValueError("base ring must be QQ") + + f = self.newform('a') + N = f.level() + + c = CremonaDatabase() + if N > c.largest_conductor(): + raise RuntimeError("Elliptic curve not found" + + " in installed database") + + isogeny_classes = c.isogeny_classes(N) + curves = [EllipticCurve(x[0][0]) for x in isogeny_classes] + + if len(curves) == 1: + return curves[0] + for p in Primes(): + for E in curves: + if E.ap(p) != f.coefficient(p): + curves.remove(E) + if len(curves) == 1: + return curves[0] + def _isogeny_to_newform_abelian_variety(self): r""" Return an isogeny from self to an abelian variety `A_f` @@ -1539,6 +1720,44 @@ def dimension(self): """ return self.lattice().rank() // 2 + def conductor(self): + """ + Return the conductor of this abelian variety. + + EXAMPLES:: + + sage: A = J0(23) + sage: A.conductor().factor() + 23^2 + + sage: A = J1(25) + sage: A.conductor().factor() + 5^24 + + sage: A = J0(11^2); A.decomposition() + [ + Simple abelian subvariety 11a(1,121) of dimension 1 of J0(121), + Simple abelian subvariety 11a(11,121) of dimension 1 of J0(121), + Simple abelian subvariety 121a(1,121) of dimension 1 of J0(121), + Simple abelian subvariety 121b(1,121) of dimension 1 of J0(121), + Simple abelian subvariety 121c(1,121) of dimension 1 of J0(121), + Simple abelian subvariety 121d(1,121) of dimension 1 of J0(121) + ] + sage: A.conductor().factor() + 11^10 + + sage: A = J0(33)[0]; A + Simple abelian subvariety 11a(1,33) of dimension 1 of J0(33) + sage: A.conductor() + 11 + sage: A.elliptic_curve().conductor() + 11 + """ + if not self.base_ring() == QQ: + raise ValueError("base ring must be QQ") + return prod(f.level() ** f.base_ring().degree() + for f in self.newform_decomposition('a')) + def rank(self): """ Return the rank of the underlying lattice of self. @@ -1773,9 +1992,9 @@ def newform_level(self, none_if_not_known=False): except AttributeError: if none_if_not_known: return None - N = [A.newform_level() for A in self.decomposition()] - level = LCM([z[0] for z in N]) - groups = sorted(set([z[1] for z in N])) + level = LCM([f.level() for f in self.newform_decomposition('a')]) + groups = sorted(set([f.group() for f in + self.newform_decomposition('a')])) if len(groups) == 1: groups = groups[0] self.__newform_level = level, groups @@ -1993,6 +2212,152 @@ def _ambient_hecke_matrix_on_modular_symbols(self, n): self.__ambient_hecke_matrix_on_modular_symbols[n] = T return T + def number_of_rational_points(self): + """ + Return the number of rational points of this modular abelian variety. + + It is not always possible to compute the order of the torsion + subgroup. The BSD conjecture is assumed to compute the algebraic rank. + + OUTPUT: a positive integer or infinity + + EXAMPLES:: + + sage: from sage_modabvar import J0, J1 + sage: J0(23).number_of_rational_points() + 11 + sage: J0(29).number_of_rational_points() + 7 + sage: J0(37).number_of_rational_points() + +Infinity + + sage: J0(12); J0(12).number_of_rational_points() + Abelian variety J0(12) of dimension 0 + 1 + + sage: J1(17).number_of_rational_points() + 584 + + sage: J1(16).number_of_rational_points() + Traceback (most recent call last): + ... + RuntimeError: Unable to compute order of torsion subgroup (it is in [1, 2, 4, 5, 10, 20]) + """ + # Check easy dimension zero case + if self.dimension() == 0: + return ZZ(1) + + positive_rank = self.lseries().vanishes_at_1() + + if positive_rank: + return infinity + else: + return self.rational_torsion_subgroup().order() + + def frobenius_polynomial(self, p, var='x'): + """ + Computes the frobenius polynomial at `p`. + + INPUT: + + - ``p`` -- prime number + + OUTPUT: + + a monic integral polynomial + + EXAMPLES:: + + sage: f = Newform('39b','a') + sage: A=AbelianVariety(f) + sage: A.frobenius_polynomial(5) + x^4 + 2*x^2 + 25 + + sage: J=J0(23) + sage: J.frobenius_polynomial(997) + x^4 + 20*x^3 + 1374*x^2 + 19940*x + 994009 + + sage: J = J0(33) + sage: J.frobenius_polynomial(7) + x^6 + 9*x^4 - 16*x^3 + 63*x^2 + 343 + + sage: J = J0(19) + sage: J.frobenius_polynomial(3, var='y') + y^2 + 2*y + 3 + + sage: J = J0(3); J + Abelian variety J0(3) of dimension 0 + sage: J.frobenius_polynomial(11) + 1 + + sage: A = J1(27)[1]; A + Simple abelian subvariety 27bG1(1,27) of dimension 12 of J1(27) + sage: A.frobenius_polynomial(11) + x^24 - 3*x^23 - 15*x^22 + 126*x^21 - 201*x^20 - 1488*x^19 + 7145*x^18 - 1530*x^17 - 61974*x^16 + 202716*x^15 - 19692*x^14 - 1304451*x^13 + 4526883*x^12 - 14348961*x^11 - 2382732*x^10 + 269814996*x^9 - 907361334*x^8 - 246408030*x^7 + 12657803345*x^6 - 28996910448*x^5 - 43086135081*x^4 + 297101409066*x^3 - 389061369015*x^2 - 855935011833*x + 3138428376721 + + sage: J = J1(33) + sage: J.frobenius_polynomial(11) + Traceback (most recent call last): + ... + ValueError: p must not divide the level of self + sage: J.frobenius_polynomial(4) + Traceback (most recent call last): + ... + ValueError: p must be prime + """ + if self.dimension() == 0: + return ZZ(1) + if self.level() % p == 0: + raise ValueError("p must not divide the level of self") + if not is_prime(p): + raise ValueError("p must be prime") + if not self.is_simple(): + from .constructor import AbelianVariety + decomp = [AbelianVariety(f) for f in + self.newform_decomposition('a')] + return prod((s.frobenius_polynomial(p) for s in + decomp)) + f = self.newform('a') + Kf = f.base_ring() + eps = f.character() + Qe = eps.base_ring() + d = Kf.degree() / Qe.degree() + + ZZpoly = PolynomialRing(ZZ, var) + z = ZZpoly.gens()[0] + + if Qe == QQ: + Gp = self.hecke_polynomial(p) + ans = ZZpoly(z**d * Gp(z+p/z)) + return ans + + # relativize number fields to compute charpoly of + # left multiplication of ap on Kf as a Qe-vector + # space. + Lf = Kf.relativize(Qe.gen(), 'a') + to_Lf = Lf.structure()[1] + + # Re is just Qe with a different name + Re = Lf.base_field() + to_Re = Qe.hom(Re.gens()[0], Re) + + name = Kf._names[0] + ap = to_Lf(f.modular_symbols(1).eigenvalue(p, name)) + Gp = ap.charpoly(var='x') + + S = PolynomialRing(Re, 'x') + x = S.gens()[0] + h = S(x**d * Gp(x=x+to_Re(eps(p))*p/x)) + + # take Qe norm + R = PolynomialRing(QQ, ['x', 'y']) + x, y = R.gens() + g = Qe.defining_polynomial() + H = sum(h[i].lift() * y**i for i in range(h.degree()+1)) + ans = g.resultant(H) + + return ZZpoly(ans.univariate_polynomial()) + ############################################################################### # Rational and Integral Homology ############################################################################### @@ -2351,7 +2716,8 @@ def rational_torsion_subgroup(self): sage: A = J.new_subvariety() sage: A Abelian subvariety of dimension 1 of J0(33) - sage: t = A.rational_torsion_subgroup() + sage: t = A.rational_torsion_subgroup(); t + Torsion subgroup of Abelian subvariety of dimension 1 of J0(33) sage: t.multiple_of_order() 4 sage: t.divisor_of_order() @@ -2360,8 +2726,6 @@ def rational_torsion_subgroup(self): 4 sage: t.gens() [[(1/2, 0, 0, -1/2, 0, 0)], [(0, 0, 1/2, 0, 1/2, -1/2)]] - sage: t - Torsion subgroup of Abelian subvariety of dimension 1 of J0(33) """ try: return self.__rational_torsion_subgroup diff --git a/src/sage/modular/abvar/abvar_ambient_jacobian.py b/src/sage/modular/abvar/abvar_ambient_jacobian.py index 796900caa86..02cd66e61ae 100644 --- a/src/sage/modular/abvar/abvar_ambient_jacobian.py +++ b/src/sage/modular/abvar/abvar_ambient_jacobian.py @@ -16,9 +16,11 @@ from .abvar import (ModularAbelianVariety_modsym_abstract, ModularAbelianVariety, simple_factorization_of_modsym_space, modsym_lattices, ModularAbelianVariety_modsym) -from sage.rings.all import QQ +from sage.rings.all import QQ, Integer from sage.modular.modsym.modsym import ModularSymbols +from sage.modular.modform.constructor import Newforms +from sage.modular.arithgroup.all import is_Gamma0, is_Gamma1 from . import morphism @@ -391,6 +393,32 @@ def decomposition(self, simple=True, bound=None): self.__decomposition[simple] = D return D + def newform_decomposition(self, names=None): + """ + Return the newforms of the simple subvarieties in the decomposition of + self as a product of simple subvarieties, up to isogeny. + + OUTPUT: + - an array of newforms + EXAMPLES:: + sage: J0(81).newform_decomposition('a') + [q - 2*q^4 + O(q^6), q - 2*q^4 + O(q^6), q + a0*q^2 + q^4 - a0*q^5 + O(q^6)] + + sage: J1(19).newform_decomposition('a') + [q - 2*q^3 - 2*q^4 + 3*q^5 + O(q^6), + q + a1*q^2 + (-1/9*a1^5 - 1/3*a1^4 - 1/3*a1^3 + 1/3*a1^2 - a1 - 1)*q^3 + (4/9*a1^5 + 2*a1^4 + 14/3*a1^3 + 17/3*a1^2 + 6*a1 + 2)*q^4 + (-2/3*a1^5 - 11/3*a1^4 - 10*a1^3 - 14*a1^2 - 15*a1 - 9)*q^5 + O(q^6)] + """ + if self.dimension() == 0: + return [] + G = self.group() + if not (is_Gamma0(G) or is_Gamma1(G)): + return [S.newform(names=names) for S in self.decomposition()] + Gtype = G.parent() + N = G.level() + preans = [Newforms(Gtype(d), names=names) * + len(Integer(N/d).divisors()) for d in N.divisors()] + ans = [newform for l in preans for newform in l] + return ans diff --git a/src/sage/modular/abvar/lseries.py b/src/sage/modular/abvar/lseries.py index 3e3f88be1f9..62174d443b0 100644 --- a/src/sage/modular/abvar/lseries.py +++ b/src/sage/modular/abvar/lseries.py @@ -1,9 +1,6 @@ """ `L`-series of modular abelian varieties -At the moment very little functionality is implemented -- this is mostly a -placeholder for future planned work. - AUTHOR: - William Stein (2007-03) @@ -25,7 +22,12 @@ ########################################################################### from sage.structure.sage_object import SageObject -from sage.rings.all import Integer +from sage.rings.all import Integer, infinity, ZZ, QQ, CC +from sage.modules.free_module import span +from sage.modular.modform.constructor import Newform, CuspForms +from sage.modular.arithgroup.congroup_gamma0 import is_Gamma0 +from sage.misc.misc_c import prod + class Lseries(SageObject): """ @@ -72,7 +74,7 @@ class Lseries_complex(Lseries): sage: A.lseries() Complex L-series attached to Abelian variety J0(37) of dimension 2 """ - def __call__(self, s): + def __call__(self, s, prec=53): """ Evaluate this complex `L`-series at `s`. @@ -80,19 +82,61 @@ def __call__(self, s): - ``s`` -- complex number + - ``prec`` -- integer (default: 53) the number of bits of precision + used in computing the lseries of the newforms. + OUTPUT: a complex number L(A, s). - EXAMPLES: - This is not yet implemented:: + EXAMPLES:: + + sage: L = J0(23).lseries() + sage: L(1) + 0.248431866590600 + sage: L(1, prec=100) + 0.24843186659059968120725033931 + + sage: L = J0(389)[0].lseries() + sage: L(1) # long time (2s) + -1.33139759782370e-19 + sage: L(1, prec=100) # long time (2s) + 6.0129758648142797032650287762e-39 + sage: L.rational_part() + 0 + + sage: L = J1(23)[0].lseries() + sage: L(1) + 0.248431866590600 + + sage: J = J0(11) * J1(11) + sage: J.lseries()(1) + 0.0644356903227915 + + sage: L = JH(17,[2]).lseries() + sage: L(1) + 0.386769938387780 - sage: L = J0(37).lseries() - sage: L(2) - Traceback (most recent call last): - ... - NotImplementedError """ - raise NotImplementedError + abelian_variety = self.abelian_variety() + # Check for easy dimension zero case + if abelian_variety.dimension() == 0: + return CC(1) + try: + factors = self.__factors[prec] + return prod(L(s) for L in factors) + except AttributeError: + self.__factors = {} + except KeyError: + pass + abelian_variety = self.abelian_variety() + newforms = abelian_variety.newform_decomposition('a') + + factors = [newform.lseries(embedding=i, prec=prec) + for newform in newforms + for i in range(newform.base_ring().degree())] + self.__factors[prec] = factors + + return prod(L(s) for L in factors) def __cmp__(self, other): """ @@ -135,22 +179,90 @@ def _repr_(self): """ return "Complex L-series attached to %s"%self.abelian_variety() + def vanishes_at_1(self): + """ + Return True if `L(1)=0` and return False otherwise. + + OUTPUT: + a boolean + + EXAMPLES: + + Numerically, it appears that the `L`-series for J0(389) vanishes at 1. + This is confirmed by this algebraic computation. :: + + sage: L = J0(389)[0].lseries(); L + Complex L-series attached to Simple abelian subvariety 389a(1,389) of dimension 1 of J0(389) + sage: L(1) # long time (2s) + -1.33139759782370e-19 + sage: L.vanishes_at_1() + True + + Numerically, it appears that the `L`-series for J1(23) vanishes at 1. + But this algebraic computation shows otherwise. :: + + sage: L = J1(23).lseries(); L + Complex L-series attached to Abelian variety J1(23) of dimension 12 + sage: L(1) + 1.71571957480487e-7 + sage: L.vanishes_at_1() + False + sage: L(1, prec=100) + 1.7157195748048518516191946658e-7 + """ + abelian_variety = self.abelian_variety() + # Check for easy dimension zero case + if abelian_variety.dimension() == 0: + return False + if not abelian_variety.is_simple(): + from .constructor import AbelianVariety + decomp = (AbelianVariety(f) for f in + abelian_variety.newform_decomposition('a')) + return any(S.lseries().vanishes_at_1() for S in decomp) + modular_symbols = abelian_variety.modular_symbols() + Phi = modular_symbols.rational_period_mapping() + ambient_module = modular_symbols.ambient_module() + + e = ambient_module([0, infinity]) + return Phi(e).is_zero() + def rational_part(self): """ Return the rational part of this `L`-function at the central critical value 1. - NOTE: This is not yet implemented. + OUTPUT: + a rational number EXAMPLES:: - sage: J0(37).lseries().rational_part() - Traceback (most recent call last): - ... - NotImplementedError + sage: A, B = J0(43).decomposition() + sage: A.lseries().rational_part() + 0 + sage: B.lseries().rational_part() + 2/7 """ - raise NotImplementedError - + abelian_variety = self.abelian_variety() + modular_symbols = abelian_variety.modular_symbols() + Phi = modular_symbols.rational_period_mapping() + ambient_module = modular_symbols.ambient_module() + + if self.vanishes_at_1(): + return QQ(0) + else: + s = ambient_module.sturm_bound() + I = ambient_module.hecke_images(0, range(1, s+1)) + PhiTe = span([Phi(ambient_module(I[n])) + for n in range(I.nrows())], ZZ) + + ambient_plus = ambient_module.sign_submodule(1) + ambient_plus_cusp = ambient_plus.cuspidal_submodule() + PhiH1plus = span([Phi(x) for + x in ambient_plus_cusp.integral_basis()], ZZ) + + return PhiTe.index_in(PhiH1plus) + + lratio = rational_part class Lseries_padic(Lseries): """ diff --git a/src/sage/modular/abvar/torsion_subgroup.py b/src/sage/modular/abvar/torsion_subgroup.py index 41b331561e8..6d4f8c42996 100644 --- a/src/sage/modular/abvar/torsion_subgroup.py +++ b/src/sage/modular/abvar/torsion_subgroup.py @@ -60,7 +60,7 @@ 26 1 7 7 27 1 3 3 29 2 7 7 - 30 1 6 12 + 30 1 6 6 31 2 5 5 32 1 4 4 33 1 4 4 @@ -93,10 +93,13 @@ from sage.modular.abvar.torsion_point import TorsionPoint from sage.modules.module import Module from .finite_subgroup import FiniteSubgroup -from sage.rings.all import ZZ +from sage.rings.all import ZZ, QQ from sage.sets.primes import Primes -from sage.modular.arithgroup.all import is_Gamma0 -from sage.arith.all import divisors, gcd, prime_range +from sage.modular.arithgroup.all import is_Gamma0, is_Gamma1 +from sage.all import divisors, gcd, prime_range +from sage.modular.dirichlet import DirichletGroup +from sage.misc.misc_c import prod + class RationalTorsionSubgroup(FiniteSubgroup): """ @@ -164,32 +167,61 @@ def __cmp__(self, other): return cmp(self.abelian_variety(), other.abelian_variety()) return FiniteSubgroup.__cmp__(self, other) - def order(self): + def order(self, proof=True): """ Return the order of the torsion subgroup of this modular abelian variety. - This may fail if the multiple obtained by counting points modulo - `p` exceeds the divisor obtained from the rational cuspidal + This function may fail if the multiple obtained by counting points + modulo `p` exceeds the divisor obtained from the rational cuspidal subgroup. + The computation of the rational torsion order of J1(p) is conjectural + and will only be used if proof=False. See Section 6.2.3 of [CES2003]_. + + INPUT: + + - ``proof`` -- a boolean (default: True) + + OUTPUT: + + The order of this torsion subgroup. + + REFERENCE: + + .. [CES2003] Brian Conrad, Bas Edixhoven, William Stein + `J_1(p)` Has Connected Fibers + Documenta Math. 8 (2003) 331--408 + EXAMPLES:: - sage: a = J0(11) - sage: a.rational_torsion_subgroup().order() + sage: A = J0(11) + sage: A.rational_torsion_subgroup().order() 5 - sage: a = J0(23) - sage: a.rational_torsion_subgroup().order() + sage: A = J0(23) + sage: A.rational_torsion_subgroup().order() 11 - sage: t = J0(37)[1].rational_torsion_subgroup() - sage: t.order() + sage: T = J0(37)[1].rational_torsion_subgroup() + sage: T.order() 3 + + sage: J = J1(13) + sage: J.rational_torsion_subgroup().order() + 19 + + Sometimes the order can only be computed with proof=False. :: + + sage: J = J1(23) + sage: J.rational_torsion_subgroup().order() + Traceback (most recent call last): + ... + RuntimeError: Unable to compute order of torsion subgroup (it is in [408991, 9406793]) + + sage: J.rational_torsion_subgroup().order(proof=False) + 408991 + """ - try: - return self._order - except AttributeError: - pass - O = self.possible_orders() + O = self.possible_orders(proof=proof) if len(O) == 1: n = O[0] self._order = n @@ -246,10 +278,21 @@ def lattice(self): else: raise NotImplementedError("unable to compute the rational torsion subgroup in this case (there is no known general algorithm yet)") - def possible_orders(self): + def possible_orders(self, proof=True): """ - Return the possible orders of this torsion subgroup, computed from - a known divisor and multiple of the order. + Return the possible orders of this torsion subgroup. Outside of special + cases, this is done by computing a divisor and multiple of the order. + + INPUT: + + - ``proof`` -- a boolean (default: True) + + OUTPUT: + + - an array of positive integers + + The computation of the rational torsion order of J1(p) is conjectural + and will only be used if proof=False. See Section 6.2.3 of [CES2003]_. EXAMPLES:: @@ -258,24 +301,49 @@ def possible_orders(self): sage: J0(33).rational_torsion_subgroup().possible_orders() [100, 200] - Note that this function has not been implemented for `J_1(N)`, - though it should be reasonably easy to do so soon (see Conrad, - Edixhoven, Stein):: - sage: J1(13).rational_torsion_subgroup().possible_orders() - Traceback (most recent call last): - ... - NotImplementedError: torsion multiple only implemented for Gamma0 + [19] + sage: J1(16).rational_torsion_subgroup().possible_orders() + [1, 2, 4, 5, 10, 20] """ try: - return self._possible_orders + if proof: + return self._possible_orders + else: + return self._possible_orders_proof_false except AttributeError: pass + + A = self.abelian_variety() + N = A.level() + # return the order of the cuspidal subgroup in the J0(p) case + if A.is_J0() and N.is_prime(): + self._possible_orders = [QQ((A.level()-1)/12).numerator()] + self._possible_orders_proof_false = self._possible_orders + return self._possible_orders + + # the elliptic curve case + if A.dimension() == 1: + self._possible_orders = [A.elliptic_curve().torsion_order()] + self._possible_orders_proof_false = self._possible_orders + return self._possible_orders + + # the conjectural J1(p) case + if not proof and A.is_J1() and N.is_prime(): + epsilons = [epsilon for epsilon in DirichletGroup(N) + if not epsilon.is_trivial() and epsilon.is_even()] + bernoullis = [epsilon.bernoulli(2) for epsilon in epsilons] + self._possible_orders_proof_false = [ZZ(N/(2**(N-3))*prod(bernoullis))] + return self._possible_orders_proof_false + u = self.multiple_of_order() l = self.divisor_of_order() + assert u % l == 0 O = [l * d for d in divisors(u//l)] self._possible_orders = O + if u == l: + self._possible_orders_proof_false = O return O def divisor_of_order(self): @@ -283,24 +351,146 @@ def divisor_of_order(self): Return a divisor of the order of this torsion subgroup of a modular abelian variety. + OUTPUT: + + A divisor of this torsion subgroup. + EXAMPLES:: sage: t = J0(37)[1].rational_torsion_subgroup() sage: t.divisor_of_order() 3 + + sage: J = J1(19) + sage: J.rational_torsion_subgroup().divisor_of_order() + 4383 + """ + try: + return self._divisor_of_order + except: + pass + A = self.abelian_variety() + N = A.level() + if A.dimension() == 0: - return ZZ(1) - R = A.rational_cusp_subgroup() - return R.order() + self._divisor_of_order = ZZ(1) + return self._divisor_of_order + + # return the order of the cuspidal subgroup in the J0(p) case + if A.is_J0() and N.is_prime(): + self._divisor_of_order = QQ((A.level()-1)/12).numerator() + return self._divisor_of_order + + # The elliptic curve case + if A.dimension() == 1: + self._divisor_of_order = A.elliptic_curve().torsion_order() + return self._divisor_of_order + + # The J1(p) case + if A.is_J1() and N.is_prime(): + epsilons = [epsilon for epsilon in DirichletGroup(N) + if not epsilon.is_trivial() and epsilon.is_even()] + bernoullis = [epsilon.bernoulli(2) for epsilon in epsilons] + self._divisor_of_order = ZZ(N/(2**(N-3))*prod(bernoullis)) + return self._divisor_of_order + + # The Gamma0 case + if all(is_Gamma0(G) for G in A.groups()): + self._divisor_of_order = A.rational_cusp_subgroup().order() + return self._divisor_of_order + + # Unhandled case + self._divisor_of_order = ZZ(1) + return self._divisor_of_order + + def multiple_of_order(self, maxp=None, proof=True): + """ + Return a multiple of the order. + + INPUT: + + - ``proof`` -- a boolean (default: True) + + The computation of the rational torsion order of J1(p) is conjectural + and will only be used if proof=False. See Section 6.2.3 of [CES2003]_. + + EXAMPLES:: + + sage: from sage_modabvar import J0, J1 + sage: J = J1(11); J + Abelian variety J1(11) of dimension 1 + sage: J.rational_torsion_subgroup().multiple_of_order() + 5 + + sage: J = J0(17) + sage: J.rational_torsion_subgroup().order() + 4 + + This is an example where proof=False leads to a better bound and better + performance. :: + + sage: J = J1(23) + sage: J.rational_torsion_subgroup().multiple_of_order() # long time (2s) + 9406793 + sage: J.rational_torsion_subgroup().multiple_of_order(proof=False) + 408991 + """ + + try: + if proof: + return self._multiple_of_order + else: + return self._multiple_of_order_proof_false + except: + pass + + A = self.abelian_variety() + N = A.level() - def multiple_of_order(self, maxp=None): + if A.dimension() == 0: + self._multiple_of_order = ZZ(1) + self._multiple_of_order_proof_false = self._multiple_of_order + return self._multiple_of_order + + # return the order of the cuspidal subgroup in the J0(p) case + if A.is_J0() and N.is_prime(): + self._multiple_of_order = QQ((A.level()-1)/12).numerator() + self._multiple_of_order_proof_false = self._multiple_of_order + return self._multiple_of_order + + # The elliptic curve case + if A.dimension() == 1: + self._multiple_of_order = A.elliptic_curve().torsion_order() + self._multiple_of_order_proof_false = self._multiple_of_order + return self._multiple_of_order + + # The conjectural J1(p) case + if not proof and A.is_J1() and N.is_prime(): + epsilons = [epsilon for epsilon in DirichletGroup(N) + if not epsilon.is_trivial() and epsilon.is_even()] + bernoullis = [epsilon.bernoulli(2) for epsilon in epsilons] + self._multiple_of_order_proof_false = ZZ(N/(2**(N-3))*prod(bernoullis)) + return self._multiple_of_order_proof_false + + # The Gamma0 and Gamma1 case + if all((is_Gamma0(G) or is_Gamma1(G) for G in A.groups())): + self._multiple_of_order = self.multiple_of_order_using_frobp() + return self._multiple_of_order + + # Unhandled case + raise NotImplementedError("No implemented algorithm") + + def multiple_of_order_using_frobp(self, maxp=None): """ Return a multiple of the order of this torsion group. - The multiple is computed using characteristic polynomials of Hecke - operators of odd index not dividing the level. + In the `Gamma_0` case, the multiple is computed using characteristic + polynomials of Hecke operators of odd index not dividing the level. In + the `Gamma_1` case, the multiple is computed by expressing the + frobenius polynomial in terms of the characteristic polynomial of left + multiplication by `a_p` for odd primes p not dividing the level. INPUT: @@ -316,61 +506,84 @@ def multiple_of_order(self, maxp=None): sage: J = J0(11) sage: G = J.rational_torsion_subgroup() - sage: G.multiple_of_order(11) + sage: G.multiple_of_order_using_frobp(11) 5 + + Increasing maxp may yield a tighter bound. If maxp=None, then Sage + will use more primes until the multiple stabilizes for 3 successive + primes. :: + sage: J = J0(389) sage: G = J.rational_torsion_subgroup(); G Torsion subgroup of Abelian variety J0(389) of dimension 32 - sage: G.multiple_of_order() + sage: G.multiple_of_order_using_frobp() 97 - sage: [G.multiple_of_order(p) for p in prime_range(3,11)] + sage: [G.multiple_of_order_using_frobp(p) for p in prime_range(3,11)] [92645296242160800, 7275, 291] - sage: [G.multiple_of_order(p) for p in prime_range(3,13)] + sage: [G.multiple_of_order_using_frobp(p) for p in prime_range(3,13)] [92645296242160800, 7275, 291, 97] - sage: [G.multiple_of_order(p) for p in prime_range(3,19)] + sage: [G.multiple_of_order_using_frobp(p) for p in prime_range(3,19)] [92645296242160800, 7275, 291, 97, 97, 97] - :: + We can compute the multiple of order of the torsion subgroup for Gamma0 + and Gamma1 varieties, and their products. :: - sage: J = J0(33) * J0(11) ; J.rational_torsion_subgroup().order() - Traceback (most recent call last): - ... - NotImplementedError: torsion multiple only implemented for Gamma0 + sage: A = J0(11) * J0(33) + sage: A.rational_torsion_subgroup().multiple_of_order_using_frobp() + 1000 + + sage: A = J1(23) + sage: A.rational_torsion_subgroup().multiple_of_order_using_frobp() + 9406793 + sage: A.rational_torsion_subgroup().multiple_of_order_using_frobp(maxp=50) + 408991 + + sage: A = J1(19) * J0(21) + sage: A.rational_torsion_subgroup().multiple_of_order_using_frobp() + 35064 The next example illustrates calling this function with a larger input and how the result may be cached when maxp is None:: sage: T = J0(43)[1].rational_torsion_subgroup() - sage: T.multiple_of_order() + sage: T.multiple_of_order_using_frobp() 14 - sage: T.multiple_of_order(50) + sage: T.multiple_of_order_using_frobp(50) 7 - sage: T.multiple_of_order() + sage: T.multiple_of_order_using_frobp() 7 + + This function is not implemented for general congruence subgroups + unless the dimension is zero. :: + + sage: from sage_modabvar import JH + sage: A = JH(13,[2]); A + Abelian variety J0(13) of dimension 0 + sage: A.rational_torsion_subgroup().multiple_of_order_using_frobp() + 1 + + sage: A = JH(15, [2]); A + Abelian variety JH(15,[2]) of dimension 1 + sage: A.rational_torsion_subgroup().multiple_of_order_using_frobp() + Traceback (most recent call last): + ... + NotImplementedError: torsion multiple only implemented for Gamma0 and Gamma1 """ if maxp is None: try: - return self.__multiple_of_order + return self.__multiple_of_order_using_frobp except AttributeError: pass - bnd = ZZ(0) A = self.abelian_variety() if A.dimension() == 0: T = ZZ(1) - self.__multiple_of_order = T + self.__multiple_of_order_using_frobp = T return T + if not all((is_Gamma0(G) or is_Gamma1(G) for G in A.groups())): + raise NotImplementedError("torsion multiple only implemented for Gamma0 and Gamma1") + + bnd = ZZ(0) N = A.level() - if not (len(A.groups()) == 1 and is_Gamma0(A.groups()[0])): - # to generalize to this case, you'll need to - # (1) define a charpoly_of_frob function: - # this is tricky because I don't know a simple - # way to do this for Gamma1 and GammaH. Will - # probably have to compute explicit matrix for - #

operator (add to modular symbols code), - # then compute some charpoly involving - # that directly... - # (2) use (1) -- see my MAGMA code. - raise NotImplementedError("torsion multiple only implemented for Gamma0") cnt = 0 if maxp is None: X = Primes() @@ -380,8 +593,40 @@ def multiple_of_order(self, maxp=None): if (2*N) % p == 0: continue - f = A.hecke_polynomial(p) - b = ZZ(f(p+1)) + if (len(A.groups()) == 1 and is_Gamma0(A.groups()[0])): + f = A.hecke_polynomial(p) + b = ZZ(f(p+1)) + else: + from .constructor import AbelianVariety + D = [AbelianVariety(f) for f in + A.newform_decomposition('a')] + b = 1 + for simple in D: + G = simple.newform_level()[1] + if is_Gamma0(G): + f = simple.hecke_polynomial(p) + b *= ZZ(f(p+1)) + else: + f = simple.newform('a') + Kf = f.base_ring() + eps = f.character() + Qe = eps.base_ring() + + if Kf != QQ: + # relativize number fields to compute charpoly of + # left multiplication of ap on Kf as a Qe-vector + # space. + Lf = Kf.relativize(Qe.gen(), 'a') + to_Lf = Lf.structure()[1] + + name = Kf._names[0] + ap = to_Lf(f.modular_symbols(1).eigenvalue(p, name)) + + G_ps = ap.matrix().charpoly() + b *= ZZ(Qe(G_ps(1 + to_Lf(eps(p))*p)).norm()) + else: + ap = f.modular_symbols(1).eigenvalue(p) + b *= ZZ(1 + eps(p)*p - ap) if bnd == 0: bnd = b @@ -399,21 +644,23 @@ def multiple_of_order(self, maxp=None): # will be used if this function is called # again with maxp equal to None (the default). if maxp is None: - # maxp is None but self.__multiple_of_order is + # maxp is None but self.__multiple_of_order_using_frobp is # not set, since otherwise we would have immediately # returned at the top of this function - self.__multiple_of_order = bnd + self.__multiple_of_order_using_frobp = bnd else: # maxp is given -- record new info we get as # a gcd... try: - self.__multiple_of_order = gcd(self.__multiple_of_order, bnd) + self.__multiple_of_order_using_frobp = \ + gcd(self.__multiple_of_order_using_frobp, bnd) except AttributeError: - # ... except in the case when self.__multiple_of_order - # was never set. In this case, we just set - # it as long as the gcd stabilized for 3 in a row. + # ... except in the case when + # self.__multiple_of_order_using_frobp was never set. In this + # case, we just set it as long as the gcd stabilized for 3 in a + # row. if cnt >= 2: - self.__multiple_of_order = bnd + self.__multiple_of_order_using_frobp = bnd return bnd