Ticket #8335: trac_8335-pseudo_conway-5.8.b0.patch

File trac_8335-pseudo_conway-5.8.b0.patch, 50.6 KB (added by jpflori, 9 years ago)
  • sage/rings/finite_rings/constructor.py

    # HG changeset patch
    # User David Roe <roed@math.harvard.edu>
    # Date 1266946097 18000
    # Node ID cf291487ba675ef650010e9bdef907bc9853a4fa
    # Parent  54981fe37e4accc93ee69ca514dcef9b3ecb63e5
    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 --git a/sage/rings/finite_rings/constructor.py b/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.
     
    206213    -  ``check_irreducible`` - verify that the polynomial
    207214       modulus is irreducible
    208215
     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).
     219
    209220    - ``proof`` -- bool (default: True); if True use provable
    210221      primality test; otherwise only use pseudoprimality test.
    211222   
     
    227238        Finite Field in a of size 3^2
    228239        sage: charpoly(a, 'y')
    229240        y^2 + 2*y + 2
     241        sage: L = GF(27); L
     242        Finite Field in z3 of size 3^3
    230243   
    231244    We illustrate the proof flag.  The following example would hang
    232245    for a very long time if we didn't use proof=False.  (NOTE: Magma
     
    337350
    338351    """
    339352    def create_key_and_extra_args(self, order, name=None, modulus=None, names=None,
    340                                   impl=None, proof=None, **kwds):
     353                                  impl=None, prefix='z', proof=None, **kwds):
    341354        """
    342355        EXAMPLES::
    343356       
    344357            sage: GF.create_key_and_extra_args(9, 'a')
    345             ((9, ('a',), 'conway', None, '{}', 3, 2, True), {})
     358            ((9, ('a',), 'conwayz', None, '{}', 3, 2, True), {})
    346359            sage: GF.create_key_and_extra_args(9, 'a', foo='value')
    347             ((9, ('a',), 'conway', None, "{'foo': 'value'}", 3, 2, True), {'foo': 'value'})
     360            ((9, ('a',), 'conwayz', None, "{'foo': 'value'}", 3, 2, True), {'foo': 'value'})
    348361        """
    349362        from sage.structure.proof.all import WithProof, arithmetic
    350363        if proof is None: proof = arithmetic()
     
    360373                n = integer.Integer(1)
    361374            elif arith.is_prime_power(order):
    362375                if not names is None: name = names
    363                 name = normalize_names(1,name)
    364376
    365377                p, n = arith.factor(order)[0]
    366378
    367379                if modulus is None or modulus == "default":
    368                     if exists_conway_polynomial(p, n):
     380                    if exists_conway_polynomial(p, n) or name is None:
    369381                        modulus = "conway"
    370382                    else:
    371383                        if p == 2:
    372384                            modulus = "minimal_weight"
    373385                        else:
    374386                            modulus = "random"
     387                elif p != 2 and modulus == 'first_lexicographic':
     388                    R = FiniteField(p)['x']
     389                    x = R.gen()
     390                    for f in R.polynomials(of_degree=n):
     391                        if f.is_irreducible():
     392                            f = f.monic() # probably not necessary...
     393                            break
     394                    modulus = f
     395
     396                if modulus == "conway":
     397                    modulus += prefix
     398                    if name is None:
     399                        name = prefix + str(n)
    375400                elif modulus == "random":
    376401                    modulus += str(random.randint(0, 1<<128))
     402                name = normalize_names(1,name)
    377403
    378404                if isinstance(modulus, (list, tuple)):
    379405                    modulus = FiniteField(p)['x'](modulus)
     
    445471                    if modulus.degree() != n:
    446472                        raise ValueError("The degree of the modulus does not correspond to the cardinality of the field.")
    447473                if name is None:
    448                     raise TypeError("you must specify the generator name.")
    449                 if order < zech_log_bound
     474                    raise TypeError("you must specify the generator name or use Conway polynomials.")
     475                if order < zech_log_bound and (impl is None or impl =='givaro'):
    450476                    # DO *NOT* use for prime subfield, since that would lead to
    451477                    # a circular reference in the call to ParentWithGens in the
    452478                    # __init__ method.
     
    466492        EXAMPLES::
    467493       
    468494            sage: key, extra = GF.create_key_and_extra_args(9, 'a'); key
    469             (9, ('a',), 'conway', None, '{}', 3, 2, True)
     495            (9, ('a',), 'conwayz', None, '{}', 3, 2, True)
    470496            sage: K = GF.create_object(0, key); K
    471497            Finite Field in a of size 3^2
    472498            sage: GF.other_keys(key, K)
     
    604630    """
    605631    return sage.databases.conway.ConwayPolynomials().has_polynomial(p,n)
    606632
     633pseudo_conway_poly = {}
     634
     635class PseudoConwayPolyTree(SageObject):
     636    """
     637    An object holding references to pseudo-Conway polynomials for divisors of the given degree, so that they aren't garbage collected.
     638    """
     639    def __init__(self, p, n, nodes_dict, f):
     640        """
     641        INPUT::
     642
     643        - p -- The prime for which this defines an extension of GF(p).
     644        - n -- The degree of the extension.
     645        - 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.
     646                        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).
     647        - f -- The polynomial defining this extension.
     648
     649        EXAMPLES::
     650
     651            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     652            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     653            sage: PCPT.get_pseudo_conway_poly(3)
     654            x^3 + x + 1
     655        """
     656        self.p = p
     657        self.n = n
     658        self.nodes_dict = nodes_dict
     659        self.f = f
     660       
     661    def get_pseudo_conway_poly(self, m):
     662        """
     663        Returns the pseudo-Conway polynomial associated to a divisor of the degree of this tree.
     664
     665        EXAMPLES::
     666       
     667            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     668            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     669            sage: PCPT.get_pseudo_conway_poly(3)
     670            x^3 + x + 1
     671            sage: PCPT.get_pseudo_conway_poly(4)
     672            Traceback (most recent call last):
     673            ...
     674            RuntimeError: 4 does not divide 6
     675        """
     676        if m == self.n:
     677            return self.f
     678        if not m.divides(self.n):
     679            raise RuntimeError, "%s does not divide %s"%(m, self.n)
     680        try:
     681            return pseudo_conway_poly[self.p][m]().f
     682        except (KeyError, AttributeError):
     683            return conway_polynomial(self.p, m)
     684           
     685    def check_consistency(self):
     686        """
     687        Checks that the pseudo-Conway polynomials dividing the degree of this tree satisfy the required compatibility conditions.
     688
     689        EXAMPLES::
     690
     691            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     692            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     693            sage: PCPT.check_consistency()
     694            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 60, False) # long
     695            sage: PCPT.check_consistency() # long
     696        """
     697        p, n = self.p, self.n
     698        K = GF(p**n, modulus = self.f, names='a')
     699        a = K.gen()
     700        for m in n.divisors():
     701            assert (a**((p**n-1)//(p**m-1))).minimal_polynomial() == self.get_pseudo_conway_poly(m)
     702
     703def find_pseudo_conway_polynomial_tree(p, n, use_database=True):
     704    """
     705    Returns an object holding references to a set of consistent pseudo-Conway polynomials for degrees dividing n.
     706
     707    A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
     708   
     709    - `f` is irreducible.
     710    - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
     711    - 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`.
     712    - `f` is lexicographically least among all such polynomials, under a certain ordering.
     713
     714    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.
     715
     716    INPUT:
     717
     718    - `p` -- a prime.
     719    - `n` -- an integer greater than 1.
     720    - `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)
     721
     722    EXAMPLES::
     723
     724        sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     725        sage: PCPT = find_pseudo_conway_polynomial_tree(2, 12, False)
     726        sage: PCPT.f
     727        x^12 + x^8 + x^7 + x^6 + x^4 + x^3 + 1
     728    """
     729    if not pseudo_conway_poly.has_key(p):
     730        pseudo_conway_poly[p] = {}
     731    pdict = pseudo_conway_poly[p]
     732    n = sage.rings.integer.Integer(n)
     733    if pdict.has_key(n):
     734        pcp = pdict[n]()
     735        if pcp is not None:
     736            return pcp
     737    if use_database and exists_conway_polynomial(p, n):
     738        ans = PseudoConwayPolyTree(p, n, True, GF(p)['x'](conway_polynomial(p, n)))
     739        pdict[n] = weakref.ref(ans)
     740        return ans
     741    if n == 1 or n.is_prime():
     742        ans = PseudoConwayPolyTree(p, n, False, compute_pseudo_conway_polynomial(p, n, None))
     743    nodes = {}
     744    for q in n.prime_factors():
     745        nodes[q] = find_pseudo_conway_polynomial_tree(p, n//q, use_database)
     746    ans = PseudoConwayPolyTree(p, n, nodes, compute_pseudo_conway_polynomial(p, n, nodes))
     747    pdict[n] = weakref.ref(ans)
     748    return ans
     749
     750def _find_pow_of_frobenius(p, x, y, field_degree):
     751    """
     752    Finds the power of Frobenius which yields `x` when applied to `y`.
     753
     754    INPUT:
     755
     756    - `p` -- a prime.
     757    - ``field_degree`` -- a positive integer
     758    - `x` -- an element of a field `K` of characeristic p so that the multiplicative order of `x` is ``p^field_degree``.
     759    - `y` -- an element of `K` with the same minimal polynomial.
     760   
     761    OUTPUT:
     762   
     763    An element `i` of ``Integers(field_degree)`` so that `x = y^{p^i}`
     764
     765    EXAMPLES::
     766
     767        sage: from sage.rings.finite_rings.constructor import _find_pow_of_frobenius
     768        sage: K.<a> = GF(3^14)
     769        sage: x = K.multiplicative_generator()
     770        sage: y = x^27
     771        sage: _find_pow_of_frobenius(3, x, y, 14)
     772        11
     773    """
     774    assert x.minimal_polynomial() == y.minimal_polynomial()
     775    assert x.multiplicative_order() == p**field_degree - 1
     776    assert y.multiplicative_order() == p**field_degree - 1
     777    from integer_mod import mod
     778    for i in range(field_degree):
     779        if x == y: break
     780        y = y**p
     781    else:
     782        raise RuntimeError, "No appropriate power of Frobenius found"
     783    return mod(i, field_degree)
     784
     785def _crt_non_coprime(running, a):
     786    """
     787    Extends the ``crt`` method of IntegerMod to the case of non-relatively prime modulus.
     788
     789    EXAMPLES::
     790
     791        sage: from sage.rings.finite_rings.constructor import _crt_non_coprime
     792        sage: a = _crt_non_coprime(mod(14, 18), mod(20,30)); a
     793        50
     794        sage: a.modulus()
     795        90
     796        sage: _crt_non_coprime(mod(13, 18), mod(20,30))
     797        Traceback (most recent call last):
     798        ...
     799        AssertionError
     800    """
     801    g = running.modulus().gcd(a.modulus())
     802    if g == 1:
     803        return running.crt(a)
     804    else:
     805        assert running % g == a % g
     806        running_modulus = running.modulus()
     807        a_modulus = a.modulus()
     808        for qq in g.prime_divisors():
     809            a_val_unit = a_modulus.val_unit(qq)
     810            running_val_unit = running_modulus.val_unit(qq)
     811            if a_val_unit[0] > running_val_unit[0]:
     812                running_modulus = running_val_unit[1]
     813            else:
     814                a_modulus = a_val_unit[1]
     815        return (running % running_modulus).crt(a % a_modulus)
     816
     817def _frobenius_shift(K, generators, check_only=False):
     818    """
     819    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.
     820
     821    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
     822
     823    ``generators[q1]^((p^(n/q1)-1)/(p^m-1)) == generators[q2]^((p^(n/q2)-1)/(p^m-1))``
     824
     825    holds.
     826
     827    INPUT:
     828
     829    - K -- a finite field of degree n   
     830    - generators -- a dictionary, indexed by prime divisors `q` of `n`, whose entries are elements of `K` satisfying the `n/q` pseudo-Conway polynomial.
     831    - check_only -- if True, just asserts that the generators given form a compatible system.
     832
     833    EXAMPLES::
     834
     835        sage: R.<x> = GF(2)[]
     836        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
     837        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
     838        sage: f12 = x^12 + x^10 + x^9 + x^8 + x^4 + x^2 + 1
     839        sage: K.<a> = GF(2^60, modulus='first_lexicographic')
     840        sage: x30 = f30.any_root(K)
     841        sage: x20 = f20.any_root(K)
     842        sage: x12 = f12.any_root(K)
     843        sage: generators = {2: x30, 3: x20, 5: x12}
     844        sage: from sage.rings.finite_rings.constructor import _frobenius_shift, _find_pow_of_frobenius
     845        sage: _frobenius_shift(K, generators)
     846        sage: _find_pow_of_frobenius(2, x30, generators[2], 30)
     847        0
     848        sage: _find_pow_of_frobenius(2, x20, generators[3], 20)
     849        13
     850        sage: _find_pow_of_frobenius(2, x12, generators[5], 12)
     851        8
     852    """
     853    p = K.characteristic()
     854    n = K.degree()
     855    if len(generators) == 1:
     856        return generators
     857    compatible = {}
     858    from integer_mod import mod
     859    for m in n.divisors():
     860        compatible[m] = {}
     861    for q, x in generators.iteritems():
     862        for m in (n//q).divisors():
     863            compatible[m][q] = x**((p**(n//q)-1)//(p**m-1))
     864    if check_only:
     865        for m in n.divisors():
     866            try:
     867                q, x = compatible[m].popitem()
     868            except KeyError:
     869                break
     870            for qq, xx in compatible[m].iteritems():
     871                assert x == xx
     872        return
     873    crt = {}
     874    qlist = sorted(generators.keys())
     875    for j in range(1, len(qlist)):
     876        for i in range(j):
     877            crt[(i, j)] = []
     878    for m in n.divisors():
     879        mqlist = sorted(compatible[m].keys())
     880        for k in range(1,len(mqlist)):
     881            j = qlist.index(mqlist[k])
     882            i = qlist.index(mqlist[k-1])
     883            crt[(i,j)].append(_find_pow_of_frobenius(p, compatible[m][qlist[j]], compatible[m][qlist[i]], m))
     884    from integer_mod import mod
     885    pairs = crt.keys()
     886    for i, j in pairs:
     887        L = crt[(i,j)]
     888        #print qlist[i], qlist[j]
     889        #print " ".join(["(%s, %s)"%(mm, mm.modulus()) for mm in L])
     890        running = mod(0,1)
     891        for a in L:
     892            running = _crt_non_coprime(running, a)
     893        #print "(%s, %s)"%(running, running.modulus())
     894        crt[(i,j)] = [(mod(running, q**(running.modulus().valuation(q))), running.modulus().valuation(q)) for q in qlist]
     895        crt[(j,i)] = [(-a, level) for a, level in crt[(i,j)]]
     896    # Let x_j be the power of Frobenius we apply to generators[qlist[j]], for 0 < j < len(qlist)
     897    # We have some direct conditions on the x_j: x_j reduces to each entry in crt[(0,j)].
     898    # But we also have the equations x_j - x_i reduces to each entry in crt[(i,j)].
     899    # We solve for x_j one prime at a time.  For each prime, we have an equations of the form
     900    # 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
     901    # (possibly 0, possibly different) of the same prime.
     902
     903    # We can set x_0=0 everywhere, can get an initial setting of x_j from the c_0j.
     904    # We go through prime by prime.
     905    import bisect
     906    frob_powers=[mod(0,1) for q in qlist]
     907    def find_leveller(qindex, level, x, xleveled, searched, i):
     908        searched[i] = True
     909        crt_possibles = []
     910        for j in range(1,len(qlist)):
     911            if i==j: continue
     912            if crt[(i,j)][qindex][1] >= level:
     913                if xleveled[j]:
     914                    return [j]
     915                elif not searched.has_key(j):
     916                    crt_possibles.append(j)
     917        for j in crt_possibles:
     918            path = find_leveller(qindex, level, x, xleveled, searched, j)
     919            if path is not None:
     920                path.append(j)
     921                return path
     922        return None
     923    def propogate_levelling(qindex, level, x, xleveled, i):
     924        for j in range(1, len(qlist)):
     925            if i==j: continue
     926            if not xleveled[j] and crt[(i,j)][qindex][1] >= level:
     927                newxj = x[i][0] + crt[(i,j)][qindex][0]
     928                x[j] = (newxj, min(x[i][1], crt[(i,j)][qindex][1]))
     929                xleveled[j] = True
     930                propogate_levelling(qindex, level, x, xleveled, j)
     931               
     932    for qindex in range(len(qlist)):
     933        q = qlist[qindex]
     934        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.               
     935        # We first check that our equations are consistent and determine which powers of q occur as moduli.
     936        levels = []
     937        for j in range(2, len(qlist)):
     938            for i in range(j):
     939                # we need crt[(0,j)] = crt[(0,i)] + crt[(i,j)]
     940                if i != 0: assert x[j][0] == x[i][0] + crt[(i,j)][qindex][0]
     941                level = crt[(i,j)][qindex][1]
     942                if level > 0:
     943                    ins = bisect.bisect_left(levels,level)
     944                    if ins == len(levels):
     945                        levels.append(level)
     946                    elif levels[ins] != level:
     947                        levels.insert(ins, level)
     948        for level in levels:
     949            xleveled = [0] + [x[i][1] >= level for i in range(1,len(qlist))]
     950            while True:
     951                try:
     952                    i = xleveled.index(False, 1)
     953                    searched = {}
     954                    levelling_path = find_leveller(qindex, level, x, xleveled, searched, i)
     955                    if levelling_path is None:
     956                        # Any lift will work, since there are no constraints
     957                        x[i] = (mod(x[i][0].lift(), q**level), level)
     958                        xleveled[i] = True
     959                        propogate_levelling(qindex, level, x, xleveled, i)
     960                    else:
     961                        levelling_path.append(i)
     962                        for m in range(1,len(path)):
     963                            # this point on the path may have already been leveled in a previous propogation
     964                            if not xleveled[path[m]]:
     965                                newx = x[path[m-1]][0] + crt[(path[m-1],path[m])][qindex][0]
     966                                x[path[m]] = (newx, min(x[path[m-1]][1], crt[(path[m-1],path[m])][qindex][1]))
     967                                xleveled[path[m]] = True
     968                                propogate_levelling(qindex, level, x, xleveled, path[m])
     969                except ValueError:
     970                    break
     971        for j in range(1,len(qlist)):
     972            frob_powers[j] = frob_powers[j].crt(x[j][0])
     973    for j in range(1, len(qlist)):
     974        #print "frob_power[%s] = mod(%s, %s)"%(qlist[j], frob_powers[j], frob_powers[j].modulus())
     975        generators[qlist[j]] = generators[qlist[j]]**(p**(-frob_powers[j]).lift())
     976    _frobenius_shift(K, generators, check_only=True)
     977
     978def compute_pseudo_conway_polynomial(p, n, nodes):
     979    """
     980    Computes a pseudo-Conway polynomial over the prime `p` of degree `n`.
     981
     982    A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
     983   
     984    - `f` is irreducible.
     985    - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
     986    - 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`.
     987    - `f` is lexicographically least among all such polynomials, under a certain ordering.
     988
     989    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.
     990
     991    INPUT:
     992   
     993    - `p` -- a prime
     994    - `n` -- a positive integer
     995    - ``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`.
     996
     997    OUTPUT:
     998
     999    A pseudo-Conway polynomial over the prime `p` of degree `n`.
     1000
     1001    ALGORITHM:
     1002   
     1003    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.
     1004
     1005    REFERENCE:
     1006
     1007    .. [HL99] Heath L. and Loehr N. (1999).
     1008       New algorithms for generating Conway polynomials over finite fields.
     1009       Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms, pp. 429-437.
     1010
     1011    EXAMPLES::
     1012
     1013        sage: from sage.rings.finite_rings.constructor import compute_pseudo_conway_polynomial, find_pseudo_conway_polynomial_tree
     1014        sage: compute_pseudo_conway_polynomial(next_prime(10000),11,None) # long because the arithmetic in FiniteField_ext_pari fields is slow.
     1015        x^11 + x + 7
     1016        sage: PCPT30 = find_pseudo_conway_polynomial_tree(2, 30, False)
     1017        sage: PCPT20 = find_pseudo_conway_polynomial_tree(2, 20, False)
     1018        sage: PCPT12 = find_pseudo_conway_polynomial_tree(2, 12, False)
     1019        sage: compute_pseudo_conway_polynomial(2, 60, {2:PCPT30,3:PCPT20,5:PCPT12})
     1020        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
     1021    """
     1022    k = GF(p)
     1023    from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     1024    R = PolynomialRing(k, names='x')
     1025    if n == 1:
     1026        return R([-k.multiplicative_generator(), 1])
     1027    if nodes is None:
     1028        # n is prime.  We require only that the chosen generator (x) be a multiplicative generator for the units.
     1029        if p in [2,3,5,7,11]: # here we have Cunningham data
     1030            from sage.rings.factorint import factor_cunningham
     1031            F = [(p**n-1)//a for a, e in factor_cunningham(n)]
     1032        else:
     1033            F = (p**n-1).prime_factors() # could be VERY expensive
     1034        x = R.gen()
     1035        from sage.rings.arith import power_mod
     1036        for f in R.polynomials(of_degree=n):
     1037            # By using this iterator our chosen polynomial will be fairly sparse.
     1038            if f.is_irreducible():
     1039                if not f.is_monic():
     1040                    f = f.monic()
     1041                for m in F:
     1042                    if power_mod(x, m, f) == 1:
     1043                        break
     1044                else:
     1045                    return f
     1046    # Our strategy is as follows:
     1047    # Work in an arbitrary field K of order p**n.
     1048    # We're looking for a multiplicative generator of K whose
     1049    K = GF(p**n, modulus="first_lexicographic", names='a')
     1050    r = p**n - 1
     1051    xi = {}
     1052    for q in nodes.iterkeys():
     1053        #print nodes[q].f.degree(), n//q
     1054        xi[q] = nodes[q].f.any_root(K, -n//q, True)
     1055    _frobenius_shift(K, xi)
     1056           
     1057    q, x = xi.popitem()
     1058    v = p**(n//q) - 1
     1059    for q, xitem in xi.iteritems():
     1060        fi = p**(n//q) - 1
     1061        g, alpha, beta = v.xgcd(fi)
     1062        x = x**beta * xitem**alpha
     1063        v = v.lcm(fi)
     1064    g = r // v
     1065    #print g, n, g.divides(p**n-1)
     1066    z = x.nth_root(g, cunningham=(p in [2,3,5,7,11]))
     1067    groot1 = K.multiplicative_generator()**v
     1068    while z.multiplicative_order() != r:
     1069        z *= groot1
     1070    #print z
     1071    return z.minimal_polynomial()
    6071072
    6081073zech_log_bound = 2**16
  • sage/rings/finite_rings/element_ext_pari.py

    diff --git a/sage/rings/finite_rings/element_ext_pari.py b/sage/rings/finite_rings/element_ext_pari.py
    a b  
    862862       
    863863        EXAMPLES::
    864864       
    865             sage: from sage.rings.finite_rings.finite_field_ext_pari import FiniteField_ext_pari
    866865            sage: k = GF(next_prime(10**10))
    867866            sage: a = k(17)/k(19)
    868867            sage: b = a.lift(); b
  • sage/rings/finite_rings/finite_field_ext_pari.py

    diff --git a/sage/rings/finite_rings/finite_field_ext_pari.py b/sage/rings/finite_rings/finite_field_ext_pari.py
    a b  
    5050    - ``modulus`` -- you may provide a polynomial to use for reduction or
    5151      a string:
    5252
    53       - ``'conway'`` -- force the use of a Conway polynomial, will
    54         raise a ``RuntimeError`` if none is found in the database
     53      - ``'conway'`` -- force the use of a Conway polynomial; if
     54        none is found in the database, will use a pseudo-Conway
     55        polynomial instead;
    5556      - ``'random'`` -- use a random irreducible polynomial
    5657      - ``'default'`` -- a Conway polynomial is used if found. Otherwise
    5758        a random polynomial is used
     
    165166        We do NOT yet define natural consistent inclusion maps
    166167        between different finite fields.
    167168    """
    168     def __init__(self, q, name, modulus=None):
     169    def __init__(self, q, name, modulus='conway'):
    169170        """
    170171        Create finite field of order `q` with variable printed as name.
    171172           
     
    205206        self.__order = q
    206207        self.__is_field = True
    207208
    208         if modulus is None or modulus == "default":
    209             from constructor import exists_conway_polynomial
    210             if exists_conway_polynomial(self.__char, self.__degree):
    211                 modulus = "conway"
    212             else:
    213                 modulus = "random"
    214 
    215         if isinstance(modulus,str):
    216             if modulus == "conway":
    217                 from constructor import conway_polynomial
    218                 modulus = conway_polynomial(self.__char, self.__degree)
     209        if isinstance(modulus, str):
     210            if modulus.startswith('conway'):
     211                from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     212                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.
     213                self._prefix = modulus[6:]
     214                if len(self._prefix) == 0: self._prefix = 'z'
     215                modulus = self._PCPT.get_pseudo_conway_poly(self.__degree)
    219216            elif modulus == "random":
    220217                # The following is fast/deterministic, but has serious problems since
    221218                # it crashes on 64-bit machines, and I can't figure out why:
  • sage/rings/finite_rings/finite_field_givaro.py

    diff --git a/sage/rings/finite_rings/finite_field_givaro.py b/sage/rings/finite_rings/finite_field_givaro.py
    a b  
    2727    - ``name`` -- (default: ``'a'``) variable used for
    2828      :meth:`~sage.rings.finite_rings.element_givaro.FiniteField_givaroElement.poly_repr()`
    2929
    30     - ``modulus`` -- (default: ``None``, a conway polynomial is used if found.
    31       Otherwise a random polynomial is used) A minimal polynomial to use for
    32       reduction or ``'random'`` to force a random irreducible polynomial.
     30    - ``modulus`` -- (default: ``conway``, a conway polynomial is used if in
     31      the database, otherwise a pseudo-Conway polynomial is generated)
     32      A minimal polynomial to use for reduction or ``'random'`` to force a
     33      random irreducible polynomial.
    3334
    3435    - ``repr`` -- (default: ``'poly'``) controls the way elements are printed
    3536      to the user:
     
    8687        sage: sage.rings.finite_rings.finite_field_givaro.FiniteField_givaro(9,repr='log').gen()
    8788        1
    8889    """
    89     def __init__(self, q, name="a",  modulus=None, repr="poly", cache=False):
     90    def __init__(self, q, name="a",  modulus="conway", repr="poly", cache=False):
    9091        """
    9192        Initialize ``self``.
    9293
     
    125126        self._kwargs['cache'] = cache
    126127
    127128        self._is_conway = False
    128         if modulus is None or modulus == 'conway':
     129        if isinstance(modulus, str) and modulus.startswith('conway'):
    129130            if k == 1:
    130131                modulus = 'random' # this will use the gfq_factory_pk function.
    131             elif ConwayPolynomials().has_polynomial(p, k):
    132                 from sage.rings.finite_rings.constructor import conway_polynomial
    133                 modulus = conway_polynomial(p, k)
    134                 self._is_conway = True
    135             elif modulus is None:
    136                 modulus = 'random'
    137132            else:
    138                 raise ValueError, "Conway polynomial not found"
     133                from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     134                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.
     135                self._prefix = modulus[6:]
     136                if len(self._prefix) == 0: self._prefix = 'z'
     137                modulus = self._PCPT.get_pseudo_conway_poly(k)
     138                # The following is only used for _gap_init_, so we check to see if the defining polynomial is actually a conway polynomial.
     139                self._is_conway = ConwayPolynomials().has_polynomial(p, k)
    139140               
    140141        from sage.rings.polynomial.all import is_Polynomial
    141142        if is_Polynomial(modulus):
  • sage/rings/finite_rings/finite_field_ntl_gf2e.py

    diff --git a/sage/rings/finite_rings/finite_field_ntl_gf2e.py b/sage/rings/finite_rings/finite_field_ntl_gf2e.py
    a b  
    22Finite Fields of Characteristic 2
    33"""
    44
    5 
    65from sage.rings.finite_rings.finite_field_base import FiniteField
    76from sage.libs.pari.all import pari
    87from finite_field_ext_pari import FiniteField_ext_pari
     
    2221    """
    2322    if globals().has_key("GF2"):
    2423        return
    25     global ResidueField_generic, is_FiniteField, exists_conway_polynomial, conway_polynomial, Cache_ntl_gf2e, GF, GF2, is_Polynomial
     24    global ResidueField_generic, is_FiniteField, exists_conway_polynomial, conway_polynomial, Cache_ntl_gf2e, GF, GF2, is_Polynomial, ConwayPolynomials
    2625    import sage.rings.residue_field
    2726    ResidueField_generic = sage.rings.residue_field.ResidueField_generic
    2827
     
    3938    import sage.rings.finite_rings.element_ntl_gf2e
    4039    Cache_ntl_gf2e = sage.rings.finite_rings.element_ntl_gf2e.Cache_ntl_gf2e
    4140
     41    import sage.databases.conway
     42    ConwayPolynomials = sage.databases.conway.ConwayPolynomials
     43
    4244    import sage.rings.finite_rings.constructor
    4345    GF = sage.rings.finite_rings.constructor.GF
    4446    GF2 = GF(2)
     
    5961    - ``modulus`` -- you may provide a polynomial to use for reduction or
    6062      a string:
    6163
    62       - ``'conway'`` -- force the use of a Conway polynomial, will
    63         raise a ``RuntimeError`` if ``None`` is found in the database;
     64      - ``'conway'`` -- force the use of a Conway polynomial; if
     65        none is found in the database use a generated pseudo-Conway
     66        polynomial instead;
    6467      - ``'minimal_weight'`` -- use a minimal weight polynomial, should
    6568        result in faster arithmetic;
     69      - ``'first_lexicographic'`` -- uses the first irreducible polynomial;
    6670      - ``'random'`` -- use a random irreducible polynomial.
    67       - ``'default'`` -- a Conway polynomial is used if found. Otherwise
    68         a sparse polynomial is used.
    6971
    7072    - ``repr`` -- controls the way elements are printed to the user:
    7173                 (default: ``'poly'``)
     
    9698        sage: k.<a> = GF(2^211, modulus='conway')
    9799        sage: k.modulus()
    98100        x^211 + x^9 + x^6 + x^5 + x^3 + x + 1
    99         sage: k.<a> = GF(2^411, modulus='conway')
    100         Traceback (most recent call last):
    101         ...
    102         RuntimeError: requested conway polynomial not in database.
     101        sage: k.<a> = GF(2^23, modulus='conway')
     102        sage: a.multiplicative_order() == k.order()-1
     103        True
    103104    """
    104105
    105106    def __init__(self, q, names="a",  modulus=None, repr="poly"):
     
    139140        self._kwargs = {'repr':repr}
    140141        self._is_conway = False
    141142
    142         if modulus is None or modulus == 'default':
    143             if exists_conway_polynomial(p, k):
    144                 modulus = "conway"
    145             else:
    146                 modulus = "minimal_weight"
    147         if modulus == "conway":
    148             modulus = conway_polynomial(p, k)
    149             self._is_conway = True
     143        if isinstance(modulus, str) and modulus.startswith('conway'):
     144             from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     145             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.
     146             self._prefix = modulus[6:]
     147             if len(self._prefix) == 0: self._prefix = 'z'
     148             modulus = self._PCPT.get_pseudo_conway_poly(k)
     149             self._is_conway = ConwayPolynomials().has_polynomial(p, k)
    150150        if is_Polynomial(modulus):
    151151            modulus = modulus.list()
    152152        self._cache = Cache_ntl_gf2e(self, k, modulus)
  • sage/rings/polynomial/polynomial_element.pyx

    diff --git a/sage/rings/polynomial/polynomial_element.pyx b/sage/rings/polynomial/polynomial_element.pyx
    a b  
    12151215            sage: f = QQbar['x'](1)
    12161216            sage: f.squarefree_decomposition()
    12171217            1
    1218         """
     1218
     1219        TESTS::
     1220
     1221            sage: K.<a> = GF(3^2)
     1222            sage: R.<x> = K[]
     1223            sage: f = x^243+2*x^81+x^9+1
     1224            sage: f.squarefree_decomposition()
     1225            (x^27 + 2*x^9 + x + 1)^9
     1226            sage: f = x^243+a*x^27+1
     1227            sage: f.squarefree_decomposition()
     1228            (x^9 + (2*a + 1)*x + 1)^27
     1229
     1230            sage: for K in [GF(2^18,'a'), GF(3^2,'a'), GF(47^3,'a')]:
     1231            ...       R.<x> = K[]
     1232            ...       if K.characteristic() < 5: m = 4
     1233            ...       else: m = 1
     1234            ...       for _ in range(m):
     1235            ...           f = (R.random_element(4)^3*R.random_element(m)^(m+1))(x^6)
     1236            ...           F = f.squarefree_decomposition()
     1237            ...           assert F.prod() == f
     1238            ...           for i in range(len(F)):
     1239            ...               assert gcd(F[i][0], F[i][0].derivative()) == 1
     1240            ...               for j in range(len(F)):
     1241            ...                   if i == j: continue
     1242            ...                   assert gcd(F[i][0], F[j][0]) == 1
     1243            ...
     1244        """
     1245        if not self.base_ring().is_unique_factorization_domain():
     1246            raise NotImplementedError, "Squarefree decomposition not implemented for " + str(self.parent())
     1247
    12191248        if self.degree() == 0:
    12201249            return Factorization([], unit=self[0])
    12211250
    1222         # Wikipedia says this works for arbitrary fields of
    1223         # characteristic 0.
    1224 
    1225         if self.base_ring().characteristic() != 0:
    1226             raise NotImplementedError("Squarefree decomposition not implemented for " + str(self.parent()))
    1227 
    1228         f = [self]
    1229         cur = self
    1230         while cur.degree() > 0:
    1231             cur = cur.gcd(cur.derivative())
    1232             f.append(cur)
    1233 
    1234         g = []
    1235         for i in range(len(f) - 1):
    1236             g.append(f[i] // f[i+1])
    1237 
    1238         a = []
    1239         for i in range(len(g) - 1):
    1240             a.append(g[i] // g[i+1])
    1241         if g == []: # do not factor the unit part
    1242             return Factorization([(self,1)])
    1243         a.append(g[-1])
    1244 
     1251        p = self.base_ring().characteristic()
    12451252        factors = []
    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)
     1253        if p == 0:
     1254            # This is Yun's algorithm which works for arbitrary ringss of
     1255            # characteristic 0.
     1256            f = [self]
     1257            cur = self
     1258            while cur.degree() > 0:
     1259                cur = cur.gcd(cur.derivative())
     1260                f.append(cur)
     1261
     1262            g = []
     1263            for i in range(len(f) - 1):
     1264                g.append(f[i] // f[i+1])
     1265
     1266            a = []
     1267            for i in range(len(g) - 1):
     1268                a.append(g[i] // g[i+1])
     1269            a.append(g[-1])
     1270
     1271            unit = f[-1]
     1272            for i in range(len(a)):
     1273                if a[i].degree() > 0:
     1274                    factors.append((a[i], i+1))
     1275                else:
     1276                    unit = unit * a[i].constant_coefficient() ** (i + 1)
     1277        else:
     1278            # In characteristic p, we use algorithm 3.4.2 in Cohen's Course in
     1279            # Computational Algebraic Number Theory
     1280            # It is basically Yun's algorithm with special treatment for powers
     1281            # divisible by p
     1282            # Beware that p-th roots might not exist
     1283            unit = self.leading_coefficient()
     1284            T0 = self.monic()
     1285            e = 1
     1286            if T0.degree() > 0:
     1287                der = T0.derivative()
     1288                while der.is_zero():
     1289                    T0 = T0.parent()([T0[p*i].pth_root() for i in range(T0.degree()//p + 1)])
     1290                    if T0 == 1: raise RuntimeError
     1291                    der = T0.derivative()
     1292                    e = e*p
     1293                T = T0.gcd(der)
     1294                V = T0 // T
     1295                k = 0
     1296                while T0.degree() > 0:
     1297                    k += 1
     1298                    if p.divides(k):
     1299                        T = T // V
     1300                        k += 1
     1301                    W = V.gcd(T)
     1302                    if W.degree() < V.degree():
     1303                        factors.append((V // W, e*k))
     1304                        V = W
     1305                        T = T // V
     1306                        if V.degree() == 0:
     1307                            if T.degree() == 0:
     1308                                break
     1309                            # T is of the form sum_{i=0}^n t_i X^{pi}
     1310                            T0 = T0.parent()([T[p*i].pth_root() for i in range(T.degree()//p + 1)])
     1311                            der = T0.derivative()
     1312                            e = p*e
     1313                            while der.is_zero():
     1314                                T0 = T0.parent()([T0[p*i].pth_root() for i in range(T0.degree()//p + 1)])
     1315                                der = T0.derivative()
     1316                                e = p*e
     1317                            T = T0.gcd(der)
     1318                            V = T0 // T
     1319                            k = 0
     1320                    else:
     1321                        T = T//V
    12521322
    12531323        return Factorization(factors, unit=unit, sort=False)
    12541324
     
    52925362
    52935363        return self.roots(ring=CC, multiplicities=False)
    52945364
     5365    def any_root(self, ring=None, degree=None, assume_squarefree=False, extend=False):
     5366        """
     5367        Returns a root of this polynomial in the given ring.
     5368
     5369        INPUT:
     5370       
     5371        - ``ring`` -- The ring in which a root is sought.  By default this is the coefficient ring.
     5372
     5373        - ``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.
     5374
     5375        - ``assume_squarefree`` (bool) -- Used for polynomials over finite fields.  If True, this polynomial is assumed to be squarefree.
     5376
     5377        - ``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.
     5378
     5379        EXAMPLES::
     5380
     5381            sage: R.<x> = GF(11)[]
     5382            sage: f = 7*x^7 + 8*x^6 + 4*x^5 + x^4 + 6*x^3 + 10*x^2 + 8*x + 5
     5383            sage: f.any_root()
     5384            2
     5385            sage: f.factor()
     5386            (7) * (x + 9) * (x^6 + 10*x^4 + 6*x^3 + 5*x^2 + 2*x + 2)
     5387            sage: f = x^6 + 10*x^4 + 6*x^3 + 5*x^2 + 2*x + 2
     5388            sage: f.any_root(GF(11^6))
     5389            10*z6^5 + 3*z6^4 + 8*z6^3 + z6^2 + 10*z6 + 4
     5390            sage: sorted(f.roots(GF(11^6)))
     5391            [(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)]
     5392            sage: f.any_root(extend=True)
     5393            10*z6^5 + 3*z6^4 + 8*z6^3 + z6^2 + 10*z6 + 4
     5394
     5395            sage: g = (x-1)*(x^2 + 3*x + 9) * (x^5 + 5*x^4 + 8*x^3 + 5*x^2 + 3*x + 5)
     5396            sage: g.any_root(ring=GF(11^10), degree=1)
     5397            1
     5398            sage: g.any_root(ring=GF(11^10), degree=2)
     5399            6*z10^9 + 7*z10^7 + 7*z10^6 + 3*z10^5 + z10^2 + z10 + 3
     5400            sage: g.any_root(ring=GF(11^10), degree=5)
     5401            7*z10^9 + 2*z10^8 + 6*z10^7 + 3*z10^6 + 2*z10^5 + 7*z10^3 + 7*z10^2 + 2
     5402
     5403        TESTS::
     5404
     5405            sage: R.<x> = GF(5)[]
     5406            sage: K.<a> = GF(5^12)
     5407            sage: for _ in range(40):
     5408            ...       f = R.random_element(degree=4)
     5409            ...       assert f(f.any_root(K)) == 0
     5410        """
     5411        if self.base_ring().is_finite() and self.base_ring().is_field():
     5412            if self.degree() < 0:
     5413                return ring(0)
     5414            if self.degree() == 0:
     5415                raise ValueError, "no roots A %s"%self
     5416            if not assume_squarefree:
     5417                SFD = self.squarefree_decomposition()
     5418                SFD.sort()
     5419                for f, e in SFD:
     5420                    try:
     5421                        return f.any_root(ring, degree, True, extend)
     5422                    except ValueError:
     5423                        pass
     5424            if self.degree() == 1 and (degree is None or degree == 1):
     5425                if ring is None:
     5426                    return -self[0]/self[1]
     5427                else:
     5428                    return ring(-self[0]/self[1])
     5429            q = self.base_ring().order()
     5430            if ring is None:
     5431                allowed_deg_mult = Integer(1)
     5432            else:
     5433                if not self.base_ring().is_prime_field(): raise NotImplementedError
     5434                if ring.characteristic() != self.base_ring().characteristic(): raise ValueError, "ring must be an extension of the base ring"
     5435                if not (ring.is_field() and ring.is_finite()): raise NotImplementedError
     5436                allowed_deg_mult = Integer(ring.factored_order()[0][1]) # generally it will be the quotient of this by the degree of the base ring.
     5437            if degree is None:
     5438                x = self.parent().gen()
     5439                if allowed_deg_mult == 1 and not extend:
     5440                    self = self.gcd(x**q-x)
     5441                    degree = -1
     5442                    if self.degree() == 0: raise ValueError, "no roots B %s"%self
     5443                else:
     5444                    xq = x
     5445                    d = Integer(0)
     5446                    while True: # one pass for roots that actually lie within ring.
     5447                        e = self.degree()
     5448                        if 2*d+2 > e: # this polynomial has no factors dividing allowed_deg_mult
     5449                            if allowed_deg_mult % e == 0:
     5450                                degree = -e
     5451                            break
     5452                        while d < allowed_deg_mult:
     5453                            d = d+1
     5454                            xq = xq**q % self
     5455                            if d.divides(allowed_deg_mult): break
     5456                        A = self.gcd(xq-x)
     5457                        if A != 1:
     5458                            self = A
     5459                            degree = -d
     5460                            break
     5461                        if d == allowed_deg_mult:
     5462                            break
     5463                    if degree is None:
     5464                        if allowed_deg_mult==1 and not extend: raise ValueError, "no roots C %s"%self
     5465                        xq = x
     5466                        d = Integer(0)
     5467                        while True: # now another for roots that will lie in an extension.
     5468                            e = self.degree()
     5469                            if 2*d+2 > e: # this polynomial is irreducible.
     5470                                degree = -e
     5471                                break
     5472                            while True: # we waste a little effort here in computing the xq again.
     5473                                d = d+1
     5474                                xq = xq**q % self
     5475                                if allowed_deg_mult.divides(d): break
     5476                            A = self.gcd(xq-x)
     5477                            if A != 1:
     5478                                self = A
     5479                                degree = -d
     5480                                break                           
     5481            if degree == 0: raise ValueError, "degree should be nonzero"
     5482            R = self.parent()
     5483            x = R.gen()
     5484            if degree > 0:
     5485                xq = x
     5486                d = 0
     5487                while True:
     5488                    e = self.degree()
     5489                    if 2*d > e:
     5490                        if degree != e:
     5491                            raise ValueError, "no roots D %s"%self
     5492                        break
     5493                    d = d+1
     5494                    xq = xq**q % self
     5495                    if d == degree:
     5496                        break
     5497                    A = self.gcd(xq-x)
     5498                    if A != 1:
     5499                        self = self // A
     5500                if d == degree:
     5501                    self = self.gcd(xq-x)
     5502                    if self.degree() == 0:
     5503                        raise ValueError, "no roots E %s"%self
     5504            else:
     5505                degree = -degree
     5506            if ring is None:
     5507                if degree == 1:
     5508                    ring = self.base_ring()
     5509                else:
     5510                    ring = self.base_ring().extension(degree) # this won't work yet.
     5511            # now self has only roots of degree ``degree``.
     5512            # for now, we only implement the Cantor-Zassenhaus split
     5513            k = self.degree() // degree
     5514            if k == 1:
     5515                try:
     5516                    return self.roots(ring, multiplicities=False)[0] # is there something better to do here?
     5517                except IndexError:
     5518                    raise ValueError, "no roots F %s"%self
     5519            if q % 2 == 0:
     5520                T = x
     5521                while True:
     5522                    C = T
     5523                    for i in range(degree-1):
     5524                        C = T + C**q % self
     5525                    h = self.gcd(C)
     5526                    hd = h.degree()
     5527                    if hd == 0 or hd == self.degree():
     5528                        T = T * x**q
     5529                    else:
     5530                        if 2*hd <= self.degree():
     5531                            return h.any_root(ring, -degree, True, extend)
     5532                        else:
     5533                            return (self//h).any_root(ring, -degree, True, extend)
     5534            else:
     5535                from sage.rings.arith import power_mod
     5536                while True:
     5537                    T = R.random_element(2*degree-1)
     5538                    if T == 0: continue
     5539                    T = T.monic()
     5540                    h = self.gcd(power_mod(T, Integer((q**degree-1)/2), self)-1)
     5541                    hd = h.degree()
     5542                    if hd != 0 and hd != self.degree():
     5543                        if 2*hd <= self.degree():
     5544                            return h.any_root(ring, -degree, True, extend)
     5545                        else:
     5546                            return (self//h).any_root(ring, -degree, True, extend)
     5547        else:
     5548            return self.roots(ring=ring, multiplicities=False)[0]
     5549       
    52955550
    52965551    def variable_name(self):
    52975552        """
     
    57686023        # over rings of positive characteristic, we rely on the square-free decomposition if available
    57696024        try:
    57706025            F = self.squarefree_decomposition()
    5771         except NotImplementedError:
     6026        # We catch:
     6027        # - NotImplementedError in case squarefree decomposition is not implemented
     6028        # - AttributeError in case p-th roots are not (or do not exist)
     6029        except (NotImplementedError, AttributeError):
    57726030            F = self.factor()
    57736031        return all([e<=1 for (f,e) in F])
    57746032