"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
  188     H_add_one = H_new.copy()
  189     if H_add_one.column == n:
  190         raise ValueError
  191     H_add_one.extendColumn(vector.Vector([F_p.one] + [F_p.zero] * (n - 1)))
  192     H_add_one.toSubspace(True)
  193     full_base = H_add_one.supplementBasis()
  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         )