NZMATH  1.2.0
About: NZMATH is a Python based number theory oriented calculation system.
  Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

round2.py
Go to the documentation of this file.
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 
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 
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
nzmath.finitefield.FinitePrimeFieldElement
Definition: finitefield.py:159
nzmath.round2._prepare_squarefactors
def _prepare_squarefactors(disc)
Definition: round2.py:55
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.finitefield.FinitePrimeField
Definition: finitefield.py:205
nzmath.round2._rational_polynomial
def _rational_polynomial(coeffs)
Definition: round2.py:385
nzmath.round2._coeff_list
def _coeff_list(upoly, size)
Definition: round2.py:334
nzmath.round2.ModuleWithDenominator.dimension
dimension
Definition: round2.py:410
nzmath.round2._null_linear_combination
def _null_linear_combination(zeta, alpha, j, p, theminpoly)
Definition: round2.py:236
nzmath.factor.misc
Definition: misc.py:1
nzmath.round2._mod_p
def _mod_p(poly, p)
Definition: round2.py:315
nzmath.round2.ModuleWithDenominator.rank
rank
Definition: round2.py:408
nzmath.round2.ModuleWithDenominator.determinant
def determinant(self)
Definition: round2.py:483
nzmath.matrix
Definition: matrix.py:1
nzmath.gcd
Definition: gcd.py:1
nzmath.squarefree
Definition: squarefree.py:1
nzmath.rational
Definition: rational.py:1
nzmath.round2.ModuleWithDenominator.poly_basis
poly_basis
Definition: round2.py:414
nzmath.rational.Rational
Definition: rational.py:10
nzmath.intresidue.IntegerResidueClass
Definition: intresidue.py:12
nzmath.vector.Vector
Definition: vector.py:3
nzmath.round2._min_abs_poly
def _min_abs_poly(poly_p)
Definition: round2.py:324
nzmath.round2.ModuleWithDenominator.__mul__
def __mul__(self, scale)
Definition: round2.py:468
nzmath.round2.ModuleWithDenominator.denominator
denominator
Definition: round2.py:407
nzmath.round2._standard_base
def _standard_base(degree, i)
Definition: round2.py:377
nzmath.gcd.gcd
def gcd(a, b)
Definition: gcd.py:7
nzmath.intresidue
Definition: intresidue.py:1
nzmath.vector
Definition: vector.py:1
nzmath.arith1
Definition: arith1.py:1
nzmath.poly.uniutil
Definition: uniutil.py:1
nzmath.round2._default_omega
def _default_omega(degree)
Definition: round2.py:370
nzmath.round2.ModuleWithDenominator.__truediv__
def __truediv__(self, divisor)
Definition: round2.py:480
nzmath.round2.ModuleWithDenominator.__add__
def __add__(self, other)
Definition: round2.py:439
nzmath.round2._p_maximal
def _p_maximal(p, e, minpoly_coeff)
Definition: round2.py:93
nzmath.round2.ModuleWithDenominator.get_polynomials
def get_polynomials(self)
Definition: round2.py:427
nzmath.round2._pull_back
def _pull_back(elem, p)
Definition: round2.py:341
nzmath.round2._p_radical
def _p_radical(omega, p, q, minpoly, n)
Definition: round2.py:137
nzmath.round2.ModuleWithDenominator
Definition: round2.py:393
nzmath.round2.round2
def round2(minpoly_coeff)
Definition: round2.py:32
nzmath.round2._p_module
def _p_module(alpha, l, p, theminpoly)
Definition: round2.py:205
nzmath.round2._normalize_int
def _normalize_int(elem)
Definition: round2.py:358
nzmath.round2.Dedekind
def Dedekind(minpoly_coeff, p, e)
Definition: round2.py:264
nzmath.round2.ModuleWithDenominator.__init__
def __init__(self, basis, denominator, **hints)
Definition: round2.py:397
nzmath.round2.ModuleWithDenominator.rational_basis
rational_basis
Definition: round2.py:413
nzmath.round2.ModuleWithDenominator.get_rationals
def get_rationals(self)
Definition: round2.py:416
nzmath.round2._factor_minpoly_modp
def _factor_minpoly_modp(minpoly_coeff, p)
Definition: round2.py:292
nzmath.round2.ModuleWithDenominator.basis
basis
Definition: round2.py:406
nzmath.round2._kernel_of_qpow
def _kernel_of_qpow(omega, q, p, minpoly, n)
Definition: round2.py:166
nzmath.round2.ModuleWithDenominator.linear_combination
def linear_combination(self, coefficients)
Definition: round2.py:490
nzmath.finitefield
Definition: finitefield.py:1