Ticket #7931: 7931_nth_root.patch

File 7931_nth_root.patch, 19.1 KB (added by roed, 10 years ago)

Should fix the problem observed earlier

  • sage/rings/finite_rings/element_base.pyx

    # HG changeset patch
    # User David Roe <roed@math.harvard.edu>
    # Date 1266938910 28800
    # Node ID 6072d15e7b6e28cdd51a403c7c0be9f6495a1912
    # Parent  349e0c4b8392d69f9530c1fa1924ef82b29ab4d3
    #7931: implements a better algorithm for nth roots in finite fields
    and Z/nZ.
    
    diff -r 349e0c4b8392 -r 6072d15e7b6e sage/rings/finite_rings/element_base.pyx
    a b  
    396396
    397397        return order
    398398
    399     def nth_root(self, int n, extend = False, all = False):
     399    def nth_root(self, n, extend = False, all = False, algorithm=None):
    400400        r"""
    401401        Returns an nth root of self.
    402402
    403403        INPUT:
    404404
    405         - ``n`` - integer >= 1 (must fit in C int type)
     405        - ``n`` - integer >= 1
    406406
    407407        - ``extend`` - bool (default: True); if True, return an nth
    408408          root in an extension ring, if necessary. Otherwise, raise a
     
    412412        - ``all`` - bool (default: False); if True, return all nth
    413413          roots of self, instead of just one.
    414414
     415        - ``algorithm`` - string (default: None); 'Johnston' is the only
     416          currently supported option.
     417
    415418        OUTPUT:
    416419
    417420        If self has an nth root, returns one (if all == False) or a
     
    423426       
    424427           The 'extend' option is not implemented (yet).
    425428
     429        EXAMPLES::
     430
     431            sage: K = GF(31)
     432            sage: a = K(22)
     433            sage: K(22).nth_root(7)
     434            13
     435            sage: K(25).nth_root(5)
     436            5
     437            sage: K(23).nth_root(3)
     438            29
     439            sage: K.<a> = GF(625)
     440            sage: (3*a^2+a+1).nth_root(13)**13
     441            3*a^2 + a + 1
     442
     443            sage: k.<a> = GF(29^2)
     444            sage: b = a^2 + 5*a + 1
     445            sage: b.nth_root(11)
     446            3*a + 20
     447            sage: b.nth_root(5)
     448            Traceback (most recent call last):
     449            ...
     450            ValueError: no nth root
     451            sage: b.nth_root(5, all = True)
     452            []
     453            sage: b.nth_root(3, all = True)
     454            [14*a + 18, 10*a + 13, 5*a + 27]
     455
     456            sage: k.<a> = GF(29^5)
     457            sage: b = a^2 + 5*a + 1
     458            sage: b.nth_root(5)
     459            19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
     460            sage: b.nth_root(7)
     461            Traceback (most recent call last):
     462            ...
     463            ValueError: no nth root
     464            sage: b.nth_root(4, all=True)
     465            []
     466
     467        TESTS::
     468
     469            sage: for p in [2,3,5,7,11]: # long
     470            ...       for n in [2,5,10]:
     471            ...           q = p^n
     472            ...           K.<a> = GF(q)
     473            ...           for r in (q-1).divisors():
     474            ...               if r == 1: continue
     475            ...               x = K.random_element()
     476            ...               y = x^r
     477            ...               if y.nth_root(r)**r != y: raise RuntimeError
     478            ...               if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     479            ...               if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     480
     481        ALGORITHMS:
     482
     483        - The default is currently an algorithm described in the following paper:
     484
     485        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     486
    426487        AUTHOR:
    427488
    428         - David Roe (2007-10-3)
    429 
     489        - David Roe (2010-02-13)
    430490        """
    431         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    432         from sage.rings.integer import Integer
     491        if self.is_zero():
     492            if n < 0:
     493                raise ZeroDivisionError
     494            if all: return [self]
     495            else: return self
     496        if algorithm is None:
     497            algorithm = 'Johnston'
    433498        if extend:
    434499            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]
     500        from sage.rings.integer import Integer
     501        n = Integer(n)
     502        K = self.parent()
     503        q = K.order()
     504        if self.is_one():
     505            gcd = n.gcd(q-1)
     506            if gcd == 1:
     507                if all: return [self]
     508                else: return self
     509            else:
     510                g = K.multiplicative_generator()
     511                q1overn = (q-1)//gcd
     512                nthroot = g**q1overn
     513                if all:
     514                    return [nthroot**a for a in range(gcd)]
     515                else:
     516                    return nthroot
     517        n = n % (q-1)
     518        gcd, alpha, beta = n.xgcd(q-1)
     519        if gcd == 1:
     520            if all: return [self**alpha]
     521            else: return self**alpha
     522        n = gcd
     523        self = self**alpha
     524        F = n.factor()
     525        from sage.groups.generic import discrete_log
     526        if algorithm == 'Johnston':
     527            q1overn = (q-1)//n
     528            if self**q1overn != 1:
     529                if all:
     530                    return []
     531                else:
     532                    raise ValueError, "no nth root"
     533            g = K.multiplicative_generator()
     534            for r, v in F:
     535                k, h = (q-1).val_unit(r)
     536                z = h * (-h).inverse_mod(r**v)
     537                x = (1 + z) // r**v
     538                if k == 1:
     539                    self = self**x
     540                else:
     541                    t = discrete_log(self**h, g**(r**v*h), r**(k-v), operation='*')
     542                    self = self**x * g**(-z*t)
     543            if all:
     544                nthroot = g**q1overn
     545                L = [self]
     546                for i in range(1,n):
     547                    self *= nthroot
     548                    L.append(self)
     549                return L
     550            else:
     551                return self
    440552        else:
    441             if len(L) == 0:
    442                 raise ValueError, "no nth root"
    443             else:
    444                 return L[0][0]
     553            raise ValueError, "unknown algorithm"
    445554
    446555    def additive_order(self):
    447556        """
  • sage/rings/finite_rings/element_ext_pari.py

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

    diff -r 349e0c4b8392 -r 6072d15e7b6e sage/rings/finite_rings/integer_mod.pyx
    a b  
    950950               
    951951    square_root = sqrt
    952952
    953     def nth_root(self, int n, extend = False, all = False):
     953    def nth_root(self, n, extend = False, all = False, algorithm = None):
    954954        r"""
    955955        Returns an `n`\th root of ``self``.
    956956       
    957957        INPUT:
    958        
    959        
    960         -  ``n`` - integer `\geq 1` (must fit in C
    961            ``int`` type)
    962        
    963         -  ``all`` - bool (default: ``False``); if
    964            ``True``, return all `n`\th roots of
    965            ``self``, instead of just one.
    966        
    967        
    968         OUTPUT: If self has an `n`\th root, returns one (if
    969         ``all`` is false) or a list of all of them (if
    970         ``all`` is true). Otherwise, raises a
    971         ``ValueError``.
    972        
    973         AUTHORS:
    974958
    975         - David Roe (2007-10-3)
     959        - ``n`` - integer `\geq 1`
     960
     961        - ``extend`` - bool (default: True); if True, return an nth
     962          root in an extension ring, if necessary. Otherwise, raise a
     963          ValueError if the root is not in the base ring.  Warning:
     964          this option is not implemented!
     965
     966        - ``all`` - bool (default: ``False``); if ``True``, return all `n`\th
     967          roots of ``self``, instead of just one.
     968
     969        - ``algorithm`` - string (default: None); The algorithm for the prime modulus case.
     970          CRT and p-adic log techniques are used to reduce to this case.
     971          'Johnston' is the only currently supported option.
     972
     973        OUTPUT:
     974
     975        If self has an `n`\th root, returns one (if ``all`` is ``False``) or a
     976        list of all of them (if ``all`` is ``True``).  Otherwise, raises a
     977        ``ValueError`` (if ``extend`` is ``False``) or a ``NotImplementedError`` (if
     978        ``extend`` is ``True``).
     979
     980        .. warning::
     981           The 'extend' option is not implemented (yet).
    976982       
    977983        EXAMPLES::
    978984       
    979             sage: k.<a> = GF(29)
    980             sage: b = a^2 + 5*a + 1
    981             sage: b.nth_root(5)
    982             24
    983             sage: b.nth_root(7)
    984             Traceback (most recent call last):
    985             ...
    986             ValueError: no nth root
    987             sage: b.nth_root(4, all=True)
    988             [21, 20, 9, 8]
     985
     986            sage: K = GF(31)
     987            sage: a = K(22)
     988            sage: K(22).nth_root(7)
     989            13
     990            sage: K(25).nth_root(5)
     991            5
     992            sage: K(23).nth_root(3)
     993            29
     994            sage: mod(225,2^5*3^2).nth_root(4, all=True)
     995            [225, 129, 33, 63, 255, 159, 9, 201, 105, 279, 183, 87, 81, 273, 177, 207, 111, 15, 153, 57, 249, 135, 39, 231]
     996            sage: mod(275,2^5*7^4).nth_root(7, all=True)
     997            [58235, 25307, 69211, 36283, 3355, 47259, 14331]
     998            sage: mod(1,8).nth_root(2,all=True)
     999            [1, 7, 5, 3]
     1000            sage: mod(4,8).nth_root(2,all=True)
     1001            [2, 6]
     1002            sage: mod(1,16).nth_root(4,all=True)
     1003            [1, 15, 13, 3, 9, 7, 5, 11]
     1004            sage: (mod(22,31)^200).nth_root(200)
     1005            5
     1006
     1007        TESTS::
     1008
     1009            sage: for p in [1009,2003,10007,100003]: # long
     1010            ...       K = GF(p)
     1011            ...       for r in (p-1).divisors():
     1012            ...           if r == 1: continue
     1013            ...           x = K.random_element()
     1014            ...           y = x^r
     1015            ...           if y.nth_root(r)**r != y: raise RuntimeError
     1016            ...           if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     1017            ...           if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     1018
     1019        ALGORITHMS:
     1020
     1021        - The default for prime modulus is currently an algorithm described in the following paper:
     1022
     1023        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     1024
     1025        AUTHORS:
     1026
     1027        - David Roe (2010-2-13)
    9891028        """
    990 
    991         # I removed the following text from the docstring, because
    992         # this behavior is not implemented:
    993 #             extend -- bool (default: True); if True, return a square
    994 #                  root in an extension ring, if necessary. Otherwise,
    995 #                  raise a \class{ValueError} if the square is not in the base
    996 #                  ring.
    997 # ...
    998 #                                                           (if
    999 #            extend = False) or a NotImplementedError (if extend = True).
    1000 
    1001 
    10021029        if extend:
    10031030            raise NotImplementedError
    1004         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    1005         R = PolynomialRing(self.parent(), "x")
    1006         f = R([-self] + [self.parent()(0)] * (n - 1) + [self.parent()(1)])
    1007         L = f.roots()
    1008         if all:
    1009             return [x[0] for x in L]
     1031        K = self.parent()
     1032        n = Integer(n)
     1033        F = K.factored_order()
     1034        if len(F) == 0:
     1035            if all:
     1036                return [self]
     1037            else:
     1038                return self
     1039        if len(F) != 1:
     1040            if all:
     1041                # 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...
     1042                L = []
     1043                for p, k in F:
     1044                    L.append(mod(self, p**k).nth_root(n, all=True, algorithm=algorithm))
     1045                ans = L[0]
     1046                for i in range(1, len(L)):
     1047                    ans = [a.crt(b) for a in ans for b in L[i]]
     1048            else:
     1049                ans = mod(0,1)
     1050                for p, k in F:
     1051                    ans = ans.crt(mod(self, p**k).nth_root(n, algorithm=algorithm))
     1052            return ans
     1053        p, k = F[0]
     1054        if self == 0:
     1055            if n < 0:
     1056                raise ZeroDivisionError
     1057            if all:
     1058                if k == 1:
     1059                    return [self]
     1060                else:
     1061                    minval = max(1, (k/n).ceil())
     1062                    return [K(a*p**minval) for a in range(p**(k-minval))]
     1063            else:
     1064                return self
     1065        if p == 2 and k == 1:
     1066            if all: return [self]
     1067            else: return self
     1068        if k > 1:
     1069            pval, upart = self.lift().val_unit(p)
     1070            if not n.divides(pval):
     1071                if all:
     1072                    return []
     1073                else:
     1074                    raise ValueError, "no nth root"
     1075            if pval > 0:
     1076                if all:
     1077                    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))]
     1078                else:
     1079                    return K(p**(pval // n) * mod(upart.lift(), p**(k-pval)).nth_root(n, algorithm=algorithm).lift())
     1080            from sage.rings.padics.all import ZpFM
     1081            if p == 2:
     1082                if self % 4 == 3:
     1083                    if n % 2 == 0:
     1084                        if all: return []
     1085                        else: raise ValueError, "no nth root"
     1086                    else:
     1087                        sign = [-1]
     1088                        self = -self
     1089                elif n % 2 == 0:
     1090                    sign = [1, -1]
     1091                else:
     1092                    sign = [1]
     1093                if k == 2:
     1094                    if all: return [K(s) for s in sign]
     1095                    else: return K(sign[0])
     1096                if all: modp = [mod(1,2)]
     1097                else: modp = mod(1,2)
     1098            else:
     1099                sign = [1]
     1100                modp = self % p
     1101                self = self / K(ZpFM(p,k).teichmuller(modp))
     1102                modp = modp.nth_root(n, all=all, algorithm=algorithm)
     1103            # now self is congruent to 1 mod 4 or 1 mod p (for odd p), so the power series for p-adic log converges.
     1104            # Hensel lifting is probably better, but this is easier at the moment.
     1105            R = ZpFM(p,k)
     1106            plog = R(self).log()
     1107            nval = n.valuation(p)
     1108            if nval >= plog.valuation() + (-1 if p == 2 else 0):
     1109                if self == 1:
     1110                    if all: return [self]
     1111                    else: return self
     1112                else:
     1113                    if all: return []
     1114                    else: raise ValueError, "no nth root"
     1115            if all:
     1116                ans = [plog // n + p**(k - nval) * i for i in range(p**nval)]
     1117                ans = [s*K(R.teichmuller(m) * a.exp()) for a in ans for m in modp for s in sign]
     1118                return ans
     1119            else:
     1120                return sign[0] * K(R.teichmuller(modp) * (plog // n).exp())
     1121        if n < 0:
     1122            self = ~self
     1123            n = -n
     1124        if algorithm is None:
     1125            algorithm = 'Johnston'
     1126        from sage.groups.generic import discrete_log
     1127        p = K.order()
     1128        if self == 1:
     1129            if all:
     1130                gcd = n.gcd(p-1)
     1131                # the following may eventually be improved to not need a multiplicative generator.
     1132                base = K.multiplicative_generator()**((p-1)//gcd)
     1133                return [base**a for a in range(gcd)]
     1134            else:
     1135                return self
     1136        n = n % (p-1)
     1137        gcd, alpha, beta = n.xgcd(p-1) # gcd = alpha*n + beta*(p-1), so 1/n = alpha/gcd (mod p-1)
     1138        if gcd == 1:
     1139            if all:
     1140                return [self**alpha]
     1141            else:
     1142                return self**alpha
     1143        n = Integer(gcd)
     1144        self = self**alpha
     1145        F = n.factor()
     1146        if algorithm == 'Johnston':
     1147            p1overn = (p-1)//n
     1148            if self**p1overn != 1:
     1149                if all:
     1150                    return []
     1151                else:
     1152                    raise ValueError, "no nth root"
     1153            g = K.multiplicative_generator()
     1154            for r, v in F:
     1155                k, h = (p-1).val_unit(r)
     1156                z = h * (-h).inverse_mod(r**v)
     1157                x = (1 + z) // r**v
     1158                if k == 1:
     1159                    self = self**x
     1160                else:
     1161                    t = discrete_log(self**h, g**(r**v*h), r**(k-v), operation='*')
     1162                    self = self**x * g**(-z*t)
     1163            if all:
     1164                nthroot = g**p1overn
     1165                L = [self]
     1166                for i in range(1,n):
     1167                    self *= nthroot
     1168                    L.append(self)
     1169                return L
     1170            else:
     1171                return self
    10101172        else:
    1011             if len(L) == 0:
    1012                 raise ValueError, "no nth root"
    1013             else:
    1014                 return L[0][0]
     1173            raise ValueError, "unknown algorithm"
    10151174
    10161175       
    10171176    def _balanced_abs(self):