Ticket #8335: trac_8335-pseudo_conway-fixes.patch

File trac_8335-pseudo_conway-fixes.patch, 48.1 KB (added by jpflori, 8 years ago)
  • sage/rings/finite_rings/constructor.py

    # HG changeset patch
    # User Jean-Pierre Flori <jean-pierre.flori@ssi.gouv.fr>
    # Date 1374228320 -7200
    #      Fri Jul 19 12:05:20 2013 +0200
    # Node ID 6ed8118fc53f0996e11274d8c3bb7a1b379c1bef
    # Parent  9ca55ab215a991709a147cf3e305d74a3a2fcf07
    bla
    
    diff --git a/sage/rings/finite_rings/constructor.py b/sage/rings/finite_rings/constructor.py
    a b  
    168168from finite_field_givaro import FiniteField_givaro
    169169
    170170import sage.interfaces.gap
    171 import sage.databases.conway
    172171
    173172from sage.structure.factory import UniqueFactory
    174173from sage.structure.sage_object import SageObject
     
    185184   
    186185    -  ``order`` - int
    187186   
    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.
     187    -  ``name`` - string; must be specified unless ``order`` is prime.
    192188   
    193189    -  ``modulus`` - (optional) either a defining polynomial for the
    194190       field, i.e., generator of the field will be a root of this
    195191       polynomial; or a string:
    196192
    197            - 'conway': force the use of a pseudo-Conway polynomial.
    198             If a conway polynomial is found in the database, that is used.
     193           - 'conway': use a Conway polynomial (if available).
     194           - 'pseudo-conway': If a conway polynomial is found in the database,
     195            that is used.
    199196            Otherwise, a polynomial is used with the same algebraic properties
    200197            but without the lexicographic constraint implying uniqueness.
    201198           - 'random': use a random irreducible polynomial;
     
    350347
    351348    """
    352349    def create_key_and_extra_args(self, order, name=None, modulus=None, names=None,
    353                                   impl=None, proof=None, **kwds):
     350                                  impl=None, prefix='z', proof=None, **kwds):
    354351        """
    355352        EXAMPLES::
    356353       
     
    372369                p = integer.Integer(order)
    373370                n = integer.Integer(1)
    374371            elif arith.is_prime_power(order):
    375                 if not names is None: name = names
    376                 name = normalize_names(1,name)
    377 
    378372                p, n = arith.factor(order)[0]
    379373
     374                # compute modulus with chosen algo
    380375                if modulus is None or isinstance(modulus, str):
    381376                    # A string specifies an algorithm to find a suitable modulus.
    382377                    if modulus == "default":    # for backward compatibility
     
    386381                    modulus = GF(p)['x'](modulus)
    387382                elif sage.rings.polynomial.polynomial_element.is_Polynomial(modulus):
    388383                    modulus = modulus.change_variable_name('x')
     384
     385                # tweak generator name
     386                if not names is None: name = names
     387
     388                if modulus == "pseudo-conway" and name is None:
     389                    name = prefix + str(n)
     390
     391                name = normalize_names(1,name)
     392
    389393                else:
    390394                    raise TypeError("wrong type for modulus parameter")
    391395            else:
     
    541545
    542546    return isinstance(x, FiniteField_prime_modn) or \
    543547           (isinstance(x, FiniteField_generic) and x.degree() == 1)
    544    
    545 ##################################################################
    546    
    547 def conway_polynomial(p, n):
    548     r"""
    549     Return the Conway polynomial of degree n over GF(p), which is
    550     loaded from a table.
    551    
    552     If the requested polynomial is not known, this function raises a
    553     RuntimeError exception.
    554    
    555     INPUT:
    556    
    557    
    558     -  ``p`` - int
    559    
    560     -  ``n`` - int
    561    
    562    
    563     OUTPUT:
    564    
    565    
    566     -  ``Polynomial`` - a polynomial over the prime finite
    567        field GF(p).
    568    
    569    
    570     .. note::
    571 
    572        The first time this function is called a table is read from
    573        disk, which takes a fraction of a second. Subsequent calls do
    574        not require reloading the table.
    575    
    576     See also the ``ConwayPolynomials()`` object, which is a
    577     table of Conway polynomials. For example, if
    578     ``c=ConwayPolynomials()``, then
    579     ``c.primes()`` is a list of all primes for which the
    580     polynomials are known, and for a given prime `p`,
    581     ``c.degree(p)`` is a list of all degrees for which the
    582     Conway polynomials are known.
    583    
    584     EXAMPLES::
    585    
    586         sage: conway_polynomial(2,5)
    587         x^5 + x^2 + 1
    588         sage: conway_polynomial(101,5)
    589         x^5 + 2*x + 99
    590         sage: conway_polynomial(97,101)
    591         Traceback (most recent call last):
    592         ...
    593         RuntimeError: requested conway polynomial not in database.
    594     """
    595     (p,n)=(int(p),int(n))
    596     R = FiniteField(p)['x']
    597     try:
    598         return R(sage.databases.conway.ConwayPolynomials()[p][n])
    599     except KeyError:
    600         raise RuntimeError("requested conway polynomial not in database.")
    601 
    602 def exists_conway_polynomial(p, n):
    603     r"""
    604     Return True if the Conway polynomial over `F_p` of degree
    605     `n` is in the database and False otherwise.
    606    
    607     If the Conway polynomial is in the database, to obtain it use the
    608     command ``conway_polynomial(p,n)``.
    609    
    610     EXAMPLES::
    611    
    612         sage: exists_conway_polynomial(2,3)
    613         True
    614         sage: exists_conway_polynomial(2,-1)
    615         False
    616         sage: exists_conway_polynomial(97,200)
    617         False
    618         sage: exists_conway_polynomial(6,6)
    619         False
    620     """
    621     return sage.databases.conway.ConwayPolynomials().has_polynomial(p,n)
    622 
    623 pseudo_conway_poly = {}
    624 
    625 class PseudoConwayPolyTree(SageObject):
    626     """
    627     An object holding references to pseudo-Conway polynomials for divisors of the given degree, so that they aren't garbage collected.
    628     """
    629     def __init__(self, p, n, nodes_dict, f):
    630         """
    631         INPUT::
    632 
    633         - p -- The prime for which this defines an extension of GF(p).
    634         - n -- The degree of the extension.
    635         - 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.
    636                         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).
    637         - f -- The polynomial defining this extension.
    638 
    639         EXAMPLES::
    640 
    641             sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
    642             sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
    643             sage: PCPT.get_pseudo_conway_poly(3)
    644             x^3 + x + 1
    645         """
    646         self.p = p
    647         self.n = n
    648         self.nodes_dict = nodes_dict
    649         self.f = f
    650        
    651     def get_pseudo_conway_poly(self, m):
    652         """
    653         Returns the pseudo-Conway polynomial associated to a divisor of the degree of this tree.
    654 
    655         EXAMPLES::
    656        
    657             sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
    658             sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
    659             sage: PCPT.get_pseudo_conway_poly(3)
    660             x^3 + x + 1
    661             sage: PCPT.get_pseudo_conway_poly(4)
    662             Traceback (most recent call last):
    663             ...
    664             RuntimeError: 4 does not divide 6
    665         """
    666         if m == self.n:
    667             return self.f
    668         if not m.divides(self.n):
    669             raise RuntimeError, "%s does not divide %s"%(m, self.n)
    670         try:
    671             return pseudo_conway_poly[self.p][m]().f
    672         except (KeyError, AttributeError):
    673             return conway_polynomial(self.p, m)
    674            
    675     def check_consistency(self):
    676         """
    677         Checks that the pseudo-Conway polynomials dividing the degree of this tree satisfy the required compatibility conditions.
    678 
    679         EXAMPLES::
    680 
    681             sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
    682             sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
    683             sage: PCPT.check_consistency()
    684             sage: PCPT = find_pseudo_conway_polynomial_tree(2, 60, False) # long
    685             sage: PCPT.check_consistency() # long
    686         """
    687         p, n = self.p, self.n
    688         K = GF(p**n, modulus = self.f, names='a')
    689         a = K.gen()
    690         for m in n.divisors():
    691             assert (a**((p**n-1)//(p**m-1))).minimal_polynomial() == self.get_pseudo_conway_poly(m)
    692 
    693 def find_pseudo_conway_polynomial_tree(p, n, use_database=True):
    694     """
    695     Returns an object holding references to a set of consistent pseudo-Conway polynomials for degrees dividing n.
    696 
    697     A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
    698    
    699     - `f` is irreducible.
    700     - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
    701     - 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`.
    702     - `f` is lexicographically least among all such polynomials, under a certain ordering.
    703 
    704     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.
    705 
    706     INPUT:
    707 
    708     - `p` -- a prime.
    709     - `n` -- an integer greater than 1.
    710     - `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)
    711 
    712     EXAMPLES::
    713 
    714         sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
    715         sage: PCPT = find_pseudo_conway_polynomial_tree(2, 12, False)
    716         sage: PCPT.f
    717         x^12 + x^8 + x^7 + x^6 + x^4 + x^3 + 1
    718     """
    719     if not pseudo_conway_poly.has_key(p):
    720         pseudo_conway_poly[p] = {}
    721     pdict = pseudo_conway_poly[p]
    722     n = sage.rings.integer.Integer(n)
    723     if pdict.has_key(n):
    724         pcp = pdict[n]()
    725         if pcp is not None:
    726             return pcp
    727     if use_database and exists_conway_polynomial(p, n):
    728         ans = PseudoConwayPolyTree(p, n, True, GF(p)['x'](conway_polynomial(p, n)))
    729         pdict[n] = weakref.ref(ans)
    730         return ans
    731     if n == 1 or n.is_prime():
    732         ans = PseudoConwayPolyTree(p, n, False, compute_pseudo_conway_polynomial(p, n, None))
    733     nodes = {}
    734     for q in n.prime_factors():
    735         nodes[q] = find_pseudo_conway_polynomial_tree(p, n//q, use_database)
    736     ans = PseudoConwayPolyTree(p, n, nodes, compute_pseudo_conway_polynomial(p, n, nodes))
    737     pdict[n] = weakref.ref(ans)
    738     return ans
    739 
    740 def _find_pow_of_frobenius(p, x, y, field_degree):
    741     """
    742     Finds the power of Frobenius which yields `x` when applied to `y`.
    743 
    744     INPUT:
    745 
    746     - `p` -- a prime.
    747     - ``field_degree`` -- a positive integer
    748     - `x` -- an element of a field `K` of characeristic p so that the multiplicative order of `x` is ``p^field_degree``.
    749     - `y` -- an element of `K` with the same minimal polynomial.
    750    
    751     OUTPUT:
    752    
    753     An element `i` of ``Integers(field_degree)`` so that `x = y^{p^i}`
    754 
    755     EXAMPLES::
    756 
    757         sage: from sage.rings.finite_rings.constructor import _find_pow_of_frobenius
    758         sage: K.<a> = GF(3^14)
    759         sage: x = K.multiplicative_generator()
    760         sage: y = x^27
    761         sage: _find_pow_of_frobenius(3, x, y, 14)
    762         11
    763     """
    764     assert x.minimal_polynomial() == y.minimal_polynomial()
    765     assert x.multiplicative_order() == p**field_degree - 1
    766     assert y.multiplicative_order() == p**field_degree - 1
    767     from integer_mod import mod
    768     for i in range(field_degree):
    769         if x == y: break
    770         y = y**p
    771     else:
    772         raise RuntimeError, "No appropriate power of Frobenius found"
    773     return mod(i, field_degree)
    774 
    775 def _crt_non_coprime(running, a):
    776     """
    777     Extends the ``crt`` method of IntegerMod to the case of non-relatively prime modulus.
    778 
    779     EXAMPLES::
    780 
    781         sage: from sage.rings.finite_rings.constructor import _crt_non_coprime
    782         sage: a = _crt_non_coprime(mod(14, 18), mod(20,30)); a
    783         50
    784         sage: a.modulus()
    785         90
    786         sage: _crt_non_coprime(mod(13, 18), mod(20,30))
    787         Traceback (most recent call last):
    788         ...
    789         AssertionError
    790     """
    791     g = running.modulus().gcd(a.modulus())
    792     if g == 1:
    793         return running.crt(a)
    794     else:
    795         assert running % g == a % g
    796         running_modulus = running.modulus()
    797         a_modulus = a.modulus()
    798         for qq in g.prime_divisors():
    799             a_val_unit = a_modulus.val_unit(qq)
    800             running_val_unit = running_modulus.val_unit(qq)
    801             if a_val_unit[0] > running_val_unit[0]:
    802                 running_modulus = running_val_unit[1]
    803             else:
    804                 a_modulus = a_val_unit[1]
    805         return (running % running_modulus).crt(a % a_modulus)
    806 
    807 def _frobenius_shift(K, generators, check_only=False):
    808     """
    809     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.
    810 
    811     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
    812 
    813     ``generators[q1]^((p^(n/q1)-1)/(p^m-1)) == generators[q2]^((p^(n/q2)-1)/(p^m-1))``
    814 
    815     holds.
    816 
    817     INPUT:
    818 
    819     - K -- a finite field of degree n   
    820     - generators -- a dictionary, indexed by prime divisors `q` of `n`, whose entries are elements of `K` satisfying the `n/q` pseudo-Conway polynomial.
    821     - check_only -- if True, just asserts that the generators given form a compatible system.
    822 
    823     EXAMPLES::
    824 
    825         sage: R.<x> = GF(2)[]
    826         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
    827         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
    828         sage: f12 = x^12 + x^10 + x^9 + x^8 + x^4 + x^2 + 1
    829         sage: K.<a> = GF(2^60, modulus='first_lexicographic')
    830         sage: x30 = f30.any_root(K)
    831         sage: x20 = f20.any_root(K)
    832         sage: x12 = f12.any_root(K)
    833         sage: generators = {2: x30, 3: x20, 5: x12}
    834         sage: from sage.rings.finite_rings.constructor import _frobenius_shift, _find_pow_of_frobenius
    835         sage: _frobenius_shift(K, generators)
    836         sage: _find_pow_of_frobenius(2, x30, generators[2], 30)
    837         0
    838         sage: _find_pow_of_frobenius(2, x20, generators[3], 20)
    839         13
    840         sage: _find_pow_of_frobenius(2, x12, generators[5], 12)
    841         8
    842     """
    843     p = K.characteristic()
    844     n = K.degree()
    845     if len(generators) == 1:
    846         return generators
    847     compatible = {}
    848     from integer_mod import mod
    849     for m in n.divisors():
    850         compatible[m] = {}
    851     for q, x in generators.iteritems():
    852         for m in (n//q).divisors():
    853             compatible[m][q] = x**((p**(n//q)-1)//(p**m-1))
    854     if check_only:
    855         for m in n.divisors():
    856             try:
    857                 q, x = compatible[m].popitem()
    858             except KeyError:
    859                 break
    860             for qq, xx in compatible[m].iteritems():
    861                 assert x == xx
    862         return
    863     crt = {}
    864     qlist = sorted(generators.keys())
    865     for j in range(1, len(qlist)):
    866         for i in range(j):
    867             crt[(i, j)] = []
    868     for m in n.divisors():
    869         mqlist = sorted(compatible[m].keys())
    870         for k in range(1,len(mqlist)):
    871             j = qlist.index(mqlist[k])
    872             i = qlist.index(mqlist[k-1])
    873             crt[(i,j)].append(_find_pow_of_frobenius(p, compatible[m][qlist[j]], compatible[m][qlist[i]], m))
    874     from integer_mod import mod
    875     pairs = crt.keys()
    876     for i, j in pairs:
    877         L = crt[(i,j)]
    878         #print qlist[i], qlist[j]
    879         #print " ".join(["(%s, %s)"%(mm, mm.modulus()) for mm in L])
    880         running = mod(0,1)
    881         for a in L:
    882             running = _crt_non_coprime(running, a)
    883         #print "(%s, %s)"%(running, running.modulus())
    884         crt[(i,j)] = [(mod(running, q**(running.modulus().valuation(q))), running.modulus().valuation(q)) for q in qlist]
    885         crt[(j,i)] = [(-a, level) for a, level in crt[(i,j)]]
    886     # Let x_j be the power of Frobenius we apply to generators[qlist[j]], for 0 < j < len(qlist)
    887     # We have some direct conditions on the x_j: x_j reduces to each entry in crt[(0,j)].
    888     # But we also have the equations x_j - x_i reduces to each entry in crt[(i,j)].
    889     # We solve for x_j one prime at a time.  For each prime, we have an equations of the form
    890     # 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
    891     # (possibly 0, possibly different) of the same prime.
    892 
    893     # We can set x_0=0 everywhere, can get an initial setting of x_j from the c_0j.
    894     # We go through prime by prime.
    895     import bisect
    896     frob_powers=[mod(0,1) for q in qlist]
    897     def find_leveller(qindex, level, x, xleveled, searched, i):
    898         searched[i] = True
    899         crt_possibles = []
    900         for j in range(1,len(qlist)):
    901             if i==j: continue
    902             if crt[(i,j)][qindex][1] >= level:
    903                 if xleveled[j]:
    904                     return [j]
    905                 elif not searched.has_key(j):
    906                     crt_possibles.append(j)
    907         for j in crt_possibles:
    908             path = find_leveller(qindex, level, x, xleveled, searched, j)
    909             if path is not None:
    910                 path.append(j)
    911                 return path
    912         return None
    913     def propogate_levelling(qindex, level, x, xleveled, i):
    914         for j in range(1, len(qlist)):
    915             if i==j: continue
    916             if not xleveled[j] and crt[(i,j)][qindex][1] >= level:
    917                 newxj = x[i][0] + crt[(i,j)][qindex][0]
    918                 x[j] = (newxj, min(x[i][1], crt[(i,j)][qindex][1]))
    919                 xleveled[j] = True
    920                 propogate_levelling(qindex, level, x, xleveled, j)
    921                
    922     for qindex in range(len(qlist)):
    923         q = qlist[qindex]
    924         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.               
    925         # We first check that our equations are consistent and determine which powers of q occur as moduli.
    926         levels = []
    927         for j in range(2, len(qlist)):
    928             for i in range(j):
    929                 # we need crt[(0,j)] = crt[(0,i)] + crt[(i,j)]
    930                 if i != 0: assert x[j][0] == x[i][0] + crt[(i,j)][qindex][0]
    931                 level = crt[(i,j)][qindex][1]
    932                 if level > 0:
    933                     ins = bisect.bisect_left(levels,level)
    934                     if ins == len(levels):
    935                         levels.append(level)
    936                     elif levels[ins] != level:
    937                         levels.insert(ins, level)
    938         for level in levels:
    939             xleveled = [0] + [x[i][1] >= level for i in range(1,len(qlist))]
    940             while True:
    941                 try:
    942                     i = xleveled.index(False, 1)
    943                     searched = {}
    944                     levelling_path = find_leveller(qindex, level, x, xleveled, searched, i)
    945                     if levelling_path is None:
    946                         # Any lift will work, since there are no constraints
    947                         x[i] = (mod(x[i][0].lift(), q**level), level)
    948                         xleveled[i] = True
    949                         propogate_levelling(qindex, level, x, xleveled, i)
    950                     else:
    951                         levelling_path.append(i)
    952                         for m in range(1,len(path)):
    953                             # this point on the path may have already been leveled in a previous propogation
    954                             if not xleveled[path[m]]:
    955                                 newx = x[path[m-1]][0] + crt[(path[m-1],path[m])][qindex][0]
    956                                 x[path[m]] = (newx, min(x[path[m-1]][1], crt[(path[m-1],path[m])][qindex][1]))
    957                                 xleveled[path[m]] = True
    958                                 propogate_levelling(qindex, level, x, xleveled, path[m])
    959                 except ValueError:
    960                     break
    961         for j in range(1,len(qlist)):
    962             frob_powers[j] = frob_powers[j].crt(x[j][0])
    963     for j in range(1, len(qlist)):
    964         #print "frob_power[%s] = mod(%s, %s)"%(qlist[j], frob_powers[j], frob_powers[j].modulus())
    965         generators[qlist[j]] = generators[qlist[j]]**(p**(-frob_powers[j]).lift())
    966     _frobenius_shift(K, generators, check_only=True)
    967 
    968 def compute_pseudo_conway_polynomial(p, n, nodes):
    969     """
    970     Computes a pseudo-Conway polynomial over the prime `p` of degree `n`.
    971 
    972     A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
    973    
    974     - `f` is irreducible.
    975     - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
    976     - 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`.
    977     - `f` is lexicographically least among all such polynomials, under a certain ordering.
    978 
    979     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.
    980 
    981     INPUT:
    982    
    983     - `p` -- a prime
    984     - `n` -- a positive integer
    985     - ``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`.
    986 
    987     OUTPUT:
    988 
    989     A pseudo-Conway polynomial over the prime `p` of degree `n`.
    990 
    991     ALGORITHM:
    992    
    993     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.
    994 
    995     REFERENCE:
    996 
    997     .. [HL99] Heath L. and Loehr N. (1999).
    998        New algorithms for generating Conway polynomials over finite fields.
    999        Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms, pp. 429-437.
    1000 
    1001     EXAMPLES::
    1002 
    1003         sage: from sage.rings.finite_rings.constructor import compute_pseudo_conway_polynomial, find_pseudo_conway_polynomial_tree
    1004         sage: compute_pseudo_conway_polynomial(next_prime(10000),11,None) # long because the arithmetic in FiniteField_ext_pari fields is slow.
    1005         x^11 + x + 7
    1006         sage: PCPT30 = find_pseudo_conway_polynomial_tree(2, 30, False)
    1007         sage: PCPT20 = find_pseudo_conway_polynomial_tree(2, 20, False)
    1008         sage: PCPT12 = find_pseudo_conway_polynomial_tree(2, 12, False)
    1009         sage: compute_pseudo_conway_polynomial(2, 60, {2:PCPT30,3:PCPT20,5:PCPT12})
    1010         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
    1011     """
    1012     k = GF(p)
    1013     from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    1014     R = PolynomialRing(k, names='x')
    1015     if n == 1:
    1016         return R([-k.multiplicative_generator(), 1])
    1017     if nodes is None:
    1018         # n is prime.  We require only that the chosen generator (x) be a multiplicative generator for the units.
    1019         if p in [2,3,5,7,11]: # here we have Cunningham data
    1020             from sage.rings.factorint import factor_cunningham
    1021             F = [(p**n-1)//a for a, e in factor_cunningham(n)]
    1022         else:
    1023             F = (p**n-1).prime_factors() # could be VERY expensive
    1024         x = R.gen()
    1025         from sage.rings.arith import power_mod
    1026         for f in R.polynomials(of_degree=n):
    1027             # By using this iterator our chosen polynomial will be fairly sparse.
    1028             if f.is_irreducible():
    1029                 if not f.is_monic():
    1030                     f = f.monic()
    1031                 for m in F:
    1032                     if power_mod(x, m, f) == 1:
    1033                         break
    1034                 else:
    1035                     return f
    1036     # Our strategy is as follows:
    1037     # Work in an arbitrary field K of order p**n.
    1038     # We're looking for a multiplicative generator of K whose
    1039     K = GF(p**n, modulus="first_lexicographic", names='a')
    1040     r = p**n - 1
    1041     xi = {}
    1042     for q in nodes.iterkeys():
    1043         #print nodes[q].f.degree(), n//q
    1044         xi[q] = nodes[q].f.any_root(K, -n//q, True)
    1045     _frobenius_shift(K, xi)
    1046            
    1047     q, x = xi.popitem()
    1048     v = p**(n//q) - 1
    1049     for q, xitem in xi.iteritems():
    1050         fi = p**(n//q) - 1
    1051         g, alpha, beta = v.xgcd(fi)
    1052         x = x**beta * xitem**alpha
    1053         v = v.lcm(fi)
    1054     g = r // v
    1055     #print g, n, g.divides(p**n-1)
    1056     z = x.nth_root(g, cunningham=(p in [2,3,5,7,11]))
    1057     groot1 = K.multiplicative_generator()**v
    1058     while z.multiplicative_order() != r:
    1059         z *= groot1
    1060     #print z
    1061     return z.minimal_polynomial()
    1062548
    1063549zech_log_bound = 2**16
  • new file sage/rings/finite_rings/conway_polynomials.py

    diff --git a/sage/rings/finite_rings/conway_polynomials.py b/sage/rings/finite_rings/conway_polynomials.py
    new file mode 100644
    - +  
     1"""
     2Routines for Conway and pseudo-Conway polynomials.
     3
     4AUTHORS:
     5
     6- David Roe
     7
     8- Jean-Pierre Flori
     9
     10"""
     11import sage.databases.conway
     12
     13def conway_polynomial(p, n):
     14    r"""
     15    Return the Conway polynomial of degree n over GF(p), which is
     16    loaded from a table.
     17   
     18    If the requested polynomial is not known, this function raises a
     19    RuntimeError exception.
     20   
     21    INPUT:
     22   
     23   
     24    -  ``p`` - int
     25   
     26    -  ``n`` - int
     27   
     28   
     29    OUTPUT:
     30   
     31   
     32    -  ``Polynomial`` - a polynomial over the prime finite
     33       field GF(p).
     34   
     35   
     36    .. note::
     37
     38       The first time this function is called a table is read from
     39       disk, which takes a fraction of a second. Subsequent calls do
     40       not require reloading the table.
     41   
     42    See also the ``ConwayPolynomials()`` object, which is a
     43    table of Conway polynomials. For example, if
     44    ``c=ConwayPolynomials()``, then
     45    ``c.primes()`` is a list of all primes for which the
     46    polynomials are known, and for a given prime `p`,
     47    ``c.degree(p)`` is a list of all degrees for which the
     48    Conway polynomials are known.
     49   
     50    EXAMPLES::
     51   
     52        sage: conway_polynomial(2,5)
     53        x^5 + x^2 + 1
     54        sage: conway_polynomial(101,5)
     55        x^5 + 2*x + 99
     56        sage: conway_polynomial(97,101)
     57        Traceback (most recent call last):
     58        ...
     59        RuntimeError: requested conway polynomial not in database.
     60    """
     61    (p,n)=(int(p),int(n))
     62    R = FiniteField(p)['x']
     63    try:
     64        return R(sage.databases.conway.ConwayPolynomials()[p][n])
     65    except KeyError:
     66        raise RuntimeError("requested conway polynomial not in database.")
     67
     68def exists_conway_polynomial(p, n):
     69    r"""
     70    Return True if the Conway polynomial over `F_p` of degree
     71    `n` is in the database and False otherwise.
     72   
     73    If the Conway polynomial is in the database, to obtain it use the
     74    command ``conway_polynomial(p,n)``.
     75   
     76    EXAMPLES::
     77   
     78        sage: exists_conway_polynomial(2,3)
     79        True
     80        sage: exists_conway_polynomial(2,-1)
     81        False
     82        sage: exists_conway_polynomial(97,200)
     83        False
     84        sage: exists_conway_polynomial(6,6)
     85        False
     86    """
     87    return sage.databases.conway.ConwayPolynomials().has_polynomial(p,n)
     88
     89pseudo_conway_poly = {}
     90
     91class PseudoConwayPolyTree(SageObject):
     92    """
     93    An object holding references to pseudo-Conway polynomials for divisors of the given degree, so that they aren't garbage collected.
     94    """
     95    def __init__(self, p, n, nodes_dict, f):
     96        """
     97        INPUT::
     98
     99        - p -- The prime for which this defines an extension of GF(p).
     100        - n -- The degree of the extension.
     101        - 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.
     102                        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).
     103        - f -- The polynomial defining this extension.
     104
     105        EXAMPLES::
     106
     107            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     108            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     109            sage: PCPT.get_pseudo_conway_poly(3)
     110            x^3 + x + 1
     111        """
     112        self.p = p
     113        self.n = n
     114        self.nodes_dict = nodes_dict
     115        self.f = f
     116       
     117    def get_pseudo_conway_poly(self, m):
     118        """
     119        Returns the pseudo-Conway polynomial associated to a divisor of the degree of this tree.
     120
     121        EXAMPLES::
     122       
     123            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     124            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     125            sage: PCPT.get_pseudo_conway_poly(3)
     126            x^3 + x + 1
     127            sage: PCPT.get_pseudo_conway_poly(4)
     128            Traceback (most recent call last):
     129            ...
     130            RuntimeError: 4 does not divide 6
     131        """
     132        if m == self.n:
     133            return self.f
     134        if not m.divides(self.n):
     135            raise RuntimeError, "%s does not divide %s"%(m, self.n)
     136        try:
     137            return pseudo_conway_poly[self.p][m]().f
     138        except (KeyError, AttributeError):
     139            return conway_polynomial(self.p, m)
     140           
     141    def check_consistency(self):
     142        """
     143        Checks that the pseudo-Conway polynomials dividing the degree of this tree satisfy the required compatibility conditions.
     144
     145        EXAMPLES::
     146
     147            sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     148            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 6, False)
     149            sage: PCPT.check_consistency()
     150            sage: PCPT = find_pseudo_conway_polynomial_tree(2, 60, False) # long
     151            sage: PCPT.check_consistency() # long
     152        """
     153        p, n = self.p, self.n
     154        K = GF(p**n, modulus = self.f, names='a')
     155        a = K.gen()
     156        for m in n.divisors():
     157            assert (a**((p**n-1)//(p**m-1))).minimal_polynomial() == self.get_pseudo_conway_poly(m)
     158
     159def find_pseudo_conway_polynomial_tree(p, n, use_database=True):
     160    """
     161    Returns an object holding references to a set of consistent pseudo-Conway polynomials for degrees dividing n.
     162
     163    A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
     164   
     165    - `f` is irreducible.
     166    - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
     167    - 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`.
     168    - `f` is lexicographically least among all such polynomials, under a certain ordering.
     169
     170    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.
     171
     172    INPUT:
     173
     174    - `p` -- a prime.
     175    - `n` -- an integer greater than 1.
     176    - `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)
     177
     178    EXAMPLES::
     179
     180        sage: from sage.rings.finite_rings.constructor import find_pseudo_conway_polynomial_tree
     181        sage: PCPT = find_pseudo_conway_polynomial_tree(2, 12, False)
     182        sage: PCPT.f
     183        x^12 + x^8 + x^7 + x^6 + x^4 + x^3 + 1
     184    """
     185    if not pseudo_conway_poly.has_key(p):
     186        pseudo_conway_poly[p] = {}
     187    pdict = pseudo_conway_poly[p]
     188    n = sage.rings.integer.Integer(n)
     189    if pdict.has_key(n):
     190        pcp = pdict[n]()
     191        if pcp is not None:
     192            return pcp
     193    if use_database and exists_conway_polynomial(p, n):
     194        ans = PseudoConwayPolyTree(p, n, True, GF(p)['x'](conway_polynomial(p, n)))
     195        pdict[n] = weakref.ref(ans)
     196        return ans
     197    if n == 1 or n.is_prime():
     198        ans = PseudoConwayPolyTree(p, n, False, compute_pseudo_conway_polynomial(p, n, None))
     199    nodes = {}
     200    for q in n.prime_factors():
     201        nodes[q] = find_pseudo_conway_polynomial_tree(p, n//q, use_database)
     202    ans = PseudoConwayPolyTree(p, n, nodes, compute_pseudo_conway_polynomial(p, n, nodes))
     203    pdict[n] = weakref.ref(ans)
     204    return ans
     205
     206def _find_pow_of_frobenius(p, x, y, field_degree):
     207    """
     208    Finds the power of Frobenius which yields `x` when applied to `y`.
     209
     210    INPUT:
     211
     212    - `p` -- a prime.
     213    - ``field_degree`` -- a positive integer
     214    - `x` -- an element of a field `K` of characeristic p so that the multiplicative order of `x` is ``p^field_degree``.
     215    - `y` -- an element of `K` with the same minimal polynomial.
     216   
     217    OUTPUT:
     218   
     219    An element `i` of ``Integers(field_degree)`` so that `x = y^{p^i}`
     220
     221    EXAMPLES::
     222
     223        sage: from sage.rings.finite_rings.constructor import _find_pow_of_frobenius
     224        sage: K.<a> = GF(3^14)
     225        sage: x = K.multiplicative_generator()
     226        sage: y = x^27
     227        sage: _find_pow_of_frobenius(3, x, y, 14)
     228        11
     229    """
     230    assert x.minimal_polynomial() == y.minimal_polynomial()
     231    assert x.multiplicative_order() == p**field_degree - 1
     232    assert y.multiplicative_order() == p**field_degree - 1
     233    from integer_mod import mod
     234    for i in range(field_degree):
     235        if x == y: break
     236        y = y**p
     237    else:
     238        raise RuntimeError, "No appropriate power of Frobenius found"
     239    return mod(i, field_degree)
     240
     241def _crt_non_coprime(running, a):
     242    """
     243    Extends the ``crt`` method of IntegerMod to the case of non-relatively prime modulus.
     244
     245    EXAMPLES::
     246
     247        sage: from sage.rings.finite_rings.constructor import _crt_non_coprime
     248        sage: a = _crt_non_coprime(mod(14, 18), mod(20,30)); a
     249        50
     250        sage: a.modulus()
     251        90
     252        sage: _crt_non_coprime(mod(13, 18), mod(20,30))
     253        Traceback (most recent call last):
     254        ...
     255        AssertionError
     256    """
     257    g = running.modulus().gcd(a.modulus())
     258    if g == 1:
     259        return running.crt(a)
     260    else:
     261        assert running % g == a % g
     262        running_modulus = running.modulus()
     263        a_modulus = a.modulus()
     264        for qq in g.prime_divisors():
     265            a_val_unit = a_modulus.val_unit(qq)
     266            running_val_unit = running_modulus.val_unit(qq)
     267            if a_val_unit[0] > running_val_unit[0]:
     268                running_modulus = running_val_unit[1]
     269            else:
     270                a_modulus = a_val_unit[1]
     271        return (running % running_modulus).crt(a % a_modulus)
     272
     273def _frobenius_shift(K, generators, check_only=False):
     274    """
     275    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.
     276
     277    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
     278
     279    ``generators[q1]^((p^(n/q1)-1)/(p^m-1)) == generators[q2]^((p^(n/q2)-1)/(p^m-1))``
     280
     281    holds.
     282
     283    INPUT:
     284
     285    - K -- a finite field of degree n   
     286    - generators -- a dictionary, indexed by prime divisors `q` of `n`, whose entries are elements of `K` satisfying the `n/q` pseudo-Conway polynomial.
     287    - check_only -- if True, just asserts that the generators given form a compatible system.
     288
     289    EXAMPLES::
     290
     291        sage: R.<x> = GF(2)[]
     292        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
     293        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
     294        sage: f12 = x^12 + x^10 + x^9 + x^8 + x^4 + x^2 + 1
     295        sage: K.<a> = GF(2^60, modulus='first_lexicographic')
     296        sage: x30 = f30.any_root(K)
     297        sage: x20 = f20.any_root(K)
     298        sage: x12 = f12.any_root(K)
     299        sage: generators = {2: x30, 3: x20, 5: x12}
     300        sage: from sage.rings.finite_rings.constructor import _frobenius_shift, _find_pow_of_frobenius
     301        sage: _frobenius_shift(K, generators)
     302        sage: _find_pow_of_frobenius(2, x30, generators[2], 30)
     303        0
     304        sage: _find_pow_of_frobenius(2, x20, generators[3], 20)
     305        13
     306        sage: _find_pow_of_frobenius(2, x12, generators[5], 12)
     307        8
     308    """
     309    p = K.characteristic()
     310    n = K.degree()
     311    if len(generators) == 1:
     312        return generators
     313    compatible = {}
     314    from integer_mod import mod
     315    for m in n.divisors():
     316        compatible[m] = {}
     317    for q, x in generators.iteritems():
     318        for m in (n//q).divisors():
     319            compatible[m][q] = x**((p**(n//q)-1)//(p**m-1))
     320    if check_only:
     321        for m in n.divisors():
     322            try:
     323                q, x = compatible[m].popitem()
     324            except KeyError:
     325                break
     326            for qq, xx in compatible[m].iteritems():
     327                assert x == xx
     328        return
     329    crt = {}
     330    qlist = sorted(generators.keys())
     331    for j in range(1, len(qlist)):
     332        for i in range(j):
     333            crt[(i, j)] = []
     334    for m in n.divisors():
     335        mqlist = sorted(compatible[m].keys())
     336        for k in range(1,len(mqlist)):
     337            j = qlist.index(mqlist[k])
     338            i = qlist.index(mqlist[k-1])
     339            crt[(i,j)].append(_find_pow_of_frobenius(p, compatible[m][qlist[j]], compatible[m][qlist[i]], m))
     340    from integer_mod import mod
     341    pairs = crt.keys()
     342    for i, j in pairs:
     343        L = crt[(i,j)]
     344        #print qlist[i], qlist[j]
     345        #print " ".join(["(%s, %s)"%(mm, mm.modulus()) for mm in L])
     346        running = mod(0,1)
     347        for a in L:
     348            running = _crt_non_coprime(running, a)
     349        #print "(%s, %s)"%(running, running.modulus())
     350        crt[(i,j)] = [(mod(running, q**(running.modulus().valuation(q))), running.modulus().valuation(q)) for q in qlist]
     351        crt[(j,i)] = [(-a, level) for a, level in crt[(i,j)]]
     352    # Let x_j be the power of Frobenius we apply to generators[qlist[j]], for 0 < j < len(qlist)
     353    # We have some direct conditions on the x_j: x_j reduces to each entry in crt[(0,j)].
     354    # But we also have the equations x_j - x_i reduces to each entry in crt[(i,j)].
     355    # We solve for x_j one prime at a time.  For each prime, we have an equations of the form
     356    # 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
     357    # (possibly 0, possibly different) of the same prime.
     358
     359    # We can set x_0=0 everywhere, can get an initial setting of x_j from the c_0j.
     360    # We go through prime by prime.
     361    import bisect
     362    frob_powers=[mod(0,1) for q in qlist]
     363    def find_leveller(qindex, level, x, xleveled, searched, i):
     364        searched[i] = True
     365        crt_possibles = []
     366        for j in range(1,len(qlist)):
     367            if i==j: continue
     368            if crt[(i,j)][qindex][1] >= level:
     369                if xleveled[j]:
     370                    return [j]
     371                elif not searched.has_key(j):
     372                    crt_possibles.append(j)
     373        for j in crt_possibles:
     374            path = find_leveller(qindex, level, x, xleveled, searched, j)
     375            if path is not None:
     376                path.append(j)
     377                return path
     378        return None
     379    def propogate_levelling(qindex, level, x, xleveled, i):
     380        for j in range(1, len(qlist)):
     381            if i==j: continue
     382            if not xleveled[j] and crt[(i,j)][qindex][1] >= level:
     383                newxj = x[i][0] + crt[(i,j)][qindex][0]
     384                x[j] = (newxj, min(x[i][1], crt[(i,j)][qindex][1]))
     385                xleveled[j] = True
     386                propogate_levelling(qindex, level, x, xleveled, j)
     387               
     388    for qindex in range(len(qlist)):
     389        q = qlist[qindex]
     390        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.               
     391        # We first check that our equations are consistent and determine which powers of q occur as moduli.
     392        levels = []
     393        for j in range(2, len(qlist)):
     394            for i in range(j):
     395                # we need crt[(0,j)] = crt[(0,i)] + crt[(i,j)]
     396                if i != 0: assert x[j][0] == x[i][0] + crt[(i,j)][qindex][0]
     397                level = crt[(i,j)][qindex][1]
     398                if level > 0:
     399                    ins = bisect.bisect_left(levels,level)
     400                    if ins == len(levels):
     401                        levels.append(level)
     402                    elif levels[ins] != level:
     403                        levels.insert(ins, level)
     404        for level in levels:
     405            xleveled = [0] + [x[i][1] >= level for i in range(1,len(qlist))]
     406            while True:
     407                try:
     408                    i = xleveled.index(False, 1)
     409                    searched = {}
     410                    levelling_path = find_leveller(qindex, level, x, xleveled, searched, i)
     411                    if levelling_path is None:
     412                        # Any lift will work, since there are no constraints
     413                        x[i] = (mod(x[i][0].lift(), q**level), level)
     414                        xleveled[i] = True
     415                        propogate_levelling(qindex, level, x, xleveled, i)
     416                    else:
     417                        levelling_path.append(i)
     418                        for m in range(1,len(path)):
     419                            # this point on the path may have already been leveled in a previous propogation
     420                            if not xleveled[path[m]]:
     421                                newx = x[path[m-1]][0] + crt[(path[m-1],path[m])][qindex][0]
     422                                x[path[m]] = (newx, min(x[path[m-1]][1], crt[(path[m-1],path[m])][qindex][1]))
     423                                xleveled[path[m]] = True
     424                                propogate_levelling(qindex, level, x, xleveled, path[m])
     425                except ValueError:
     426                    break
     427        for j in range(1,len(qlist)):
     428            frob_powers[j] = frob_powers[j].crt(x[j][0])
     429    for j in range(1, len(qlist)):
     430        #print "frob_power[%s] = mod(%s, %s)"%(qlist[j], frob_powers[j], frob_powers[j].modulus())
     431        generators[qlist[j]] = generators[qlist[j]]**(p**(-frob_powers[j]).lift())
     432    _frobenius_shift(K, generators, check_only=True)
     433
     434def compute_pseudo_conway_polynomial(p, n, nodes):
     435    """
     436    Computes a pseudo-Conway polynomial over the prime `p` of degree `n`.
     437
     438    A Conway polynomial `f \in \Bold{F}_p` of degree `n` satisfies the following three conditions:
     439   
     440    - `f` is irreducible.
     441    - In the quotient field `\Bold{F}_p[x]/(f)`, the indeterminant `x` is a multiplicative generator.
     442    - 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`.
     443    - `f` is lexicographically least among all such polynomials, under a certain ordering.
     444
     445    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.
     446
     447    INPUT:
     448   
     449    - `p` -- a prime
     450    - `n` -- a positive integer
     451    - ``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`.
     452
     453    OUTPUT:
     454
     455    A pseudo-Conway polynomial over the prime `p` of degree `n`.
     456
     457    ALGORITHM:
     458   
     459    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.
     460
     461    REFERENCE:
     462
     463    .. [HL99] Heath L. and Loehr N. (1999).
     464       New algorithms for generating Conway polynomials over finite fields.
     465       Proceedings of the tenth annual ACM-SIAM symposium on Discrete algorithms, pp. 429-437.
     466
     467    EXAMPLES::
     468
     469        sage: from sage.rings.finite_rings.constructor import compute_pseudo_conway_polynomial, find_pseudo_conway_polynomial_tree
     470        sage: compute_pseudo_conway_polynomial(next_prime(10000),11,None) # long because the arithmetic in FiniteField_ext_pari fields is slow.
     471        x^11 + x + 7
     472        sage: PCPT30 = find_pseudo_conway_polynomial_tree(2, 30, False)
     473        sage: PCPT20 = find_pseudo_conway_polynomial_tree(2, 20, False)
     474        sage: PCPT12 = find_pseudo_conway_polynomial_tree(2, 12, False)
     475        sage: compute_pseudo_conway_polynomial(2, 60, {2:PCPT30,3:PCPT20,5:PCPT12})
     476        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
     477    """
     478    k = GF(p)
     479    from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     480    R = PolynomialRing(k, names='x')
     481    if n == 1:
     482        return R([-k.multiplicative_generator(), 1])
     483    if nodes is None:
     484        # n is prime.  We require only that the chosen generator (x) be a multiplicative generator for the units.
     485        if p in [2,3,5,7,11]: # here we have Cunningham data
     486            from sage.rings.factorint import factor_cunningham
     487            F = [(p**n-1)//a for a, e in factor_cunningham(n)]
     488        else:
     489            F = (p**n-1).prime_factors() # could be VERY expensive
     490        x = R.gen()
     491        from sage.rings.arith import power_mod
     492        for f in R.polynomials(of_degree=n):
     493            # By using this iterator our chosen polynomial will be fairly sparse.
     494            if f.is_irreducible():
     495                if not f.is_monic():
     496                    f = f.monic()
     497                for m in F:
     498                    if power_mod(x, m, f) == 1:
     499                        break
     500                else:
     501                    return f
     502    # Our strategy is as follows:
     503    # Work in an arbitrary field K of order p**n.
     504    # We're looking for a multiplicative generator of K whose
     505    K = GF(p**n, modulus="first_lexicographic", names='a')
     506    r = p**n - 1
     507    xi = {}
     508    for q in nodes.iterkeys():
     509        #print nodes[q].f.degree(), n//q
     510        xi[q] = nodes[q].f.any_root(K, -n//q, True)
     511    _frobenius_shift(K, xi)
     512           
     513    q, x = xi.popitem()
     514    v = p**(n//q) - 1
     515    for q, xitem in xi.iteritems():
     516        fi = p**(n//q) - 1
     517        g, alpha, beta = v.xgcd(fi)
     518        x = x**beta * xitem**alpha
     519        v = v.lcm(fi)
     520    g = r // v
     521    #print g, n, g.divides(p**n-1)
     522    z = x.nth_root(g, cunningham=(p in [2,3,5,7,11]))
     523    groot1 = K.multiplicative_generator()**v
     524    while z.multiplicative_order() != r:
     525        z *= groot1
     526    #print z
     527    return z.minimal_polynomial()
  • sage/rings/polynomial/polynomial_ring.py

    diff --git a/sage/rings/polynomial/polynomial_ring.py b/sage/rings/polynomial/polynomial_ring.py
    a b  
    21722172        - Peter Bruin (June 2013)
    21732173        """
    21742174        from sage.libs.pari.all import pari
    2175         from sage.rings.finite_rings.constructor import (conway_polynomial,
    2176                                                          exists_conway_polynomial)
     2175        from sage.rings.finite_rings.conway_polynomials import (conway_polynomial,
     2176                  exists_conway_polynomial, find_pseudo_conway_polynomial_tree)
    21772177        from polynomial_gf2x import (GF2X_BuildIrred_list, GF2X_BuildSparseIrred_list,
    21782178                                     GF2X_BuildRandomIrred_list)
    21792179
     
    21932193            return self(pari(p).ffinit(n))
    21942194        elif algorithm == "conway":
    21952195            return self(conway_polynomial(p, n))
     2196        elif algorithm == "pseudo-conway":
     2197            _PCPT = find_pseudo_conway_polynomial_tree(p, k)
     2198            return self(_PCPT.get_pseudo_conway_poly(k))
    21962199        elif algorithm == "first_lexicographic":
    21972200            if p == 2:
    21982201                return self(GF2X_BuildIrred_list(n))