Ticket #7931: 7931_nth_root.2.patch

File 7931_nth_root.2.patch, 22.0 KB (added by roed, 9 years ago)

Apply just this patch

  • sage/rings/finite_rings/element_base.pyx

    # HG changeset patch
    # User David Roe <roed@math.harvard.edu>
    # Date 1266938910 28800
    # Node ID 6d6c6a580546f87b06e55064ccd269a9e001abb6
    # Parent  a3d717d075ecd6533ee7ceb0bd149ac4f36a422a
    #7931: implements a better algorithm for nth roots in finite fields and Z/nZ.
    
    diff -r a3d717d075ec -r 6d6c6a580546 sage/rings/finite_rings/element_base.pyx
    a b  
    386386        if self.is_zero():
    387387            raise ArithmeticError, "Multiplicative order of 0 not defined."
    388388        n = self._parent.order() - 1
     389        F = self._parent.factored_unit_order()[0]
    389390        order = 1
    390         for p, e in sage.rings.arith.factor(n):
     391        for p, e in F:
    391392            # Determine the power of p that divides the order.
    392             a = self**(n/(p**e))
     393            a = self**(n//(p**e))
    393394            while a != 1:
    394395                order = order * p
    395396                a = a**p
    396397
    397398        return order
    398399
    399     def nth_root(self, int n, extend = False, all = False):
     400    def nth_root(self, n, extend = False, all = False, algorithm=None, cunningham=False):
    400401        r"""
    401402        Returns an nth root of self.
    402403
    403404        INPUT:
    404405
    405         - ``n`` - integer >= 1 (must fit in C int type)
     406        - ``n`` - integer >= 1
    406407
    407408        - ``extend`` - bool (default: True); if True, return an nth
    408409          root in an extension ring, if necessary. Otherwise, raise a
     
    412413        - ``all`` - bool (default: False); if True, return all nth
    413414          roots of self, instead of just one.
    414415
     416        - ``algorithm`` - string (default: None); 'Johnston' is the only
     417          currently supported option.
     418
    415419        OUTPUT:
    416420
    417421        If self has an nth root, returns one (if all == False) or a
     
    423427       
    424428           The 'extend' option is not implemented (yet).
    425429
     430        EXAMPLES::
     431
     432            sage: K = GF(31)
     433            sage: a = K(22)
     434            sage: K(22).nth_root(7)
     435            13
     436            sage: K(25).nth_root(5)
     437            5
     438            sage: K(23).nth_root(3)
     439            29
     440            sage: K.<a> = GF(625)
     441            sage: (3*a^2+a+1).nth_root(13)**13
     442            3*a^2 + a + 1
     443
     444            sage: k.<a> = GF(29^2)
     445            sage: b = a^2 + 5*a + 1
     446            sage: b.nth_root(11)
     447            3*a + 20
     448            sage: b.nth_root(5)
     449            Traceback (most recent call last):
     450            ...
     451            ValueError: no nth root
     452            sage: b.nth_root(5, all = True)
     453            []
     454            sage: b.nth_root(3, all = True)
     455            [14*a + 18, 10*a + 13, 5*a + 27]
     456
     457            sage: k.<a> = GF(29^5)
     458            sage: b = a^2 + 5*a + 1
     459            sage: b.nth_root(5)
     460            19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
     461            sage: b.nth_root(7)
     462            Traceback (most recent call last):
     463            ...
     464            ValueError: no nth root
     465            sage: b.nth_root(4, all=True)
     466            []
     467
     468        TESTS::
     469
     470            sage: for p in [2,3,5,7,11]: # long
     471            ...       for n in [2,5,10]:
     472            ...           q = p^n
     473            ...           K.<a> = GF(q)
     474            ...           for r in (q-1).divisors():
     475            ...               if r == 1: continue
     476            ...               x = K.random_element()
     477            ...               y = x^r
     478            ...               if y.nth_root(r)**r != y: raise RuntimeError
     479            ...               if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     480            ...               if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     481
     482        ALGORITHMS:
     483
     484        - The default is currently an algorithm described in the following paper:
     485
     486        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     487
    426488        AUTHOR:
    427489
    428         - David Roe (2007-10-3)
    429 
     490        - David Roe (2010-02-13)
    430491        """
    431         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    432         from sage.rings.integer import Integer
     492        if self.is_zero():
     493            if n < 0:
     494                raise ZeroDivisionError
     495            if all: return [self]
     496            else: return self
    433497        if extend:
    434498            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]
     499        from sage.rings.integer import Integer
     500        n = Integer(n)
     501        return nth_root_common(self, n, all, algorithm, cunningham)
    445502
    446503    def additive_order(self):
    447504        """
     
    516573            b^11 + b^10 + b^9 + b^7 + b^5 + b^4 + b^2 + b
    517574        """
    518575        return self.pth_power(-k)
     576
     577def nth_root_common(self, n, all, algorithm, cunningham):
     578    """
     579    This function exists to reduce code duplication between finite field nth roots and integer_mod nth roots.
     580
     581    The inputs are described there.
     582
     583    TESTS::
     584
     585        sage: a = Zmod(17)(13)
     586        sage: from sage.rings.finite_rings.element_base import nth_root_common
     587        sage: nth_root_common(a, 4, True, "Johnston", False)
     588        [3, 5, 14, 12]
     589    """
     590    K = self.parent()
     591    q = K.order()
     592    if self.is_one():
     593        gcd = n.gcd(q-1)
     594        if gcd == 1:
     595            if all: return [self]
     596            else: return self
     597        else:
     598            # the following may eventually be improved to not need a multiplicative generator.
     599            g = K.multiplicative_generator()
     600            q1overn = (q-1)//gcd
     601            nthroot = g**q1overn
     602            return [nthroot**a for a in range(gcd)] if all else nthroot
     603    n = n % (q-1)
     604    gcd, alpha, beta = n.xgcd(q-1) # gcd = alpha*n + beta*(p-1), so 1/n = alpha/gcd (mod p-1)
     605    if gcd == 1:
     606        return [self**alpha] if all else self**alpha
     607    n = gcd
     608    self = self**alpha
     609    if cunningham:
     610        F = n._factor_cunningham()
     611    else:
     612        F = n.factor()
     613    from sage.groups.generic import discrete_log
     614    if algorithm is None or algorithm == 'Johnston':
     615        q1overn = (q-1)//n
     616        if self**q1overn != 1:
     617            if all: return []
     618            else: raise ValueError, "no nth root"
     619        g = K.multiplicative_generator()
     620        for r, v in F:
     621            k, h = (q-1).val_unit(r)
     622            z = h * (-h).inverse_mod(r**v)
     623            x = (1 + z) // r**v
     624            if k == 1:
     625                self = self**x
     626            else:
     627                t = discrete_log(self**h, g**(r**v*h), r**(k-v), operation='*')
     628                self = self**x * g**(-z*t)
     629        if all:
     630            nthroot = g**q1overn
     631            L = [self]
     632            for i in range(1,n):
     633                self *= nthroot
     634                L.append(self)
     635            return L
     636        else:
     637            return self
     638    else:
     639        raise ValueError, "unknown algorithm"
  • sage/rings/finite_rings/element_ext_pari.py

    diff -r a3d717d075ec -r 6d6c6a580546 sage/rings/finite_rings/element_ext_pari.py
    a b  
    286286        """
    287287        return self.square_root(extend=extend, all=all)
    288288
    289     def nth_root(self, n, extend = False, all = False):
    290         r"""
    291         Returns an nth root of self.
    292        
    293         INPUT:
    294        
    295        
    296         -  ``n`` - integer = 1 (must fit in C int type)
    297        
    298         -  ``extend`` - bool (default: True); if True, return
    299            an nth root in an extension ring, if necessary. Otherwise, raise a
    300            ValueError if the root is not in the base ring. Warning: this
    301            option is not implemented!
    302        
    303         -  ``all`` - bool (default: False); if True, return all
    304            nth roots of self, instead of just one.
    305        
    306        
    307         OUTPUT: If self has an nth root, returns one (if all == False) or a
    308         list of all of them (if all == True). Otherwise, raises a
    309         ValueError (if extend = False) or a NotImplementedError (if extend
    310         = True).
    311        
    312         .. warning::
    313 
    314            The 'extend' option is not implemented (yet).
    315        
    316         AUTHORS:
    317 
    318         - David Roe (2007-10-3)
    319        
    320         EXAMPLES::
    321        
    322             sage: k.<a> = GF(29^5)
    323             sage: b = a^2 + 5*a + 1
    324             sage: b.nth_root(5)
    325             19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
    326             sage: b.nth_root(7)
    327             Traceback (most recent call last):
    328             ...
    329             ValueError: no nth root
    330             sage: b.nth_root(4, all=True)
    331             []
    332         """
    333         if extend:
    334             raise NotImplementedError
    335         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    336         R = PolynomialRing(self.parent(), "x")
    337         f = R([-self] + [self.parent()(0)] * (n - 1) + [self.parent()(1)])
    338         L = f.roots()
    339         if all:
    340             return [x[0] for x in L]
    341         else:
    342             if len(L) == 0:
    343                 raise ValueError, "no nth root"
    344             else:
    345                 return L[0][0]
    346        
    347289    def rational_reconstruction(self):
    348290        """
    349291        If the parent field is a prime field, uses rational reconstruction
  • sage/rings/finite_rings/finite_field_base.pxd

    diff -r a3d717d075ec -r 6d6c6a580546 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 a3d717d075ec -r 6d6c6a580546 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.pyx

    diff -r a3d717d075ec -r 6d6c6a580546 sage/rings/finite_rings/integer_mod.pyx
    a b  
    974974               
    975975    square_root = sqrt
    976976
    977     def nth_root(self, int n, extend = False, all = False):
     977    def nth_root(self, n, extend = False, all = False, algorithm = None, cunningham = False):
    978978        r"""
    979979        Returns an `n`\th root of ``self``.
    980980       
    981981        INPUT:
    982        
    983        
    984         -  ``n`` - integer `\geq 1` (must fit in C
    985            ``int`` type)
    986        
    987         -  ``all`` - bool (default: ``False``); if
    988            ``True``, return all `n`\th roots of
    989            ``self``, instead of just one.
    990        
    991        
    992         OUTPUT: If self has an `n`\th root, returns one (if
    993         ``all`` is false) or a list of all of them (if
    994         ``all`` is true). Otherwise, raises a
    995         ``ValueError``.
    996        
    997         AUTHORS:
    998982
    999         - David Roe (2007-10-3)
     983        - ``n`` - integer `\geq 1`
     984
     985        - ``extend`` - bool (default: True); if True, return an nth
     986          root in an extension ring, if necessary. Otherwise, raise a
     987          ValueError if the root is not in the base ring.  Warning:
     988          this option is not implemented!
     989
     990        - ``all`` - bool (default: ``False``); if ``True``, return all `n`\th
     991          roots of ``self``, instead of just one.
     992
     993        - ``algorithm`` - string (default: None); The algorithm for the prime modulus case.
     994          CRT and p-adic log techniques are used to reduce to this case.
     995          'Johnston' is the only currently supported option.
     996
     997        OUTPUT:
     998
     999        If self has an `n`\th root, returns one (if ``all`` is ``False``) or a
     1000        list of all of them (if ``all`` is ``True``).  Otherwise, raises a
     1001        ``ValueError`` (if ``extend`` is ``False``) or a ``NotImplementedError`` (if
     1002        ``extend`` is ``True``).
     1003
     1004        .. warning::
     1005           The 'extend' option is not implemented (yet).
    10001006       
    10011007        EXAMPLES::
    10021008       
    1003             sage: k.<a> = GF(29)
    1004             sage: b = a^2 + 5*a + 1
    1005             sage: b.nth_root(5)
    1006             24
    1007             sage: b.nth_root(7)
    1008             Traceback (most recent call last):
    1009             ...
    1010             ValueError: no nth root
    1011             sage: b.nth_root(4, all=True)
    1012             [21, 20, 9, 8]
     1009
     1010            sage: K = GF(31)
     1011            sage: a = K(22)
     1012            sage: K(22).nth_root(7)
     1013            13
     1014            sage: K(25).nth_root(5)
     1015            5
     1016            sage: K(23).nth_root(3)
     1017            29
     1018            sage: mod(225,2^5*3^2).nth_root(4, all=True)
     1019            [225, 129, 33, 63, 255, 159, 9, 201, 105, 279, 183, 87, 81, 273, 177, 207, 111, 15, 153, 57, 249, 135, 39, 231]
     1020            sage: mod(275,2^5*7^4).nth_root(7, all=True)
     1021            [58235, 25307, 69211, 36283, 3355, 47259, 14331]
     1022            sage: mod(1,8).nth_root(2,all=True)
     1023            [1, 7, 5, 3]
     1024            sage: mod(4,8).nth_root(2,all=True)
     1025            [2, 6]
     1026            sage: mod(1,16).nth_root(4,all=True)
     1027            [1, 15, 13, 3, 9, 7, 5, 11]
     1028            sage: (mod(22,31)^200).nth_root(200)
     1029            5
     1030
     1031        TESTS::
     1032
     1033            sage: for p in [1009,2003,10007,100003]: # long
     1034            ...       K = GF(p)
     1035            ...       for r in (p-1).divisors():
     1036            ...           if r == 1: continue
     1037            ...           x = K.random_element()
     1038            ...           y = x^r
     1039            ...           if y.nth_root(r)**r != y: raise RuntimeError
     1040            ...           if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     1041            ...           if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     1042
     1043        ALGORITHMS:
     1044
     1045        - The default for prime modulus is currently an algorithm described in the following paper:
     1046
     1047        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     1048
     1049        AUTHORS:
     1050
     1051        - David Roe (2010-2-13)
    10131052        """
    1014 
    1015         # I removed the following text from the docstring, because
    1016         # this behavior is not implemented:
    1017 #             extend -- bool (default: True); if True, return a square
    1018 #                  root in an extension ring, if necessary. Otherwise,
    1019 #                  raise a \class{ValueError} if the square is not in the base
    1020 #                  ring.
    1021 # ...
    1022 #                                                           (if
    1023 #            extend = False) or a NotImplementedError (if extend = True).
    1024 
    1025 
    10261053        if extend:
    10271054            raise NotImplementedError
    1028         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    1029         R = PolynomialRing(self.parent(), "x")
    1030         f = R([-self] + [self.parent()(0)] * (n - 1) + [self.parent()(1)])
    1031         L = f.roots()
    1032         if all:
    1033             return [x[0] for x in L]
    1034         else:
    1035             if len(L) == 0:
    1036                 raise ValueError, "no nth root"
     1055        K = self.parent()
     1056        n = Integer(n)
     1057        F = K.factored_order()
     1058        if len(F) == 0:
     1059            if all:
     1060                return [self]
    10371061            else:
    1038                 return L[0][0]
    1039 
     1062                return self
     1063        if len(F) != 1:
     1064            if all:
     1065                # 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...
     1066                L = []
     1067                for p, k in F:
     1068                    L.append(mod(self, p**k).nth_root(n, all=True, algorithm=algorithm))
     1069                ans = L[0]
     1070                for i in range(1, len(L)):
     1071                    ans = [a.crt(b) for a in ans for b in L[i]]
     1072            else:
     1073                ans = mod(0,1)
     1074                for p, k in F:
     1075                    ans = ans.crt(mod(self, p**k).nth_root(n, algorithm=algorithm))
     1076            return ans
     1077        p, k = F[0]
     1078        if self.is_zero():
     1079            if n < 0:
     1080                raise ZeroDivisionError
     1081            if all:
     1082                if k == 1:
     1083                    return [self]
     1084                else:
     1085                    minval = max(1, (k/n).ceil())
     1086                    return [K(a*p**minval) for a in range(p**(k-minval))]
     1087            else:
     1088                return self
     1089        if n < 0:
     1090            self = ~self
     1091            n = -n
     1092        if p == 2 and k == 1:
     1093            if all: return [self]
     1094            else: return self
     1095        if k > 1:
     1096            pval, upart = self.lift().val_unit(p)
     1097            if not n.divides(pval):
     1098                if all:
     1099                    return []
     1100                else:
     1101                    raise ValueError, "no nth root"
     1102            if pval > 0:
     1103                if all:
     1104                    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))]
     1105                else:
     1106                    return K(p**(pval // n) * mod(upart.lift(), p**(k-pval)).nth_root(n, algorithm=algorithm).lift())
     1107            from sage.rings.padics.all import ZpFM
     1108            if p == 2:
     1109                if self % 4 == 3:
     1110                    if n % 2 == 0:
     1111                        if all: return []
     1112                        else: raise ValueError, "no nth root"
     1113                    else:
     1114                        sign = [-1]
     1115                        self = -self
     1116                elif n % 2 == 0:
     1117                    sign = [1, -1]
     1118                else:
     1119                    sign = [1]
     1120                if k == 2:
     1121                    if all: return [K(s) for s in sign]
     1122                    else: return K(sign[0])
     1123                if all: modp = [mod(1,2)]
     1124                else: modp = mod(1,2)
     1125            else:
     1126                sign = [1]
     1127                modp = self % p
     1128                self = self / K(ZpFM(p,k).teichmuller(modp))
     1129                modp = modp.nth_root(n, all=all, algorithm=algorithm)
     1130            # now self is congruent to 1 mod 4 or 1 mod p (for odd p), so the power series for p-adic log converges.
     1131            # Hensel lifting is probably better, but this is easier at the moment.
     1132            R = ZpFM(p,k)
     1133            plog = R(self).log()
     1134            nval = n.valuation(p)
     1135            if nval >= plog.valuation() + (-1 if p == 2 else 0):
     1136                if self == 1:
     1137                    if all: return [self]
     1138                    else: return self
     1139                else:
     1140                    if all: return []
     1141                    else: raise ValueError, "no nth root"
     1142            if all:
     1143                ans = [plog // n + p**(k - nval) * i for i in range(p**nval)]
     1144                ans = [s*K(R.teichmuller(m) * a.exp()) for a in ans for m in modp for s in sign]
     1145                return ans
     1146            else:
     1147                return sign[0] * K(R.teichmuller(modp) * (plog // n).exp())
     1148        from element_base import nth_root_common
     1149        return nth_root_common(self, n, all, algorithm, cunningham)
    10401150       
    10411151    def _balanced_abs(self):
    10421152        """
  • sage/rings/finite_rings/integer_mod_ring.py

    diff -r a3d717d075ec -r 6d6c6a580546 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        """