"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
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
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