"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/round2.py" (19 Nov 2012, 16572 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "round2.py" see the Fossies "Dox" file reference documentation.

    1 """
    2 Round 2 method
    3 
    4 The method is for obtaining the maximal order of a number field from a
    5 subring generated by a root of a defining polynomial of the field.
    6 
    7 - H.Cohen CCANT Algorithm 6.1.8
    8 - Kida Y. LN Chapter 3
    9 """
   10 
   11 from __future__ import division
   12 import logging
   13 import nzmath.arith1 as arith1
   14 import nzmath.finitefield as finitefield
   15 import nzmath.factor.misc as factor_misc
   16 import nzmath.gcd as gcd
   17 import nzmath.intresidue as intresidue
   18 import nzmath.matrix as matrix
   19 import nzmath.poly.uniutil as uniutil
   20 import nzmath.rational as rational
   21 import nzmath.vector as vector
   22 import nzmath.squarefree as squarefree
   23 
   24 
   25 _log = logging.getLogger('nzmath.round2')
   26 
   27 uniutil.special_ring_table[finitefield.FinitePrimeField] = uniutil.FinitePrimeFieldPolynomial
   28 
   29 Z = rational.theIntegerRing
   30 Q = rational.theRationalField
   31 
   32 def round2(minpoly_coeff):
   33     """
   34     Return integral basis of the ring of integers of a field with its
   35     discriminant.  The field is given by a list of integers, which is
   36     a polynomial of generating element theta.  The polynomial ought to
   37     be monic, in other word, the generating element ought to be an
   38     algebraic integer.
   39 
   40     The integral basis will be given as a list of rational vectors
   41     with respect to theta.  (In other functions, bases are returned in
   42     the same fashion.)
   43     """
   44     minpoly_int = uniutil.polynomial(enumerate(minpoly_coeff), Z)
   45     d = minpoly_int.discriminant()
   46     squarefactors = _prepare_squarefactors(d)
   47     omega = _default_omega(minpoly_int.degree())
   48     for p, e in squarefactors:
   49         _log.debug("p = %d" % p)
   50         omega = omega + _p_maximal(p, e, minpoly_coeff)
   51 
   52     G = omega.determinant()
   53     return omega.get_rationals(), d * G**2
   54 
   55 def _prepare_squarefactors(disc):
   56     """
   57     Return a list of square factors of disc (=discriminant).
   58 
   59     PRECOND: d is integer
   60     """
   61     squarefactors = []
   62     if disc < 0:
   63         fund_disc, absd = -1, -disc
   64     else:
   65         fund_disc, absd = 1, disc
   66     v2, absd = arith1.vp(absd, 2)
   67     if squarefree.trivial_test_ternary(absd):
   68         fund_disc *= absd
   69     else:
   70         for p, e in factor_misc.FactoredInteger(absd):
   71             if e > 1:
   72                 squareness, oddity = divmod(e, 2)
   73                 squarefactors.append((p, squareness))
   74                 if oddity:
   75                     fund_disc *= p
   76             else:
   77                 fund_disc *= p
   78     if fund_disc & 3 == 1:
   79         if v2 & 1:
   80             squareness, oddity = divmod(v2, 2)
   81             squarefactors.append((2, squareness))
   82             if oddity:
   83                 fund_disc *= 2
   84         else:
   85             squarefactors.append((2, v2 >> 1))
   86     else: # fund_disc & 3 == 3
   87         assert v2 >= 2
   88         fund_disc *= 4
   89         if v2 > 2:
   90             squarefactors.append((2, (v2 - 2) >> 1))
   91     return squarefactors
   92 
   93 def _p_maximal(p, e, minpoly_coeff):
   94     """
   95     Return p-maximal basis with some related informations.
   96 
   97     The arguments:
   98       p: the prime
   99       e: the exponent
  100       minpoly_coeff: (intefer) list of coefficients of the minimal
  101         polynomial of theta
  102     """
  103     # Apply the Dedekind criterion
  104     finished, omega = Dedekind(minpoly_coeff, p, e)
  105     if finished:
  106         _log.info("Dedekind(%d)" % p)
  107         return omega
  108 
  109     # main loop to construct p-maximal order
  110     minpoly = uniutil.polynomial(enumerate(minpoly_coeff), Z)
  111     theminpoly = minpoly.to_field_polynomial()
  112     n = theminpoly.degree()
  113     q = p ** (arith1.log(n, p) + 1)
  114     while True:
  115         # Ip: radical of pO
  116         # Ip = <alpha>, l = dim Ip/pO (essential part of Ip)
  117         alpha, l = _p_radical(omega, p, q, minpoly, n)
  118 
  119         # instead of step 9 big matrix,
  120         # Kida's LN section 2.2
  121         # Up = {x in Ip | xIp \subset pIp} = <zeta> + p<omega>
  122         zeta = _p_module(alpha, l, p, theminpoly)
  123         if zeta.rank == 0:
  124             # no new basis is found
  125             break
  126 
  127         # new order
  128         # 1/p Up = 1/p<zeta> + <omega>
  129         omega2 = zeta / p + omega
  130         if all(o1 == o2 for o1, o2 in zip(omega.basis, omega2.basis)):
  131             break
  132         omega = omega2
  133 
  134     # now <omega> is p-maximal.
  135     return omega
  136 
  137 def _p_radical(omega, p, q, minpoly, n):
  138     """
  139     Return module Ip with dimension of Ip/pO.
  140 
  141     Ip is the radical of pO, or
  142     Ip = {x in O | ~x in kernel f},
  143     where ~x is x mod pO, f is q(=p^e > n)-th powering, which in fact an
  144     Fp-linear map.
  145     """
  146     # Ip/pO = kernel of q-th powering.
  147     # omega_j^q = \sum a_i_j omega_i
  148     base_p = _kernel_of_qpow(omega, q, p, minpoly, n)
  149     l = base_p.column
  150 
  151     # expand basis of Ip/pO to O/pO
  152     base_p.toSubspace()
  153     beta_p = base_p.supplementBasis()
  154 
  155     # basis of Ip
  156     # pulled back from bases of Ip/pO and O/pO
  157     omega_poly = omega.get_polynomials()
  158     alpha_basis = []
  159     for j in range(1, l + 1):
  160         alpha_basis.append(omega.linear_combination([_pull_back(c, p) for c in beta_p[j]]))
  161     for j in range(l + 1, n + 1):
  162         alpha_basis.append(omega.linear_combination([_pull_back(c, p) * p for c in beta_p[j]]))
  163     alpha = ModuleWithDenominator(alpha_basis, omega.denominator)
  164     return alpha, l
  165 
  166 def _kernel_of_qpow(omega, q, p, minpoly, n):
  167     """
  168     Return the kernel of q-th powering, which is a linear map over Fp.
  169     q is a power of p which exceeds n.
  170 
  171     (omega_j^q (mod theminpoly) = \sum a_i_j omega_i   a_i_j in Fp)
  172     """
  173     omega_poly = omega.get_polynomials()
  174     denom = omega.denominator
  175     theminpoly = minpoly.to_field_polynomial()
  176     field_p = finitefield.FinitePrimeField.getInstance(p)
  177     zero = field_p.zero
  178     qpow = matrix.zeroMatrix(n, n, field_p) # Fp matrix
  179     for j in range(n):
  180         a_j = [zero] * n
  181         omega_poly_j = uniutil.polynomial(enumerate(omega.basis[j]), Z)
  182         omega_j_qpow = pow(omega_poly_j, q, minpoly)
  183         redundancy = gcd.gcd(omega_j_qpow.content(), denom ** q)
  184         omega_j_qpow = omega_j_qpow.coefficients_map(lambda c: c // redundancy)
  185         essential = denom ** q // redundancy
  186         while omega_j_qpow:
  187             i = omega_j_qpow.degree()
  188             a_ji = int(omega_j_qpow[i] / (omega_poly[i][i] * essential))
  189             omega_j_qpow -= a_ji * (omega_poly[i] * essential)
  190             if omega_j_qpow.degree() < i:
  191                 a_j[i] = field_p.createElement(a_ji)
  192             else:
  193                 _log.debug("%s / %d" % (str(omega_j_qpow), essential))
  194                 _log.debug("j = %d, a_ji = %s" % (j, a_ji))
  195                 raise ValueError("omega is not a basis")
  196         qpow.setColumn(j + 1, a_j)
  197 
  198     result = qpow.kernel()
  199     if result is None:
  200         _log.debug("(_k_q): trivial kernel")
  201         return matrix.zeroMatrix(n, 1, field_p)
  202     else:
  203         return result
  204 
  205 def _p_module(alpha, l, p, theminpoly):
  206     """
  207     Return basis of Up/pO, where Up = {x in Ip | xIp \subset pIp}.
  208     """
  209     zeta = alpha.get_polynomials()[:l]
  210     for j in range(l):
  211         # refine zeta so that zeta * alpha[j] (mod theminpoly) = 0 (mod pIp)
  212         kernel = _null_linear_combination(zeta, alpha, j, p, theminpoly)
  213         if kernel is None:
  214             zeta = []
  215             break
  216         newzeta = []
  217         for k in range(1, kernel.column + 1):
  218             newzeta.append(sum(_pull_back(c, p) * z for (c, z) in zip(kernel[k], zeta)))
  219         zeta = newzeta
  220     n = theminpoly.degree()
  221     denominator = 1
  222     zeta_list = []
  223     for z in zeta:
  224         while True:
  225             try:
  226                 zeta_list.append([_normalize_int(c * denominator) for c in _coeff_list(z, n)])
  227                 break
  228             except ValueError:
  229                 for i in range(len(zeta_list)):
  230                     zeta_list[i] = [c * p for c in zeta_list[i]]
  231                 denominator *= p
  232                 _log.debug("denominator = %d" % denominator)
  233     # since zeta may be empty, a hint 'dimension' is necessary
  234     return ModuleWithDenominator(zeta_list, denominator, dimension=n)
  235 
  236 def _null_linear_combination(zeta, alpha, j, p, theminpoly):
  237     """
  238     Return linear combination coefficients of tau_i = z_i * alpha[j],
  239     which is congruent to 0 modulo theminpoly and pIp.
  240 
  241     alpha is a module.
  242 
  243     zeta_{j+1} = {z in zeta_j | z * alpha[j] (mod theminpoly) = 0 (mod pIp)}
  244     """
  245     n = theminpoly.degree()
  246     l = len(zeta)
  247     assert n == len(alpha.basis)
  248     alpha_basis = [tuple(b) for b in alpha.get_rationals()]
  249     alpha_mat = matrix.createMatrix(n, alpha_basis, coeff_ring=Q)
  250     alpha_poly = alpha.get_polynomials()
  251 
  252     taus = []
  253     for i in range(l):
  254         tau_i = zeta[i] * alpha_poly[j] % theminpoly
  255         taus.append(tuple(_coeff_list(tau_i, n)))
  256     tau_mat = matrix.createMatrix(n, l, taus, coeff_ring=Q)
  257 
  258     xi = alpha_mat.inverseImage(tau_mat)
  259 
  260     field_p = finitefield.FinitePrimeField.getInstance(p)
  261     xi_p = xi.map(field_p.createElement)
  262     return xi_p.kernel() # linear combination of tau_i's
  263 
  264 def Dedekind(minpoly_coeff, p, e):
  265     """
  266     Return (finished or not, an order)
  267 
  268     the Dedekind criterion
  269 
  270     Arguments:
  271     - minpoly_coeff: (integer) list of the minimal polynomial of theta.
  272     - p, e: p**e divides the discriminant of the minimal polynomial.
  273     """
  274     n = len(minpoly_coeff) - 1  # degree of the minimal polynomial
  275 
  276     m, uniq = _factor_minpoly_modp(minpoly_coeff, p)
  277     omega = _default_omega(n)
  278 
  279     if m == 0:
  280         return True, omega
  281 
  282     minpoly = uniutil.polynomial(enumerate(minpoly_coeff), Z)
  283     v = [_coeff_list(uniq, n)]
  284     shift = uniq
  285     for i in range(1, m):
  286         shift = shift.upshift_degree(1).pseudo_mod(minpoly)
  287         v.append(_coeff_list(shift, n))
  288     updater = ModuleWithDenominator(v, p)
  289 
  290     return (m + 1 > e), updater + omega
  291 
  292 def _factor_minpoly_modp(minpoly_coeff, p):
  293     """
  294     Factor theminpoly modulo p, and return two values in a tuple.
  295     We call gcd(square factors mod p, difference of minpoly and its modp) Z.
  296     1) degree of Z
  297     2) (minpoly mod p) / Z
  298     """
  299     Fp = finitefield.FinitePrimeField.getInstance(p)
  300     theminpoly_p = uniutil.polynomial([(d, Fp.createElement(c)) for (d, c) in enumerate(minpoly_coeff)], Fp)
  301     modpfactors = theminpoly_p.factor()
  302     mini_p = arith1.product([t for (t, e) in modpfactors])
  303     quot_p = theminpoly_p.exact_division(mini_p)
  304     mini = _min_abs_poly(mini_p)
  305     quot = _min_abs_poly(quot_p)
  306     minpoly = uniutil.polynomial(enumerate(minpoly_coeff), Z)
  307     f_p = _mod_p((mini * quot - minpoly).scalar_exact_division(p), p)
  308     gcd = f_p.getRing().gcd
  309     common_p = gcd(gcd(mini_p, quot_p), f_p) # called Z
  310     uniq_p = theminpoly_p // common_p
  311     uniq = _min_abs_poly(uniq_p)
  312 
  313     return common_p.degree(), uniq
  314 
  315 def _mod_p(poly, p):
  316     """
  317     Return modulo p reduction of given integer coefficient polynomial.
  318     """
  319     coeff = {}
  320     for d, c in poly:
  321         coeff[d] = finitefield.FinitePrimeFieldElement(c, p)
  322     return uniutil.polynomial(poly, finitefield.FinitePrimeField.getInstance(p))
  323 
  324 def _min_abs_poly(poly_p):
  325     """
  326     Return minimal absolute mapping of given F_p coefficient polynomial.
  327     """
  328     coeff = {}
  329     p = poly_p.getCoefficientRing().char
  330     for d, c in poly_p:
  331         coeff[d] = _pull_back(c, p)
  332     return uniutil.polynomial(coeff, Z)
  333 
  334 def _coeff_list(upoly, size):
  335     """
  336     Return a list of given size consisting of coefficients of upoly
  337     and possibly zeros padded.
  338     """
  339     return [upoly[i] for i in range(size)]
  340 
  341 def _pull_back(elem, p):
  342     """
  343     Return an integer which is a pull back of elem in Fp.
  344     """
  345     if not isinstance(elem, finitefield.FinitePrimeFieldElement):
  346         if isinstance(elem, intresidue.IntegerResidueClass):
  347             # expecting Z/(p^2 Z)
  348             result = finitefield.FinitePrimeFieldElement(elem.n, p).n
  349         else:
  350             # expecting Z or Q
  351             result = finitefield.FinitePrimeFieldElement(elem, p).n
  352     else:
  353         result = elem.n
  354     if result > (p >> 1): # minimum absolute
  355         result -= p
  356     return result
  357 
  358 def _normalize_int(elem):
  359     """
  360     Return integer object, which is equal to given elem, whose type is
  361     either integer or rational number.
  362     """
  363     if isinstance(elem, (int, long)):
  364         return elem
  365     elif elem.denominator == 1:
  366         return elem.numerator
  367     else:
  368         raise ValueError("non integer: %s" % str(elem))
  369 
  370 def _default_omega(degree):
  371     """
  372     Return the default omega
  373     """
  374     # omega is initialized to basis of Z[theta]
  375     return ModuleWithDenominator([_standard_base(degree, i) for i in range(degree)], 1)
  376 
  377 def _standard_base(degree, i):
  378     """
  379     Return i-th standard unit base
  380     """
  381     base = [0] * degree
  382     base[i] = 1
  383     return base
  384 
  385 def _rational_polynomial(coeffs):
  386     """
  387     Return rational polynomial with given coefficients in ascending
  388     order.
  389     """
  390     return uniutil.polynomial(enumerate(coeffs), Q)
  391 
  392 
  393 class ModuleWithDenominator (object):
  394     """
  395     Represent basis of Z-module with denominator
  396     """
  397     def __init__(self, basis, denominator, **hints):
  398         """
  399         basis is a list of integer sequences.
  400         denominator is a common denominator of all bases.
  401 
  402         Example:
  403         ModuleWithDenominator([[1, 2], [4, 6]], 2) represents a module
  404         <[1/2, 1], [2, 3]>.
  405         """
  406         self.basis = basis
  407         self.denominator = denominator
  408         self.rank = len(self.basis)
  409         if self.rank:
  410             self.dimension = len(self.basis[0])
  411         else:
  412             self.dimension = hints['dimension']
  413         self.rational_basis = None
  414         self.poly_basis = None
  415 
  416     def get_rationals(self):
  417         """
  418         Return a list of rational list, which is basis divided by
  419         denominator.
  420         """
  421         if self.rational_basis is None:
  422             self.rational_basis = []
  423             for base in self.basis:
  424                 self.rational_basis.append([rational.Rational(c, self.denominator) for c in base])
  425         return self.rational_basis
  426 
  427     def get_polynomials(self):
  428         """
  429         Return a list of rational polynomials, which is made from basis
  430         divided by denominator.
  431         """
  432         if self.poly_basis is None:
  433             self.poly_basis = []
  434             for base in self.basis:
  435                 self.poly_basis.append(_rational_polynomial([rational.Rational(c, self.denominator) for c in base]))
  436 
  437         return self.poly_basis
  438 
  439     def __add__(self, other):
  440         """
  441         Return sum of two modules.
  442         """
  443         denominator = gcd.lcm(self.denominator, other.denominator)
  444         row = self.dimension
  445         assert row == other.dimension
  446         column = self.rank + other.rank
  447         mat = matrix.createMatrix(row, column)
  448         adjust = denominator // self.denominator
  449         for i, base in enumerate(self.basis):
  450             mat.setColumn(i + 1, [c * adjust for c in base])
  451         adjust = denominator // other.denominator
  452         for i, base in enumerate(other.basis):
  453             mat.setColumn(self.rank + i + 1, [c * adjust for c in base])
  454 
  455         # Hermite normal form
  456         hnf = mat.hermiteNormalForm()
  457         # The hnf returned by the hermiteNormalForm method may have columns
  458         # of zeros, and they are not needed.
  459         zerovec = vector.Vector([0] * hnf.row)
  460         while hnf.row < hnf.column or hnf[1] == zerovec:
  461             hnf.deleteColumn(1)
  462 
  463         omega = []
  464         for j in range(1, hnf.column + 1):
  465             omega.append(list(hnf[j]))
  466         return ModuleWithDenominator(omega, denominator)
  467 
  468     def __mul__(self, scale):
  469         """
  470         scalar multiplication.
  471         """
  472         if not (self.denominator % scale):
  473             return ModuleWithDenominator(self.basis, self.denominator // scale)
  474         else:
  475             muled = [[c * scale for c in base] for base in self.basis]
  476             return ModuleWithDenominator(muled, self.denominator)
  477 
  478     __rmul__ = __mul__
  479 
  480     def __truediv__(self, divisor):
  481         return ModuleWithDenominator(self.basis, self.denominator * divisor)
  482 
  483     def determinant(self):
  484         """
  485         Return determinant of the basis (basis ought to be of full
  486         rank and in Hermite normal form).
  487         """
  488         return arith1.product([rational.Rational(self.basis[i][i], self.denominator) for i in range(self.rank)])
  489 
  490     def linear_combination(self, coefficients):
  491         """
  492         Return a list corresponding a linear combination of basis with
  493         given coefficients.  The denominator is ignored.
  494         """
  495         new_basis = [0] * self.dimension
  496         for c, base in zip(coefficients, self.basis):
  497             for i in range(self.dimension):
  498                 new_basis[i] += c * base[i]
  499         return new_basis