Ticket #14963: trac_14963-Sunits.patch

File trac_14963-Sunits.patch, 18.8 KB (added by cremona, 8 years ago)

applies to 5.12.beta1 + #14519 + #14746

  • sage/groups/abelian_gps/element_base.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1374874890 -3600
    # Node ID e2494b9eb61f4efdf3e0ed00bbaac02bd89ecd82
    # Parent  71b7bbac30a37f8ad2a80052f7d08ef5c57b50b9
    trac 14963: extend unit group class to handle S-units of number fields
    
    diff --git a/sage/groups/abelian_gps/element_base.py b/sage/groups/abelian_gps/element_base.py
    a b  
    7171        else:
    7272            self._exponents = tuple( ZZ(e) for e in exponents )
    7373            if len(self._exponents) != n:
    74                 raise IndexError('argument length (= %s) must be %s.'%(len(X), n))
     74                raise IndexError('argument length (= %s) must be %s.'%(len(exponents), n))
    7575
    7676    def exponents(self):
    7777        """
  • sage/libs/pari/gen.pyx

    diff --git a/sage/libs/pari/gen.pyx b/sage/libs/pari/gen.pyx
    a b  
    69086908        sig_on()
    69096909        return self.new_gen(bnfisunit(self.g, t0))
    69106910
     6911    def bnfissunit(self, sunit_data, x):
     6912        t0GEN(x)
     6913        t1GEN(sunit_data)
     6914        sig_on()
     6915        return self.new_gen(bnfissunit(self.g, t1, t0))
     6916
    69116917    def dirzetak(self, n):
    69126918        t0GEN(n)
    69136919        sig_on()
  • sage/rings/number_field/number_field.py

    diff --git a/sage/rings/number_field/number_field.py b/sage/rings/number_field/number_field.py
    a b  
    34533453       
    34543454        A list of generators of the unit group.
    34553455
     3456       .. note::
     3457
     3458            For more functionality see the S_unit_group() function.
     3459       
    34563460        EXAMPLE::
    34573461
    34583462            sage: K.<a> = QuadraticField(-3)
     
    53295333            sage: U.gens_values()  # result not independently verified
    53305334            [-1, a^9 + a - 1, a^16 - a^15 + a^14 - a^12 + a^11 - a^10 - a^8 + a^7 - 2*a^6 + a^4 - 3*a^3 + 2*a^2 - 2*a + 1, 2*a^16 - a^14 - a^13 + 3*a^12 - 2*a^10 + a^9 + 3*a^8 - 3*a^6 + 3*a^5 + 3*a^4 - 2*a^3 - 2*a^2 + 3*a + 4, a^15 + a^14 + 2*a^11 + a^10 - a^9 + a^8 + 2*a^7 - a^5 + 2*a^3 - a^2 - 3*a + 1, a^16 + a^15 + a^14 + a^13 + a^12 + a^11 + a^10 + a^9 + a^8 + a^7 + a^6 + a^5 + a^4 + a^3 + a^2 - 2, 2*a^16 - 3*a^15 + 3*a^14 - 3*a^13 + 3*a^12 - a^11 + a^9 - 3*a^8 + 4*a^7 - 5*a^6 + 6*a^5 - 4*a^4 + 3*a^3 - 2*a^2 - 2*a + 4, a^15 - a^12 + a^10 - a^9 - 2*a^8 + 3*a^7 + a^6 - 3*a^5 + a^4 + 4*a^3 - 3*a^2 - 2*a + 2, 2*a^16 + a^15 - a^11 - 3*a^10 - 4*a^9 - 4*a^8 - 4*a^7 - 5*a^6 - 7*a^5 - 8*a^4 - 6*a^3 - 5*a^2 - 6*a - 7]
    53315335        """
     5336        proof = proof_flag(proof)
     5337
    53325338        try:
    53335339            return self._unit_group
    53345340        except AttributeError:
     
    53475353            self._unit_group_no_proof = U
    53485354        return U
    53495355
     5356    def S_unit_group(self, proof=None, S=None):
     5357        """
     5358        Return the S-unit group (including torsion) of this number field.
     5359
     5360        ALGORITHM: Uses PARI's bnfsunit command.
     5361
     5362        INPUT:
     5363
     5364        - ``proof`` (bool, default True) flag passed to ``pari``.
     5365        - ``S`` - list or tuple of prime ideals, or an ideal, or a single
     5366          ideal or element from which an ideal can be constructed, in
     5367          which case the support is used.  If None, the global unit
     5368          group is constructed; otherwise, the S-unit group is
     5369          constructed.
     5370
     5371        .. note::
     5372
     5373           The group is cached.
     5374       
     5375        EXAMPLES::
     5376
     5377            sage: x = polygen(QQ)                                                 
     5378            sage: K.<a> = NumberField(x^4 - 10*x^3 + 20*5*x^2 - 15*5^2*x + 11*5^3)
     5379            sage: U = K.S_unit_group(S=a); U                                     
     5380            S-unit group with structure C10 x Z x Z x Z of Number Field in a with defining polynomial x^4 - 10*x^3 + 100*x^2 - 375*x + 1375 with S = (Fractional ideal (11, 1/275*a^3 + 4/55*a^2 - 5/11*a + 9), Fractional ideal (5, 1/275*a^3 + 4/55*a^2 - 5/11*a + 5))
     5381            sage: U.gens()
     5382            (u0, u1, u2, u3)
     5383            sage: U.gens_values()
     5384            [-7/275*a^3 + 1/11*a^2 - 9/11*a - 1, 6/275*a^3 - 9/55*a^2 + 14/11*a - 2, 14/275*a^3 - 21/55*a^2 + 29/11*a - 6, 1/275*a^3 + 4/55*a^2 - 5/11*a + 5]
     5385            sage: U.invariants()
     5386            (10, 0, 0, 0)
     5387            sage: [u.multiplicative_order() for u in U.gens()]
     5388            [10, +Infinity, +Infinity, +Infinity]
     5389            sage: U.primes()
     5390            (Fractional ideal (11, 1/275*a^3 + 4/55*a^2 - 5/11*a + 9), Fractional ideal (5, 1/275*a^3 + 4/55*a^2 - 5/11*a + 5))
     5391
     5392        With the default value of `S`, the S-unit group is the same as
     5393        the global unit group::
     5394
     5395            sage: x = polygen(QQ)
     5396            sage: K.<a> = NumberField(x^3 + 3)
     5397            sage: U = K.unit_group(proof=False)   
     5398            sage: U == K.S_unit_group(proof=False)
     5399            True
     5400
     5401        The value of `S` may be specified as a list of prime ideals,
     5402        or an ideal, or an element of the field::
     5403
     5404            sage: K.<a> = NumberField(x^3 + 3)
     5405            sage: U = K.S_unit_group(proof=False, S=K.ideal(6).prime_factors()); U
     5406            S-unit group with structure C2 x Z x Z x Z x Z of Number Field in a with defining polynomial x^3 + 3 with S = (Fractional ideal (-a^2 + a - 1), Fractional ideal (a + 1), Fractional ideal (a))
     5407            sage: K.<a> = NumberField(x^3 + 3)
     5408            sage: U = K.S_unit_group(proof=False, S=K.ideal(6)); U
     5409            S-unit group with structure C2 x Z x Z x Z x Z of Number Field in a with defining polynomial x^3 + 3 with S = (Fractional ideal (-a^2 + a - 1), Fractional ideal (a + 1), Fractional ideal (a))
     5410            sage: K.<a> = NumberField(x^3 + 3)
     5411            sage: U = K.S_unit_group(proof=False, S=6); U
     5412            S-unit group with structure C2 x Z x Z x Z x Z of Number Field in a with defining polynomial x^3 + 3 with S = (Fractional ideal (-a^2 + a - 1), Fractional ideal (a + 1), Fractional ideal (a))
     5413
     5414            sage: U
     5415            S-unit group with structure C2 x Z x Z x Z x Z of Number Field in a with defining polynomial x^3 + 3 with S = (Fractional ideal (-a^2 + a - 1), Fractional ideal (a + 1), Fractional ideal (a))
     5416            sage: U.primes()
     5417            (Fractional ideal (-a^2 + a - 1),
     5418            Fractional ideal (a + 1),
     5419            Fractional ideal (a))
     5420            sage: U.gens()
     5421            (u0, u1, u2, u3, u4)
     5422            sage: U.gens_values()
     5423            [-1, a^2 - 2, -a^2 + a - 1, a + 1, a]
     5424
     5425        The exp and log methods can be used to create `S`-units from
     5426        sequences of exponents, and recover the exponents::
     5427
     5428            sage: U.gens_orders()
     5429            (2, 0, 0, 0, 0)
     5430            sage: u = U.exp((3,1,4,1,5)); u
     5431            -6*a^2 + 18*a - 54
     5432            sage: u.norm().factor()       
     5433            -1 * 2^9 * 3^5
     5434            sage: U.log(u)                 
     5435            (1, 1, 4, 1, 5)
     5436           
     5437
     5438        """
     5439        proof = proof_flag(proof)
     5440
     5441        # process the parameter S:
     5442        if not S:
     5443            S = ()
     5444        else:
     5445            if type(S)==list:
     5446                S = tuple(S)
     5447            if not type(S)==tuple:
     5448                try:
     5449                    S = tuple(self.ideal(S).prime_factors())
     5450                except (NameError, TypeError, ValueError):
     5451                    raise ValueError("Cannot make a set of primes from %s"%(S,))
     5452            else:
     5453                try:
     5454                    S = tuple(self.ideal(P) for P in S)
     5455                except (NameError, TypeError, ValueError):
     5456                    raise ValueError("Cannot make a set of primes from %s"%(S,))
     5457                if not all([P.is_prime() for P in S]):
     5458                    raise ValueError("Not all elements of %s are prime ideals"%(S,))
     5459
     5460        try:
     5461            return self._S_unit_group_cache[S]
     5462        except AttributeError:
     5463            self._S_unit_group_cache = {}
     5464        except KeyError:
     5465            pass
     5466
     5467        if proof == False:
     5468            try:
     5469                return self._S_unit_group_no_proof_cache[S]
     5470            except AttributeError:
     5471                self._S_unit_group_no_proof_cache = {}
     5472            except KeyError:
     5473                pass
     5474
     5475        U = UnitGroup(self,proof,S=S)
     5476        if proof:
     5477            self._S_unit_group_cache[S] = U
     5478        else:
     5479            self._S_unit_group_no_proof_cache[S] = U
     5480        return U
     5481
    53505482    def zeta(self, n=2, all=False):
    53515483        """
    53525484        Return one, or a list of all, primitive n-th root of unity in this field.
  • sage/rings/number_field/unit_group.py

    diff --git a/sage/rings/number_field/unit_group.py b/sage/rings/number_field/unit_group.py
    a b  
    11r"""
    2 Unit Groups of Number Fields
     2Unit and S-unit groups of Number Fields
    33
    44EXAMPLES::
    55
     
    7070    sage: UK.fundamental_units()
    7171    [a^3 + a^2 - 1, a - 1]   
    7272
     73S-unit groups may be constructed, where S is a set of primes::
     74
     75    sage: K.<a> = NumberField(x^6+2)
     76    sage: S = K.ideal(3).prime_factors(); S
     77    [Fractional ideal (3, a + 1), Fractional ideal (3, a - 1)]
     78    sage: SUK = UnitGroup(K,S=tuple(S)); SUK
     79    S-unit group with structure C2 x Z x Z x Z x Z of Number Field in a with defining polynomial x^6 + 2 with S = (Fractional ideal (3, a + 1), Fractional ideal (3, a - 1))
     80    sage: SUK.primes()
     81    (Fractional ideal (3, a + 1), Fractional ideal (3, a - 1))
     82    sage: SUK.rank() 
     83    4
     84    sage: SUK.gens_values()
     85    [-1, a^2 + 1, a^5 + a^4 - a^2 - a - 1, a + 1, -a + 1]
     86    sage: u = 9*prod(SUK.gens_values()); u
     87    -18*a^5 - 18*a^4 - 18*a^3 - 9*a^2 + 9*a + 27
     88    sage: SUK.log(u)
     89    (1, 3, 1, 7, 7)
     90    sage: u == SUK.exp((1,3,1,7,7))
     91    True
     92
    7393A relative number field example::
    7494
    7595    sage: L.<a, b> = NumberField([x^2 + x + 1, x^4 + 1])
     
    126146
    127147class UnitGroup(AbelianGroupWithValues_class):
    128148    """
    129     The unit group of a number field.
     149    The unit group or an S-unit group of a number field.
    130150
    131151    TESTS::
    132152
     
    166186        z^9 + z^5 + z^4 + 1
    167187        sage: UK.gen(5) # random
    168188        z^5 + z
     189
     190    An S-unit group::
     191
     192        sage: SUK = UnitGroup(K,S=21); SUK
     193        S-unit group with structure C26 x Z x Z x Z x Z x Z x Z x Z x Z x Z x Z of Cyclotomic Field of order 13 and degree 12 with S = (Fractional ideal (3, z^3 - z - 1), Fractional ideal (3, z^3 + z^2 + z - 1), Fractional ideal (3, z^3 + z^2 - 1), Fractional ideal (3, z^3 - z^2 - z - 1), Fractional ideal (7))
     194        sage: SUK.rank()     
     195        10
     196        sage: SUK.zeta_order()
     197        26
     198        sage: SUK.log(21*z)
     199        (6, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
    169200    """
    170201    # This structure is not a parent in the usual sense. The
    171202    # "elements" are NumberFieldElement_absolute. Instead, they should
    172203    # derive from AbelianGroupElement and coerce into
    173204    # NumberFieldElement_absolute.
    174205
    175     def __init__(self, number_field, proof=True):
     206    def __init__(self, number_field, proof=True, S=None):
    176207        """
    177208        Create a unit group of a number field.
    178209
     
    180211
    181212        - ``number_field`` - a number field
    182213        - ``proof`` - boolean (default True): proof flag
     214        - ``S`` - tuple of prime ideals, or an ideal, or a single
     215          ideal or element from which an ideal can be constructed, in
     216          which case the support is used.  If None, the global unit
     217          group is constructed; otherwise, the S-unit group is
     218          constructed.
    183219
    184220        The proof flag is passed to pari via the ``pari_bnf()`` function
    185221        which computes the unit group.  See the documentation for the
     
    211247            (u0, u1, u2, u3, u4, u5)
    212248            sage: UK.gens_values() # random
    213249            [-z^11, z^5 + z^3, z^6 + z^5, z^9 + z^7 + z^5, z^9 + z^5 + z^4 + 1, z^5 + z]
     250            sage: SUK = UnitGroup(K,S=2); SUK
     251            S-unit group with structure C26 x Z x Z x Z x Z x Z x Z of Cyclotomic Field of order 13 and degree 12 with S = (Fractional ideal (2),)
     252
    214253            """
    215254        proof = get_flag(proof, "number_field")
    216255        K = number_field
    217256        pK = K.pari_bnf(proof)
    218257        self.__number_field = K
     258        self.__pari_number_field = pK
    219259
    220         # compute the units via pari:
     260        # process the parameter S:
     261        if not S:
     262            S = self.__S = ()
     263        else:
     264            if type(S)==list:
     265                S = tuple(S)
     266            if not type(S)==tuple:
     267                try:
     268                    S = tuple(K.ideal(S).prime_factors())
     269                except (NameError, TypeError, ValueError):
     270                    raise ValueError("Cannot make a set of primes from %s"%(S,))
     271            else:
     272                try:
     273                    S = tuple(K.ideal(P) for P in S)
     274                except (NameError, TypeError, ValueError):
     275                    raise ValueError("Cannot make a set of primes from %s"%(S,))
     276                if not all([P.is_prime() for P in S]):
     277                    raise ValueError("Not all elements of %s are prime ideals"%(S,))
     278            self.__S = S
     279            self.__pS = pS = [P.pari_prime() for P in S]
     280
     281        # compute the fundemental units via pari:
    221282        fu = [K(u) for u in pK.bnfunit()]
     283        self.__nfu = len(fu)
     284
     285        # compute the additional S-unit generators:
     286        if S:
     287            self.__S_unit_data = pK.bnfsunit(pS)
     288            su = [K(u) for u in self.__S_unit_data[0]]
     289        else:
     290            su = []
     291        self.__nsu = len(su)
     292        self.__rank = self.__nfu + self.__nsu
    222293
    223294        # compute a torsion generator and pick the 'simplest' one:
    224295        n, z = pK.nfrootsof1()
     
    230301        # to allow for this in the dlog function!  So we do not.
    231302
    232303        # Store the actual generators (torsion first):
    233         gens = [z] + fu
     304        gens = [z] + fu + su
    234305        values = Sequence(gens, immutable=True, universe=self, check=False)
    235306        # Construct the abtract group:       
    236         gens_orders = tuple([ZZ(n)]+[ZZ(0)]*len(fu))
     307        gens_orders = tuple([ZZ(n)]+[ZZ(0)]*(self.__rank))
    237308        AbelianGroupWithValues_class.__init__(self, gens_orders, 'u', values, number_field)
    238309
    239310    def _element_constructor_(self, u):
     
    272343            ValueError: a is not a unit
    273344        """
    274345        K = self.__number_field
     346        pK = self.__pari_number_field
    275347        try:
    276348            u = K(u)
    277349        except TypeError:
    278350            raise ValueError, "%s is not an element of %s"%(u,K)
    279         if not u.is_integral() or u.norm().abs() != 1:
    280             raise ValueError, "%s is not a unit"%u
    281         m = K.pari_bnf().bnfisunit(pari(u)).mattranspose()
     351        if self.__S:
     352            m = pK.bnfissunit(self.__S_unit_data, pari(u)).mattranspose()
     353            if m.ncols()==0:
     354                raise ValueError, "%s is not an S-unit"%u               
     355        else:
     356            if not u.is_integral() or u.norm().abs() != 1:
     357                raise ValueError, "%s is not a unit"%u
     358            m = pK.bnfisunit(pari(u)).mattranspose()
     359
    282360        # convert column matrix to a list:
    283361        m = [ZZ(m[0,i].python()) for i in range(m.ncols())]
    284         # NB pari puts the torsion at the end!
    285         m.insert(0,m.pop())
     362
     363        # NB pari puts the torsion after the fundamental units, before
     364        # the extra S-units but we have the torsion first:
     365        m = [m[self.__nfu]] + m[:self.__nfu] + m[self.__nfu+1:]
     366
    286367        return self.element_class(self, m)
    287368
    288369    def rank(self):
     
    294375            sage: K.<z> = CyclotomicField(13)
    295376            sage: UnitGroup(K).rank()
    296377            5
     378            sage: SUK = UnitGroup(K,S=2); SUK.rank()
     379            6
    297380        """
    298381        return self.ngens()-1
    299382
     
    309392            Unit group with structure C2 x Z of Number Field in a with defining polynomial x^3 - 2
    310393            sage: U._repr_()
    311394            'Unit group with structure C2 x Z of Number Field in a with defining polynomial x^3 - 2'
     395            sage: UnitGroup(NumberField(x^3 - 2, 'a'),S=2)   
     396            S-unit group with structure C2 x Z x Z of Number Field in a with defining polynomial x^3 - 2 with S = (Fractional ideal (a),)
    312397        """
     398        if self.__S:
     399            return 'S-unit group with structure %s of %s with S = %s'%(
     400                self._group_notation(self.gens_orders()),
     401                self.number_field(),
     402                self.primes())
    313403        return 'Unit group with structure %s of %s'%(
    314404            self._group_notation(self.gens_orders()),
    315405            self.number_field())
     
    448538       
    449539        EXAMPLES::
    450540
    451             sage: x = polygen(QQ)
    452             sage: U = UnitGroup(NumberField(x^2 + 23, 'w')); U
     541            sage: U = UnitGroup(QuadraticField(-23, 'w')); U
    453542            Unit group with structure C2 of Number Field in w with defining polynomial x^2 + 23
    454543            sage: U.number_field()
    455544            Number Field in w with defining polynomial x^2 + 23
     
    457546        return self.__number_field
    458547
    459548
     549    def primes(self):
     550        """
     551        Return the (possibly empty) list of primes associated with this S-unit group.
     552       
     553        EXAMPLES::
     554
     555            sage: K.<a> = QuadraticField(-23)     
     556            sage: S = tuple(K.ideal(3).prime_factors()); S
     557            (Fractional ideal (3, 1/2*a - 1/2), Fractional ideal (3, 1/2*a + 1/2))
     558            sage: U = UnitGroup(K,S=tuple(S)); U
     559            S-unit group with structure C2 x Z x Z of Number Field in a with defining polynomial x^2 + 23 with S = (Fractional ideal (3, 1/2*a - 1/2), Fractional ideal (3, 1/2*a + 1/2))
     560            sage: U.primes() == S                     
     561            True
     562        """
     563        return self.__S
     564
     565
    460566    def log(self, u):
    461567        r"""
    462568        Return the exponents of the unit ``u`` with respect to group generators.
     
    488594            -253576*z^11 + 7003*z^10 - 395532*z^9 - 35275*z^8 - 500326*z^7 - 35275*z^6 - 395532*z^5 + 7003*z^4 - 253576*z^3 - 59925*z - 59925
    489595            sage: UK.log(unit)
    490596            (13, 6, 7, 8, 9, 10)
     597
     598        An S-unit example::
     599       
     600           sage: SUK = UnitGroup(K,S=2)           
     601           sage: v = (3,1,4,1,5,9,2)
     602           sage: u = SUK.exp(v); u
     603           -997204*z^11 - 2419728*z^10 - 413812*z^9 - 413812*z^8 - 2419728*z^7 - 997204*z^6 - 2129888*z^4 - 1616524*z^3 + 149364*z^2 - 1616524*z - 2129888
     604           sage: SUK.log(u)         
     605           (3, 1, 4, 1, 5, 9, 2)
     606           sage: SUK.log(u) == v   
     607           True
    491608        """
    492609        return self(u).exponents()
    493610
     
    523640            (13, 6, 7, 8, 9, 10)
    524641            sage: UK.exp(UK.log(u)) == u.value()
    525642            True
     643
     644        An S-unit example::
     645       
     646           sage: SUK = UnitGroup(K,S=2)           
     647           sage: v = (3,1,4,1,5,9,2)
     648           sage: u = SUK.exp(v); u
     649           -997204*z^11 - 2419728*z^10 - 413812*z^9 - 413812*z^8 - 2419728*z^7 - 997204*z^6 - 2129888*z^4 - 1616524*z^3 + 149364*z^2 - 1616524*z - 2129888
     650           sage: SUK.log(u)         
     651           (3, 1, 4, 1, 5, 9, 2)
     652           sage: SUK.log(u) == v   
     653           True
    526654        """
    527655        return prod([u**e for u,e in zip(self.gens_values(),exponents)], self.number_field().one_element())
    528656