## "Fossies" - the Fresh Open Source Software Archive

### Member "NZMATH-1.2.0/nzmath/prime_decomp.py" (19 Nov 2012, 14285 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 "prime_decomp.py" see the Fossies "Dox" file reference documentation.

```    1 """
2 prime ideal decomposition over an (absolute) number field
3 """
4
5 import nzmath.arith1 as arith1
6 import nzmath.rational as rational
7 import nzmath.matrix as matrix
8 import nzmath.vector as vector
9 import nzmath.algfield as algfield
10 import nzmath.finitefield as finitefield
11 import nzmath.poly.uniutil as uniutil
12 import nzmath.module as module
13 import nzmath.algorithm as algorithm
14
15
16 def prime_decomp(p, polynomial):
17     """
18     Return prime decomposition of (p) over Q[x]/(polynomial).
19     This method returns a list of (P_k, e_k, f_k),
20     where P_k is an instance of Ideal_with_generator,
21           e_k is the ramification index of P_k,
22           f_k is the residue degree of P_k.
23     """
24     field = algfield.NumberField(polynomial)
25     field_disc = field.field_discriminant()
26     if (field.disc() // field_disc) % p != 0:
27         return _easy_prime_decomp(p, polynomial)
28     else:
29         return _main_prime_decomp(p, polynomial)
30
31 def _easy_prime_decomp(p, polynomial):
32     """
33     prime decomposition by factoring polynomial
34     """
35     f = algfield.fppoly(polynomial, p)
36     d = f.degree()
37     factor = f.factor()
38     p_alg = algfield.BasicAlgNumber([[p] + [0] * (d - 1), 1], polynomial)
39     if len(factor) == 1 and factor[0][1] == 1:
40         return [(module.Ideal_with_generator([p_alg]),
41             1, factor[0][0].degree())]
42     else:
43         dlist = []
44         for i in range(len(factor)):
45             basis_list = []
46             for j in range(d):
47                 if factor[i][0][j] == 0:
48                     basis_list.append(0)
49                 else:
50                     basis_list.append(factor[i][0][j].toInteger())
51             prime_ideal = module.Ideal_with_generator([p_alg,
52                           algfield.BasicAlgNumber([basis_list, 1], polynomial)])
53             dlist.append((prime_ideal, factor[i][1], factor[i][0].degree()))
54         return dlist
55
56 def _main_prime_decomp(p, polynomial):
57     """
58     main step of prime decomposition
59     """
60     field = algfield.NumberField(polynomial)
61     n = field.degree
62     int_basis = field.integer_ring()
63     base_multiply = module._base_multiplication(int_basis, field)
64     H_lst = _squarefree_decomp(p, field, base_multiply) # step 2-5
65     done_list = []
66     for j in range(len(H_lst)):
67         H_j = H_lst[j]
68         # step 7
69         if H_j != None and H_j.column == n: #rank(H_j)
70             continue
71         L = [H_j]
72         while L:
73             H = L.pop()
74             A, gamma_multiply = _separable_algebra(
75                 p, H, base_multiply) # step 8, 9
76             sol = _check_H(p, H, field, gamma_multiply) # step 10-11
77             if isinstance(sol, tuple):
78                 done_list.append((sol[0], j + 1, sol[1]))
79             else:
80                 L += _splitting_squarefree(p, H, base_multiply,
81                     A, gamma_multiply, sol) # step 12-15
82     return done_list
83
84 def _squarefree_decomp(p, field, base_multiply):
85     """
86     decompose (p) as (p)= A_1 (A_2)^2 (A_3)^2 ..., where each A_i is squarefree
87     """
88     #CCANT Algo 6.2.9 step 2-5
89     n = field.degree
90     int_basis = field.integer_ring()
91     # step 2
92     log_n_p = arith1.log(n, p)
93     q = p ** log_n_p
94     if  q < n:
95         q *= p
96     base_p = _kernel_of_qpow(int_basis, q, p, field)
97     I_p_over_pO = base_p
98     # step 3
99     K_lst = [I_p_over_pO]
100     while K_lst[-1]: # K[-1] is zero matrix in F_p
101         K_lst.append(_mul_mod_pO(K_lst[0], K_lst[-1], base_multiply))
102     i = len(K_lst)
103     # step 4
104     J_lst = [K_lst[0]]
105     for j in range(1, i):
106         J_lst.append(_div_mod_pO(K_lst[j], K_lst[j - 1], base_multiply, p))
107     # step 5
108     H_lst = []
109     for j in range(i - 1):
110         H_lst.append(_div_mod_pO(J_lst[j], J_lst[j + 1], base_multiply, p))
111     H_lst.append(J_lst[-1])
112     return H_lst
113
114 def _kernel_of_qpow(basis_matrix, q, p, field):
115     """
116     compute kernel of w^q.
117     (this method is same as round2._kernel_of_qpow)
118     """
119     omega = _matrix_to_algnumber_list(basis_matrix, field)
120     omega_pow = []
121     for omega_i in omega:
122         omega_pow.append(omega_i ** q)
123     omega_pow_mat = _algnumber_list_to_matrix(omega_pow, field)
124     A_over_Z = basis_matrix.inverse(omega_pow_mat) #theta repr -> omega repr
125     A = A_over_Z.map(
126          lambda x: finitefield.FinitePrimeFieldElement(x.numerator, p))
127     return A.kernel()
128
129 def _mul_mod_pO(self, other, base_multiply):
130     """
131     ideal multiplication modulo pO using base_multiply
132     """
133     #CCANT Algorithm 6.2.5
134     n = base_multiply.column
135     new_self = _matrix_mul_by_base_mul(self, other, base_multiply)
136     return new_self.image()
137
138 def _div_mod_pO(self, other, base_multiply, p):
139     """
140     ideal division modulo pO by using base_multiply
141     """
142     n = other.row
143     m = other.column
144     F_p = finitefield.FinitePrimeField(p)
145     #Algo 2.3.7
146     if self == None:
147         self_new = matrix.zeroMatrix(n, F_p)
148         r = 0
149     else:
150         self_new = self
151         r = self.column
152     if r == m:
153         return matrix.unitMatrix(n, F_p)
154     # if r > m, raise NoInverseImage error
155     X = matrix.Subspace.fromMatrix(other.inverseImage(self_new))
156     B = X.supplementBasis()
157     other_self = other * B # first r columns are self, others are supplement
158     other_self.toFieldMatrix()
159
160     gamma_part = other_self.subMatrix(
161         range(1, n + 1), range(r + 1, m + 1))
162     omega_gamma = other_self.inverseImage(_matrix_mul_by_base_mul(
163         matrix.unitMatrix(n, F_p), gamma_part, base_multiply))
164     vect_list = []
165     for k in range(1, n + 1):
166         vect = vector.Vector([F_p.zero] * ((m - r) ** 2))
167         for i in range(1, m - r + 1):
168             for j in range(1, m - r + 1):
169                 vect[(m - r) * (i - 1) + j] = _mul_place(
170                     k, i, omega_gamma, m - r)[j + r]
171         vect_list.append(vect)
172     return matrix.FieldMatrix( (m - r) ** 2, n, vect_list).kernel()
173
174 def _separable_algebra(p, H, base_multiply):
175     """
176     return A=O/H (as matrix representation)
177     """
178     # CCANT Algo 6.2.9 step 8-9
179     F_p = finitefield.FinitePrimeField(p)
180     n = base_multiply.row
181     if H == None:
182         H_new = matrix.zeroMatrix(n, 1, F_p)
183         r = 0
184     else:
185         H_new = H
186         r = H.column
187     # step 8
190         raise ValueError
191     H_add_one.extendColumn(vector.Vector([F_p.one] + [F_p.zero] * (n - 1)))
194     full_base.toFieldMatrix()
195     # step 9
196     A = full_base.subMatrix(
197         range(1, n + 1), range(r + 1, n + 1))
198     gamma_mul = full_base.inverseImage(
199         _matrix_mul_by_base_mul(A, A, base_multiply))
200     gamma_mul = gamma_mul.subMatrix(
201         range(r + 1, n + 1), range(1, gamma_mul.column + 1))
202     return A, gamma_mul
203
204 def _check_H(p, H, field, gamma_mul):
205     """
206     If H express the prime ideal, return (H, f),
207     where f: residual degree.
208     Else return some column of M_1 (not (1,0,..,0)^t).
209     """
210     F_p = finitefield.FinitePrimeField(p)
211     if H == None:
212         f = field.degree
213     else:
214         f = H.row - H.column # rank(A), A.column
215     # CCANT Algo 6.2.9 step 10-11
216     # step 10
217     B = matrix.unitMatrix(f, F_p)
218     M_basis = []
219     for j in range(1, f + 1):
220         alpha_pow = _pow_by_base_mul(B[j], p, gamma_mul, f)
221         alpha_pow -= B[j]
222         M_basis.append(alpha_pow)
223     M_1 = matrix.FieldMatrix(f, f, M_basis).kernel()
224     # step 11
225     if M_1.column > 1: # dim(M_1) > 1
226         return M_1[M_1.column]
227     else:
228         H_simple = _two_element_repr_prime_ideal(
229                    p, H.map(lambda x: x.getResidue()), field, f)
230         p_alg = algfield.BasicAlgNumber([[p] + [0] * (field.degree - 1), 1],
231                 field.polynomial)
232         sol = module.Ideal_with_generator([p_alg, H_simple])
233         return (sol, f)
234
235 def _two_element_repr_prime_ideal(p, hnf, field, f):
236     """
237     return alpha such that (p, alpha) is same ideal to
238     the ideal represented by hnf (matrix) w.r.t. the integer ring of field.
239     we assume the ideal is a prime ideal whose norm is p^f, where p is prime.
240     """
241     # assume that first basis of the integral basis equals 1
242     int_basis = field.integer_ring()
243
244     # step 1
245     n = field.degree
246     Z = rational.theIntegerRing
247     beta = (p * matrix.unitMatrix(n, Z).subMatrix(
248         range(1, n + 1), range(2, n + 1))).copy()
249     beta.extendColumn(hnf)
250     m = beta.column
251     # step 2
252     R = 0
253     P_norm = p ** f
254     while True:
255         R += 1
256         lmd = [R] * m
257         while True:
258             # step 3
259             alpha = lmd[0] * beta[1]
260             for j in range(1, m):
261                 alpha += lmd[j] * beta[j + 1]
262             alpha = _vector_to_algnumber(int_basis * alpha, field)
263             n = alpha.norm() / P_norm
264             r = divmod(int(n), p)[1]
265             if r:
266                 return alpha
267             # step 4
268             for j in range(m - 1, 0, -1):
269                 if lmd[j] != -R:
270                     break
271             lmd[j] -= 1
272             for i in range(j + 1, m):
273                 lmd[i] = R
274             # step 5
275             for j in range(m):
276                 if lmd[j] != 0:
277                     break
278             else:
279                 break
280
281 def _splitting_squarefree(p, H, base_multiply, A, gamma_mul, alpha):
282     """
283     split squarefree part
284     """
285     F_p = finitefield.FinitePrimeField(p)
286     if H == None:
287         n = base_multiply.row
288         f = n
289         H_new = matrix.zeroMatrix(n, 1, F_p)
290     else:
291         n = H.row
292         f = n - H.column
293         H_new = H.copy()
294     # step 12
295     one_gamma = vector.Vector([F_p.one] + [F_p.zero] * (f - 1)) # w.r.t. gamma_1
296     minpoly_matrix_alpha = matrix.FieldMatrix(f, 2, [one_gamma, alpha])
297     ker = minpoly_matrix_alpha.kernel()
298     alpha_pow = alpha
299     while ker == None:
300         alpha_pow = _vect_mul_by_base_mul(alpha_pow, alpha, gamma_mul, f)
301         minpoly_matrix_alpha.extendColumn(alpha_pow)
302         ker = minpoly_matrix_alpha.kernel()
303     minpoly_alpha = uniutil.FinitePrimeFieldPolynomial(
304         enumerate(ker[1].compo), F_p)
305     # step 13
306     m_list = minpoly_alpha.factor()
307     # step 14
308     L = []
309     for (m_s, i) in m_list:
310         M_s = H_new.copy()
311         beta_s_gamma = _substitution_by_base_mul(
312             m_s, alpha, gamma_mul, one_gamma, f)
313         # beta_s_gamma (w.r.t. gamma) -> beta_s (w.r.t. omega)
314         beta_s = A * beta_s_gamma
315         omega_beta_s = _matrix_mul_by_base_mul(
316             matrix.unitMatrix(n, F_p), beta_s.toMatrix(True), base_multiply)
317         M_s.extendColumn(omega_beta_s)
318         L.append(M_s.image())
319     return L
320
321 #################################################
322 # Utility function
323 #################################################
324 def _vector_to_algnumber(vector_repr, field):
325     """
326     vector (w.r.t. theta) -> omega
327     """
328     int_vector, denom = module._toIntegerMatrix(vector_repr)
329     int_vector = int_vector[1].compo # Submodule -> list
330     omega = algfield.BasicAlgNumber(
331             [int_vector, denom], field.polynomial)
332     return omega
333
334 def _algnumber_to_vector(omega, field):
335     """
336     omega -> vector (w.r.t. theta)
337     """
338     n = field.degree
339     vect = [rational.Rational(
340         omega.coeff[j], omega.denom) for j in range(n)]
341     return vector.Vector(vect)
342
343 def _matrix_to_algnumber_list(matrix_repr, field):
344     """
345     matrix (w.r.t. theta) -> [omega_i]
346     """
347     int_matrix, denom = module._toIntegerMatrix(matrix_repr)
348     omega = []
349     for i in range(1, int_matrix.column + 1):
350         omega_i = algfield.BasicAlgNumber(
351             [int_matrix[i], denom], field.polynomial)
352         omega.append(omega_i)
353     return omega
354
355 def _algnumber_list_to_matrix(omega, field):
356     """
357     [omega_i] -> matrix (w.r.t. theta)
358     """
359     vect = []
360     for omega_i in omega:
361         vect_i = _algnumber_to_vector(omega_i, field)
362         vect.append(vector.Vector(vect_i))
363     return matrix.Matrix(omega[0].degree, len(omega), vect)
364
365 def _vect_mul_by_base_mul(self, other, base_mul, max_j=False):
366     """
367     return self * other with respect to basis (gamma_i),
368     where self, other are vector.
369     (base_mul)_ij express gamma_i gamma_j
370     """
371     n = base_mul.row
372     sol_val = vector.Vector([0] * n)
373     for i1 in range(1, n + 1):
374         for i2 in range(1, n + 1):
375             if not max_j:
376                 sol_val += self[i1] * other[
377                     i2] * module._symmetric_element(i1, i2, base_mul)
378             else:
379                 sol_val += self[i1] * other[
380                     i2] * _mul_place(i1, i2, base_mul, max_j)
381     return sol_val
382
383 def _matrix_mul_by_base_mul(self, other, base_mul, max_j=False):
384     """
385     return self * other with respect to basis (gamma_i),
386     where self, other are matrix.
387     (base_mul)_ij express gamma_i gamma_j
388     """
389     n = base_mul.row
390     sol = []
391     for k1 in range(1, self.column + 1):
392         for k2 in range(1, other.column + 1):
393             sol_ele = _vect_mul_by_base_mul(
394                 self[k1], other[k2], base_mul, max_j)
395             sol.append(sol_ele)
396     return matrix.FieldMatrix(n, len(sol), sol)
397
398 def _mul_place(i, j, mul_matrix, max_j):
399     """
400     find the vector in mul_matrix corresponding (i, j),
401     where 1 <= j <= max_j.
402     """
403     return mul_matrix[(i - 1) * max_j + j]
404
405 def _pow_by_base_mul(self, other, base_mul, max_j=False):
406     """
407     return self * other with respect to basis (gamma_i),
408     where self, other are matrix.
409     (base_mul)_ij express gamma_i * gamma_j
410     This function uses right-left binary method.
411     """
412     func = algorithm.powering_func(
413         lambda x,y:_vect_mul_by_base_mul(x, y, base_mul, max_j))
414     return func(self, other)
415
416 def _substitution_by_base_mul(P, ele, base_mul, one_gamma, max_j=False):
417     """
418     return P(ele) with respect to basis (gamma_i),
419     where P is a polynomial, ele is a vector w.r.t. gamma_i.
420     (base_mul)_ij express gamma_i * gamma_j
421     This function uses digital method (Horner method).
422     """
423     coeff = reversed(P.sorted)
424     mul = lambda x, y : _vect_mul_by_base_mul(x, y, base_mul, max_j)
425     return algorithm.digital_method(
426         coeff, ele, lambda x, y : x + y, mul,
427         lambda a, x : a * x,
428         algorithm.powering_func(one_gamma, mul),
429         vector.Vector([base_mul.coeff_ring.zero] * base_mul.row),
430         one_gamma
431         )
```