Ticket #7575: trac_7575.patch

File trac_7575.patch, 65.8 KB (added by Robert Miller, 13 years ago)

depends on #8184, #8155 (and possibly #8124)

  • new file sage/schemes/elliptic_curves/BSD.py

    # HG changeset patch
    # User Robert L. Miller <rlm@rlmiller.org>
    # Date 1263943454 28800
    # Node ID 8a90b87fca0ce33dded2bd986d413441ac42ff0f
    # Parent  1ca8aa94810aa7b13cd1671e651b3ed6a365db76
    #7575: introduce BSD.py module, catch errors, fix height bounds in heegner_index_bound, mwrank interface & doc fixes
    
    diff -r 1ca8aa94810a -r 8a90b87fca0c sage/schemes/elliptic_curves/BSD.py
    - +  
     1# -*- coding: utf-8 -*-
     2#import ell_point
     3#import formal_group
     4#import ell_torsion
     5#from ell_generic import EllipticCurve_generic, is_EllipticCurve
     6#from ell_number_field import EllipticCurve_number_field
     7
     8#import sage.groups.all
     9import sage.rings.arith as arith
     10import sage.rings.all as rings
     11from sage.rings.all import ZZ
     12
     13def prove_BSD(self, verbosity=0, simon=False, proof=None, secs_hi=30):
     14    r"""
     15    Attempts to prove the Birch and Swinnerton-Dyer conjectural
     16    formula for `E`, returning a list of primes `p` for which this
     17    function fails to prove BSD(E,p).  Here, BSD(E,p) is the
     18    statement: "the Birch and Swinnerton-Dyer formula holds up to a
     19    rational number coprime to `p`."
     20   
     21    INPUT:
     22   
     23        - ``verbosity`` - int, how much information about the proof to print.
     24       
     25            - 0 - print nothing
     26            - 1 - print sketch of proof
     27            - 2 - print information about remaining primes
     28       
     29        - ``simon`` - bool (default False), whether to use two_descent or
     30          simon_two_descent at p=2.
     31   
     32        - ``proof`` - bool or None (default: None, see
     33          proof.elliptic_curve or sage.structure.proof). If False, this
     34          function just immediately returns the empty list.
     35       
     36        - ``secs_hi`` - maximum number of seconds to try to compute the
     37          Heegner index before switching over to trying to compute the
     38          Heegner index bound. (Rank 0 only!)
     39
     40    NOTE:
     41   
     42    When printing verbose output, phrases such as "by Mazur" are referring
     43    to the following list of papers:
     44
     45    REFERENCES:
     46   
     47    .. [Cha] B. Cha. Vanishing of some cohomology goups and bounds for the
     48       Shafarevich-Tate groups of elliptic curves. J. Number Theory, 111:154-
     49       178, 2005.
     50    .. [Jetchev] D. Jetchev. Global divisibility of Heegner points and
     51       Tamagawa numbers. Compos. Math. 144 (2008), no. 4, 811--826.
     52    .. [Kato] K. Kato. p-adic Hodge theory and values of zeta functions of
     53       modular forms. Astérisque, (295):ix, 117-290, 2004.
     54    .. [Kolyvagin] V. A. Kolyvagin. On the structure of Shafarevich-Tate
     55       groups. Algebraic geometry, 94--121, Lecture Notes in Math., 1479,
     56       Springer, Berlin, 1991.
     57    .. [LumStein] A. Lum, W. Stein. Verification of the Birch and
     58       Swinnerton-Dyer Conjecture for Elliptic Curves with Complex
     59       Multiplication (unpublished)
     60    .. [Mazur] B. Mazur. Modular curves and the Eisenstein ideal. Inst.
     61       Hautes Études Sci. Publ. Math. No. 47 (1977), 33--186 (1978).
     62    .. [Rubin] K. Rubin. The "main conjectures" of Iwasawa theory for
     63       imaginary quadratic fields. Invent. Math. 103 (1991), no. 1, 25--68.
     64    .. [SteinWuthrich] W. Stein and C. Wuthrich. Computations about
     65       Tate-Shafarevich groups using Iwasawa theory.
     66       http://wstein.org/papers/shark, February 2008.
     67    .. [SteinEtAl] G. Grigorov, A. Jorza, S. Patrikis, W. Stein,
     68       C. Tarniţǎ. Computational verification of the Birch and
     69       Swinnerton-Dyer conjecture for individual elliptic curves.
     70       Math. Comp. 78 (2009), no. 268, 2397--2425.
     71
     72
     73    EXAMPLES::
     74   
     75        sage: EllipticCurve('11a').prove_BSD(verbosity=2)
     76        p = 2: True by 2-descent
     77        True for p not in {2, 5} by Kolyvagin.
     78        True for p=5 by Mazur
     79        []
     80       
     81        sage: EllipticCurve('14a').prove_BSD(verbosity=2)
     82        p = 2: True by 2-descent
     83        True for p not in {2, 3} by Kolyvagin.
     84        Remaining primes:
     85        p = 3: reducible, not surjective, good ordinary, divides a Tamagawa number
     86        [3]
     87        sage: EllipticCurve('14a').prove_BSD(simon=True)
     88        [3]
     89
     90    A rank two curve::
     91   
     92        sage: E = EllipticCurve('389a')
     93
     94    We know nothing with proof=True::
     95   
     96        sage: E.prove_BSD()
     97        Set of all prime numbers: 2, 3, 5, 7, ...
     98
     99    We (think we) know everything with proof=False::
     100   
     101        sage: E.prove_BSD(proof=False)
     102        []
     103
     104    A curve of rank 0 and prime conductor::
     105   
     106        sage: E = EllipticCurve('19a')
     107        sage: E.prove_BSD(verbosity=2)
     108        p = 2: True by 2-descent
     109        True for p not in {2, 3} by Kolyvagin.
     110        True for p=3 by Mazur
     111        []
     112
     113        sage: E = EllipticCurve('37a')
     114        sage: E.rank()
     115        1
     116        sage: E._EllipticCurve_rational_field__rank
     117        {True: 1}
     118        sage: E.analytic_rank = lambda : 0
     119        sage: E.prove_BSD()
     120        Traceback (most recent call last):
     121        ...
     122        RuntimeError: It seems that the rank conjecture does not hold for this curve (Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field)! This may be a counterexample to BSD, but is more likely a bug.
     123
     124    We test the consistency check for the 2-part of Sha::
     125
     126        sage: E = EllipticCurve('37a')
     127        sage: S = E.sha(); S
     128        Shafarevich-Tate group for the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
     129        sage: S.an = lambda : 4
     130        sage: E.prove_BSD()
     131        Traceback (most recent call last):
     132        ...
     133        RuntimeError: Two descent => ord_2(#Sha[2]) = 0, but ord_2(#sha_an) = 2.
     134
     135    An example with a Tamagawa number at 5::
     136
     137        sage: E = EllipticCurve('123a1')
     138        sage: E.prove_BSD(verbosity=2)
     139        p = 2: True by 2-descent
     140        True for p not in {2, 5} by Kolyvagin.
     141        Remaining primes:
     142        p = 5: reducible, not surjective, good ordinary, divides a Tamagawa number
     143        [5]
     144
     145    A curve for which 3 divides the order of the Shafarevich-Tate group::
     146
     147        sage: E = EllipticCurve('681b')                       
     148        sage: E.prove_BSD(verbosity=2)               # long time
     149        p = 2: True by 2-descent
     150        True for p not in {2, 3} by Kolyvagin.
     151        ALERT: p = 3 left in Kolyvagin bound
     152            0 <= ord_p(#Sha) <= 2
     153            ord_p(#Sha_an) = 2
     154        Remaining primes:
     155        p = 3: irreducible, surjective, non-split multiplicative
     156        [3]
     157   
     158    A curve for which we need to use ``heegner_index_bound``::
     159
     160        sage: E = EllipticCurve('198b')
     161        sage: E.prove_BSD(verbosity=1, secs_hi=1)
     162        p = 2: True by 2-descent
     163        Timeout stopped Heegner index computation...
     164        Proceeding to use heegner_index_bound instead.
     165        True for p not in {2, 3} by Kolyvagin.
     166        [3]
     167
     168    TESTS:
     169   
     170    This was fixed by trac #8184 and #7575::
     171   
     172        sage: EllipticCurve('438e1').prove_BSD(verbosity=1)
     173        p = 2: True by 2-descent
     174        True for p not in {2} by Kolyvagin.
     175        []
     176   
     177    ::
     178   
     179        sage: E = EllipticCurve('960d1')
     180        sage: E.prove_BSD(verbosity=1)
     181        p = 2: True by 2-descent
     182        Timeout stopped Heegner index computation...
     183        Proceeding to use heegner_index_bound instead.
     184        True for p not in {2} by Kolyvagin.
     185        []
     186   
     187    """
     188    if proof is None:
     189        from sage.structure.proof.proof import get_flag
     190        proof = get_flag(proof, "elliptic_curve")
     191    else:
     192        proof = bool(proof)
     193    if not proof:
     194        return []
     195    two_tor_rk = self.two_torsion_rank()
     196    Sha = self.sha()
     197    sha_an = Sha.an()
     198    if simon:
     199        rank_lower_bd, two_sel_rk_bd = self.simon_two_descent()[:2]
     200        two_sel_rk_bd -= two_tor_rk
     201        sha2_lower_bd = 0
     202    else:
     203        MWRC = self.mwrank_curve()
     204        two_sel_rk_bd = MWRC.selmer_rank()
     205        two_sel_rk_bd -= two_tor_rk
     206        rank_lower_bd = MWRC.rank()
     207        sha2_lower_bd = two_sel_rk_bd - MWRC.rank_bound()
     208    sha2_upper_bd = two_sel_rk_bd - rank_lower_bd
     209    # note: two_sel_rk_bd will include sha
     210    rank = None
     211    if rank_lower_bd > 1:
     212        # We do not know BSD(E,p) for even a single p, since it's
     213        # an open problem to show that L^r(E,1)/(Reg*Omega) is
     214        # rational for any curve with r >= 2.
     215        from sage.sets.all import Primes
     216        return Primes()
     217    if sha2_upper_bd == sha2_lower_bd:
     218        rank = rank_lower_bd
     219        if sha_an.ord(2) != sha2_lower_bd:
     220            raise RuntimeError("Two descent => ord_2(#Sha[2]) = %d, but ord_2(#sha_an) = %d."%(sha2_lower_bd,sha_an.ord(2)))
     221        if verbosity > 0:
     222            print 'p = 2: True by 2-descent'
     223        two_proven = True
     224    else:
     225        if sha2_upper_bd < sha2_lower_bd:
     226            raise RuntimeError("Apparent contradiction: ord_2(#Sha[2]_an) == %d, rank(2-Sel)-rank(2-tor) = %d, rank >= %d, ord_2(#Sha[2]) >= %d"%(sha_an.ord(2),two_sel_rk_bd,rank_lower_bd,sha2_lower_bd))
     227        if two_sel_rk_bd - rank_lower_bd == sha_an.ord(2):
     228            if verbosity > 0:
     229                print 'p = 2: ord_2(#Sha[2]_an) == %d >= ord_2(#Sha[2]).'%sha_an.ord(2)
     230        elif two_sel_rk_bd - rank_lower_bd > sha_an.ord(2):
     231            if verbosity > 0:
     232                print 'p = 2: ord_2(#Sha[2]_an) == %d, and ord_2(#Sha[2]) <= %d.'%(sha_an.ord(2),two_sel_rk_bd - rank_lower_bd)
     233        else:
     234            raise RuntimeError("Apparent contradiction: ord_2(#Sha[2]_an) == %d, rank(2-Sel)-rank(2-tor) = %d, rank >= %d"%(sha_an.ord(2),two_sel_rk_bd,rank_lower_bd))
     235        two_proven = False
     236        if verbosity > 0:
     237            print 'Looking rank up in database...'
     238        rank = self.rank(use_database=True, only_use_mwrank=False)
     239
     240    if rank != self.analytic_rank():
     241        raise RuntimeError("It seems that the rank conjecture does not hold for this curve (%s)! This may be a counterexample to BSD, but is more likely a bug."%(self))
     242
     243    # We replace self by the optimal curve, which we can do since
     244    # truth of BSD(E,p) is invariant under isogeny.
     245    self = self.optimal_curve()
     246   
     247    N = self.conductor()
     248   
     249    # reduce set of remaining primes to a finite set
     250    import signal
     251    remaining_primes = []
     252    kolyvagin_primes = []
     253    heegner_index = None
     254    if self.rank() == 0:
     255        try:
     256            old_alarm = signal.alarm(secs_hi)
     257            old_alarm_set = (old_alarm != 0)
     258            for D in self.heegner_discriminants_list(10):
     259                I = None
     260                while I is None:
     261                    dsl=15
     262                    try:
     263                        I = self.heegner_index(D, descent_second_limit=dsl)
     264                    except RuntimeError as err:
     265                        if err.args[0][-33:] == 'Generators not provably computed.':
     266                            dsl += 1
     267                        else: raise RuntimeError(err)
     268                J = I.is_int()
     269                if J[0] and J[1]>0:
     270                    I = J[1]
     271                else:
     272                    J = (2*I).is_int()
     273                    if J[0] and J[1]>0:
     274                        I = J[1]
     275                    else:
     276                        I = None
     277                if I is not None:
     278                    if heegner_index is None:
     279                        heegner_index = I
     280                        break # no big long loops just yet...
     281            old_alarm_sub = signal.alarm(0)
     282            if old_alarm_set:
     283                old_alarm -= old_alarm_sub
     284        except KeyboardInterrupt:
     285            if signal.alarm(0)==0:
     286                print 'Timeout stopped Heegner index computation...'
     287                print 'Proceeding to use heegner_index_bound instead.'
     288            else:
     289                raise KeyboardInterrupt
     290        except RuntimeError as err:
     291            if err.args[0][:37] == 'End Of File (EOF) in read_nonblocking':
     292                print 'Computing Heegner index failed due to mwrank interface: see #7575.'
     293            else: raise RuntimeError(err)
     294            old_alarm_sub = signal.alarm(0)
     295            if old_alarm_set:
     296                old_alarm -= old_alarm_sub
     297                if old_alarm <= 0:
     298                    raise KeyboardInterrupt
     299                signal.alarm(old_alarm)
     300        if old_alarm_set: # in case alarm was already set...
     301            if old_alarm <= 0:
     302                raise KeyboardInterrupt
     303            signal.alarm(old_alarm)
     304        if heegner_index is None:
     305            for D in self.heegner_discriminants_list(100):
     306                max_height = 12
     307                heegner_primes = -1
     308                while heegner_primes == -1:
     309                    max_height += 1
     310                    heegner_primes, _ = self.heegner_index_bound(D, max_height=max_height)
     311                if isinstance(heegner_primes, list):
     312                    break
     313            if not isinstance(heegner_primes, list):
     314                raise RuntimeError("Tried 100 Heegner discriminants, and heegner_index_bound failed each time.")
     315            if 2 in heegner_primes:
     316                heegner_primes.remove(2)
     317        else:
     318            heegner_primes = [p for p in arith.prime_divisors(heegner_index) if p!=2]
     319    else: # rank 1
     320        for D in self.heegner_discriminants_list(10):
     321            I = self.heegner_index(D)
     322            J = I.is_int()
     323            if J[0] and J[1]>0:
     324                I = J[1]
     325            else:
     326                J = (2*I).is_int()
     327                if J[0] and J[1]>0:
     328                    I = J[1]
     329                else:
     330                    continue
     331            heegner_index = I
     332            break
     333        heegner_primes = [p for p in arith.prime_divisors(heegner_index) if p!=2]
     334   
     335    if self.has_cm():
     336        # ensure that CM is by a maximal order
     337        non_max_j_invs = [-12288000, 54000, 287496, 16581375]
     338        if self.j_invariant() in non_max_j_invs:
     339            for E in self.isogeny_class()[0]:
     340                if E.j_invariant() not in non_max_j_invs:
     341                    Sha = E.sha()
     342                    sha_an = Sha.an()
     343                    if verbosity > 0:
     344                        print 'CM by non maximal order: switching curves'
     345                    break
     346        else:
     347            E = self
     348        if E.analytic_rank() == 0:
     349            if verbosity > 0:
     350                print 'p >= 5: true by Rubin'
     351            remaining_primes.append(3)
     352        else:
     353            K = rings.QuadraticField(E.cm_discriminant(), 'a')
     354            D_K = K.disc()
     355            D_E = E.discriminant()
     356            if len(K.factor(3)) == 1: # 3 does not split in K
     357                remaining_primes.append(3)
     358            for p in arith.prime_divisors(D_K):
     359                if p >= 5:
     360                    remaining_primes.append(p)
     361            for p in arith.prime_divisors(D_E):
     362                if p >= 5 and D_K%p and len(K.factor(p)) == 1: # p is inert in K
     363                    remaining_primes.append(p)
     364            for p in heegner_primes:
     365                if p >= 5 and D_E%p != 0 and D_K%p != 0 and len(K.factor(p)) == 1: # p is good for E and inert in K
     366                    kolyvagin_primes.append(p)
     367            assert sha_an in ZZ and sha_an > 0
     368            for p in arith.prime_divisors(sha_an):
     369                if p >= 5 and D_K%p != 0 and len(K.factor(p)) == 1:
     370                    if E.is_good(p):
     371                        if verbosity > 2 and p in heegner_primes and heegner_index is None:
     372                            print 'ALERT: Prime p (%d) >= 5 dividing sha_an, good for E, inert in K, in heegner_primes, should not divide the actual Heegner index'
     373                        # Note that the following check is not entirely
     374                        # exhaustive, in case there is a p not dividing
     375                        # the Heegner index in heegner_primes,
     376                        # for which only an outer bound was computed
     377                        if p not in heegner_primes:
     378                            raise RuntimeError("p = %d divides sha_an, is of good reduction for E, inert in K, and does not divide the Heegner index. This may be a counterexample to BSD, but is more likely a bug. %s"%(p,self))
     379            if verbosity > 0:
     380                print 'True for p not in {%s} by Kolyvagin (via Stein & Lum -- unpublished) and Rubin.'%str(list(set(remaining_primes).union(set(kolyvagin_primes))))[1:-1]
     381    else: # no CM
     382        E = self
     383        # do some tricks to get to a finite set without calling bound_kolyvagin
     384        remaining_primes = [p for p, reason in E.non_surjective()]
     385        for p in heegner_primes:
     386            if p not in remaining_primes:
     387                remaining_primes.append(p)
     388        assert sha_an in ZZ and sha_an > 0
     389        for p in arith.prime_divisors(sha_an):
     390            if p not in remaining_primes:
     391                remaining_primes.append(p)
     392        if 2 in remaining_primes: remaining_primes.remove(2)
     393        if verbosity > 0:
     394            print 'True for p not in {' + str([2]+list(remaining_primes))[1:-1] + '} by Kolyvagin.'
     395        primes_to_remove = []
     396        for p in remaining_primes:
     397            if E.is_surjective(p)[0] and not E.has_additive_reduction(p):
     398                if E.has_nonsplit_multiplicative_reduction(p):
     399                    if E.rank() > 0:
     400                        continue
     401                if p==3:
     402                    if (not (E.is_ordinary(p) and E.is_good(p))) and (not E.has_split_multiplicative_reduction(p)):
     403                        continue
     404                    if E.rank() > 0:
     405                        continue
     406                if verbosity > 1:
     407                    print 'p = %d: Trying p_primary_bound'%p
     408                p_bound = Sha.p_primary_bound(p)
     409                if sha_an.ord(p) == 0 and p_bound == 0:
     410                    if verbosity > 0:
     411                        print 'True for p=%d by Stein-Wuthrich.'%p
     412                    primes_to_remove.append(p)
     413                else:
     414                    print 'Analytic %d-rank is '%p + str(sha_an.ord(p)) + ', actual %d-rank is at most %d.'%(p, p_bound)
     415                    print '    by Stein-Wuthrich.\n'
     416        for p in primes_to_remove:
     417            remaining_primes.remove(p)
     418        kolyvagin_primes = []
     419        for p in remaining_primes:
     420            if E.is_surjective(p)[0]:
     421                kolyvagin_primes.append(p)
     422        for p in kolyvagin_primes:
     423            remaining_primes.remove(p)
     424    # apply other hypotheses which imply Kolyvagin's bound holds
     425    bounded_primes = []
     426    D_K = rings.QuadraticField(D, 'a').disc()
     427    assert 2 not in remaining_primes
     428    # Cha's hypothesis
     429    for p in remaining_primes:
     430        if D_K%p != 0 and N%(p**2) != 0 and E.is_irreducible(p):
     431            if verbosity > 0:
     432                print 'Kolyvagin\'s bound for p = %d applies by Cha.'%p
     433            kolyvagin_primes.append(p)
     434    # Stein et al.
     435    if not E.has_cm():
     436        L = arith.lcm([F.torsion_order() for F in E.isogeny_class()[0]])
     437        for p in remaining_primes:
     438            if p in kolyvagin_primes: continue
     439            if L%p != 0:
     440                if len(arith.prime_divisors(D_K)) == 1:
     441                    if D_K%p != 0:
     442                        if verbosity > 0:
     443                            print 'Kolyvagin\'s bound for p = %d applies by Stein et al.'%p
     444                        kolyvagin_primes.append(p)
     445                else:
     446                    if verbosity > 0:
     447                        print 'Kolyvagin\'s bound for p = %d applies by Stein et al.'%p
     448                    kolyvagin_primes.append(p)
     449    for p in kolyvagin_primes:
     450        if p in remaining_primes:
     451            remaining_primes.remove(p)
     452
     453    prime_bounds = []
     454    # apply Kolyvagin's bound
     455    primes_to_remove = []
     456    for p in kolyvagin_primes:
     457        if sha_an.ord(p) == 0 and p not in heegner_primes:
     458                if verbosity > 0:
     459                    print 'True for p = %d by Kolyvagin bound.'%p
     460                primes_to_remove.append(p)
     461                continue
     462        if heegner_index is not None: # p must divide heegner_index
     463            ord_p_bound = 2*heegner_index.ord(p)
     464            # Here Jetchev's results apply.
     465            m_max = max([E.tamagawa_number(q).ord(p) for q in N.prime_divisors()])
     466            if m_max > 0 and verbosity > 0:
     467                print 'Jetchev\'s results apply (at p = %d) with m_max ='%p, m_max
     468            ord_p_bound -= 2*m_max
     469            if ord_p_bound == 0:
     470                if sha_an.ord(p) != 0:
     471                    raise RuntimeError("p = %d: ord_p_bound == 0, but sha_an.ord(p) == %d. This appears to be a counterexample to BSD, but is more likely a bug."%(p,sha_an.ord(p)))
     472                if verbosity > 0:
     473                    print 'True for p = %d by Kolyvagin bound.'%p
     474                primes_to_remove.append(p)
     475                continue
     476        elif p not in heegner_primes:
     477            ord_p_bound = 0
     478        else:
     479            from sage.rings.infinity import Infinity
     480            ord_p_bound = Infinity
     481            if verbosity > 0:
     482                print 'p = %d may divide the Heegner index, for which only a bound was computed.'%p
     483        if verbosity > 0:
     484            print 'ALERT: p = %d left in Kolyvagin bound'%p
     485            print '    0 <= ord_p(#Sha) <=', ord_p_bound
     486            print '    ord_p(#Sha_an) =', sha_an.ord(p)
     487    for p in primes_to_remove:
     488        kolyvagin_primes.remove(p)
     489    remaining_primes = list( set(remaining_primes).union(set(kolyvagin_primes)) )
     490   
     491    # Kato's bound
     492    if rank == 0 and not E.has_cm():
     493        L_over_Omega = E.lseries().L_ratio()
     494        kato_primes = Sha.bound_kato()
     495        primes_to_remove = []
     496        for p in remaining_primes:
     497            if p not in kato_primes:
     498                if verbosity > 0:
     499                    print 'Kato further implies that #Sha[%d] is trivial.'%p
     500                primes_to_remove.append(p)
     501            if p not in [2,3] and N%p != 0:
     502                if E.is_surjective(p)[0]:
     503                    if verbosity > 1:
     504                        print 'Kato might apply nontrivially for %d'%p
     505                    # ordp(sha) <= ordp(L_over_omega)
     506        for p in primes_to_remove:
     507            remaining_primes.remove(p)
     508   
     509    # Mazur
     510    if N.is_prime():
     511        for p in remaining_primes:
     512            if E.is_reducible(p):
     513                remaining_primes.remove(p)                   
     514                if verbosity > 0:
     515                    print 'True for p=%s by Mazur'%p
     516   
     517    if two_proven is False:
     518        remaining_primes.append(2)
     519   
     520    # print some extra information
     521    remaining_primes.sort()
     522    if verbosity > 1:
     523        if len(remaining_primes) > 0:
     524            print 'Remaining primes:'
     525        for p in remaining_primes:
     526            s = 'p = ' + str(p) + ': '
     527            if not E.is_reducible(p):
     528                s += 'ir'
     529            s += 'reducible, '
     530            if not E.is_surjective(p)[0]:
     531                s += 'not '
     532            s += 'surjective, '
     533            a_p = E.an(p)
     534            if E.is_good(p):
     535                if a_p%p != 0:
     536                    s += 'good ordinary'
     537                else:
     538                    s += 'good, non-ordinary'
     539            else:
     540                assert E.is_minimal()
     541                if a_p == 0:
     542                    s += 'additive'
     543                elif a_p == 1:
     544                    s += 'split multiplicative'
     545                elif a_p == -1:
     546                    s += 'non-split multiplicative'
     547            if E.tamagawa_product()%p==0:
     548                s += ', divides a Tamagawa number'
     549            print s
     550
     551    return remaining_primes
     552
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r 1ca8aa94810a -r 8a90b87fca0c sage/schemes/elliptic_curves/ell_rational_field.py
    a b  
    7979from sage.rings.padics.precision_error import PrecisionError
    8080
    8181import heegner
     82import BSD
    8283
    8384from lseries_ell import Lseries_ell
    8485
     
    897898       
    898899        OUTPUT:
    899900       
    900         Nothing - nothing is returned (though much is printed unless
    901         verbose=False)
    902        
    903         EXAMPLES::
    904        
    905             sage: E=EllipticCurve('37a1')
    906             sage: E.two_descent(verbose=False) # no output
    907         """
    908         self.mwrank_curve().two_descent(verbose, selmer_only,
    909                                         first_limit, second_limit,
    910                                         n_aux, second_descent)
    911 
     901        Returns ``True`` if the descent succeeded, i.e. if the lower bound and
     902        the upper bound for the rank are the same. In this case, generators and
     903        the rank are cached. A return value of ``False`` indicates that either
     904        rational points were not found, or that Sha[2] is nontrivial and mwrank
     905        was unable to determine this for sure.
     906       
     907        EXAMPLES::
     908       
     909            sage: E=EllipticCurve('37a1')
     910            sage: E.two_descent(verbose=False)
     911            True
     912           
     913        """
     914        misc.verbose("Calling mwrank C++ library.")
     915        C = self.mwrank_curve()
     916        C.two_descent(verbose, selmer_only,
     917                        first_limit, second_limit,
     918                        n_aux, second_descent)
     919        if C.certain():
     920            self.__gens[True] = [self.point(x, check=True) for x in C.gens()]
     921            self.__gens[True].sort()
     922            self.__rank[True] = len(self.__gens[True])
     923        return C.certain()
    912924
    913925    ####################################################################
    914926    #  Etc.
     
    14371449
    14381450    def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30):
    14391451        r"""
    1440         Given a curve with no 2-torsion, computes (probably) the rank of
    1441         the Mordell-Weil group, with certainty the rank of the 2-Selmer
    1442         group, and a list of independent points on the curve.
     1452        Computes a lower bound for the rank of the Mordell-Weil group `E(Q)`,
     1453        the rank of the 2-Selmer group, and a list of independent points on
     1454        `E(Q)/2E(Q)`.
    14431455       
    14441456        INPUT:
    14451457       
     
    14661478        OUTPUT:
    14671479       
    14681480       
    1469         -  ``integer`` - "probably" the rank of self
     1481        -  ``integer`` - lower bound on the rank of self
    14701482       
    14711483        -  ``integer`` - the 2-rank of the Selmer group
    14721484       
    14731485        -  ``list`` - list of independent points on the
    1474            curve.
     1486           quotient `E(Q)/2E(Q)`.
    14751487       
    14761488       
    14771489        IMPLEMENTATION: Uses Denis Simon's GP/PARI scripts from
     
    15411553        """
    15421554        t = simon_two_descent(self, verbose=verbose, lim1=lim1, lim3=lim3, limtriv=limtriv,
    15431555                              maxprob=maxprob, limbigprime=limbigprime)
    1544         prob_rank = rings.Integer(t[0])
     1556        rank_low_bd = rings.Integer(t[0])
    15451557        two_selmer_rank = rings.Integer(t[1])
    1546         prob_gens = [self(P) for P in t[2]]
    1547         return prob_rank, two_selmer_rank, prob_gens
     1558        gens_mod_two = [self(P) for P in t[2]]
     1559        if rank_low_bd == two_selmer_rank - self.two_torsion_rank():
     1560            gens = [P for P in gens_mod_two if P.additive_order() != 2]
     1561            self.__gens[True] = gens
     1562            self.__gens[True].sort()
     1563            self.__rank[True] = len(self.__gens[True])
     1564        return rank_low_bd, two_selmer_rank, gens_mod_two
    15481565
    15491566    two_descent_simon = simon_two_descent
    15501567
     
    15991616
    16001617    def rank(self, use_database=False, verbose=False,
    16011618                   only_use_mwrank=True,
    1602                    algorithm='mwrank_shell',
     1619                   algorithm='mwrank_lib',
    16031620                   proof=None):
    16041621        """
    16051622        Return the rank of this elliptic curve, assuming no conjectures.
     
    17071724
    17081725        if algorithm == 'mwrank_lib':
    17091726            misc.verbose("using mwrank lib")
    1710             C = self.mwrank_curve()
     1727            if self.is_integral(): E = self
     1728            else: E = self.integral_model()
     1729            C = E.mwrank_curve()
    17111730            C.set_verbose(verbose)
    17121731            r = C.rank()
    1713             if not C.certain():
    1714                 del self.__mwrank_curve
    1715                 raise RuntimeError, "Unable to compute the rank with certainty (lower bound=%s).  This could be because Sha(E/Q)[2] is nontrivial."%C.rank() + "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again.  You could also try rank with only_use_mwrank=False."
     1732            if C.certain():
     1733                proof = True
     1734            else:
     1735                if proof:
     1736                    print "Unable to compute the rank with certainty (lower bound=%s)."%C.rank()
     1737                    print "This could be because Sha(E/Q)[2] is nontrivial."
     1738                    print "Try calling something like two_descent(second_limit=13) on the"
     1739                    print "curve then trying this command again.  You could also try rank"
     1740                    print "with only_use_mwrank=False."
     1741                    del E.__mwrank_curve
     1742                    raise RuntimeError, 'Rank not provably correct.'
     1743                else:
     1744                    misc.verbose("Warning -- rank not proven correct", level=1)
    17161745            self.__rank[proof] = r
    17171746        elif algorithm == 'mwrank_shell':
    17181747            misc.verbose("using mwrank shell")
     
    17451774        return self.__rank[proof]
    17461775
    17471776    def gens(self, verbose=False, rank1_search=10,
    1748              algorithm='mwrank_shell',
     1777             algorithm='mwrank_lib',
    17491778             only_use_mwrank=True,
    17501779             proof = None,
    1751              use_database = True):
     1780             use_database = True,
     1781             descent_second_limit=12):
    17521782        """
    17531783        Compute and return generators for the Mordell-Weil group E(Q)
    17541784        *modulo* torsion.
    17551785       
    1756         HINT: If you would like to control the height bounds used in the
    1757         2-descent, first call the two_descent function with those height
    1758         bounds. However that function, while it displays a lot of output,
    1759         returns no values.
    1760        
    1761         TODO: Allow passing of command-line parameters to mwrank.
    1762        
    17631786        .. warning::
    17641787
    17651788           If the program fails to give a provably correct result, it
     
    17931816        -  ``use_database`` - bool (default True) if True,
    17941817           attempts to find curve and gens in the (optional) database
    17951818       
     1819        -  ``descent_second_limit`` - (default: 16)- used in 2-descent
    17961820       
    17971821        OUTPUT:
    17981822       
     
    18851909        # end if (not_use_mwrank)
    18861910        if algorithm == "mwrank_lib":
    18871911            misc.verbose("Calling mwrank C++ library.")
    1888             C = self.mwrank_curve(verbose)
     1912            if not self.is_integral():
     1913                xterm = 1; yterm = 1
     1914                ai = self.a_invariants()
     1915                for a in ai:
     1916                    if not a.is_integral():
     1917                       for p, _ in a.denom().factor():
     1918                          e  = min([(ai[i].valuation(p)/[1,2,3,4,6][i]) for i in range(5)]).floor()
     1919                          ai = [ai[i]/p**(e*[1,2,3,4,6][i]) for i in range(5)]
     1920                          xterm *= p**(2*e)
     1921                          yterm *= p**(3*e)
     1922                E = constructor.EllipticCurve(list(ai))
     1923            else:
     1924                E = self; xterm = 1; yterm = 1
     1925            C = E.mwrank_curve(verbose)
    18891926            if not (verbose is None):
    18901927                C.set_verbose(verbose)
     1928            C.two_descent(verbose=verbose, second_limit=descent_second_limit)
    18911929            G = C.gens()
    18921930            if proof is True and C.certain() is False:
    18931931                del self.__mwrank_curve
    18941932                raise RuntimeError, "Unable to compute the rank, hence generators, with certainty (lower bound=%s, generators found=%s).  This could be because Sha(E/Q)[2] is nontrivial."%(C.rank(),G) + \
    1895                       "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again."
     1933                      "\nTry increasing descent_second_limit then trying this command again."
    18961934            else:
    18971935                proof = C.certain()
     1936            G = [[x*xterm,y*yterm,z] for x,y,z in G]
    18981937        else:
    18991938            # when gens() calls mwrank it passes the command-line
    19001939            # parameter "-p 100" which helps curves with large
     
    19782017        """
    19792018        return len(self.gens(proof = proof))
    19802019
    1981     def regulator(self, use_database=True, proof=None, precision=None):
     2020    def regulator(self, use_database=True, proof=None, precision=None,
     2021                        descent_second_limit=12, verbose=False):
    19822022        """
    19832023        Returns the regulator of this curve, which must be defined over Q.
    19842024       
     
    19952035        -  ``precision`` - int or None (default: None): the
    19962036           precision in bits of the result (default real precision if None)
    19972037       
     2038        -  ``descent_second_limit`` - (default: 16)- used in 2-descent
     2039       
     2040        -  ``verbose`` - whether to print mwrank's verbose output
    19982041       
    19992042        EXAMPLES::
    20002043       
     
    20442087        # Next we find the gens, taking them from the database if they
    20452088        # are there and use_database is True, else computing them:   
    20462089
    2047         G = self.gens(proof=proof, use_database=use_database)
     2090        G = self.gens(proof=proof, use_database=use_database, descent_second_limit=descent_second_limit, verbose=verbose)
    20482091
    20492092        # Finally compute the regulator of the generators found:
    20502093
     
    23482391       
    23492392       
    23502393        -  ``height_limit (float)`` - bound on naive height
    2351            (at most 21,
    2352        
    2353         -  ``or mwrank overflows`` - see below)
    23542394       
    23552395        -  ``verbose (bool)`` - (default: False)
    23562396       
     
    23622402       
    23632403        OUTPUT: points (list) - list of independent points which generate
    23642404        the subgroup of the Mordell-Weil group generated by the points
    2365         found and then p-saturated for p20.
     2405        found and then p-saturated for `p \leq 20`.
    23662406       
    23672407        .. warning::
    23682408
    23692409           height_limit is logarithmic, so increasing by 1 will cause
    23702410           the running time to increase by a factor of approximately
    2371            4.5 (=exp(1.5)). The limit of 21 is to prevent overflow,
    2372            but in any case using height_limit=20 takes rather a long
    2373            time!
    2374        
    2375         IMPLEMENTATION: Uses Cremona's mwrank package. At the heart of this
    2376         function is Cremona's port of Stoll's ratpoints program (version
    2377         1.4).
     2411           4.5 (=exp(1.5)).
     2412       
     2413        IMPLEMENTATION: Uses Michael Stoll's ratpoints library.
    23782414       
    23792415        EXAMPLES::
    23802416       
    23812417            sage: E=EllipticCurve('389a1')
    23822418            sage: E.point_search(5, verbose=False)
    2383             [(0 : -1 : 1), (-1 : 1 : 1)]
     2419            [(-1 : 1 : 1), (-3/4 : 7/8 : 1)]
    23842420       
    23852421        Increasing the height_limit takes longer, but finds no more
    23862422        points::
    23872423       
    23882424            sage: E.point_search(10, verbose=False)
    2389             [(0 : -1 : 1), (-1 : 1 : 1)]
     2425            [(-1 : 1 : 1), (-3/4 : 7/8 : 1)]
    23902426       
    23912427        In fact this curve has rank 2 so no more than 2 points will ever be
    23922428        output, but we are not using this fact.
     
    23942430        ::
    23952431       
    23962432            sage: E.saturation(_)
    2397             ([(0 : -1 : 1), (-1 : 1 : 1)], '1', 0.152460172772408)
     2433            ([(-1 : 1 : 1), (-3/4 : 7/8 : 1)], '1', 0.152460172772408)
    23982434       
    23992435        What this shows is that if the rank is 2 then the points listed do
    24002436        generate the Mordell-Weil group (mod torsion). Finally,
     
    24042440            sage: E.rank()
    24052441            2
    24062442        """
    2407         height_limit = float(height_limit)
    2408         if height_limit > _MAX_HEIGHT:
    2409             raise OverflowError, "height_limit (=%s) must be at most %s."%(height_limit,_MAX_HEIGHT)
    2410         c = self.mwrank_curve()
    2411         mw = mwrank.mwrank_MordellWeil(c, verbose)
    2412         mw.search(height_limit, verbose=verbose)
    2413         v = mw.points()
    2414         return [self(P) for P in v]
     2443        from sage.libs.ratpoints import ratpoints
     2444        from sage.functions.all import exp
     2445        from sage.rings.arith import GCD
     2446        H = exp(float(height_limit)) # max(|p|,|q|) <= H, if x = p/q coprime
     2447        coeffs = [16*self.b6(), 8*self.b4(), self.b2(), 1]
     2448        points = []
     2449        a1 = self.a1()
     2450        a3 = self.a3()
     2451        new_H = H*2 # since we change the x-coord by 2 below
     2452        for X,Y,Z in ratpoints(coeffs, new_H, verbose):
     2453            if Z == 0: continue
     2454            z = 2*Z
     2455            x = X/2
     2456            y = (Y/z - a1*x - a3*z)/2
     2457            d = GCD((x,y,z))
     2458            x = x/d
     2459            if max(x.numerator().abs(), x.denominator().abs()) <= H:
     2460                y = y/d
     2461                z = z/d
     2462                points.append(self((x,y,z)))
     2463        return self.saturation(points, verbose=verbose, max_prime=20)[0]
    24152464   
    24162465    def two_torsion_rank(self):
    24172466        r"""
     
    54565505            x+=1
    54575506        return ans
    54585507
    5459     def prove_BSD(self, verbosity=0, simon=False, proof=None, secs_hi=30):
    5460         r"""
    5461         Attempts to prove the Birch and Swinnerton-Dyer conjectural
    5462         formula for `E`, returning a list of primes `p` for which this
    5463         function fails to prove BSD(E,p).  Here, BSD(E,p) is the
    5464         statement: "the Birch and Swinnerton-Dyer formula holds up to a
    5465         rational number coprime to `p`."
    5466        
    5467         INPUT:
    5468        
    5469             - ``verbosity`` - int, how much information about the proof to print.
    5470            
    5471                 - 0 - print nothing
    5472                 - 1 - print sketch of proof
    5473                 - 2 - print information about remaining primes
    5474            
    5475             - ``simon`` - bool (default False), whether to use two_descent or
    5476               simon_two_descent at p=2.
    5477        
    5478             - ``proof`` - bool or None (default: None, see
    5479               proof.elliptic_curve or sage.structure.proof). If False, this
    5480               function just immediately returns the empty list.
    5481            
    5482             - ``secs_hi`` - maximum number of seconds to try to compute the
    5483               Heegner index before switching over to trying to compute the
    5484               Heegner index bound. (Rank 0 only!)
    5485 
    5486         NOTE:
    5487        
    5488         When printing verbose output, phrases such as "by Mazur" are referring
    5489         to the following list of papers:
    5490 
    5491         REFERENCES:
    5492        
    5493         .. [Cha] B. Cha. Vanishing of some cohomology goups and bounds for the
    5494            Shafarevich-Tate groups of elliptic curves. J. Number Theory, 111:154-
    5495            178, 2005.
    5496         .. [Jetchev] D. Jetchev. Global divisibility of Heegner points and
    5497            Tamagawa numbers. Compos. Math. 144 (2008), no. 4, 811--826.
    5498         .. [Kato] K. Kato. p-adic Hodge theory and values of zeta functions of
    5499            modular forms. Astérisque, (295):ix, 117-290, 2004.
    5500         .. [Kolyvagin] V. A. Kolyvagin. On the structure of Shafarevich-Tate
    5501            groups. Algebraic geometry, 94--121, Lecture Notes in Math., 1479,
    5502            Springer, Berlin, 1991.
    5503         .. [LumStein] A. Lum, W. Stein. Verification of the Birch and
    5504            Swinnerton-Dyer Conjecture for Elliptic Curves with Complex
    5505            Multiplication (unpublished)
    5506         .. [Mazur] B. Mazur. Modular curves and the Eisenstein ideal. Inst.
    5507            Hautes Études Sci. Publ. Math. No. 47 (1977), 33--186 (1978).
    5508         .. [Rubin] K. Rubin. The "main conjectures" of Iwasawa theory for
    5509            imaginary quadratic fields. Invent. Math. 103 (1991), no. 1, 25--68.
    5510         .. [SteinWuthrich] W. Stein and C. Wuthrich. Computations about
    5511            Tate-Shafarevich groups using Iwasawa theory.
    5512            http://wstein.org/papers/shark, February 2008.
    5513         .. [SteinEtAl] G. Grigorov, A. Jorza, S. Patrikis, W. Stein,
    5514            C. Tarniţǎ. Computational verification of the Birch and
    5515            Swinnerton-Dyer conjecture for individual elliptic curves.
    5516            Math. Comp. 78 (2009), no. 268, 2397--2425.
    5517 
    5518 
    5519         EXAMPLES::
    5520        
    5521             sage: EllipticCurve('11a').prove_BSD(verbosity=2)   
    5522             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5523             True for p not in {2, 5} by Kolyvagin.
    5524             True for p=5 by Mazur
    5525             Remaining primes:
    5526             p = 2: irreducible, surjective, good, non-ordinary
    5527             [2]
    5528            
    5529             sage: EllipticCurve('14a').prove_BSD(verbosity=2)   
    5530             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5531             True for p not in {2, 3} by Kolyvagin.
    5532             Remaining primes:
    5533             p = 2: reducible, surjective, non-split multiplicative, divides a Tamagawa number
    5534             p = 3: reducible, surjective, good ordinary, divides a Tamagawa number
    5535             [2, 3]
    5536 
    5537         A rank two curve::
    5538        
    5539             sage: E = EllipticCurve('389a')
    5540 
    5541         We know nothing with proof=True::
    5542        
    5543             sage: E.prove_BSD()
    5544             Set of all prime numbers: 2, 3, 5, 7, ...
    5545 
    5546         We (think we) know everything with proof=False::
    5547        
    5548             sage: E.prove_BSD(proof=False)
    5549             []
    5550 
    5551         A curve of rank 0 and prime conductor::
    5552        
    5553             sage: E = EllipticCurve('19a')
    5554             sage: E.prove_BSD(verbosity=2)
    5555             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5556             True for p not in {2, 3} by Kolyvagin.
    5557             True for p=3 by Mazur
    5558             Remaining primes:
    5559             p = 2: irreducible, surjective, good, non-ordinary
    5560             [2]
    5561 
    5562             sage: E = EllipticCurve('37a')
    5563             sage: E.rank()
    5564             1
    5565             sage: E._EllipticCurve_rational_field__rank
    5566             {True: 1}
    5567             sage: E._EllipticCurve_rational_field__rank = {True:0}
    5568             sage: E.prove_BSD()
    5569             Traceback (most recent call last):
    5570             ...
    5571             RuntimeError: It seems that the rank conjecture does not hold for this curve (Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field)! This may be a counterexample to BSD, but is more likely a bug.
    5572 
    5573         We check the error message indicating that this code doesn't
    5574         yet use Simon 2-descent in case of rational 2-torsion (though
    5575         it could with a little more work)::
    5576        
    5577             sage: E = EllipticCurve('14a')
    5578             sage: E.prove_BSD(simon=True)
    5579             Traceback (most recent call last):
    5580             ...
    5581             RuntimeError: Simon two-descent only valid for curves without two torsion.
    5582 
    5583         We test the consistency check for the 2-part of Sha::
    5584 
    5585             sage: E = EllipticCurve('37a')
    5586             sage: S = E.sha(); S
    5587             Shafarevich-Tate group for the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
    5588             sage: S.an = lambda : 4
    5589             sage: E.prove_BSD()
    5590             [2]
    5591 
    5592         An example with a Tamagawa number at 5::
    5593 
    5594             sage: E = EllipticCurve('123a1')
    5595             sage: E.prove_BSD(verbosity=2)
    5596             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5597             True for p not in {2, 5} by Kolyvagin.
    5598             Remaining primes:
    5599             p = 2: irreducible, surjective, good, non-ordinary
    5600             p = 5: reducible, surjective, good ordinary, divides a Tamagawa number
    5601             [2, 5]
    5602 
    5603         A curve for which 3 divides the order of the Shafarevich-Tate group::
    5604 
    5605             sage: E = EllipticCurve('681b')                       
    5606             sage: E.prove_BSD(verbosity=2)
    5607             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5608             True for p not in {2, 3} by Kolyvagin.
    5609             ALERT: p = 3 left in Kolyvagin bound
    5610                 0 <= ord_p(#Sha) <= 2
    5611                 ord_p(#Sha_an) = 2
    5612             Remaining primes:
    5613             p = 2: reducible, surjective, good ordinary, divides a Tamagawa number
    5614             p = 3: irreducible, surjective, non-split multiplicative
    5615             [2, 3]
    5616        
    5617         A curve for which we need to use ``heegner_index_bound``::
    5618 
    5619             sage: E = EllipticCurve('198b')
    5620             sage: E.prove_BSD(verbosity=1, secs_hi=1)
    5621             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5622             Timeout stopped Heegner index computation...
    5623             Proceeding to use heegner_index_bound instead.
    5624             True for p not in {2, 3} by Kolyvagin.
    5625             [2, 3]
    5626 
    5627         TESTS::
    5628        
    5629             sage: E = EllipticCurve('438e1')
    5630             sage: E.prove_BSD(verbosity=1)
    5631             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5632             True for p not in {2} by Kolyvagin.
    5633             [2]
    5634        
    5635         ::
    5636        
    5637             sage: E = EllipticCurve('960d1')
    5638             sage: E.prove_BSD(verbosity=1) # long time
    5639             p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank
    5640             Timeout stopped Heegner index computation...
    5641             Proceeding to use heegner_index_bound instead.
    5642             True for p not in {2} by Kolyvagin.
    5643             [2]
    5644        
    5645         """
    5646         if proof is None:
    5647             from sage.structure.proof.proof import get_flag
    5648             proof = get_flag(proof, "elliptic_curve")
    5649         else:
    5650             proof = bool(proof)
    5651         if not proof:
    5652             return []
    5653         two_tor_rk = self.two_torsion_rank()
    5654         if simon:
    5655             if two_tor_rk > 0:
    5656                 raise RuntimeError("Simon two-descent only valid for curves without two torsion.")
    5657             rank_lower_bd, two_sel_rk = self.simon_two_descent()[:2]
    5658             if rank_lower_bd == two_sel_rk:
    5659                 rank = rank_lower_bd
    5660             else:
    5661                 raise RuntimeError("Rank can't be computed precisely using Simon's program.")
    5662         else:
    5663             two_sel_rk_bd = self.rank_bound()
    5664             rank = self.rank()
    5665         if rank > 1:
    5666             # We do not know BSD(E,p) for even a single p, since it's
    5667             # an open problem to show that L^r(E,1)/(Reg*Omega) is
    5668             # rational for any curve with r >= 2.
    5669             from sage.sets.all import Primes
    5670             return Primes()
    5671         if rank != self.analytic_rank():
    5672             raise RuntimeError("It seems that the rank conjecture does not hold for this curve (%s)! This may be a counterexample to BSD, but is more likely a bug."%(self))
    5673 
    5674         # We replace self by the optimal curve, which we can do since
    5675         # truth of BSD(E,p) is invariant under isogeny.
    5676         self = self.optimal_curve()
    5677        
    5678         Sha = self.sha()
    5679         sha_an = Sha.an()
    5680         N = self.conductor()
    5681        
    5682         # p = 2
    5683         if two_sel_rk_bd > rank:
    5684             if verbosity > 0:
    5685                 print 'p = 2: mwrank did not achieve a tight bound on the Selmer rank.'
    5686             two_proven = False
    5687         elif two_sel_rk_bd < rank:
    5688             raise RuntimeError("MWrank seems to have computed an incorrect upper bound of %d on the rank."%two_sel_rk_bd)
    5689         else:
    5690             # until we can easily access the computed rank of Sha[2]:
    5691             two_proven = False
    5692             if verbosity > 0:
    5693                 print 'p = 2: Unverified since it is difficult to access the rank bound for Sha[2] computed by MWrank'
    5694            
    5695         # reduce set of remaining primes to a finite set
    5696         import signal
    5697         remaining_primes = []
    5698         kolyvagin_primes = []
    5699         heegner_index = None
    5700         if self.rank() == 0:
    5701             try:
    5702                 old_alarm = signal.alarm(secs_hi)
    5703                 old_alarm_set = (old_alarm != 0)
    5704                 for D in self.heegner_discriminants_list(10):
    5705                     I = self.heegner_index(D)
    5706                     J = I.is_int()
    5707                     if J[0] and J[1]>0:
    5708                         I = J[1]
    5709                     else:
    5710                         J = (2*I).is_int()
    5711                         if J[0] and J[1]>0:
    5712                             I = J[1]
    5713                         else:
    5714                             I = None
    5715                     if I is not None:
    5716                         if heegner_index is None:
    5717                             heegner_index = I
    5718                             break # no big long loops just yet...
    5719                 old_alarm_sub = signal.alarm(0)
    5720                 if old_alarm_set:
    5721                     old_alarm -= old_alarm_sub
    5722             except KeyboardInterrupt:
    5723                 if signal.alarm(0)==0:
    5724                     print 'Timeout stopped Heegner index computation...'
    5725                     print 'Proceeding to use heegner_index_bound instead.'
    5726                 else:
    5727                     raise KeyboardInterrupt
    5728             if old_alarm_set: # in case alarm was already set...
    5729                 if old_alarm <= 0:
    5730                     raise KeyboardInterrupt
    5731                 signal.alarm(old_alarm)
    5732             if heegner_index is None:
    5733                 for D in self.heegner_discriminants_list(100):
    5734                     heegner_primes, _ = self.heegner_index_bound(D)
    5735                     if isinstance(heegner_primes, list):
    5736                         break
    5737                 if not isinstance(heegner_primes, list):
    5738                     raise RuntimeError("Tried 100 Heegner discriminants, and heegner_index_bound failed each time.")
    5739                 if 2 in heegner_primes:
    5740                     heegner_primes.remove(2)
    5741             else:
    5742                 heegner_primes = [p for p in arith.prime_divisors(heegner_index) if p!=2]
    5743         else: # rank 1
    5744             for D in self.heegner_discriminants_list(10):
    5745                 I = self.heegner_index(D)
    5746                 J = I.is_int()
    5747                 if J[0] and J[1]>0:
    5748                     I = J[1]
    5749                 else:
    5750                     J = (2*I).is_int()
    5751                     if J[0] and J[1]>0:
    5752                         I = J[1]
    5753                     else:
    5754                         continue
    5755                 heegner_index = I
    5756                 break
    5757             heegner_primes = [p for p in arith.prime_divisors(heegner_index) if p!=2]
    5758        
    5759         if self.has_cm():
    5760             # ensure that CM is by a maximal order
    5761             non_max_j_invs = [-12288000, 54000, 287496, 16581375]
    5762             if self.j_invariant() in non_max_j_invs:
    5763                 for E in self.isogeny_class()[0]:
    5764                     if E.j_invariant() not in non_max_j_invs:
    5765                         Sha = E.sha()
    5766                         sha_an = Sha.an()
    5767                         if verbosity > 0:
    5768                             print 'CM by non maximal order: switching curves'
    5769                         break
    5770             else:
    5771                 E = self
    5772             if E.analytic_rank() == 0:
    5773                 if verbosity > 0:
    5774                     print 'p >= 5: true by Rubin'
    5775                 remaining_primes.append(3)
    5776             else:
    5777                 K = rings.QuadraticField(E.cm_discriminant(), 'a')
    5778                 D_K = K.disc()
    5779                 D_E = E.discriminant()
    5780                 if len(K.factor(3)) == 1: # 3 does not split in K
    5781                     remaining_primes.append(3)
    5782                 for p in arith.prime_divisors(D_K):
    5783                     if p >= 5:
    5784                         remaining_primes.append(p)
    5785                 for p in arith.prime_divisors(D_E):
    5786                     if p >= 5 and D_K%p and len(K.factor(p)) == 1: # p is inert in K
    5787                         remaining_primes.append(p)
    5788                 for p in heegner_primes:
    5789                     if p >= 5 and D_E%p != 0 and D_K%p != 0 and len(K.factor(p)) == 1: # p is good for E and inert in K
    5790                         kolyvagin_primes.append(p)
    5791                 assert sha_an in ZZ and sha_an > 0
    5792                 for p in arith.prime_divisors(sha_an):
    5793                     if p >= 5 and D_K%p != 0 and len(K.factor(p)) == 1:
    5794                         if E.is_good(p):
    5795                             if verbosity > 2 and p in heegner_primes and heegner_index is None:
    5796                                 print 'ALERT: Prime p (%d) >= 5 dividing sha_an, good for E, inert in K, in heegner_primes, should not divide the actual Heegner index'
    5797                             # Note that the following check is not entirely
    5798                             # exhaustive, in case there is a p not dividing
    5799                             # the Heegner index in heegner_primes,
    5800                             # for which only an outer bound was computed
    5801                             if p not in heegner_primes:
    5802                                 raise RuntimeError("p = %d divides sha_an, is of good reduction for E, inert in K, and does not divide the Heegner index. This may be a counterexample to BSD, but is more likely a bug. %s"%(p,self))
    5803                 if verbosity > 0:
    5804                     print 'True for p not in {%s} by Kolyvagin (via Stein & Lum -- unpublished) and Rubin.'%str(list(set(remaining_primes).union(set(kolyvagin_primes))))[1:-1]
    5805         else: # no CM
    5806             E = self
    5807             # do some tricks to get to a finite set without calling bound_kolyvagin
    5808             remaining_primes = [p for p, reason in E.non_surjective()]
    5809             for p in heegner_primes:
    5810                 if p not in remaining_primes:
    5811                     remaining_primes.append(p)
    5812             assert sha_an in ZZ and sha_an > 0
    5813             for p in arith.prime_divisors(sha_an):
    5814                 if p not in remaining_primes:
    5815                     remaining_primes.append(p)
    5816             if 2 in remaining_primes: remaining_primes.remove(2)
    5817             if verbosity > 0:
    5818                 print 'True for p not in {' + str([2]+list(remaining_primes))[1:-1] + '} by Kolyvagin.'
    5819             primes_to_remove = []
    5820             for p in remaining_primes:
    5821                 if E.is_surjective(p)[0] and not E.has_additive_reduction(p):
    5822                     if E.has_nonsplit_multiplicative_reduction(p):
    5823                         if E.rank() > 0:
    5824                             continue
    5825                     if p==3:
    5826                         if (not (E.is_ordinary(p) and E.is_good(p))) and (not E.has_split_multiplicative_reduction(p)):
    5827                             continue
    5828                         if E.rank() > 0:
    5829                             continue
    5830                     p_bound = Sha.p_primary_bound(p)
    5831                     if sha_an.ord(p) == 0 and p_bound == 0:
    5832                         if verbosity > 0:
    5833                             print 'True for p=%d by Stein-Wuthrich.'%p
    5834                         primes_to_remove.append(p)
    5835                     else:
    5836                         print 'Analytic %d-rank is '%p + str(sha_an.ord(p)) + ', actual %d-rank is at most %d.'%(p, p_bound)
    5837                         print '    by Stein-Wuthrich.\n'
    5838             for p in primes_to_remove:
    5839                 remaining_primes.remove(p)
    5840             kolyvagin_primes = []
    5841             for p in remaining_primes:
    5842                 if E.is_surjective(p)[0]:
    5843                     kolyvagin_primes.append(p)
    5844             for p in kolyvagin_primes:
    5845                 remaining_primes.remove(p)
    5846         # apply other hypotheses which imply Kolyvagin's bound holds
    5847         bounded_primes = []
    5848         D_K = rings.QuadraticField(D, 'a').disc()
    5849         assert 2 not in remaining_primes
    5850         # Cha's hypothesis
    5851         for p in remaining_primes:
    5852             if D_K%p != 0 and N%(p**2) != 0 and E.is_irreducible(p):
    5853                 if verbosity > 0:
    5854                     print 'Kolyvagin\'s bound for p = %d applies by Cha.'%p
    5855                 kolyvagin_primes.append(p)
    5856         # Stein et al.
    5857         if not E.has_cm():
    5858             L = arith.lcm([F.torsion_order() for F in E.isogeny_class()[0]])
    5859             for p in remaining_primes:
    5860                 if p in kolyvagin_primes: continue
    5861                 if L%p != 0:
    5862                     if len(arith.prime_divisors(D_K)) == 1:
    5863                         if D_K%p != 0:
    5864                             if verbosity > 0:
    5865                                 print 'Kolyvagin\'s bound for p = %d applies by Stein et al.'%p
    5866                             kolyvagin_primes.append(p)
    5867                     else:
    5868                         if verbosity > 0:
    5869                             print 'Kolyvagin\'s bound for p = %d applies by Stein et al.'%p
    5870                         kolyvagin_primes.append(p)
    5871         for p in kolyvagin_primes:
    5872             if p in remaining_primes:
    5873                 remaining_primes.remove(p)
    5874 
    5875         prime_bounds = []
    5876         # apply Kolyvagin's bound
    5877         primes_to_remove = []
    5878         for p in kolyvagin_primes:
    5879             if sha_an.ord(p) == 0 and p not in heegner_primes:
    5880                     if verbosity > 0:
    5881                         print 'True for p = %d by Kolyvagin bound.'%p
    5882                     primes_to_remove.append(p)
    5883                     continue
    5884             if heegner_index is not None: # p must divide heegner_index
    5885                 ord_p_bound = 2*heegner_index.ord(p)
    5886                 # Here Jetchev's results apply.
    5887                 m_max = max([E.tamagawa_number(q).ord(p) for q in N.prime_divisors()])
    5888                 if m_max > 0 and verbosity > 0:
    5889                     print 'Jetchev\'s results apply with m_max =', m_max
    5890                 ord_p_bound -= 2*m_max
    5891                 if ord_p_bound == 0:
    5892                     if sha_an.ord(p) != 0:
    5893                         raise RuntimeError("p = %d: ord_p_bound == 0, but sha_an.ord(p) == %d. This appears to be a counterexample to BSD, but is more likely a bug."%(p,sha_an.ord(p)))
    5894                     if verbosity > 0:
    5895                         print 'True for p = %d by Kolyvagin bound.'%p
    5896                     primes_to_remove.append(p)
    5897                     continue
    5898             elif p not in heegner_primes:
    5899                 ord_p_bound = 0
    5900             else:
    5901                 from sage.rings.infinity import Infinity
    5902                 ord_p_bound = Infinity
    5903                 if verbosity > 0:
    5904                     print 'p = %d may divide the Heegner index, for which only a bound was computed.'%p
    5905             if verbosity > 0:
    5906                 print 'ALERT: p = %d left in Kolyvagin bound'%p
    5907                 print '    0 <= ord_p(#Sha) <=', ord_p_bound
    5908                 print '    ord_p(#Sha_an) =', sha_an.ord(p)
    5909         for p in primes_to_remove:
    5910             kolyvagin_primes.remove(p)
    5911         remaining_primes = list( set(remaining_primes).union(set(kolyvagin_primes)) )
    5912        
    5913         # Kato's bound
    5914         if rank == 0 and not E.has_cm():
    5915             L_over_Omega = E.lseries().L_ratio()
    5916             kato_primes = Sha.bound_kato()
    5917             primes_to_remove = []
    5918             for p in remaining_primes:
    5919                 if p not in kato_primes:
    5920                     if verbosity > 0:
    5921                         print 'Kato further implies that #Sha[%d] is trivial.'%p
    5922                     primes_to_remove.append(p)
    5923                 if p not in [2,3] and N%p != 0:
    5924                     if E.is_surjective(p)[0]:
    5925                         if verbosity > 1:
    5926                             print 'Kato might apply nontrivially for %d'%p
    5927                         # ordp(sha) <= ordp(L_over_omega)
    5928             for p in primes_to_remove:
    5929                 remaining_primes.remove(p)
    5930        
    5931         # Mazur
    5932         if N.is_prime():
    5933             for p in remaining_primes:
    5934                 if E.is_reducible(p):
    5935                     remaining_primes.remove(p)                   
    5936                     if verbosity > 0:
    5937                         print 'True for p=%s by Mazur'%p
    5938        
    5939         if two_proven is False:
    5940             remaining_primes.append(2)
    5941        
    5942         # print some extra information
    5943         remaining_primes.sort()
    5944         if verbosity > 1:
    5945             if len(remaining_primes) > 0:
    5946                 print 'Remaining primes:'
    5947             for p in remaining_primes:
    5948                 s = 'p = ' + str(p) + ': '
    5949                 if not E.is_reducible(p):
    5950                     s += 'ir'
    5951                 s += 'reducible, '
    5952                 if not E.is_surjective(p):
    5953                     s += 'not '
    5954                 s += 'surjective, '
    5955                 a_p = E.an(p)
    5956                 if E.is_good(p):
    5957                     if a_p%p != 0:
    5958                         s += 'good ordinary'
    5959                     else:
    5960                         s += 'good, non-ordinary'
    5961                 else:
    5962                     assert E.is_minimal()
    5963                     if a_p == 0:
    5964                         s += 'additive'
    5965                     elif a_p == 1:
    5966                         s += 'split multiplicative'
    5967                     elif a_p == -1:
    5968                         s += 'non-split multiplicative'
    5969                 if E.tamagawa_product()%p==0:
    5970                     s += ', divides a Tamagawa number'
    5971                 print s
    5972 
    5973         return remaining_primes
     5508    prove_BSD = BSD.prove_BSD
    59745509
    59755510    def integral_points(self, mw_base='auto', both_signs=False, verbose=False):
    59765511        """
  • sage/schemes/elliptic_curves/heegner.py

    diff -r 1ca8aa94810a -r 8a90b87fca0c sage/schemes/elliptic_curves/heegner.py
    a b  
    63296329        return IR(alpha-MIN_ERR,alpha+MIN_ERR) * IR(LE1-err_E,LE1+err_E) * IR(LF1-err_F,LF1+err_F)
    63306330
    63316331
    6332 def heegner_index(self, D,  min_p=2, prec=5):
     6332def heegner_index(self, D,  min_p=2, prec=5, descent_second_limit=16, verbose_mwrank=False):
    63336333    r"""
    63346334    Return an interval that contains the index of the Heegner
    63356335    point `y_K` in the group of `K`-rational points modulo torsion
     
    63526352    -  ``min_p (int)`` - (default: 2) only rule out primes
    63536353       = min_p dividing the index.
    63546354
    6355     -  ``verbose (bool)`` - (default: False); print lots of
     6355    -  ``verbose_mwrank (bool)`` - (default: False); print lots of
    63566356       mwrank search status information when computing regulator
    63576357
    63586358    -  ``prec (int)`` - (default: 5), use prec\*sqrt(N) +
    63596359       20 terms of L-series in computations, where N is the conductor.
    63606360
     6361    -  ``descent_second_limit`` - (default: 16)- used in 2-descent
     6362       when computing regulator of the twist
    63616363
    63626364    OUTPUT: an interval that contains the index
    63636365
     
    64066408    index by `2`. Unfortunately, this is not an if and only if
    64076409    condition, i.e., sometimes the index must be multiplied by
    64086410    `2` even though the denominator is not `2`.
     6411
     6412    This example demonstrates the `descent_second_limit` option,
     6413    which can be used to fine tune the 2-descent used to compute
     6414    the regulator of the twist. If we set the parameter lower than
     6415    its usual value, then the point search is not high enough to
     6416    find what it is looking for::
     6417   
     6418        sage: E = EllipticCurve([0, 0, 1, -34874, -2506691])
     6419        sage: E.heegner_index(-8, descent_second_limit=10)
     6420        Traceback (most recent call last):
     6421        ...
     6422        RuntimeError: ...
     6423
     6424    However when we use the default values, we find the points we need::
     6425
     6426        sage: E.heegner_index(-8)
     6427        1.00000?
     6428
    64096429    """
    64106430    # First compute upper bound on height of Heegner point.
    64116431    tm = verbose("computing heegner point height...")
     
    64346454   
    64356455    if c > _MAX_HEIGHT or F is self:
    64366456        verbose("Doing direct computation of MW group.")
    6437         reg = F.regulator()
     6457        reg = F.regulator(descent_second_limit=descent_second_limit, verbose=verbose_mwrank)
    64386458        return self._adjust_heegner_index(ht/IR(reg))
    64396459
    64406460    # Do naive search to eliminate possibility that Heegner point
  • sage/schemes/elliptic_curves/sha_tate.py

    diff -r 1ca8aa94810a -r 8a90b87fca0c sage/schemes/elliptic_curves/sha_tate.py
    a b  
    161161        self.__an_numerical = Sha
    162162        return Sha
    163163       
    164     def an(self, use_database=False):
     164    def an(self, use_database=False, descent_second_limit=12):
    165165        r"""
    166166        Returns the Birch and Swinnerton-Dyer conjectural order of Sha
    167167        as a provably correct integer, unless the analytic rank is > 1,
    168168        in which case this function returns a numerical value.
    169169
    170         INPUT: ``use_database`` -- bool (default: False); if True, try to use any
    171         databases installed to lookup the analytic order of Sha, if
    172         possible.  The order of Sha is computed if it can't be looked up.
     170        INPUT:
     171       
     172            - ``use_database`` -- bool (default: False); if True, try to use any
     173              databases installed to lookup the analytic order of Sha, if
     174              possible.  The order of Sha is computed if it can't be looked up.
     175           
     176            - ``descent_second_limit`` -- int (default: 12); limit to use on
     177              point searching for the quartic twist in the hard case
    173178
    174179        This result is proved correct if the order of vanishing is 0
    175180        and the Manin constant is <= 2.
     
    179184        listed in Cremona's tables, if this curve appears in Cremona's
    180185        tables.
    181186
     187        NOTE:
     188       
     189        If you come across the following error::
     190       
     191            sage: E = EllipticCurve([0, 0, 1, -34874, -2506691])
     192            sage: E.sha().an()
     193            Traceback (most recent call last):
     194            ...
     195            RuntimeError: Unable to compute the rank, hence generators, with certainty (lower bound=0, generators found=[]).  This could be because Sha(E/Q)[2] is nontrivial.
     196            Try increasing descent_second_limit then trying this command again.
     197
     198        You can increase the `descent_second_limit` (in the above example,
     199        set to the default, 12) option to try again::
     200
     201            sage: E.sha().an(descent_second_limit=16)
     202            1
     203
    182204        EXAMPLES::
    183205
    184206            sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11)
     
    275297                self.__an = s
    276298                return s
    277299           
    278             regulator = E.regulator(use_database=use_database)   # this could take a *long* time; and could fail...?
     300            regulator = E.regulator(use_database=use_database, descent_second_limit=descent_second_limit)
    279301            T = E.torsion_subgroup().order()
    280302            omega = E.period_lattice().omega()
    281303            Sha = int(round ( (L1 * T * T) / (E.tamagawa_product() * regulator * omega) ))