Ticket #8829: 8829-ec-nf-sat.patch

File 8829-ec-nf-sat.patch, 20.5 KB (added by robertwb, 12 years ago)
  • sage/rings/residue_field.pyx

    # HG changeset patch
    # User Robert Bradshaw <robertwb@math.washington.edu>
    # Date 1272616711 25200
    # Node ID 44a69cc310af7b00c3889fa6205d39a176c7ce5d
    # Parent  48a42f04011f7ee195bb63ef221f708422bb30b0
    #8829 - Saturation in elliptic curves over number fields, gens() for E(K) where E/Q, [K:Q] = 2.
    
    diff -r 48a42f04011f -r 44a69cc310af sage/rings/residue_field.pyx
    a b  
    980980        try:
    981981            return FiniteField_givaro.__call__(self, x)
    982982        except TypeError:
    983             try:
    984                 return self._structure[0](x)
    985             except:
    986                 raise TypeError
    987 
     983            return self._structure[0](x)
  • sage/schemes/elliptic_curves/ell_number_field.py

    diff -r 48a42f04011f -r 44a69cc310af sage/schemes/elliptic_curves/ell_number_field.py
    a b  
    8989#                  http://www.gnu.org/licenses/
    9090#*****************************************************************************
    9191
     92import math
     93
    9294from ell_field import EllipticCurve_field
    9395import ell_point
    9496from sage.rings.ring import Ring
     
    100102
    101103from gp_simon import simon_two_descent
    102104from constructor import EllipticCurve
    103 from sage.rings.all import PolynomialRing, QQ, ZZ, is_Ideal, is_NumberFieldElement, is_NumberFieldFractionalIdeal,is_NumberField, GF, prime_range
     105from sage.rings.all import (PolynomialRing, QQ, ZZ, RR, RDF, GF,
     106    is_Ideal, is_NumberFieldElement, is_NumberFieldFractionalIdeal,
     107    is_NumberField, prime_range, primes)
    104108from sage.misc.misc import verbose, forall
    105109from sage.misc.functional import ideal
    106110from kodaira_symbol import KodairaSymbol
     
    108112from sage.structure.element import RingElement
    109113from sage.rings.infinity import Infinity # just for verbose output
    110114from sage.rings.arith import valuation
     115from sage.matrix.all import matrix
     116from sage.modules.all import vector
    111117
    112118class EllipticCurve_number_field(EllipticCurve_field):
    113119    r"""
     
    259265        two_selmer_rank = Integer(t[1])
    260266        prob_gens = [self(P) for P in t[2]]
    261267        return prob_rank, two_selmer_rank, prob_gens
     268   
     269    def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False, reg_bound=None, reg=None, debug=False):
     270        r"""
     271        Given a list of K-rational points on E, compute the saturation in
     272        E(K) of the subgroup they generate.
     273       
     274        INPUT:
     275       
     276       
     277        -  ``points (list)`` - list of points on E
     278       
     279        -  ``verbose (bool)`` - (default: False), if True, give
     280           verbose output
     281       
     282        -  ``max_prime (int)`` - (default: 0), saturation is
     283           performed for all primes up to max_prime. If max_prime==0,
     284           perform saturation at *all* primes, i.e., compute the true
     285           saturation.
     286       
     287        -  ``odd_primes_only (bool)`` - (default: False), only do saturation
     288           at odd primes
     289           
     290        The following two inputs are optional, and may be provided to speed
     291        up the computation.
     292       
     293        -  ``reg_bound (real)`` - (default: None), lower bound in the
     294           regulator E(K), if known
     295       
     296        -  ``reg (real)`` - (default: None), regulator of the span of points,
     297           if known
     298       
     299        - ``debug (int)`` - (default: 0), used for debugging and testing
     300       
     301       
     302        OUTPUT:
     303       
     304       
     305        -  ``saturation (list)`` - points that form a basis for
     306           the saturation
     307       
     308        -  ``index (int)`` - the index of the group generated
     309           by points in their saturation
     310       
     311        -  ``regulator (real with default precision or None)`` -
     312           regulator of saturated points
     313
     314        EXAMPLES::
     315       
     316            sage: K.<i> = QuadraticField(-1)
     317            sage: E = EllipticCurve('389a')
     318            sage: P,Q = E.gens()
     319            sage: EK = E.change_ring(K)
     320
     321            sage: EK.saturation([2*P], max_prime=2)
     322            ([(-1 : 1 : 1)], 2, 0.686667083306)
     323            sage: EK.saturation([12*P], max_prime=2)
     324            ([(26/361 : -5720/6859 : 1)], 4, 6.18000374975)
     325            sage: EK.saturation([12*P], reg_bound=0.1)
     326            ([(-1 : 1 : 1)], 12, 0.686667083306)
     327           
     328            sage: EK.saturation([2*P, Q], max_prime=2)
     329            ([(-1 : 1 : 1), (0 : -1 : 1)], 2, 0.152460177943)
     330            sage: EK.saturation([P+Q, P-Q], reg_bound=.1, debug=2)
     331            ([(4 : 8 : 1), (-1 : 1 : 1)], 2, 0.152460177943)
     332            sage: EK.saturation([P+Q, 17*Q], reg_bound=0.1)
     333            ([(4 : 8 : 1), (0 : -1 : 1)], 17, 0.152460177943)
     334
     335            sage: P, Q, R = EK.gens()
     336            sage: print P, Q, R
     337            (-1 : 1 : 1) (0 : -1 : 1) (i - 2 : -i - 3 : 1)
     338            sage: EK.saturation([P+Q, Q+R, R+P], reg_bound=0.1)
     339            ([(4 : 8 : 1), (-1/25*i + 18/25 : -69/125*i - 58/125 : 1), (841/1369*i - 171/1369 : 41334/50653*i - 74525/50653 : 1)], 2, 0.103174443217)
     340            sage: EK.saturation([26*R], reg_bound=0.1)
     341            ([(i - 2 : -i - 3 : 1)], 26, 1.06438645309)
     342           
     343        Another number field::
     344           
     345            sage: E = EllipticCurve('389a')
     346            sage: P, Q = E.gens()
     347            sage: K.<a> = NumberField(x^3-x+1)
     348            sage: EK = E.change_ring(K)
     349            sage: EK.saturation([P+Q, P-Q], reg_bound=0.1)
     350            ([(4 : 8 : 1), (-1 : 1 : 1)], 2, 0.152460177943)
     351            sage: EK.saturation([3*P, P+5*Q], reg_bound=0.1)
     352            ([(-1 : 1 : 1), (4 : 8 : 1)], 15, 0.152460177943)
     353
     354        A different curve::
     355
     356            sage: K.<a> = QuadraticField(3)
     357            sage: E = EllipticCurve('37a')
     358            sage: EK = EllipticCurve('37a').change_ring(K)
     359            sage: P, Q = EK.gens()
     360            sage: EK.saturation([3*P-Q, 3*P+Q], reg_bound=.01)
     361            ([(1/2*a : 1/4*a - 3/4 : 1), (0 : -1 : 1)], 6, 0.0317814530726)
     362
     363        The points must be linearly independent::
     364       
     365            sage: EK.saturation([2*P, 3*Q, P-Q])
     366            Traceback (most recent call last):
     367            ...
     368            ValueError: points not linearly independent
     369       
     370        Degenerate case::
     371           
     372            sage: EK.saturation([])
     373            ([], 1, 1.00000000000000)
     374       
     375        ALGORITHM:
     376       
     377            For rank 1 subgroups, simply do trial divison up to the maximal
     378            prime divisor. For higher rank subgroups, perform trial divison
     379            on all linear combinations for small primes, and look for
     380            projections $E(K) \rightarrow \oplus E(k) \otimes \FF_p$ which
     381            are either full rank or provide p-divisble linear combinations,
     382            where the $k$ here are residue fields of $K$.
     383           
     384        TESTS::
     385       
     386            sage: K.<i> = QuadraticField(-1)
     387            sage: E = EllipticCurve('389a')
     388            sage: P,Q = E.gens()
     389            sage: EK = E.change_ring(K)
     390
     391            sage: EK.saturation([P+Q, P-Q], reg_bound=.1, debug=2)
     392            ([(4 : 8 : 1), (-1 : 1 : 1)], 2, 0.152460177943)
     393            sage: EK.saturation([5*P+6*Q, 5*P-3*Q], reg_bound=.1)
     394            ([(-1696464703914542290/25600964182280974369 : -145311084734320503143074826065/129534210658407920857712927153 : 1), (4 : 8 : 1)], 45, 0.152460177943)
     395            sage: EK.saturation([5*P+6*Q, 5*P-3*Q], reg_bound=.1, debug=2)
     396            ([(-1696464703914542290/25600964182280974369 : -145311084734320503143074826065/129534210658407920857712927153 : 1), (4 : 8 : 1)], 45, 0.152460177943)
     397            """
     398        # This code does a lot of residue field construction and elliptic curve
     399        # group structure computations, and would benifit immensly if those
     400        # were sped up.
     401        full_saturation = max_prime == 0
     402        span = [self(P) for P in points]
     403        n = len(span)
     404       
     405        if len(span) == 0:
     406            return span, ZZ(1), RR(1)
     407
     408        if reg is None:
     409            mat = matrix(RDF, n)
     410            heights = [P.height() for P in span]
     411            for i in range(n):
     412                for j in range(n):
     413                    if i == j:
     414                        mat[i,j] = heights[i]
     415                    else:
     416                        mat[i,j] = mat[j,i] = ((span[i] + span[j]).height() - heights[i] - heights[j])/2
     417            reg = mat.det()
     418            if reg / min(heights) < 1e-6:
     419                raise ValueError, "points not linearly independent"
     420       
     421        if full_saturation:
     422            if reg_bound is None:
     423                # TODO: Hmm... verify this for rank > 1
     424                reg_bound = self.height_function().min(.1, 5) ** n
     425            bound = math.sqrt(reg/reg_bound)
     426            prime_list = prime_range(math.ceil(bound) + 1)
     427            if verbose:
     428                print "Testing primes up to", prime_list[-1]
     429        else:
     430            prime_list = prime_range(max_prime+1)
     431        if odd_primes_only and 2 in prime_list:
     432            prime_list.remove(2)
     433       
     434        index = ZZ(1)
     435       
     436        if len(span) == 1:
     437            # Rank one, trial division.
     438            P = span[0]
     439            for p in prime_list:
     440                if full_saturation and p > bound / index:
     441                    break
     442                division_points = P.division_points(p)
     443                while division_points:
     444                    if verbose:
     445                        print "    divisible by", p
     446                    P = division_points[0]
     447                    index *= p
     448                    division_points = P.division_points(p)
     449            return [P], index, reg / index**2
     450       
     451        else:
     452       
     453            def list_of_multiples(P, n):
     454                """
     455                Computes [i*P for i in range(n)] via repeated addition.
     456                """
     457                L = [P.curve()(0), P]
     458                for i in range(n-2):
     459                    L.append(L[-1] + P)
     460                return L
     461       
     462            # For small enough primes, attempt trial division on all
     463            # (normalized) linear combinations.
     464            small_primes = [p for p in prime_list if p==2 or p**n < 100]
     465            if small_primes and debug < 2:
     466                if verbose:
     467                    print "Trial divison for", small_primes
     468                max_small_prime = max(small_primes)
     469                while True:
     470                    found_divisor = False
     471                    lin_combs = {(0,) * n: self(0)}
     472                    for i, P in reversed(list(enumerate(span))):
     473                        multiples = list_of_multiples(P, max_small_prime)
     474                        for key, val in lin_combs.items():
     475                            is_first = sum(key) == 0
     476                            lin_comb = list(key)
     477                            for j, jP in enumerate(multiples):
     478                                if j == 0: continue
     479                                lin_comb[i] = j
     480                                lin_combs[tuple(lin_comb)] = val + jP
     481                                if is_first: break # The *last* non-zero entry is always 1.
     482                   
     483                    for p in small_primes:
     484                        for v in GF(p)**n:
     485                            if v == 0: continue
     486                            if v[v.nonzero_positions()[-1]] != 1: continue
     487                            divisors = lin_combs[tuple(v)].division_points(p)
     488                            if divisors:
     489                                if verbose:
     490                                    print "    divisible by", p
     491                                span[v.nonzero_positions()[-1]] = divisors[0]
     492                                found_divisor = True
     493                                index *= p
     494                                if full_saturation:
     495                                    prime_list = [p for p in prime_list if p <= bound/index]
     496                                    small_primes = [p for p in small_primes if p in prime_list]
     497                                break
     498                        if found_divisor:
     499                            # recompute lin_combs
     500                            break
     501
     502                    if not found_divisor:
     503                        # went through lin_combs for all p, exit main while loop
     504                        if verbose:
     505                            print "    saturated at", small_primes
     506                        for p in small_primes:
     507                            prime_list.remove(p)
     508                        break
     509           
     510            # Now compute maps E(K) \oplus F_p -> E(F_q) \oplus F_p.
     511            # Perhaps this could be done in parallel with the code above.
     512            K = self.base_ring()
     513            D = self.discriminant()
     514            all_projections = dict((p, ([], [], set())) for p in prime_list)
     515
     516            # First define two utility functions.
     517            def projection(Ek, p):
     518                """
     519                Project span onto distinct copies of E(k) \oplus \F_p.
     520                """
     521                G, gens = Ek.abelian_group() # cached!
     522                invs = G.invariants()
     523                c = invs[0]
     524                g = gens[0]
     525                if len(invs) == 1:
     526                    cc = c
     527                else:
     528                    cc = c // invs[1]
     529                if not p.divides(cc):
     530                    return None
     531                cofactor = c // p
     532                gg = cofactor * g
     533                logs = {}
     534                for j, jgg in enumerate(list_of_multiples(gg, p)):
     535                    logs[jgg] = j
     536                v = []
     537                for P in span:
     538                    try:
     539                        Pk = Ek(P)
     540                    except ZeroDivisionError:
     541                        Pk = Ek(0)
     542                    v.append(logs[cofactor * Pk])
     543                return vector(GF(p), v)
     544           
     545            def verify_independence(p):
     546                """
     547                Verifies that there are no mod p linear dependancies coming
     548                from the projections down to each E(k).
     549                """
     550                index = 1
     551                projections, proj_primes, faux_combs = all_projections[p]
     552                if len(projections) < n:
     553                    return 1
    262554               
     555                m = matrix(projections).transpose()
     556                if m.rank() == n:
     557                    prime_list.remove(p)
     558                    del all_projections[p]
     559                    if verbose:
     560                        if len(prime_list) < 10:
     561                            prime_list_str = str(prime_list)
     562                        else:
     563                            prime_list_str = str(prime_list[:5] + ["..."] + prime_list[-3:])
     564                        print "saturated at", p, "(remaining %s)" % prime_list_str
     565                    return 1
     566               
     567                for v in m.left_kernel().gens():
     568                    v.set_immutable()
     569                    if v in faux_combs:
     570                        continue
     571                    faux_combs.add(v)
     572                    Q = sum(int(a)*P for a, P in zip(v, span))
     573                    divisors = Q.division_points(p)
     574                    if divisors:
     575                        if verbose:
     576                            print "    divisible by", p
     577                        replaced = v.nonzero_positions()[-1]
     578                        v /= v[replaced] # normalize
     579                        Q = sum(int(a)*P for a, P in zip(v, span) if a)
     580                        span[replaced] = Q.division_points(p)[0]
     581                        # We have to fix all the other projections...
     582                        for other_p, (proj, proj_primes, faux_combs) in all_projections.items():
     583                            faux_combs.clear()
     584                            if other_p == p:
     585                                proj = []
     586                                for pp in proj_primes:
     587                                    proj.append(projection(self.reduction(pp), p))
     588                            else:
     589                                inv_p = 1/GF(other_p)(p)
     590                                for other_v in proj:
     591                                    other_v[replaced] = inv_p * sum(int(a)*b for a,b in zip(v, other_v))
     592                        return p * verify_independence(p)
     593                       
     594                return 1
     595
     596            # Now we enter the actual loop.
     597            proj_count = 0
     598            curve_count = 0
     599            for rat_p in primes(3, 1e9):
     600                for pp, e in K.factor(rat_p):
     601                    if not prime_list:
     602                        if verbose:
     603                            print "searched up to p =", rat_p, ", %s curves (%s used)" % (curve_count, proj_count)
     604                        return span, index, reg / index**2
     605                    if pp.divides(D):
     606                        continue
     607                    if verbose > 1:
     608                        print "rat_p", rat_p, "prime_list", prime_list, "count %s (%s)" % (curve_count, proj_count)
     609                    Ek = self.reduction(pp)
     610                    curve_count += 1
     611                    for p in prime_list:
     612                        v = projection(Ek, p)
     613                        if v is not None:
     614                            proj_count += 1
     615                            projections, proj_primes, faux_combs = all_projections[p]
     616                            projections.append(v)
     617                            proj_primes.append(pp)
     618                            old_index = index
     619                            index *= verify_independence(p)
     620                            if full_saturation and index != old_index:
     621                                for p in prime_list[:]:
     622                                    if p > bound / index:
     623                                        del all_projections[p]
     624                                        prime_list.remove(p)
     625   
     626    def gens(self, **kwds):
     627        """
     628        Compute and return generators for the Mordell-Weil group E(K)
     629        *modulo* torsion. Only implemented for K quadratic, E defined over Q.
     630        Any keyword arguments are passed on to E(Q).gens().
     631       
     632        EXAMPLES::
     633           
     634            sage: E = EllipticCurve('37a')
     635            sage: E.change_ring(QuadraticField(-1, 'i')).gens()
     636            [(0 : -1 : 1)]
     637            sage: E.change_ring(QuadraticField(3, 'i')).gens()
     638            [(0 : -1 : 1), (-i + 2 : 2*i - 4 : 1)]
     639           
     640            sage: K.<a> = QuadraticField(-7)
     641            sage: EllipticCurve('11a').change_ring(K).gens()
     642            [(-6 : -11/2*a - 1/2 : 1)]
     643            sage: EllipticCurve('389a').change_ring(K).gens()
     644            [(-1 : 1 : 1), (0 : -1 : 1), (1/8*a + 5/8 : -5/16*a - 9/16 : 1)]
     645
     646            sage: E = EllipticCurve([1, a])
     647            sage: E.gens()
     648            Traceback (most recent call last):
     649            ...
     650            NotImplementedError
     651        """
     652        if not kwds:
     653            try:
     654                return list(self.__gens)
     655            except AttributeError:
     656                pass
     657       
     658        K = self.base_ring()
     659        if K.absolute_degree() == 2:
     660            try:
     661                E = self.change_ring(QQ)
     662            except TypeError:
     663                pass
     664            else:
     665                ED = E.quadratic_twist(K.disc())
     666                E_gens = E.gens(**kwds)
     667                ED_gens = ED.gens(**kwds)
     668                phi = ED.change_ring(K).isomorphism_to(E)
     669                EK_span = [self(P) for P in E_gens] + [phi(P) for P in ED_gens]
     670                gens = self.saturation(EK_span, max_prime=2)[0]
     671                if not kwds:
     672                    self.__gens = gens
     673                return gens
     674        raise NotImplementedError
     675   
    263676    def is_local_integral_model(self,*P):
    264677        r"""
    265678        Tests if self is integral at the prime ideal `P`, or at all the
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r 48a42f04011f -r 44a69cc310af sage/schemes/elliptic_curves/ell_rational_field.py
    a b  
    22412241            sage: Q=5*P; Q
    22422242            (1/4 : -5/8 : 1)
    22432243            sage: E.saturation([Q])
    2244             ([(0 : 0 : 1)], '5', 0.0511114075779915)
     2244            ([(0 : 0 : 1)], 5, 0.0511114075779915)
    22452245        """
    22462246        if not isinstance(points, list):
    22472247            raise TypeError, "points (=%s) must be a list."%points
     
    22722272                break
    22732273        sat = mw.points()
    22742274        sat = [self(P) for P in sat]
    2275         return sat, index, R(reg)
     2275        return sat, ZZ(index), R(reg)
    22762276
    22772277
    22782278    def CPS_height_bound(self):