Ticket #7931: 7931_nth_root-folded.patch

File 7931_nth_root-folded.patch, 35.0 KB (added by davidloeffler, 9 years ago)

Qfolded patch. Apply only this patch.

  • sage/rings/finite_rings/all.py

    # HG changeset patch
    # User David Roe <roed@math.harvard.edu>
    # Date 1266938910 28800
    # Node ID b6a452dd35f8b70faddcf9c42bea80f56a2ccec1
    # Parent  4a9f7b41ae22e247c2a6675c7854e419f7340d09
    #7931: implements a better algorithm for nth roots in finite fields and Z/nZ.
    
    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/all.py
    a b  
    2222                          conway_polynomial, exists_conway_polynomial)
    2323GF = FiniteField
    2424
    25 from element_base import FiniteFieldElement, is_FiniteFieldElement
     25from element_base import FinitePolyExtElement as FiniteFieldElement # for backward compatibility; is this needed?
     26from element_base import is_FiniteFieldElement
    2627
  • sage/rings/finite_rings/element_base.pxd

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_base.pxd
    a b  
    1 from sage.structure.element cimport FieldElement
     1from sage.structure.element cimport CommutativeRingElement
    22
    3 cdef class FiniteFieldElement(FieldElement):
     3cdef class FiniteRingElement(CommutativeRingElement):
    44    pass
     5
     6cdef class FinitePolyExtElement(FiniteRingElement):
     7    pass
  • sage/rings/finite_rings/element_base.pyx

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_base.pyx
    a b  
    2929    from sage.rings.finite_rings.finite_field_base import is_FiniteField
    3030    return isinstance(x, Element) and is_FiniteField(x.parent())
    3131
    32 cdef class FiniteFieldElement(FieldElement):
     32cdef class FiniteRingElement(CommutativeRingElement):
     33    def _nth_root_common(self, n, all, algorithm, cunningham):
     34        """
     35        This function exists to reduce code duplication between finite field nth roots and integer_mod nth roots.
     36       
     37        The inputs are described there.
     38       
     39        TESTS::
     40       
     41            sage: a = Zmod(17)(13)
     42            sage: a._nth_root_common(4, True, "Johnston", False)
     43            [3, 5, 14, 12]
     44        """
     45        K = self.parent()
     46        q = K.order()
     47        if self.is_one():
     48            gcd = n.gcd(q-1)
     49            if gcd == 1:
     50                if all: return [self]
     51                else: return self
     52            else:
     53                # the following may eventually be improved to not need a multiplicative generator.
     54                g = K.multiplicative_generator()
     55                q1overn = (q-1)//gcd
     56                nthroot = g**q1overn
     57                return [nthroot**a for a in range(gcd)] if all else nthroot
     58        n = n % (q-1)
     59        if n == 0:
     60            if all: return []
     61            else: raise ValueError, "no nth root"
     62        gcd, alpha, beta = n.xgcd(q-1) # gcd = alpha*n + beta*(q-1), so 1/n = alpha/gcd (mod q-1)
     63        if gcd == 1:
     64            return [self**alpha] if all else self**alpha
     65        n = gcd
     66        q1overn = (q-1)//n
     67        if self**q1overn != 1:
     68            if all: return []
     69            else: raise ValueError, "no nth root"
     70        self = self**alpha
     71        if cunningham:
     72            F = n._factor_cunningham()
     73        else:
     74            F = n.factor()
     75        from sage.groups.generic import discrete_log
     76        if algorithm is None or algorithm == 'Johnston':
     77            g = K.multiplicative_generator()
     78            for r, v in F:
     79                k, h = (q-1).val_unit(r)
     80                z = h * (-h).inverse_mod(r**v)
     81                x = (1 + z) // r**v
     82                if k == 1:
     83                    self = self**x
     84                else:
     85                    t = discrete_log(self**h, g**(r**v*h), r**(k-v), operation='*')
     86                    self = self**x * g**(-z*t)
     87            if all:
     88                nthroot = g**q1overn
     89                L = [self]
     90                for i in range(1,n):
     91                    self *= nthroot
     92                    L.append(self)
     93                return L
     94            else:
     95                return self
     96        else:
     97            raise ValueError, "unknown algorithm"
     98
     99cdef class FinitePolyExtElement(FiniteRingElement):
    33100   
    34101    def _im_gens_(self, codomain, im_gens):
    35102        """
     
    386453        if self.is_zero():
    387454            raise ArithmeticError, "Multiplicative order of 0 not defined."
    388455        n = self._parent.order() - 1
     456        F = self._parent.factored_unit_order()[0]
    389457        order = 1
    390         for p, e in sage.rings.arith.factor(n):
     458        for p, e in F:
    391459            # Determine the power of p that divides the order.
    392             a = self**(n/(p**e))
     460            a = self**(n//(p**e))
    393461            while a != 1:
    394462                order = order * p
    395463                a = a**p
    396464
    397465        return order
    398466
    399     def nth_root(self, int n, extend = False, all = False):
    400         r"""
    401         Returns an nth root of self.
    402 
    403         INPUT:
    404 
    405         - ``n`` - integer >= 1 (must fit in C int type)
    406 
    407         - ``extend`` - bool (default: True); if True, return an nth
    408           root in an extension ring, if necessary. Otherwise, raise a
    409           ValueError if the root is not in the base ring.  Warning:
    410           this option is not implemented!
    411 
    412         - ``all`` - bool (default: False); if True, return all nth
    413           roots of self, instead of just one.
    414 
    415         OUTPUT:
    416 
    417         If self has an nth root, returns one (if all == False) or a
    418         list of all of them (if all == True).  Otherwise, raises a
    419         ValueError (if extend = False) or a NotImplementedError (if
    420         extend = True).
    421 
    422         .. warning::
    423        
    424            The 'extend' option is not implemented (yet).
    425 
    426         AUTHOR:
    427 
    428         - David Roe (2007-10-3)
    429 
    430         """
    431         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    432         from sage.rings.integer import Integer
    433         if extend:
    434             raise NotImplementedError
    435         R = PolynomialRing(self.parent(), "x")
    436         f = R([-self] + [self.parent()(Integer(0))] * (n - 1) + [self.parent()(1)])
    437         L = f.roots()
    438         if all:
    439             return [x[0] for x in L]
    440         else:
    441             if len(L) == 0:
    442                 raise ValueError, "no nth root"
    443             else:
    444                 return L[0][0]
    445 
    446467    def additive_order(self):
    447468        """
    448469        Return the additive order of this finite field element.
     
    461482            return Integer(1)
    462483        return self.parent().characteristic()
    463484
     485    def nth_root(self, n, extend = False, all = False, algorithm=None, cunningham=False):
     486        r"""
     487        Returns an `n`\th root of ``self``.
     488
     489        INPUT:
     490
     491        - ``n`` - integer `\geq 1`
     492
     493        - ``extend`` - bool (default: ``False``); if ``True``, return an `n`\th
     494          root in an extension ring, if necessary. Otherwise, raise a
     495          ValueError if the root is not in the base ring.  Warning:
     496          this option is not implemented!
     497
     498        - ``all`` - bool (default: ``False``); if ``True``, return all `n`\th
     499          roots of ``self``, instead of just one.
     500
     501        - ``algorithm`` - string (default: ``None``); 'Johnston' is the only
     502          currently supported option.  For IntegerMod elements, the problem
     503          is reduced to the prime modulus case using CRT and `p`-adic logs,
     504          and then this algorithm used.
     505
     506        OUTPUT:
     507
     508        If self has an `n`\th root, returns one (if ``all`` is ``False``) or a
     509        list of all of them (if ``all`` is ``True``).  Otherwise, raises a
     510        ValueError (if ``extend`` is ``False``) or a ``NotImplementedError`` (if
     511        ``extend`` is ``True``).
     512
     513        .. warning::
     514       
     515           The 'extend' option is not implemented (yet).
     516
     517        EXAMPLES::
     518
     519            sage: K = GF(31)
     520            sage: a = K(22)
     521            sage: K(22).nth_root(7)
     522            13
     523            sage: K(25).nth_root(5)
     524            5
     525            sage: K(23).nth_root(3)
     526            29
     527
     528            sage: K.<a> = GF(625)
     529            sage: (3*a^2+a+1).nth_root(13)**13
     530            3*a^2 + a + 1
     531
     532            sage: k.<a> = GF(29^2)
     533            sage: b = a^2 + 5*a + 1
     534            sage: b.nth_root(11)
     535            3*a + 20
     536            sage: b.nth_root(5)
     537            Traceback (most recent call last):
     538            ...
     539            ValueError: no nth root
     540            sage: b.nth_root(5, all = True)
     541            []
     542            sage: b.nth_root(3, all = True)
     543            [14*a + 18, 10*a + 13, 5*a + 27]
     544
     545            sage: k.<a> = GF(29^5)
     546            sage: b = a^2 + 5*a + 1
     547            sage: b.nth_root(5)
     548            19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
     549            sage: b.nth_root(7)
     550            Traceback (most recent call last):
     551            ...
     552            ValueError: no nth root
     553            sage: b.nth_root(4, all=True)
     554            []
     555
     556        TESTS::
     557
     558            sage: for p in [2,3,5,7,11]: # long
     559            ...       for n in [2,5,10]:
     560            ...           q = p^n
     561            ...           K.<a> = GF(q)
     562            ...           for r in (q-1).divisors():
     563            ...               if r == 1: continue
     564            ...               x = K.random_element()
     565            ...               y = x^r
     566            ...               if y.nth_root(r)**r != y: raise RuntimeError
     567            ...               if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     568            ...               if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     569
     570            sage: k.<a> = GF(4)
     571            sage: a.nth_root(0,all=True)
     572            []
     573            sage: k(1).nth_root(0,all=True)
     574            [a, a + 1, 1]
     575
     576        ALGORITHMS:
     577
     578        - The default is currently an algorithm described in the following paper:
     579
     580        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     581
     582        AUTHOR:
     583
     584        - David Roe (2010-02-13)
     585        """
     586        if self.is_zero():
     587            if n <= 0:
     588                if all: return []
     589                else: raise ValueError
     590            if all: return [self]
     591            else: return self
     592        if n < 0:
     593            self = ~self
     594            n = -n
     595        elif n == 0:
     596            if self == 1:
     597                if all: return [a for a in self.parent().list() if a != 0]
     598                else: return self
     599            else:
     600                if all: return []
     601                else: raise ValueError
     602        if extend:
     603            raise NotImplementedError
     604        from sage.rings.integer import Integer
     605        n = Integer(n)
     606        return self._nth_root_common(n, all, algorithm, cunningham)
     607
    464608    def pth_power(self, int k = 1):
    465609        """
    466610        Return the `(p^k)^{th}` power of self, where `p` is the
     
    516660            b^11 + b^10 + b^9 + b^7 + b^5 + b^4 + b^2 + b
    517661        """
    518662        return self.pth_power(-k)
     663
  • sage/rings/finite_rings/element_ext_pari.py

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_ext_pari.py
    a b  
    3434from sage.rings.integer import Integer
    3535import sage.rings.rational as rational
    3636from sage.libs.pari.all import pari, pari_gen
    37 from sage.rings.finite_rings.element_base import FiniteFieldElement
     37from sage.rings.finite_rings.element_base import FinitePolyExtElement
    3838import sage.rings.field_element as field_element
    3939import sage.rings.finite_rings.integer_mod as integer_mod
    4040from element_base import is_FiniteFieldElement
     
    4242from sage.structure.dynamic_class import dynamic_class
    4343from sage.categories.finite_fields import FiniteFields
    4444
    45 class FiniteField_ext_pariElement(FiniteFieldElement):
     45class FiniteField_ext_pariElement(FinitePolyExtElement):
    4646    """
    4747    An element of a finite field.
    4848   
     
    305305        """
    306306        return self.square_root(extend=extend, all=all)
    307307
    308     def nth_root(self, n, extend = False, all = False):
    309         r"""
    310         Returns an nth root of self.
    311        
    312         INPUT:
    313        
    314        
    315         -  ``n`` - integer = 1 (must fit in C int type)
    316        
    317         -  ``extend`` - bool (default: True); if True, return
    318            an nth root in an extension ring, if necessary. Otherwise, raise a
    319            ValueError if the root is not in the base ring. Warning: this
    320            option is not implemented!
    321        
    322         -  ``all`` - bool (default: False); if True, return all
    323            nth roots of self, instead of just one.
    324        
    325        
    326         OUTPUT: If self has an nth root, returns one (if all == False) or a
    327         list of all of them (if all == True). Otherwise, raises a
    328         ValueError (if extend = False) or a NotImplementedError (if extend
    329         = True).
    330        
    331         .. warning::
    332 
    333            The 'extend' option is not implemented (yet).
    334        
    335         AUTHORS:
    336 
    337         - David Roe (2007-10-3)
    338        
    339         EXAMPLES::
    340        
    341             sage: k.<a> = GF(29^5)
    342             sage: b = a^2 + 5*a + 1
    343             sage: b.nth_root(5)
    344             19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
    345             sage: b.nth_root(7)
    346             Traceback (most recent call last):
    347             ...
    348             ValueError: no nth root
    349             sage: b.nth_root(4, all=True)
    350             []
    351         """
    352         if extend:
    353             raise NotImplementedError
    354         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    355         R = PolynomialRing(self.parent(), "x")
    356         f = R([-self] + [self.parent()(0)] * (n - 1) + [self.parent()(1)])
    357         L = f.roots()
    358         if all:
    359             return [x[0] for x in L]
    360         else:
    361             if len(L) == 0:
    362                 raise ValueError, "no nth root"
    363             else:
    364                 return L[0][0]
    365        
    366308    def rational_reconstruction(self):
    367309        """
    368310        If the parent field is a prime field, uses rational reconstruction
  • sage/rings/finite_rings/element_givaro.pxd

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_givaro.pxd
    a b  
    11from sage.structure.element cimport Element, RingElement, ModuleElement
    2 from sage.rings.finite_rings.element_base cimport FiniteFieldElement
     2from sage.rings.finite_rings.element_base cimport FinitePolyExtElement
    33
    44from sage.structure.parent  cimport Parent
    55from sage.structure.sage_object cimport SageObject
     
    6161    void delete "delete "(void *o)
    6262    int gfq_element_factory "GFqDom<int>::Element"()
    6363
    64 cdef class FiniteField_givaroElement(FiniteFieldElement) #forward declaration
     64cdef class FiniteField_givaroElement(FinitePolyExtElement) #forward declaration
    6565
    6666cdef class Cache_givaro(SageObject):
    6767    cdef GivaroGfq *objectptr # C++ object
     
    8383    cdef int iterator
    8484    cdef Cache_givaro _cache
    8585
    86 cdef class FiniteField_givaroElement(FiniteFieldElement):
     86cdef class FiniteField_givaroElement(FinitePolyExtElement):
    8787    cdef int element
    8888    cdef Cache_givaro _cache
    8989    cdef object _multiplicative_order
  • sage/rings/finite_rings/element_givaro.pyx

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_givaro.pyx
    a b  
    452452            pass # handle this in next if clause
    453453
    454454        elif PY_TYPE_CHECK(e,FiniteField_ext_pariElement):
    455             # reduce FiniteFieldElements to pari
     455            # reduce FiniteField_ext_pariElements to pari
    456456            e = e._pari_()
    457457
    458458        elif sage.interfaces.gap.is_GapElement(e):
     
    838838        """
    839839        return self
    840840
    841 cdef class FiniteField_givaroElement(FiniteFieldElement):
     841cdef class FiniteField_givaroElement(FinitePolyExtElement):
    842842    """
    843843    An element of a (Givaro) finite field.
    844844    """
     
    867867            0
    868868           
    869869        """
    870         FiniteFieldElement.__init__(self, parent)
     870        FinitePolyExtElement.__init__(self, parent)
    871871        self._cache = parent._cache
    872872        self.element = 0
    873873       
  • sage/rings/finite_rings/element_ntl_gf2e.pxd

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_ntl_gf2e.pxd
    a b  
    22
    33from sage.rings.finite_rings.finite_field_base cimport FiniteField
    44from sage.structure.element cimport Element, RingElement, ModuleElement
    5 from sage.rings.finite_rings.element_base cimport FiniteFieldElement
     5from sage.rings.finite_rings.element_base cimport FinitePolyExtElement
    66
    7 cdef class FiniteField_ntl_gf2eElement(FiniteFieldElement)
     7cdef class FiniteField_ntl_gf2eElement(FinitePolyExtElement)
    88
    99cdef class FiniteField_ntl_gf2e(FiniteField):
    1010    cdef GF2EContext_c F
     
    1919    cdef int iterator
    2020    cdef FiniteField_ntl_gf2e _parent
    2121
    22 cdef class FiniteField_ntl_gf2eElement(FiniteFieldElement):
     22cdef class FiniteField_ntl_gf2eElement(FinitePolyExtElement):
    2323    cdef GF2E_c x
    2424    cdef FiniteField_ntl_gf2eElement _new(FiniteField_ntl_gf2eElement self)
  • sage/rings/finite_rings/element_ntl_gf2e.pyx

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/element_ntl_gf2e.pyx
    a b  
    482482            pass # handle this in next if clause
    483483
    484484        elif PY_TYPE_CHECK(e,FiniteField_ext_pariElement):
    485             # reduce FiniteFieldElements to pari
     485            # reduce FiniteField_ext_pariElements to pari
    486486            e = e._pari_()
    487487
    488488        elif is_GapElement(e):
     
    728728        f = pari(str(self.modulus()))
    729729        return f.subst('x', 'a') * pari("Mod(1,%s)"%self.characteristic())
    730730       
    731 cdef class FiniteField_ntl_gf2eElement(FiniteFieldElement):
     731cdef class FiniteField_ntl_gf2eElement(FinitePolyExtElement):
    732732    """
    733733    An element of an NTL:GF2E finite field.
    734734    """
  • sage/rings/finite_rings/finite_field_base.pxd

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/finite_field_base.pxd
    a b  
    11from sage.rings.ring cimport Field
    22
    33cdef class FiniteField(Field):
     4    cdef public object __factored_unit_order
    45    cdef public object __multiplicative_generator
    56    cdef public object __polynomial_ring
    67    cdef public object __vector_space
  • sage/rings/finite_rings/finite_field_base.pyx

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/finite_field_base.pyx
    a b  
    480480            997
    481481        """
    482482        raise NotImplementedError
     483
     484    def factored_order(self):
     485        """
     486        Returns the factored order of this field.  For compatibility with IntegerModRing.
     487
     488        EXAMPLES::
     489
     490            sage: GF(7^2,'a').factored_order()
     491            7^2
     492        """
     493        from sage.structure.factorization import Factorization
     494        return Factorization([(self.characteristic(), self.degree())])
     495
     496    def factored_unit_order(self):
     497        """
     498        Returns the factorization of ``self.order()-1``, as a 1-element list.  The format is for compatibility with IntegerModRing.
     499
     500        EXAMPLES::
     501
     502            sage: GF(7^2,'a').factored_unit_order()
     503            [2^4 * 3]
     504        """
     505        if self.__factored_unit_order is None:
     506            if self.characteristic() in []: # want to be [2,3,5,7,11] once #7240 is finished.
     507                self.__factored_unit_order = [(self.order()-1)._factor_cunningham()]
     508            else:
     509                self.__factored_unit_order = [(self.order()-1).factor()]
     510        return self.__factored_unit_order
    483511   
    484512    def cardinality(self):
    485513        """
  • sage/rings/finite_rings/integer_mod.pxd

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/integer_mod.pxd
    a b  
    88    int_fast32_t INTEGER_MOD_INT32_LIMIT
    99    int_fast64_t INTEGER_MOD_INT64_LIMIT
    1010
    11 cimport sage.structure.element
     11from sage.rings.finite_rings.element_base cimport FiniteRingElement
    1212from sage.rings.integer cimport Integer
    1313
    1414cdef class NativeIntStruct:
     
    1919    cdef object inverses # also a list
    2020    cdef lookup(NativeIntStruct self, Py_ssize_t value)
    2121
    22 cdef class IntegerMod_abstract(sage.structure.element.CommutativeRingElement):
     22cdef class IntegerMod_abstract(FiniteRingElement):
    2323    cdef NativeIntStruct __modulus
    2424    cdef _new_c_from_long(self, long value)
    2525    cdef void set_from_mpz(self, mpz_t value)
  • sage/rings/finite_rings/integer_mod.pyx

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/integer_mod.pyx
    a b  
    255255        return <object>PyList_GET_ITEM(self.table, value)
    256256
    257257   
    258 cdef class IntegerMod_abstract(sage.structure.element.CommutativeRingElement):
     258cdef class IntegerMod_abstract(FiniteRingElement):
    259259
    260260    def __init__(self, parent):
    261261        """
     
    972972               
    973973    square_root = sqrt
    974974
    975     def nth_root(self, int n, extend = False, all = False):
     975    def nth_root(self, n, extend = False, all = False, algorithm = None, cunningham = False):
    976976        r"""
    977977        Returns an `n`\th root of ``self``.
    978978       
    979979        INPUT:
    980        
    981        
    982         -  ``n`` - integer `\geq 1` (must fit in C
    983            ``int`` type)
    984        
    985         -  ``all`` - bool (default: ``False``); if
    986            ``True``, return all `n`\th roots of
    987            ``self``, instead of just one.
    988        
    989        
    990         OUTPUT: If self has an `n`\th root, returns one (if
    991         ``all`` is false) or a list of all of them (if
    992         ``all`` is true). Otherwise, raises a
    993         ``ValueError``.
    994        
    995         AUTHORS:
    996 
    997         - David Roe (2007-10-3)
     980
     981        - ``n`` - integer `\geq 1`
     982
     983        - ``extend`` - bool (default: True); if True, return an nth
     984          root in an extension ring, if necessary. Otherwise, raise a
     985          ValueError if the root is not in the base ring.  Warning:
     986          this option is not implemented!
     987
     988        - ``all`` - bool (default: ``False``); if ``True``, return all `n`\th
     989          roots of ``self``, instead of just one.
     990
     991        - ``algorithm`` - string (default: None); The algorithm for the prime modulus case.
     992          CRT and p-adic log techniques are used to reduce to this case.
     993          'Johnston' is the only currently supported option.
     994
     995        OUTPUT:
     996
     997        If self has an `n`\th root, returns one (if ``all`` is ``False``) or a
     998        list of all of them (if ``all`` is ``True``).  Otherwise, raises a
     999        ``ValueError`` (if ``extend`` is ``False``) or a ``NotImplementedError`` (if
     1000        ``extend`` is ``True``).
     1001
     1002        .. warning::
     1003           The 'extend' option is not implemented (yet).
     1004
     1005        NOTES:
     1006
     1007        - If `n = 0`
     1008           - ``all=True``
     1009             - ``self=1``; all nonzero elements of the parent are returned in a list.
     1010               Note that this could be very expensive for large parents.
     1011             - otherwise; an empty list is returned
     1012           - ``all=False``
     1013             - ``self=1``; ``self`` is returned
     1014             - otherwise; a ``ValueError`` is raised
     1015        - If `n < 0`
     1016          - if self is invertible, the `(-n)`\th root of the inverse of self is returned
     1017          - otherwise a ValueError is raised or empty list returned.
    9981018       
    9991019        EXAMPLES::
    10001020       
    1001             sage: k.<a> = GF(29)
    1002             sage: b = a^2 + 5*a + 1
    1003             sage: b.nth_root(5)
    1004             24
    1005             sage: b.nth_root(7)
     1021
     1022            sage: K = GF(31)
     1023            sage: a = K(22)
     1024            sage: K(22).nth_root(7)
     1025            13
     1026            sage: K(25).nth_root(5)
     1027            5
     1028            sage: K(23).nth_root(3)
     1029            29
     1030            sage: mod(225,2^5*3^2).nth_root(4, all=True)
     1031            [225, 129, 33, 63, 255, 159, 9, 201, 105, 279, 183, 87, 81, 273, 177, 207, 111, 15, 153, 57, 249, 135, 39, 231]
     1032            sage: mod(275,2^5*7^4).nth_root(7, all=True)
     1033            [58235, 25307, 69211, 36283, 3355, 47259, 14331]
     1034            sage: mod(1,8).nth_root(2,all=True)
     1035            [1, 7, 5, 3]
     1036            sage: mod(4,8).nth_root(2,all=True)
     1037            [2, 6]
     1038            sage: mod(1,16).nth_root(4,all=True)
     1039            [1, 15, 13, 3, 9, 7, 5, 11]
     1040            sage: (mod(22,31)^200).nth_root(200)
     1041            5
     1042            sage: mod(3,6).nth_root(0,all=True)
     1043            []
     1044            sage: mod(3,6).nth_root(0)
    10061045            Traceback (most recent call last):
    10071046            ...
    1008             ValueError: no nth root
    1009             sage: b.nth_root(4, all=True)
    1010             [21, 20, 9, 8]
     1047            ValueError
     1048            sage: mod(1,6).nth_root(0,all=True)
     1049            [1, 2, 3, 4, 5]
     1050
     1051        TESTS::
     1052
     1053            sage: for p in [1009,2003,10007,100003]:
     1054            ...       K = GF(p)
     1055            ...       for r in (p-1).divisors():
     1056            ...           if r == 1: continue
     1057            ...           x = K.random_element()
     1058            ...           y = x^r
     1059            ...           if y.nth_root(r)**r != y: raise RuntimeError
     1060            ...           if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     1061            ...           if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     1062
     1063            sage: for t in xrange(200):
     1064            ...       n = randint(1,2^63)
     1065            ...       K = Integers(n)
     1066            ...       b = K.random_element()
     1067            ...       e = randint(-2^62, 2^63)
     1068            ...       try:
     1069            ...           a = b.nth_root(e)
     1070            ...           if a^e != b:
     1071            ...               print n, b, e, a
     1072            ...               raise NotImplementedError
     1073            ...       except ValueError:
     1074            ...           pass
     1075
     1076        ALGORITHMS:
     1077
     1078        - The default for prime modulus is currently an algorithm described in the following paper:
     1079
     1080        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     1081
     1082        AUTHORS:
     1083
     1084        - David Roe (2010-2-13)
    10111085        """
    1012 
    1013         # I removed the following text from the docstring, because
    1014         # this behavior is not implemented:
    1015 #             extend -- bool (default: True); if True, return a square
    1016 #                  root in an extension ring, if necessary. Otherwise,
    1017 #                  raise a \class{ValueError} if the square is not in the base
    1018 #                  ring.
    1019 # ...
    1020 #                                                           (if
    1021 #            extend = False) or a NotImplementedError (if extend = True).
    1022 
    1023 
    10241086        if extend:
    10251087            raise NotImplementedError
    1026         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    1027         R = PolynomialRing(self.parent(), "x")
    1028         f = R([-self] + [self.parent()(0)] * (n - 1) + [self.parent()(1)])
    1029         L = f.roots()
    1030         if all:
    1031             return [x[0] for x in L]
    1032         else:
    1033             if len(L) == 0:
    1034                 raise ValueError, "no nth root"
     1088        K = self.parent()
     1089        n = Integer(n)
     1090        if n == 0:
     1091            if self == 1:
     1092                if all: return [K(a) for a in range(1,K.order())]
     1093                else: return self
    10351094            else:
    1036                 return L[0][0]
    1037 
     1095                if all: return []
     1096                else: raise ValueError
     1097        F = K.factored_order()
     1098        if len(F) == 0:
     1099            if all:
     1100                return [self]
     1101            else:
     1102                return self
     1103        if len(F) != 1:
     1104            if all:
     1105                # we should probably do a first pass to see if there are any solutions so that we don't get giant intermediate lists and waste time...
     1106                L = []
     1107                for p, k in F:
     1108                    L.append(mod(self, p**k).nth_root(n, all=True, algorithm=algorithm))
     1109                ans = L[0]
     1110                for i in range(1, len(L)):
     1111                    ans = [a.crt(b) for a in ans for b in L[i]]
     1112            else:
     1113                ans = mod(0,1)
     1114                for p, k in F:
     1115                    ans = ans.crt(mod(self, p**k).nth_root(n, algorithm=algorithm))
     1116            return ans
     1117        p, k = F[0]
     1118        if self.is_zero():
     1119            if n < 0:
     1120                if all: return []
     1121                else: raise ValueError
     1122            if all:
     1123                if k == 1:
     1124                    return [self]
     1125                else:
     1126                    minval = max(1, (k/n).ceil())
     1127                    return [K(a*p**minval) for a in range(p**(k-minval))]
     1128            else:
     1129                return self
     1130        if n < 0:
     1131            try:
     1132                self = ~self
     1133            except ZeroDivisionError:
     1134                if all: return []
     1135                else: raise ValueError
     1136            n = -n
     1137        if p == 2 and k == 1:
     1138            if all: return [self]
     1139            else: return self
     1140        if k > 1:
     1141            pval, upart = self.lift().val_unit(p)
     1142            if not n.divides(pval):
     1143                if all:
     1144                    return []
     1145                else:
     1146                    raise ValueError, "no nth root"
     1147            if pval > 0:
     1148                if all:
     1149                    return [K(a.lift()*p**(pval // n) + p**(k - (pval - pval//n)) * b) for a in mod(upart, p**(k-pval)).nth_root(n, all=True, algorithm=algorithm) for b in range(p**(pval - pval//n))]
     1150                else:
     1151                    return K(p**(pval // n) * mod(upart, p**(k-pval)).nth_root(n, algorithm=algorithm).lift())
     1152            from sage.rings.padics.all import ZpFM
     1153            R = ZpFM(p,k,print_mode='digits')
     1154            self_orig = self
     1155            if p == 2:
     1156                sign = [1]
     1157                if self % 4 == 3:
     1158                    if n % 2 == 0:
     1159                        if all: return []
     1160                        else: raise ValueError, "no nth root"
     1161                    else:
     1162                        sign = [-1]
     1163                        self = -self
     1164                elif n % 2 == 0:
     1165                    if k > 2 and self % 8 == 5:
     1166                        if all: return []
     1167                        else: raise ValueError, "no nth root"
     1168                    sign = [1, -1]
     1169                if k == 2:
     1170                    if all: return [K(s) for s in sign[:2]]
     1171                    else: return K(sign[0])
     1172                if all: modp = [mod(self,8)]
     1173                else: modp = mod(self,8)
     1174            else:
     1175                sign = [1]
     1176                modp = self % p
     1177                self = self / K(R.teichmuller(modp))
     1178                modp = modp.nth_root(n, all=all, algorithm=algorithm)
     1179            # now self is congruent to 1 mod 4 or 1 mod p (for odd p), so the power series for p-adic log converges.
     1180            # Hensel lifting is probably better, but this is easier at the moment.
     1181            plog = R(self).log()
     1182            nval = n.valuation(p)
     1183            if nval >= plog.valuation() + (-1 if p == 2 else 0):
     1184                if self == 1:
     1185                    if all:
     1186                        return [s*K(p*k+m.lift()) for k in range(p**(k-(2 if p==2 else 1))) for m in modp for s in sign]
     1187                    else: return self_orig
     1188                else:
     1189                    if all: return []
     1190                    else: raise ValueError, "no nth root"
     1191            if all:
     1192                ans = [plog // n + p**(k - nval) * i for i in range(p**nval)]
     1193                ans = [s*K(R.teichmuller(m) * a.exp()) for a in ans for m in modp for s in sign]
     1194                return ans
     1195            else:
     1196                return sign[0] * K(R.teichmuller(modp) * (plog // n).exp())
     1197        return self._nth_root_common(n, all, algorithm, cunningham)
     1198
     1199    def _nth_root_naive(self, n):
     1200        """
     1201        Computes all nth roots using brute force, for doc-testing.
     1202
     1203        TESTS::
     1204
     1205            sage: for n in range(2,100): # long time
     1206            ...       K=Integers(n)
     1207            ...       elist = range(1,min(2*n+2,100))
     1208            ...       for e in random_sublist(elist, 5/len(elist)):
     1209            ...           for a in random_sublist(range(1,n), min((n+2)//2,10)/(n-1)):
     1210            ...               b = K(a)
     1211            ...               try:
     1212            ...                   L = b.nth_root(e, all=True)
     1213            ...                   if len(L) > 0:
     1214            ...                       c = b.nth_root(e)
     1215            ...               except:
     1216            ...                   L = [-1]
     1217            ...               M = b._nth_root_naive(e)
     1218            ...               if sorted(L) != M:
     1219            ...                   print "mod(%s, %s).nth_root(%s,all=True), mod(%s, %s)._nth_root_naive(%s)"%(a,n,e,a,n,e)
     1220            ...                   raise ValueError
     1221            ...               if len(L) > 0 and (c not in L):
     1222            ...                   print "mod(%s, %s).nth_root(%s), mod(%s, %s).nth_root(%s,all=True)"%(a,n,e,a,n,e)
     1223            ...                   raise ValueError
     1224        """
     1225        L = []
     1226        for a in self.parent():
     1227            if a**n == self:
     1228                L.append(a)
     1229        return L
    10381230       
    10391231    def _balanced_abs(self):
    10401232        """
  • sage/rings/finite_rings/integer_mod_ring.py

    diff -r 4a9f7b41ae22 -r b6a452dd35f8 sage/rings/finite_rings/integer_mod_ring.py
    a b  
    269269        self._pyx_order = integer_mod.NativeIntStruct(order)
    270270        self.__unit_group_exponent = None
    271271        self.__factored_order = None
     272        self.__factored_unit_order = None
    272273        quotient_ring.QuotientRing_generic.__init__(self, ZZ, ZZ.ideal(order), names=None)
    273274        if category is None:
    274275            from sage.categories.commutative_rings import CommutativeRings
     
    688689        """
    689690        if self.__factored_order is not None:
    690691            return self.__factored_order
    691         self.__factored_order = factor(self.__order, int_=True)
     692        self.__factored_order = factor(self.__order, int_=(self.__order < 2**31))
    692693        return self.__factored_order
     694
     695    def factored_unit_order(self):
     696        """
     697        Returns a list of Factorization objects, each the factorization of the order of the units in a `\Z / p^n\Z` component of this group (using the Chinese Remainder Theorem).
     698
     699        EXAMPLES::
     700
     701            sage: R = Integers(8*9*25*17*29)
     702            sage: R.factored_unit_order()
     703            [2^2, 2 * 3, 2^2 * 5, 2^4, 2^2 * 7]
     704        """
     705        ans = []
     706        from sage.structure.factorization import Factorization
     707        for p, e in self.factored_order():
     708            ans.append(Factorization([(p,e-1)]) * factor(p-1, int_=(self.__order < 2**31)))
     709        return ans
    693710   
    694711    def characteristic(self):
    695712        """