Ticket #8335: 8335_pseudo_conway.patch

File 8335_pseudo_conway.patch, 50.0 KB (added by roed, 10 years ago)

Apply first

  • sage/rings/finite_rings/constructor.py

    # HG changeset patch
    # User David Roe <roed@math.harvard.edu>
    # Date 1266946097 18000
    # Node ID e05450de7d57127f25e28c1765fdb83b73e2384c
    # Parent  7484c56ea8b4bf0e9f2126f8adaca4c4d795490a
    8335: adds functions for finding pseudo-Conway polynomials, facility for creating finite fields without specifying a variable name, any_root and squarefree_decomposition to polynomials over finite fields.
    
    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/finite_rings/constructor.py
    a b  
    171171import sage.databases.conway
    172172
    173173from sage.structure.factory import UniqueFactory
     174from sage.structure.sage_object import SageObject
     175import weakref
    174176
    175177class FiniteFieldFactory(UniqueFactory):
    176178    """
     
    183185   
    184186    -  ``order`` - int
    185187   
    186     -  ``name`` - string; must be specified if not a prime
    187        field
     188    -  ``name`` - string; must be specified unless
     189       - ``order`` is prime
     190       - ``modulus`` is 'conway', in which case a name is generated
     191          from the degree and the ``prefix`` string.
    188192   
    189193    -  ``modulus`` - (optional) either a defining polynomial for the
    190194       field, i.e., generator of the field will be a root of this
    191195       polynomial; or a string:
    192196
    193           - 'conway': force the use of a Conway polynomial, will
    194             raise a RuntimeError if none is found in the database;
    195           - 'random': use a random irreducible polynomial;
    196           - 'default': a Conway polynomial is used if found. Otherwise
     197           - 'conway': force the use of a pseudo-Conway polynomial.
     198            If a conway polynomial is found in the database, that is used.
     199            Otherwise, a polynomial is used with the same algebraic properties
     200            but without the lexicographic constraint implying uniqueness.
     201           - 'random': use a random irreducible polynomial;
     202           - 'default': a Conway polynomial is used if in the database. Otherwise
    197203            a sparse polynomial is used for binary fields and a
    198204            random polynomial is used for other characteristics.
     205           - 'first_lexicographic': uses the first irreducible polynomial.
    199206
    200207       Other options might be available depending on the
    201208       implementation.
     
    205212   
    206213    -  ``check_irreducible`` - verify that the polynomial
    207214       modulus is irreducible
     215
     216    - ``prefix`` (default: 'z') - a string used to generate names
     217      for subfields and pushouts of this field (mainly for use with
     218      pseudo-Conway polynomials).
    208219   
    209220    -  ``args`` - additional parameters passed to finite
    210221       field implementations
     
    224235        Finite Field in a of size 3^2
    225236        sage: charpoly(a, 'y')
    226237        y^2 + 2*y + 2
     238        sage: L = GF(27); L
     239        Finite Field in z3 of size 3^3
    227240   
    228241    ::
    229242   
     
    309322        sage: a
    310323        2
    311324    """
    312     def create_key_and_extra_args(self, order, name=None, modulus=None, names=None, impl=None, **kwds):
     325    def create_key_and_extra_args(self, order, name=None, modulus=None, names=None, prefix='z', impl=None, **kwds):
    313326        """
    314327        EXAMPLES::
    315328       
    316329            sage: GF.create_key_and_extra_args(9, 'a')
    317             ((9, ('a',), 'conway', None, '{}'), {})
     330            ((9, ('a',), 'conwayz', None, '{}'), {})
    318331            sage: GF.create_key_and_extra_args(9, 'a', foo='value')
    319             ((9, ('a',), 'conway', None, "{'foo': 'value'}"), {'foo': 'value'})
     332            ((9, ('a',), 'conwayz', None, "{'foo': 'value'}"), {'foo': 'value'})
    320333        """
    321334        order = int(order)
    322335        if order == 1:
     
    329342            modulus = None
    330343        else:
    331344            if not names is None: name = names
    332             name = normalize_names(1,name)
    333345
    334346            p,n = arith.factor(order)[0]
    335347
    336348            if modulus is None or modulus == "default":
    337                 if exists_conway_polynomial(p,n):
     349                if exists_conway_polynomial(p,n) or name is None:
    338350                    modulus = "conway"
    339351                else:
    340352                    if p==2:
    341353                        modulus = "minimal_weight"
    342354                    else:
    343355                        modulus = "random"
     356            elif p != 2 and modulus == 'first_lexicographic':
     357                R = FiniteField(p)['x']
     358                x = R.gen()
     359                for f in R.polynomials(of_degree=n):
     360                    if f.is_irreducible():
     361                        f = f.monic() # probably not necessary...
     362                        break
     363                modulus = f
     364               
     365            if modulus == "conway":
     366                modulus += prefix
     367                if name is None:
     368                    name = prefix + str(n)
    344369            elif modulus == "random":
    345370                modulus += str(random.randint(0, 1<<128))
     371            name = normalize_names(1,name)
    346372
    347373            if isinstance(modulus, (list, tuple)):
    348374                modulus = FiniteField(p)['x'](modulus)
     
    395421                if not modulus.is_irreducible():
    396422                    raise ValueError, "finite field modulus must be irreducible but it is not"
    397423            if name is None:
    398                 raise TypeError, "you must specify the generator name"
    399             if order < zech_log_bound
     424                raise TypeError, "you must specify the generator name or use Conway polynomials"
     425            if order < zech_log_bound and (impl is None or impl =='givaro')
    400426                # DO *NOT* use for prime subfield, since that would lead to
    401427                # a circular reference in the call to ParentWithGens in the
    402428                # __init__ method.
     
    416442        EXAMPLES::
    417443       
    418444            sage: key, extra = GF.create_key_and_extra_args(9, 'a'); key
    419             (9, ('a',), 'conway', None, '{}')
     445            (9, ('a',), 'conwayz', None, '{}')
    420446            sage: K = GF.create_object(0, key); K
    421447            Finite Field in a of size 3^2
    422448            sage: GF.other_keys(key, K)
     
    546572    """
    547573    return sage.databases.conway.ConwayPolynomials().has_polynomial(p,n)
    548574
     575pseudo_conway_poly = {}
     576
     577class PseudoConwayPolyTree(SageObject):
     578    """
     579    An object holding references to pseudo-Conway polynomials for divisors of the given degree, so that they aren't garbage collected.
     580    """
     581    def __init__(self, p, n, nodes_dict, f):
     582        """
     583        INPUT::
     584
     585        - p -- The prime for which this defines an extension of GF(p).
     586        - n -- The degree of the extension.
     587        - nodes_dict -- (dict or bool).  A dictionary of PseudoConwayPolyTrees, indexed by prime divisors of n.  The entry for q corresponds to the pseudo-Conway extension of degree n//q.
     588                        OR a boolean.  If True, then this polynomial is actually in the Conway polynomials database, and no references to other PCPTs are stored.  If False, then n is prime and no references are stored (since there is no compatiblity condition).
     589        - f -- The polynomial defining this extension.
     590
     591        EXAMPLES::
     592
     593            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     594            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     595            sage: PCPT.get_pseudo_conway_poly(3)
     596            x^3 + x + 1
     597        """
     598        self.p = p
     599        self.n = n
     600        self.nodes_dict = nodes_dict
     601        self.f = f
     602       
     603    def get_pseudo_conway_poly(self, m):
     604        """
     605        Returns the pseudo-Conway polynomial associated to a divisor of the degree of this tree.
     606
     607        EXAMPLES::
     608       
     609            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     610            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     611            sage: PCPT.get_pseudo_conway_poly(3)
     612            x^3 + x + 1
     613            sage: PCPT.get_pseudo_conway_poly(4)
     614            Traceback (most recent call last):
     615            ...
     616            RuntimeError: 4 does not divide 6
     617        """
     618        if m == self.n:
     619            return self.f
     620        if not m.divides(self.n):
     621            raise RuntimeError, "%s does not divide %s"%(m, self.n)
     622        try:
     623            return pseudo_conway_poly[self.p][m]().f
     624        except (KeyError, AttributeError):
     625            return conway_polynomial(self.p, m)
     626           
     627    def check_consistency(self):
     628        """
     629        Checks that the pseudo-Conway polynomials dividing the degree of this tree satisfy the required compatibility conditions.
     630
     631        EXAMPLES::
     632
     633            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     634            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     635            sage: PCPT.check_consistency()
     636            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 60, False) # long
     637            sage: PCPT.check_consistency() # long
     638        """
     639        p, n = self.p, self.n
     640        K = GF(p**n, modulus = self.f, names='a')
     641        a = K.gen()
     642        for m in n.divisors():
     643            assert (a**((p**n-1)//(p**m-1))).minimal_polynomial() == self.get_pseudo_conway_poly(m)
     644
     645def find_pseudo_conway_polynomial_tree(p, n, use_database=True):
     646    """
     647    Returns an object holding references to a set of consistent pseudo-Conway polynomials for degrees dividing n.
     648
     649    A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
     650   
     651    - `f` is irreducible.
     652    - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
     653    - In this same quotient field, the minimal polynomial of `x^{\frac{p^n-1}{p^m-1}}` is the Conway polynomial of degree `m` for every divisor `m` of `n`.
     654    - `f` is lexicographically least among all such polynomials, under a certain ordering.
     655
     656    The final condition is needed only in order to make the Conway polynomial unique.  We define a pseudo-Conway polynomial to be one satisfying the first three conditions.
     657
     658    INPUT:
     659
     660    - `p` -- a prime.
     661    - `n` -- an integer greater than 1.
     662    - `use_database` -- whether to use the Conway polynomials database if the Conway polynomial for `p, n` exists within it (versus computing a pseudo-Conway polynomial)
     663
     664    EXAMPLES::
     665
     666        sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     667        sage: PCPT = find_pseudo_conway_polynomial_tree(2, 12, False)
     668        sage: PCPT.f
     669        x^12 + x^8 + x^7 + x^6 + x^4 + x^3 + 1
     670    """
     671    if not pseudo_conway_poly.has_key(p):
     672        pseudo_conway_poly[p] = {}
     673    pdict = pseudo_conway_poly[p]
     674    n = sage.rings.integer.Integer(n)
     675    if pdict.has_key(n):
     676        pcp = pdict[n]()
     677        if pcp is not None:
     678            return pcp
     679    if use_database and exists_conway_polynomial(p, n):
     680        ans = PseudoConwayPolyTree(p, n, True, GF(p)['x'](conway_polynomial(p, n)))
     681        pdict[n] = weakref.ref(ans)
     682        return ans
     683    if n == 1 or n.is_prime():
     684        ans = PseudoConwayPolyTree(p, n, False, compute_pseudo_conway_polynomial(p, n, None))
     685    nodes = {}
     686    for q in n.prime_factors():
     687        nodes[q] = find_pseudo_conway_polynomial_tree(p, n//q, use_database)
     688    ans = PseudoConwayPolyTree(p, n, nodes, compute_pseudo_conway_polynomial(p, n, nodes))
     689    pdict[n] = weakref.ref(ans)
     690    return ans
     691
     692def _find_pow_of_frobenius(p, x, y, field_degree):
     693    """
     694    Finds the power of Frobenius which yields `x` when applied to `y`.
     695
     696    INPUT:
     697
     698    - `p` -- a prime.
     699    - ``field_degree`` -- a positive integer
     700    - `x` -- an element of a field `K` of characeristic p so that the multiplicative order of `x` is ``p^field_degree``.
     701    - `y` -- an element of `K` with the same minimal polynomial.
     702   
     703    OUTPUT:
     704   
     705    An element `i` of ``Integers(field_degree)`` so that `x = y^{p^i}`
     706
     707    EXAMPLES::
     708
     709        sage: from sage.rings.finite_rings.constructor import _find_pow_of_frobenius
     710        sage: K.<a> = GF(3^14)
     711        sage: x = K.multiplicative_generator()
     712        sage: y = x^27
     713        sage: _find_pow_of_frobenius(3, x, y, 14)
     714        11
     715    """
     716    assert x.minimal_polynomial() == y.minimal_polynomial()
     717    assert x.multiplicative_order() == p**field_degree - 1
     718    assert y.multiplicative_order() == p**field_degree - 1
     719    from integer_mod import mod
     720    for i in range(field_degree):
     721        if x == y: break
     722        y = y**p
     723    else:
     724        raise RuntimeError, "No appropriate power of Frobenius found"
     725    return mod(i, field_degree)
     726
     727def _crt_non_coprime(running, a):
     728    """
     729    Extends the ``crt`` method of IntegerMod to the case of non-relatively prime modulus.
     730
     731    EXAMPLES::
     732
     733        sage: from sage.rings.finite_rings.constructor import _crt_non_coprime
     734        sage: a = _crt_non_coprime(mod(14, 18), mod(20,30)); a
     735        50
     736        sage: a.modulus()
     737        90
     738        sage: _crt_non_coprime(mod(13, 18), mod(20,30))
     739        Traceback (most recent call last):
     740        ...
     741        AssertionError
     742    """
     743    g = running.modulus().gcd(a.modulus())
     744    if g == 1:
     745        return running.crt(a)
     746    else:
     747        assert running % g == a % g
     748        running_modulus = running.modulus()
     749        a_modulus = a.modulus()
     750        for qq in g.prime_divisors():
     751            a_val_unit = a_modulus.val_unit(qq)
     752            running_val_unit = running_modulus.val_unit(qq)
     753            if a_val_unit[0] > running_val_unit[0]:
     754                running_modulus = running_val_unit[1]
     755            else:
     756                a_modulus = a_val_unit[1]
     757        return (running % running_modulus).crt(a % a_modulus)
     758
     759def _frobenius_shift(K, generators, check_only=False):
     760    """
     761    Given a field `K` of degree `n` over ``GF(p)`` and a dictionary holding, for each divisor `q` of `n`, an element with minimal polynomial a pseudo-Conway polynomial of degree `n/q`, modifies these generators into a compatibile system.
     762
     763    Such a system of generators is said to be compatible if for each pair of prime divisors `q_1` and `q_2` and each common divisor `m` of `n/q_1` and `n/q_2`, the equality
     764
     765    ``generators[q1]^((p^(n/q1)-1)/(p^m-1)) == generators[q2]^((p^(n/q2)-1)/(p^m-1))``
     766
     767    holds.
     768
     769    INPUT:
     770
     771    - K -- a finite field of degree n   
     772    - generators -- a dictionary, indexed by prime divisors `q` of `n`, whose entries are elements of `K` satisfying the `n/q` pseudo-Conway polynomial.
     773    - check_only -- if True, just asserts that the generators given form a compatible system.
     774
     775    EXAMPLES::
     776
     777        sage: R.<x> = GF(2)[]
     778        sage: f30 = x^30 + x^28 + x^27 + x^25 + x^24 + x^20 + x^19 + x^18 + x^16 + x^15 + x^12 + x^10 + x^7 + x^2 + 1
     779        sage: f20 = x^20 + x^19 + x^15 + x^13 + x^12 + x^11 + x^9 + x^8 + x^7 + x^4 + x^2 + x + 1
     780        sage: f12 = x^12 + x^10 + x^9 + x^8 + x^4 + x^2 + 1
     781        sage: K.<a> = GF(2^60, modulus='first_lexicographic')
     782        sage: x30 = f30.any_root(K)
     783        sage: x20 = f20.any_root(K)
     784        sage: x12 = f12.any_root(K)
     785        sage: generators = {2: x30, 3: x20, 5: x12}
     786        sage: from sage.rings.finite_rings.constructor import _frobenius_shift, _find_pow_of_frobenius
     787        sage: _frobenius_shift(K, generators)
     788        sage: _find_pow_of_frobenius(2, x30, generators[2], 30)
     789        0
     790        sage: _find_pow_of_frobenius(2, x20, generators[3], 20)
     791        13
     792        sage: _find_pow_of_frobenius(2, x12, generators[5], 12)
     793        8
     794    """
     795    p = K.characteristic()
     796    n = K.degree()
     797    if len(generators) == 1:
     798        return generators
     799    compatible = {}
     800    from integer_mod import mod
     801    for m in n.divisors():
     802        compatible[m] = {}
     803    for q, x in generators.iteritems():
     804        for m in (n//q).divisors():
     805            compatible[m][q] = x**((p**(n//q)-1)//(p**m-1))
     806    if check_only:
     807        for m in n.divisors():
     808            try:
     809                q, x = compatible[m].popitem()
     810            except KeyError:
     811                break
     812            for qq, xx in compatible[m].iteritems():
     813                assert x == xx
     814        return
     815    crt = {}
     816    qlist = sorted(generators.keys())
     817    for j in range(1, len(qlist)):
     818        for i in range(j):
     819            crt[(i, j)] = []
     820    for m in n.divisors():
     821        mqlist = sorted(compatible[m].keys())
     822        for k in range(1,len(mqlist)):
     823            j = qlist.index(mqlist[k])
     824            i = qlist.index(mqlist[k-1])
     825            crt[(i,j)].append(_find_pow_of_frobenius(p, compatible[m][qlist[j]], compatible[m][qlist[i]], m))
     826    from integer_mod import mod
     827    pairs = crt.keys()
     828    for i, j in pairs:
     829        L = crt[(i,j)]
     830        #print qlist[i], qlist[j]
     831        #print " ".join(["(%s, %s)"%(mm, mm.modulus()) for mm in L])
     832        running = mod(0,1)
     833        for a in L:
     834            running = _crt_non_coprime(running, a)
     835        #print "(%s, %s)"%(running, running.modulus())
     836        crt[(i,j)] = [(mod(running, q**(running.modulus().valuation(q))), running.modulus().valuation(q)) for q in qlist]
     837        crt[(j,i)] = [(-a, level) for a, level in crt[(i,j)]]
     838    # Let x_j be the power of Frobenius we apply to generators[qlist[j]], for 0 < j < len(qlist)
     839    # We have some direct conditions on the x_j: x_j reduces to each entry in crt[(0,j)].
     840    # But we also have the equations x_j - x_i reduces to each entry in crt[(i,j)].
     841    # We solve for x_j one prime at a time.  For each prime, we have an equations of the form
     842    # x_j - x_i = c_ij.  The modulus of the currently known value of x_j, x_i and c_ij will all be powers
     843    # (possibly 0, possibly different) of the same prime.
     844
     845    # We can set x_0=0 everywhere, can get an initial setting of x_j from the c_0j.
     846    # We go through prime by prime.
     847    import bisect
     848    frob_powers=[mod(0,1) for q in qlist]
     849    def find_leveller(qindex, level, x, xleveled, searched, i):
     850        searched[i] = True
     851        crt_possibles = []
     852        for j in range(1,len(qlist)):
     853            if i==j: continue
     854            if crt[(i,j)][qindex][1] >= level:
     855                if xleveled[j]:
     856                    return [j]
     857                elif not searched.has_key(j):
     858                    crt_possibles.append(j)
     859        for j in crt_possibles:
     860            path = find_leveller(qindex, level, x, xleveled, searched, j)
     861            if path is not None:
     862                path.append(j)
     863                return path
     864        return None
     865    def propogate_levelling(qindex, level, x, xleveled, i):
     866        for j in range(1, len(qlist)):
     867            if i==j: continue
     868            if not xleveled[j] and crt[(i,j)][qindex][1] >= level:
     869                newxj = x[i][0] + crt[(i,j)][qindex][0]
     870                x[j] = (newxj, min(x[i][1], crt[(i,j)][qindex][1]))
     871                xleveled[j] = True
     872                propogate_levelling(qindex, level, x, xleveled, j)
     873               
     874    for qindex in range(len(qlist)):
     875        q = qlist[qindex]
     876        x = [0] + [crt[(0,j)][qindex] for j in range(1,len(qlist))] # we include the initial 0 to match up our indexing with crt.               
     877        # We first check that our equations are consistent and determine which powers of q occur as moduli.
     878        levels = []
     879        for j in range(2, len(qlist)):
     880            for i in range(j):
     881                # we need crt[(0,j)] = crt[(0,i)] + crt[(i,j)]
     882                if i != 0: assert x[j][0] == x[i][0] + crt[(i,j)][qindex][0]
     883                level = crt[(i,j)][qindex][1]
     884                if level > 0:
     885                    ins = bisect.bisect_left(levels,level)
     886                    if ins == len(levels):
     887                        levels.append(level)
     888                    elif levels[ins] != level:
     889                        levels.insert(ins, level)
     890        for level in levels:
     891            xleveled = [0] + [x[i][1] >= level for i in range(1,len(qlist))]
     892            while True:
     893                try:
     894                    i = xleveled.index(False, 1)
     895                    searched = {}
     896                    levelling_path = find_leveller(qindex, level, x, xleveled, searched, i)
     897                    if levelling_path is None:
     898                        # Any lift will work, since there are no constraints
     899                        x[i] = (mod(x[i][0].lift(), q**level), level)
     900                        xleveled[i] = True
     901                        propogate_levelling(qindex, level, x, xleveled, i)
     902                    else:
     903                        levelling_path.append(i)
     904                        for m in range(1,len(path)):
     905                            # this point on the path may have already been leveled in a previous propogation
     906                            if not xleveled[path[m]]:
     907                                newx = x[path[m-1]][0] + crt[(path[m-1],path[m])][qindex][0]
     908                                x[path[m]] = (newx, min(x[path[m-1]][1], crt[(path[m-1],path[m])][qindex][1]))
     909                                xleveled[path[m]] = True
     910                                propogate_levelling(qindex, level, x, xleveled, path[m])
     911                except ValueError:
     912                    break
     913        for j in range(1,len(qlist)):
     914            frob_powers[j] = frob_powers[j].crt(x[j][0])
     915    for j in range(1, len(qlist)):
     916        #print "frob_power[%s] = mod(%s, %s)"%(qlist[j], frob_powers[j], frob_powers[j].modulus())
     917        generators[qlist[j]] = generators[qlist[j]]**(p**(-frob_powers[j]).lift())
     918    _frobenius_shift(K, generators, check_only=True)
     919
     920def compute_pseudo_conway_polynomial(p, n, nodes):
     921    """
     922    Computes a pseudo-Conway polynomial over the prime `p` of degree `n`.
     923
     924    A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
     925   
     926    - `f` is irreducible.
     927    - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
     928    - In this same quotient field, the minimal polynomial of `x^{\frac{p^n-1}{p^m-1}}` is the Conway polynomial of degree `m` for every divisor `m` of `n`.
     929    - `f` is lexicographically least among all such polynomials, under a certain ordering.
     930
     931    The final condition is needed only in order to make the Conway polynomial unique.  We define a pseudo-Conway polynomial to be one satisfying the first three conditions.
     932
     933    INPUT:
     934   
     935    - `p` -- a prime
     936    - `n` -- a positive integer
     937    - ``nodes`` -- None (in the case that `n` is prime) or a dictionary of PseudoConwayPolyTree objects, indexed by prime divisors `q` of `n` with entries corresponding to pseudo-Conway polynomials of degree `n/q`.
     938
     939    OUTPUT:
     940
     941    A pseudo-Conway polynomial over the prime `p` of degree `n`.
     942
     943    ALGORITHM:
     944   
     945    Uses an algorithm described in [HL99]_, modified to find pseudo-Conway polynomials rather than Conway polynomials.  The major modification was the addition of the function _frobenius_shift.
     946
     947    REFERENCE:
     948
     949    .. [HL99] Heath L. and Loehr N. (1999).
     950       New algorithms for generating Conway polynomials over finite fields.
     951       Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms, pp. 429-437.
     952
     953    EXAMPLES::
     954
     955        sage: from sage.rings.finite_rings.constructor import compute_pseudo_conway_polynomial, find_pseudo_conway_polynomial_tree
     956        sage: compute_pseudo_conway_polynomial(next_prime(10000),11,None) # long because the arithmetic in FiniteField_ext_pari fields is slow.
     957        x^11 + x + 7
     958        sage: PCPT30 = find_pseudo_conway_polynomial_tree(2, 30, False)
     959        sage: PCPT20 = find_pseudo_conway_polynomial_tree(2, 20, False)
     960        sage: PCPT12 = find_pseudo_conway_polynomial_tree(2, 12, False)
     961        sage: compute_pseudo_conway_polynomial(2, 60, {2:PCPT30,3:PCPT20,5:PCPT12})
     962        x^60 + x^59 + x^56 + x^53 + x^49 + x^45 + x^43 + x^42 + x^41 + x^40 + x^39 + x^37 + x^36 + x^34 + x^31 + x^30 + x^27 + x^26 + x^25 + x^23 + x^21 + x^19 + x^17 + x^15 + x^14 + x^11 + x^9 + x^8 + x^7 + x^6 + x^5 + x^2 + 1
     963    """
     964    k = GF(p)
     965    from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     966    R = PolynomialRing(k, names='x')
     967    if n == 1:
     968        return R([-k.multiplicative_generator(), 1])
     969    if nodes is None:
     970        # n is prime.  We require only that the chosen generator (x) be a multiplicative generator for the units.
     971        if p in [2,3,5,7,11]: # here we have Cunningham data
     972            from sage.rings.factorint import factor_cunningham
     973            F = [(p**n-1)//a for a, e in factor_cunningham(n)]
     974        else:
     975            F = (p**n-1).prime_factors() # could be VERY expensive
     976        x = R.gen()
     977        from sage.rings.arith import power_mod
     978        for f in R.polynomials(of_degree=n):
     979            # By using this iterator our chosen polynomial will be fairly sparse.
     980            if f.is_irreducible():
     981                if not f.is_monic():
     982                    f = f.monic()
     983                for m in F:
     984                    if power_mod(x, m, f) == 1:
     985                        break
     986                else:
     987                    return f
     988    # Our strategy is as follows:
     989    # Work in an arbitrary field K of order p**n.
     990    # We're looking for a multiplicative generator of K whose
     991    K = GF(p**n, modulus="first_lexicographic", names='a')
     992    r = p**n - 1
     993    xi = {}
     994    for q in nodes.iterkeys():
     995        #print nodes[q].f.degree(), n//q
     996        xi[q] = nodes[q].f.any_root(K, -n//q, True)
     997    _frobenius_shift(K, xi)
     998           
     999    q, x = xi.popitem()
     1000    v = p**(n//q) - 1
     1001    for q, xitem in xi.iteritems():
     1002        fi = p**(n//q) - 1
     1003        g, alpha, beta = v.xgcd(fi)
     1004        x = x**beta * xitem**alpha
     1005        v = v.lcm(fi)
     1006    g = r // v
     1007    #print g, n, g.divides(p**n-1)
     1008    z = x.nth_root(g, cunningham=(p in [2,3,5,7,11]))
     1009    groot1 = K.multiplicative_generator()**v
     1010    while z.multiplicative_order() != r:
     1011        z *= groot1
     1012    #print z
     1013    return z.minimal_polynomial()
    5491014
    5501015zech_log_bound = 2**16
  • sage/rings/finite_rings/element_ext_pari.py

    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/finite_rings/element_ext_pari.py
    a b  
    565565       
    566566        EXAMPLES::
    567567       
    568             sage: from sage.rings.finite_rings.finite_field_ext_pari import FiniteField_ext_pari
    569568            sage: k = GF(next_prime(10**10))
    570569            sage: a = k(17)/k(19)
    571570            sage: b = a.lift(); b
  • sage/rings/finite_rings/element_ntl_gf2e.pxd

    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/finite_rings/element_ntl_gf2e.pxd
    a b  
    88
    99cdef class FiniteField_ntl_gf2e(FiniteField):
    1010    cdef GF2EContext_c F
     11    cdef public object _PCPT
     12    cdef public object _prefix
    1113    cdef object _polynomial
    1214    cdef object _polynomial_ring
    1315    cdef object _is_conway
  • sage/rings/finite_rings/element_ntl_gf2e.pyx

    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/finite_rings/element_ntl_gf2e.pyx
    a b  
    145145            names   -- variable used for poly_repr (default: 'a')
    146146            modulus -- you may provide a polynomial to use for reduction or
    147147                     a string:
    148                      'conway': force the use of a Conway polynomial, will
    149                      raise a RuntimeError if none is found in the database;
     148                     'conway': force the use of a Conway polynomial; if
     149                     none is found in the database use a generated pseudo-Conway
     150                     polynomial instead;
    150151                     'minimal_weight': use a minimal weight polynomial, should
    151152                     result in faster arithmetic;
     153                     'first_lexicographic': uses the first irreducible polynomial;
    152154                     'random': use a random irreducible polynomial.
    153                      'default':a Conway polynomial is used if found. Otherwise
    154                      a sparse polynomial is used.
    155155            repr  -- controls the way elements are printed to the user:
    156156                     (default: 'poly')
    157157                     'poly': polynomial representation
     
    179179            sage: k.<a> = GF(2^211, modulus='conway')
    180180            sage: k.modulus()
    181181            x^211 + x^9 + x^6 + x^5 + x^3 + x + 1
    182             sage: k.<a> = GF(2^411, modulus='conway')
    183             Traceback (most recent call last):
    184             ...
    185             RuntimeError: requested conway polynomial not in database.
     182            sage: k.<a> = GF(2^23, modulus='conway')
     183            sage: a.multiplicative_order() == k.order()-1
     184            True
    186185
    187186        TESTS::
    188187
     
    226225
    227226        unknown_modulus = True
    228227
    229         if modulus is None or modulus == 'default':
    230             import sage.rings.finite_rings.constructor
    231             if sage.rings.finite_rings.constructor.exists_conway_polynomial(p, k):
    232                 modulus = "conway"
    233             else:
    234                 modulus = "minimal_weight"
    235 
    236         if isinstance(modulus,str):
     228        if isinstance(modulus, str):
    237229            unknown_modulus = False
    238             if modulus == "conway":
    239                 modulus = conway_polynomial(p, k)
    240                 self._is_conway = True
     230            if modulus.startswith('conway'):
     231                from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     232                self._PCPT = find_pseudo_conway_polynomial_tree(p, k) # we store the tree so that the weakrefs don't fly away.  It also provides a nice interface for creating subfields.
     233                self._prefix = modulus[6:]
     234                if len(self._prefix) == 0: self._prefix = 'z'
     235                modulus = self._PCPT.get_pseudo_conway_poly(k)
     236                self._is_conway = ConwayPolynomials().has_polynomial(p, k)
    241237            elif modulus == "minimal_weight":
    242238                GF2X_BuildSparseIrred(ntl_m, k)
    243239            elif modulus == "first_lexicographic":
     
    248244                GF2X_BuildRandomIrred(ntl_m, ntl_tmp)
    249245            else:
    250246                unknown_modulus = True
     247        else:
     248            unknown_modulus = True
    251249
    252250        if is_Polynomial(modulus):
    253251            modulus = modulus.list()
  • sage/rings/finite_rings/finite_field_ext_pari.py

    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/finite_rings/finite_field_ext_pari.py
    a b  
    157157        We do NOT yet define natural consistent inclusion maps
    158158        between different finite fields.
    159159    """
    160     def __init__(self, q, name, modulus=None):
     160    def __init__(self, q, name, modulus='conway'):
    161161        """
    162162        Create finite field of order q with variable printed as name.
    163163
     
    170170                    otherwise.
    171171        - ``modulus`` -- you may provide a polynomial to use for reduction or
    172172                     a string:
    173                      'conway': force the use of a Conway polynomial, will
    174                      raise a RuntimeError if none is found in the database;
     173                     'conway': force the use of a Conway polynomial; if
     174                     none is found in the database, will use a pseudo-Conway
     175                     polynomial instead;
    175176                     'random': use a random irreducible polynomial.
    176177                     'default': a Conway polynomial is used if found. Otherwise
    177178                     a random polynomial is used.
     
    227228        self.__order = q
    228229        self.__is_field = True
    229230
    230         if modulus is None or modulus == "default":
    231             from constructor import exists_conway_polynomial
    232             if exists_conway_polynomial(self.__char, self.__degree):
    233                 modulus = "conway"
    234             else:
    235                 modulus = "random"
    236 
    237         if isinstance(modulus,str):
    238             if modulus == "conway":
    239                 from constructor import conway_polynomial
    240                 modulus = conway_polynomial(self.__char, self.__degree)
     231        if isinstance(modulus, str):
     232            if modulus.startswith('conway'):
     233                from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     234                self._PCPT = find_pseudo_conway_polynomial_tree(self.__char, self.__degree) # we store the tree so that the weakrefs don't fly away.  It also provides a nice interface for creating subfields.
     235                self._prefix = modulus[6:]
     236                if len(self._prefix) == 0: self._prefix = 'z'
     237                modulus = self._PCPT.get_pseudo_conway_poly(self.__degree)
    241238            elif modulus == "random":
    242239                # The following is fast/deterministic, but has serious problems since
    243240                # it crashes on 64-bit machines, and I can't figure out why:
  • sage/rings/finite_rings/finite_field_givaro.py

    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/finite_rings/finite_field_givaro.py
    a b  
    99from sage.libs.pari.all import pari
    1010
    1111class FiniteField_givaro(FiniteField):
    12     def __init__(self, q, name="a",  modulus=None, repr="poly", cache=False):
     12    def __init__(self, q, name="a",  modulus='conway', repr="poly", cache=False):
    1313        """
    1414        Finite Field. These are implemented using Zech logs and the
    1515        cardinality must be < 2^16. By default conway polynomials are
     
    2020            name  -- variable used for poly_repr (default: 'a')
    2121            modulus -- you may provide a minimal polynomial to use for
    2222                     reduction or 'random' to force a random
    23                      irreducible polynomial. (default: None, a conway
    24                      polynomial is used if found. Otherwise a random
    25                      polynomial is used)
     23                     irreducible polynomial. (default: 'conway', a conway
     24                     polynomial is used if in the database, otherwise
     25                     a pseudo-Conway polynomial is generated)
    2626            repr  -- controls the way elements are printed to the user:
    2727                     (default: 'poly')
    2828                     'log': repr is element.log_repr()
     
    105105        self._kwargs['cache'] = cache
    106106
    107107        self._is_conway = False
    108         if modulus is None or modulus == 'conway':
     108        if isinstance(modulus, str) and modulus.startswith('conway'):
    109109            if k == 1:
    110110                modulus = 'random' # this will use the gfq_factory_pk function.
    111             elif ConwayPolynomials().has_polynomial(p, k):
    112                 from sage.rings.finite_rings.constructor import conway_polynomial
    113                 modulus = conway_polynomial(p, k)
    114                 self._is_conway = True
    115             elif modulus is None:
    116                 modulus = 'random'
    117111            else:
    118                 raise ValueError, "Conway polynomial not found"
     112                from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     113                self._PCPT = find_pseudo_conway_polynomial_tree(p, k) # we store the tree so that the weakrefs don't fly away.  It also provides a nice interface for creating subfields.
     114                self._prefix = modulus[6:]
     115                if len(self._prefix) == 0: self._prefix = 'z'
     116                modulus = self._PCPT.get_pseudo_conway_poly(k)
     117                # The following is only used for _gap_init_, so we check to see if the defining polynomial is actually a conway polynomial.
     118                self._is_conway = ConwayPolynomials().has_polynomial(p, k)
    119119               
    120120        from sage.rings.polynomial.all import is_Polynomial
    121121        if is_Polynomial(modulus):
  • sage/rings/polynomial/polynomial_element.pyx

    diff -r 7484c56ea8b4 -r e05450de7d57 sage/rings/polynomial/polynomial_element.pyx
    a b  
    11931193            sage: f = QQbar['x'](1)
    11941194            sage: f.squarefree_decomposition()
    11951195            1
     1196
     1197        TESTS::
     1198
     1199            sage: K.<a> = GF(3^2)
     1200            sage: R.<x> = K[]
     1201            sage: f = x^243+2*x^81+x^9+1
     1202            sage: f.squarefree_decomposition()
     1203            (x^27 + 2*x^9 + x + 1)^9
     1204            sage: f = x^243+a*x^27+1
     1205            sage: f.squarefree_decomposition()
     1206            (x^9 + (2*a + 1)*x + 1)^27
     1207
     1208            sage: for K in [GF(2^18,'a'), GF(3^2,'a'), GF(47^3,'a')]:
     1209            ...       R.<x> = K[]
     1210            ...       if K.characteristic() < 5: m = 4
     1211            ...       else: m = 1
     1212            ...       for _ in range(m):
     1213            ...           f = (R.random_element(4)^3*R.random_element(m)^(m+1))(x^6)
     1214            ...           F = f.squarefree_decomposition()
     1215            ...           assert F.prod() == f
     1216            ...           for i in range(len(F)):
     1217            ...               assert gcd(F[i][0], F[i][0].derivative()) == 1
     1218            ...               for j in range(len(F)):
     1219            ...                   if i == j: continue
     1220            ...                   assert gcd(F[i][0], F[j][0]) == 1
    11961221        """
    11971222        if self.degree() == 0:
    11981223            return Factorization([], unit=self[0])
    11991224
    1200         # Wikipedia says this works for arbitrary fields of
    1201         # characteristic 0.
    1202 
    1203         if self.base_ring().characteristic() != 0:
    1204             raise NotImplementedError, "Squarefree decomposition not implemented for " + str(self.parent())
    1205 
    1206         f = [self]
    1207         cur = self
    1208         while cur.degree() > 0:
    1209             cur = cur.gcd(cur.derivative())
    1210             f.append(cur)
    1211 
    1212         g = []
    1213         for i in range(len(f) - 1):
    1214             g.append(f[i] // f[i+1])
    1215 
    1216         a = []
    1217         for i in range(len(g) - 1):
    1218             a.append(g[i] // g[i+1])
    1219         if g == []: # do not factor the unit part
    1220             return Factorization([(self,1)])
    1221         a.append(g[-1])
    1222 
     1225        p = self.base_ring().characteristic()
    12231226        factors = []
    1224         unit = f[-1]
    1225         for i in range(len(a)):
    1226             if a[i].degree() > 0:
    1227                 factors.append((a[i], i+1))
    1228             else:
    1229                 unit = unit * a[i].constant_coefficient() ** (i + 1)
    1230 
     1227        if p == 0:
     1228            # Wikipedia says this works for arbitrary fields of
     1229            # characteristic 0.
     1230           
     1231            f = [self]
     1232            cur = self
     1233            while cur.degree() > 0:
     1234                cur = cur.gcd(cur.derivative())
     1235                f.append(cur)
     1236               
     1237            g = []
     1238            for i in range(len(f) - 1):
     1239                g.append(f[i] // f[i+1])
     1240               
     1241            a = []
     1242            for i in range(len(g) - 1):
     1243                a.append(g[i] // g[i+1])
     1244            a.append(g[-1])
     1245
     1246            unit = f[-1]
     1247            for i in range(len(a)):
     1248                if a[i].degree() > 0:
     1249                    factors.append((a[i], i+1))
     1250                else:
     1251                    unit = unit * a[i].constant_coefficient() ** (i + 1)
     1252        else:
     1253            if not self.base_ring().is_field():
     1254                raise NotImplementedError, "Squarefree decomposition not implemented for " + str(self.parent())
     1255            # In characteristic p, we use algorithm 3.4.2 in Cohen's Course in Computational Algebraic Number Theory
     1256           
     1257           
     1258            unit = self.leading_coefficient()
     1259            T0 = self.monic()
     1260            e = 1
     1261            if T0.degree() > 0:
     1262                der = T0.derivative()
     1263                while der.is_zero():
     1264                    T0 = T0.parent()([T0[p*i].pth_root() for i in range(T0.degree()//p + 1)])
     1265                    if T0 == 1: raise RuntimeError
     1266                    der = T0.derivative()
     1267                    e = e*p
     1268                T = T0.gcd(der)
     1269                V = T0 // T
     1270                k = 0
     1271                while T0.degree() > 0:
     1272                    k += 1
     1273                    if p.divides(k):
     1274                        T = T // V
     1275                        k += 1
     1276                    W = V.gcd(T)
     1277                    if W.degree() < V.degree():
     1278                        factors.append((V // W, e*k))
     1279                        V = W
     1280                        T = T // V
     1281                        if V.degree() == 0:
     1282                            if T.degree() == 0:
     1283                                break
     1284                            # T is of the form sum_{i=0}^n t_i X^{pi}
     1285                            T0 = T0.parent()([T[p*i].pth_root() for i in range(T.degree()//p + 1)])
     1286                            der = T0.derivative()
     1287                            e = p*e
     1288                            while der.is_zero():
     1289                                T0 = T0.parent()([T0[p*i].pth_root() for i in range(T0.degree()//p + 1)])
     1290                                der = T0.derivative()
     1291                                e = p*e
     1292                            T = T0.gcd(der)
     1293                            V = T0 // T
     1294                            k = 0
     1295                    else:
     1296                        T = T//V
    12311297        return Factorization(factors, unit=unit, sort=False)
    12321298
    12331299    def is_square(self, root=False):
     
    50105076
    50115077        return self.roots(ring=CC, multiplicities=False)
    50125078
     5079    def any_root(self, ring=None, degree=None, assume_squarefree=False, extend=False):
     5080        """
     5081        Returns a root of this polynomial in the given ring.
     5082
     5083        INPUT:
     5084       
     5085        - ``ring`` -- The ring in which a root is sought.  By default this is the coefficient ring.
     5086
     5087        - ``degree`` (None or nonzero integer) -- Used for polynomials over finite fields.  Returns a root of degree ``abs(degree)`` over the ground field.  If negative, also assumes that all factors of this polynomial are of degree ``abs(degree)``.  If None, returns a root of minimal degree contained within the given ring.
     5088
     5089        - ``assume_squarefree`` (bool) -- Used for polynomials over finite fields.  If True, this polynomial is assumed to be squarefree.
     5090
     5091        - ``extend`` -- Used for polynomials over finite fields.  If ``ring is None`` and ``extend`` is true, finds a root in the smallest extension of the base ring over which this polynomial has a root.  If ``ring is None`` and ``extend`` is false, only looks for a root in the base ring.  If ``ring`` is specified, then extend=True allows a root to lie in an extension of the given ring.
     5092
     5093        EXAMPLES::
     5094
     5095            sage: R.<x> = GF(11)[]
     5096            sage: f = 7*x^7 + 8*x^6 + 4*x^5 + x^4 + 6*x^3 + 10*x^2 + 8*x + 5
     5097            sage: f.any_root()
     5098            2
     5099            sage: f.factor()
     5100            (7) * (x + 9) * (x^6 + 10*x^4 + 6*x^3 + 5*x^2 + 2*x + 2)
     5101            sage: f = x^6 + 10*x^4 + 6*x^3 + 5*x^2 + 2*x + 2
     5102            sage: f.any_root(GF(11^6))
     5103            10*z6^5 + 3*z6^4 + 8*z6^3 + z6^2 + 10*z6 + 4
     5104            sage: sorted(f.roots(GF(11^6)))
     5105            [(z6^5 + z6^4 + 7*z6^3 + 2*z6^2 + 10*z6, 1), (z6^5 + 3*z6^4 + 8*z6^3 + 2*z6^2 + 3*z6 + 4, 1), (10*z6^5 + 2*z6^4 + 8*z6^3 + 9*z6^2 + z6, 1), (10*z6^5 + 3*z6^4 + 8*z6^3 + z6^2 + 10*z6 + 4, 1), (2*z6^5 + 8*z6^4 + 3*z6^3 + 6*z6 + 2, 1), (9*z6^5 + 5*z6^4 + 10*z6^3 + 8*z6^2 + 3*z6 + 1, 1)]
     5106            sage: f.any_root(extend=True)
     5107            10*z6^5 + 3*z6^4 + 8*z6^3 + z6^2 + 10*z6 + 4
     5108
     5109            sage: g = (x-1)*(x^2 + 3*x + 9) * (x^5 + 5*x^4 + 8*x^3 + 5*x^2 + 3*x + 5)
     5110            sage: g.any_root(ring=GF(11^10), degree=1)
     5111            1
     5112            sage: g.any_root(ring=GF(11^10), degree=2)
     5113            6*z10^9 + 7*z10^7 + 7*z10^6 + 3*z10^5 + z10^2 + z10 + 3
     5114            sage: g.any_root(ring=GF(11^10), degree=5)
     5115            7*z10^9 + 2*z10^8 + 6*z10^7 + 3*z10^6 + 2*z10^5 + 7*z10^3 + 7*z10^2 + 2
     5116
     5117        TESTS::
     5118
     5119            sage: R.<x> = GF(5)[]
     5120            sage: K.<a> = GF(5^12)
     5121            sage: for _ in range(40):
     5122            ...       f = R.random_element(degree=4)
     5123            ...       assert f(f.any_root(K)) == 0
     5124        """
     5125        if self.base_ring().is_finite() and self.base_ring().is_field():
     5126            if self.degree() < 0:
     5127                return ring(0)
     5128            if self.degree() == 0:
     5129                raise ValueError, "no roots A %s"%self
     5130            if not assume_squarefree:
     5131                SFD = self.squarefree_decomposition()
     5132                SFD.sort()
     5133                for f, e in SFD:
     5134                    try:
     5135                        return f.any_root(ring, degree, True, extend)
     5136                    except ValueError:
     5137                        pass
     5138            if self.degree() == 1 and (degree is None or degree == 1):
     5139                if ring is None:
     5140                    return -self[0]/self[1]
     5141                else:
     5142                    return ring(-self[0]/self[1])
     5143            q = self.base_ring().order()
     5144            if ring is None:
     5145                allowed_deg_mult = Integer(1)
     5146            else:
     5147                if not self.base_ring().is_prime_field(): raise NotImplementedError
     5148                if ring.characteristic() != self.base_ring().characteristic(): raise ValueError, "ring must be an extension of the base ring"
     5149                if not (ring.is_field() and ring.is_finite()): raise NotImplementedError
     5150                allowed_deg_mult = Integer(ring.factored_order()[0][1]) # generally it will be the quotient of this by the degree of the base ring.
     5151            if degree is None:
     5152                x = self.parent().gen()
     5153                if allowed_deg_mult == 1 and not extend:
     5154                    self = self.gcd(x**q-x)
     5155                    degree = -1
     5156                    if self.degree() == 0: raise ValueError, "no roots B %s"%self
     5157                else:
     5158                    xq = x
     5159                    d = Integer(0)
     5160                    while True: # one pass for roots that actually lie within ring.
     5161                        e = self.degree()
     5162                        if 2*d+2 > e: # this polynomial has no factors dividing allowed_deg_mult
     5163                            if allowed_deg_mult % e == 0:
     5164                                degree = -e
     5165                            break
     5166                        while d < allowed_deg_mult:
     5167                            d = d+1
     5168                            xq = xq**q % self
     5169                            if d.divides(allowed_deg_mult): break
     5170                        A = self.gcd(xq-x)
     5171                        if A != 1:
     5172                            self = A
     5173                            degree = -d
     5174                            break
     5175                        if d == allowed_deg_mult:
     5176                            break
     5177                    if degree is None:
     5178                        if allowed_deg_mult==1 and not extend: raise ValueError, "no roots C %s"%self
     5179                        xq = x
     5180                        d = Integer(0)
     5181                        while True: # now another for roots that will lie in an extension.
     5182                            e = self.degree()
     5183                            if 2*d+2 > e: # this polynomial is irreducible.
     5184                                degree = -e
     5185                                break
     5186                            while True: # we waste a little effort here in computing the xq again.
     5187                                d = d+1
     5188                                xq = xq**q % self
     5189                                if allowed_deg_mult.divides(d): break
     5190                            A = self.gcd(xq-x)
     5191                            if A != 1:
     5192                                self = A
     5193                                degree = -d
     5194                                break                           
     5195            if degree == 0: raise ValueError, "degree should be nonzero"
     5196            R = self.parent()
     5197            x = R.gen()
     5198            if degree > 0:
     5199                xq = x
     5200                d = 0
     5201                while True:
     5202                    e = self.degree()
     5203                    if 2*d > e:
     5204                        if degree != e:
     5205                            raise ValueError, "no roots D %s"%self
     5206                        break
     5207                    d = d+1
     5208                    xq = xq**q % self
     5209                    if d == degree:
     5210                        break
     5211                    A = self.gcd(xq-x)
     5212                    if A != 1:
     5213                        self = self // A
     5214                if d == degree:
     5215                    self = self.gcd(xq-x)
     5216                    if self.degree() == 0:
     5217                        raise ValueError, "no roots E %s"%self
     5218            else:
     5219                degree = -degree
     5220            if ring is None:
     5221                if degree == 1:
     5222                    ring = self.base_ring()
     5223                else:
     5224                    ring = self.base_ring().extension(degree) # this won't work yet.
     5225            # now self has only roots of degree ``degree``.
     5226            # for now, we only implement the Cantor-Zassenhaus split
     5227            k = self.degree() // degree
     5228            if k == 1:
     5229                try:
     5230                    return self.roots(ring, multiplicities=False)[0] # is there something better to do here?
     5231                except IndexError:
     5232                    raise ValueError, "no roots F %s"%self
     5233            if q % 2 == 0:
     5234                T = x
     5235                while True:
     5236                    C = T
     5237                    for i in range(degree-1):
     5238                        C = T + C**q % self
     5239                    h = self.gcd(C)
     5240                    hd = h.degree()
     5241                    if hd == 0 or hd == self.degree():
     5242                        T = T * x**q
     5243                    else:
     5244                        if 2*hd <= self.degree():
     5245                            return h.any_root(ring, -degree, True, extend)
     5246                        else:
     5247                            return (self//h).any_root(ring, -degree, True, extend)
     5248            else:
     5249                from sage.rings.arith import power_mod
     5250                while True:
     5251                    T = R.random_element(2*degree-1)
     5252                    if T == 0: continue
     5253                    T = T.monic()
     5254                    h = self.gcd(power_mod(T, Integer((q**degree-1)/2), self)-1)
     5255                    hd = h.degree()
     5256                    if hd != 0 and hd != self.degree():
     5257                        if 2*hd <= self.degree():
     5258                            return h.any_root(ring, -degree, True, extend)
     5259                        else:
     5260                            return (self//h).any_root(ring, -degree, True, extend)
     5261        else:
     5262            return self.roots(ring=ring, multiplicities=False)[0]
     5263       
    50135264
    50145265    def variable_name(self):
    50155266        """