Ticket #7931: 7931.patch

File 7931.patch, 21.2 KB (added by roed, 10 years ago)
  • sage/rings/finite_field_element.py

    # HG changeset patch
    # User David Roe <roed@math.harvard.edu>
    # Date 1263482240 18000
    # Node ID be9e35864adfafee742794e6d80e9c844aa295c8
    # Parent  7fdc5c05c22090674fa14eda3f9ec6e0f06caf6f
    Updates nth_root methods for finite fields and Integers(n) to use specialized algorithms rather than just polynomial factoring.
    
    diff -r 7fdc5c05c220 -r be9e35864adf sage/rings/finite_field_element.py
    a b  
    299299        """
    300300        return self.square_root(extend=extend, all=all)
    301301
    302     def nth_root(self, n, extend = False, all = False):
    303         r"""
    304         Returns an nth root of self.
    305        
    306         INPUT:
    307        
    308        
    309         -  ``n`` - integer = 1 (must fit in C int type)
    310        
    311         -  ``extend`` - bool (default: True); if True, return
    312            an nth root in an extension ring, if necessary. Otherwise, raise a
    313            ValueError if the root is not in the base ring. Warning: this
    314            option is not implemented!
    315        
    316         -  ``all`` - bool (default: False); if True, return all
    317            nth roots of self, instead of just one.
    318        
    319        
    320         OUTPUT: If self has an nth root, returns one (if all == False) or a
    321         list of all of them (if all == True). Otherwise, raises a
    322         ValueError (if extend = False) or a NotImplementedError (if extend
    323         = True).
    324        
    325         .. warning::
    326 
    327            The 'extend' option is not implemented (yet).
    328        
    329         AUTHORS:
    330 
    331         - David Roe (2007-10-3)
    332        
    333         EXAMPLES::
    334        
    335             sage: k.<a> = GF(29^5)
    336             sage: b = a^2 + 5*a + 1
    337             sage: b.nth_root(5)
    338             19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
    339             sage: b.nth_root(7)
    340             Traceback (most recent call last):
    341             ...
    342             ValueError: no nth root
    343             sage: b.nth_root(4, all=True)
    344             []
    345         """
    346         if extend:
    347             raise NotImplementedError
    348         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    349         R = PolynomialRing(self.parent(), "x")
    350         f = R([-self] + [self.parent()(0)] * (n - 1) + [self.parent()(1)])
    351         L = f.roots()
    352         if all:
    353             return [x[0] for x in L]
    354         else:
    355             if len(L) == 0:
    356                 raise ValueError, "no nth root"
    357             else:
    358                 return L[0][0]
    359        
    360302    def rational_reconstruction(self):
    361303        """
    362304        If the parent field is a prime field, uses rational reconstruction
  • sage/rings/finite_field_givaro.pyx

    diff -r 7fdc5c05c220 -r be9e35864adf sage/rings/finite_field_givaro.pyx
    a b  
    13881388        else:
    13891389            raise ValueError, "must be a perfect square."
    13901390
    1391     def nth_root(self, int n, extend = False, all = False):
    1392         r"""
    1393         Returns an nth root of self.
    1394 
    1395         INPUT:
    1396             n -- integer >= 1 (must fit in C int type)
    1397             extend -- bool (default: True); if True, return an nth
    1398                  root in an extension ring, if necessary. Otherwise,
    1399                  raise a ValueError if the root is not in the base
    1400                  ring.  Warning: this option is not implemented!
    1401             all -- bool (default: False); if True, return all nth
    1402                  roots of self, instead of just one.
    1403 
    1404         OUTPUT:
    1405             If self has an nth root, returns one (if all == False) or a list of
    1406             all of them (if all == True).  Otherwise, raises a ValueError (if
    1407             extend = False) or a NotImplementedError (if extend = True).
    1408 
    1409         AUTHOR:
    1410             -- David Roe (2007-10-3)
    1411 
    1412         WARNING:
    1413             The 'extend' option is not implemented (yet).
    1414 
    1415         EXAMPLES:
    1416             sage: k.<a> = GF(29^2)
    1417             sage: b = a^2 + 5*a + 1
    1418             sage: b.nth_root(11)
    1419             3*a + 20
    1420             sage: b.nth_root(5)
    1421             Traceback (most recent call last):
    1422             ...
    1423             ValueError: no nth root
    1424             sage: b.nth_root(5, all = True)
    1425             []
    1426             sage: b.nth_root(3, all = True)
    1427             [14*a + 18, 10*a + 13, 5*a + 27]
    1428         """
    1429         if extend:
    1430             raise NotImplementedError
    1431         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    1432         R = PolynomialRing(self.parent(), "x")
    1433         f = R([-self] + [self.parent()(Integer(0))] * (n - 1) + [self.parent()(1)])
    1434         L = f.roots()
    1435         if all:
    1436             return [x[0] for x in L]
    1437         else:
    1438             if len(L) == 0:
    1439                 raise ValueError, "no nth root"
    1440             else:
    1441                 return L[0][0]
    1442        
    14431391    cpdef ModuleElement _add_(self, ModuleElement right):
    14441392        """
    14451393        Add two elements.
  • sage/rings/integer_mod.pyx

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

    diff -r 7fdc5c05c220 -r be9e35864adf sage/structure/element.pyx
    a b  
    26102610
    26112611        return order
    26122612
    2613     def nth_root(self, int n, extend = False, all = False):
     2613    def nth_root(self, n, extend = False, all = False, algorithm=None):
    26142614        r"""
    26152615        Returns an nth root of self.
    26162616
    26172617        INPUT:
    26182618
    2619         - ``n`` - integer >= 1 (must fit in C int type)
     2619        - ``n`` - integer >= 1
    26202620
    26212621        - ``extend`` - bool (default: True); if True, return an nth
    26222622          root in an extension ring, if necessary. Otherwise, raise a
     
    26262626        - ``all`` - bool (default: False); if True, return all nth
    26272627          roots of self, instead of just one.
    26282628
     2629        - ``algorithm`` - string (default: None); 'Johnston' is the only
     2630          currently supported option.
     2631
    26292632        OUTPUT:
    26302633
    26312634        If self has an nth root, returns one (if all == False) or a
     
    26372640       
    26382641           The 'extend' option is not implemented (yet).
    26392642
     2643        EXAMPLES::
     2644
     2645            sage: K = GF(31)
     2646            sage: a = K(22)
     2647            sage: K(22).nth_root(7)
     2648            13
     2649            sage: K(25).nth_root(5)
     2650            5
     2651            sage: K(23).nth_root(3)
     2652            29
     2653            sage: K.<a> = GF(625)
     2654            sage: (3*a^2+a+1).nth_root(13)**13
     2655            3*a^2 + a + 1
     2656
     2657            sage: k.<a> = GF(29^2)
     2658            sage: b = a^2 + 5*a + 1
     2659            sage: b.nth_root(11)
     2660            3*a + 20
     2661            sage: b.nth_root(5)
     2662            Traceback (most recent call last):
     2663            ...
     2664            ValueError: no nth root
     2665            sage: b.nth_root(5, all = True)
     2666            []
     2667            sage: b.nth_root(3, all = True)
     2668            [14*a + 18, 10*a + 13, 5*a + 27]
     2669
     2670            sage: k.<a> = GF(29^5)
     2671            sage: b = a^2 + 5*a + 1
     2672            sage: b.nth_root(5)
     2673            19*a^4 + 2*a^3 + 2*a^2 + 16*a + 3       
     2674            sage: b.nth_root(7)
     2675            Traceback (most recent call last):
     2676            ...
     2677            ValueError: no nth root
     2678            sage: b.nth_root(4, all=True)
     2679            []
     2680
     2681        TESTS::
     2682
     2683            sage: for p in [2,3,5,7,11]: # long
     2684            ...       for n in [2,5,10]:
     2685            ...           q = p^n
     2686            ...           K.<a> = GF(q)
     2687            ...           for r in (q-1).divisors():
     2688            ...               if r == 1: continue
     2689            ...               x = K.random_element()
     2690            ...               y = x^r
     2691            ...               if y.nth_root(r)**r != y: raise RuntimeError
     2692            ...               if (y^41).nth_root(41*r)**(41*r) != y^41: raise RuntimeError
     2693            ...               if (y^307).nth_root(307*r)**(307*r) != y^307: raise RuntimeError
     2694
     2695        ALGORITHMS:
     2696
     2697        - The default is currently an algorithm described in the following paper:
     2698
     2699        Johnston, Anna M. A generalized qth root algorithm. Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms. Baltimore, 1999: pp 929-930.
     2700
    26402701        AUTHOR:
    26412702
    2642         - David Roe (2007-10-3)
    2643 
     2703        - David Roe (2010-1-13)
    26442704        """
    2645         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    2646         from sage.rings.integer import Integer
     2705        if self.is_zero():
     2706            if n < 0:
     2707                raise ZeroDivisionError
     2708            if all: return [self]
     2709            else: return self
     2710        if algorithm is None:
     2711            algorithm = 'Johnston'
    26472712        if extend:
    26482713            raise NotImplementedError
    2649         R = PolynomialRing(self.parent(), "x")
    2650         f = R([-self] + [self.parent()(Integer(0))] * (n - 1) + [self.parent()(1)])
    2651         L = f.roots()
    2652         if all:
    2653             return [x[0] for x in L]
     2714        from sage.rings.integer import Integer
     2715        n = Integer(n)
     2716        K = self.parent()
     2717        q = K.order()
     2718        gcd, alpha, beta = n.xgcd(q-1)
     2719        if self.is_one():
     2720            if gcd == 1:
     2721                if all: return [self]
     2722                else: return self
     2723            else:
     2724                g = K.multiplicative_generator()
     2725                q1overn = (q-1)//gcd
     2726                nthroot = g**q1overn
     2727                if all:
     2728                    return [nthroot**a for a in range(gcd)]
     2729                else:
     2730                    return nthroot
     2731        if gcd == 1:
     2732            if all: return [self**alpha]
     2733            else: return self**alpha
     2734        m = n // gcd
     2735        n = gcd
     2736        gcd, alpha, beta = m.xgcd(q-1)
     2737        if gcd != 1:
     2738            if all:
     2739                return []
     2740            else:
     2741                raise ValueError, "no nth root"
     2742        self = self**alpha
     2743        F = n.factor()
     2744        from sage.groups.generic import discrete_log
     2745        if algorithm == 'Johnston':
     2746            q1overn = (q-1)//n
     2747            if self**q1overn != 1:
     2748                if all:
     2749                    return []
     2750                else:
     2751                    raise ValueError, "no nth root"
     2752            g = K.multiplicative_generator()
     2753            for r, v in F:
     2754                k, h = (q-1).val_unit(r)
     2755                z = h * (-h).inverse_mod(r**v)
     2756                x = (1 + z) // r**v
     2757                if k == 1:
     2758                    self = self**x
     2759                else:
     2760                    t = discrete_log(self**h, g**(r**v*h), r**(k-v), operation='*')
     2761                    self = self**x * g**(-z*t)
     2762            if all:
     2763                nthroot = g**q1overn
     2764                L = [self]
     2765                for i in range(1,n):
     2766                    self *= nthroot
     2767                    L.append(self)
     2768                return L
     2769            else:
     2770                return self
    26542771        else:
    2655             if len(L) == 0:
    2656                 raise ValueError, "no nth root"
    2657             else:
    2658                 return L[0][0]
     2772            raise ValueError, "unknown algorithm"
    26592773
    26602774    def additive_order(self):
    26612775        """