1 #include "cado.h" // IWYU pragma: keep
2 /*
3 * Authors: Joshua Peignier and Emmanuel Thomé
4 */
5 #include "numbertheory.hpp"
6 #include <gmp.h> // for mpz_t, mpz_mul, mpq_denref, mpq_numref, mpz_c...
7 #include <algorithm> // for sort
8 #include <cstdio> // for NULL, size_t
9 #include <sstream> // std::ostringstream // IWYU pragma: keep
10 #include <memory> // for allocator_traits<>::value_type
11 #include <string> // for char_traits, operator<<, string
12 #include <type_traits> // for __strip_reference_wrapper<>::__type
13 #include "cxx_mpz.hpp" // for cxx_mpz, operator<<
14 #include "gmp_aux.h" // for mpz_p_valuation
15 #include "macros.h" // for ASSERT_ALWAYS
16
17
18 using namespace std;
19
20 /*{{{ conversion of rows and columns to polynomials*/
21 void mpz_mat_row_to_poly(mpz_poly_ptr f, mpz_mat_srcptr M, const unsigned int i);
22 void mpz_mat_row_to_poly_rev(mpz_poly_ptr f, mpz_mat_srcptr M, const unsigned int i);
23 void mpz_mat_column_to_poly(mpz_poly_ptr f, mpz_mat_srcptr M, const unsigned int j);
24 void mpq_mat_row_to_poly(mpz_poly_ptr f, mpz_ptr lcm, mpq_mat_srcptr M, const unsigned int i);
25 void mpq_poly_to_mat_row(mpq_mat_ptr M, const unsigned int i, mpz_poly_srcptr f, mpz_srcptr denom);
26 void mpq_mat_column_to_poly(mpz_poly_ptr f, mpz_ptr lcm, mpq_mat_srcptr M, const unsigned int j);
27 /* }}} */
28
29 /*{{{ conversion of rows and columns to polynomials*/
mpz_mat_row_to_poly(mpz_poly_ptr f,mpz_mat_srcptr M,const unsigned int i)30 void mpz_mat_row_to_poly(mpz_poly_ptr f, mpz_mat_srcptr M, const unsigned int i)
31 {
32 mpz_poly_realloc(f,M->n);
33 unsigned int j;
34 for (j = 0 ; j < M->n; j++){
35 mpz_poly_setcoeff(f,j,mpz_mat_entry_const(M,i,j));
36 }
37 mpz_poly_cleandeg(f, M->n - 1);
38 }
39
mpz_mat_row_to_poly_rev(mpz_poly_ptr f,mpz_mat_srcptr M,const unsigned int i)40 void mpz_mat_row_to_poly_rev(mpz_poly_ptr f, mpz_mat_srcptr M, const unsigned int i)
41 {
42 mpz_poly_realloc(f,M->n);
43 unsigned int j;
44 for (j = 0 ; j < M->n; j++){
45 mpz_poly_setcoeff(f,j,mpz_mat_entry_const(M,i,M->n-1-j));
46 }
47 mpz_poly_cleandeg(f, M->n - 1);
48 }
49
mpz_mat_column_to_poly(mpz_poly_ptr f,mpz_mat_srcptr M,const unsigned int j)50 void mpz_mat_column_to_poly(mpz_poly_ptr f, mpz_mat_srcptr M, const unsigned int j)
51 {
52 mpz_poly_realloc(f,M->m);
53 unsigned int i;
54 for (i = 0 ; i < M->m; i++){
55 mpz_poly_setcoeff(f,i,mpz_mat_entry_const(M,i,j));
56 }
57 mpz_poly_cleandeg(f, M->m - 1);
58 }
59
mpq_mat_row_to_poly(mpz_poly_ptr f,mpz_ptr lcm,mpq_mat_srcptr M,const unsigned int i)60 void mpq_mat_row_to_poly(mpz_poly_ptr f, mpz_ptr lcm, mpq_mat_srcptr M, const unsigned int i)
61 {
62 /* read element w[i] as a polynomial. beware, the element might
63 * have rational coordinates, that's life ! */
64 unsigned int n = M->n;
65 ASSERT_ALWAYS(i < M->m);
66 mpz_set_ui(lcm, 1);
67 for (unsigned int j = 0; j < n; j++) {
68 mpz_lcm(lcm, lcm, mpq_denref(mpq_mat_entry_const(M, i, j)));
69 }
70 mpz_poly_realloc(f, n);
71 for (unsigned int j = 0; j < n; j++) {
72 mpq_srcptr mij = mpq_mat_entry_const(M, i, j);
73 mpz_divexact(f->coeff[j], lcm, mpq_denref(mij));
74 mpz_mul(f->coeff[j], f->coeff[j], mpq_numref(mij));
75 }
76 mpz_poly_cleandeg(f, n-1);
77 }
78
mpq_poly_to_mat_row(mpq_mat_ptr M,const unsigned int i,mpz_poly_srcptr f,mpz_srcptr denom)79 void mpq_poly_to_mat_row(mpq_mat_ptr M, const unsigned int i, mpz_poly_srcptr f, mpz_srcptr denom)
80 {
81 ASSERT_ALWAYS(f->deg < (int) M->n);
82 for (unsigned int j = 0 ; j < M->n; j++){
83 mpq_ptr mij = mpq_mat_entry(M,i,j);
84 mpz_poly_getcoeff(mpq_numref(mij), j, f);
85 mpz_set(mpq_denref(mij), denom);
86 mpq_canonicalize(mij);
87 }
88 }
89
mpq_mat_column_to_poly(mpz_poly_ptr f,mpz_ptr lcm,mpq_mat_srcptr M,const unsigned int j)90 void mpq_mat_column_to_poly(mpz_poly_ptr f, mpz_ptr lcm, mpq_mat_srcptr M, const unsigned int j)
91 {
92 mpz_poly_realloc(f,M->m);
93 mpz_set_si(lcm,1);
94 for (unsigned int i = 0 ; i < M->m ; i++) {
95 mpz_lcm(lcm,lcm,mpq_denref(mpq_mat_entry_const(M,i,j)));
96 }
97 for (unsigned int i = 0 ; i < M->m; i++){
98 mpq_srcptr mij = mpq_mat_entry_const(M, i, j);
99 mpz_divexact(f->coeff[i], lcm, mpq_denref(mij));
100 mpz_mul(f->coeff[i], f->coeff[i], mpq_numref(mij));
101 }
102 mpz_poly_cleandeg(f, M->m - 1);
103 }
104 /*}}}*/
105
106 /*{{{ commonly used wrappers around HNF functions */
mpz_mat_hnf_backend_rev(mpz_mat_ptr M,mpz_mat_ptr T)107 int mpz_mat_hnf_backend_rev(mpz_mat_ptr M, mpz_mat_ptr T) // {{{
108 {
109 /* This is almost like hnf_backend, except that we do it in a
110 * different order, which is more suitable for displaying number
111 * field elements in a way which ends up being similar to magma's
112 *
113 * T receives the transformation matrix.
114 * M is put into HNF form.
115 */
116 mpz_mat_reverse_rows(M, M);
117 mpz_mat_reverse_columns(M, M);
118 int s = mpz_mat_hnf_backend(M, T);
119 mpz_mat_reverse_rows(M, M);
120 mpz_mat_reverse_columns(M, M);
121 if (T) mpz_mat_reverse_rows(T, T);
122 if (T) mpz_mat_reverse_columns(T, T);
123 if (M->m > M->n) {
124 /* we need some swaps... */
125 mpz_mat sM;
126 mpz_mat_init(sM, M->m, M->n);
127 mpz_mat_submat_swap(sM, 0, 0, M, M->m-M->n, 0, M->n, M->n);
128 mpz_mat_submat_swap(sM, M->n, 0, M, 0, 0, M->m-M->n, M->n);
129 mpz_mat_swap(sM, M);
130 mpz_mat_clear(sM);
131 if (T) {
132 mpz_mat sT;
133 mpz_mat_init(sT, T->m, T->n);
134 mpz_mat_submat_swap(sT, 0, 0, T, T->m-T->n, 0, T->n, T->n);
135 mpz_mat_submat_swap(sT, T->n, 0, T, 0, 0, T->m-T->n, T->n);
136 mpz_mat_swap(sT, T);
137 mpz_mat_clear(sT);
138 }
139 /* While the transformations above had no effect on s (because
140 * they compensate), this one has.
141 * we have n circular shifts on length m, plus a reversal on m-n.
142 * a circular shift on length k is exactly k-1 inversions, so
143 * that sums up to n*(m-1) inversions. Then we add
144 * (m-n)*(m-n-1)/2 inversions. This is, in total,
145 * (m(m-1)+n(n-1))/2 inversions.
146 * m*(m-1) is congruent to 2 mod 4 when m is 2 or 3 mod 4
147 */
148 int ninvs = ((M->m&2)+(M->n&2))/2;
149 if (ninvs) s=-s;
150 }
151 return s;
152 }//}}}
153
join_HNF(cxx_mpz_mat const & K,cxx_mpz const & p)154 cxx_mpz_mat join_HNF(cxx_mpz_mat const& K, cxx_mpz const& p)//{{{
155 {
156 cxx_mpz_mat J(K->n, K->n, p);
157 cxx_mpz_mat T0;
158 cxx_mpz_mat I;
159
160 mpz_mat_vertical_join(I, J, K);
161 mpz_mat_hnf_backend(I, T0);
162 mpz_mat_submat_swap(I, 0, 0, J, 0, 0, K->n, K->n);
163 return J;
164 }
165 //}}}
166
join_HNF_rev(cxx_mpz_mat const & K,cxx_mpz const & p)167 cxx_mpz_mat join_HNF_rev(cxx_mpz_mat const& K, cxx_mpz const& p)//{{{
168 {
169 // Builds the block matrix containing p*identity in the top, and K in
170 // the bottom Then computes its HNF and stores it in I
171 cxx_mpz_mat J(K->n, K->n, p);
172 cxx_mpz_mat T0;
173 cxx_mpz_mat I;
174
175 mpz_mat_vertical_join(I, J, K);
176 mpz_mat_hnf_backend_rev(I, T0);
177 mpz_mat_submat_swap(I, 0, 0, J, 0, 0, K->n, K->n);
178 return J;
179 }//}}}
180 /*}}}*/
181
182 /* {{{ multiplication_table_of_order */
multiplication_table_of_order(cxx_mpq_mat const & O,cxx_mpz_poly const & g)183 cxx_mpz_mat multiplication_table_of_order(cxx_mpq_mat const& O,
184 cxx_mpz_poly const& g)
185 {
186 /* Let O be an order, with basis written with respect to the
187 * polynomial basis defined by g (of degree denoted below by n). O is
188 * thus an n*n matrix.
189 *
190 * This function computes the integer matrix of size n*n^2 such that
191 * the n coordinates at position (i,j*n) to (i,j*n+n-1) are the
192 * coordinates of the i-th times the j-th generator of O, expressed
193 * as combinations of the generators of O.
194 */
195 unsigned int n = g->deg;
196 ASSERT_ALWAYS(O->m == n);
197 ASSERT_ALWAYS(O->n == n);
198 cxx_mpq_mat R;
199 mpq_mat_inv(R, O);
200 cxx_mpz_mat M(n, n * n);
201
202 for(unsigned int i = 0; i < n; i++) {
203 cxx_mpz_poly w;
204 cxx_mpz dw;
205 mpq_mat_row_to_poly(w, dw, O, i);
206 cxx_mpq_mat T(n, n);
207 for (unsigned int j = 0; j < n; j++) {
208 cxx_mpz_poly c;
209 cxx_mpz dc;
210 mpq_mat_row_to_poly(c, dc, O, j);
211 mpz_poly_mul_mpz(c, c, mpz_poly_lc(g));
212 mpz_poly_mul_mod_f(c, c, w, g);
213 mpz_mul(dc, dc, dw);
214 mpz_mul(dc, dc, mpz_poly_lc(g));
215 mpq_poly_to_mat_row(T, j, c, dc);
216 }
217 mpq_mat_mul(T, T, R);
218 cxx_mpz_mat Tz;
219 int rc = mpq_mat_numden(Tz, NULL, T);
220 ASSERT_ALWAYS(rc);
221 for (unsigned int j = 0; j < n; j++) {
222 mpz_mat_submat_swap(M, i, j*n, Tz, j, 0, 1, n);
223 }
224 }
225 return M;
226 }/*}}}*/
227
228 /* {{{ multiplication_table_of_ideal*/
multiplication_table_of_ideal(cxx_mpz_mat const & M,cxx_mpz_mat const & I)229 cxx_mpz_mat multiplication_table_of_ideal(cxx_mpz_mat const& M,
230 cxx_mpz_mat const& I)
231 {
232 /* Let O be an order of a degree n number field. Let M be the n*n^2
233 * multiplication matrix of O.
234 *
235 * Let I be an ideal of O, given as an n*n matrix with entries
236 * expressed as coordinate vectors with respect to the basis O.
237 *
238 * This function computes the integer matrix of size n*n^2 such that
239 * the n coordinates at position (i,j*n) to (i,j*n+n-1) are the
240 * coordinates of the i-th generator of O times the j-th generator of
241 * I, expressed as combinations of the generators of I.
242 */
243 unsigned int n = M->m;
244 ASSERT_ALWAYS(M->n == n * n);
245 ASSERT_ALWAYS(I->m == n);
246 ASSERT_ALWAYS(I->n == n);
247
248 cxx_mpq_mat R;
249 mpq_mat_inv(R, cxx_mpq_mat(I));
250 cxx_mpz_mat MI(n, n * n);
251 for(unsigned int i = 0 ; i < n ; i++) {
252 cxx_mpz_mat Tz(n, n);
253 for (unsigned int j = 0; j < n; j++) {
254 mpz_mat_submat_set(Tz, j, 0, M, i, j*n, 1, n);
255 }
256 /* We have the matrix of multiplication by w (some generator of
257 * O). Row i is w*wi.
258 *
259 * We need to compute I*T*I^-1 in order to have that represent
260 * how multiplication by w affects the generators of I.
261 */
262 mpz_mat_mul(Tz, I, Tz);
263 cxx_mpq_mat T(Tz);
264 mpq_mat_mul(T, T, R);
265 int rc = mpq_mat_numden(Tz, NULL, T);
266 ASSERT_ALWAYS(rc);
267 for (unsigned int j = 0; j < n; j++) {
268 mpz_mat_submat_swap(MI, i, j*n, Tz, j, 0, 1, n);
269 }
270 }
271 return MI;
272 }/*}}}*/
273
274 /*{{{ multiply_elements_in_order */
multiply_elements_in_order(cxx_mpz_mat const & M,cxx_mpz_mat const & E,cxx_mpz_mat const & F)275 cxx_mpz_mat multiply_elements_in_order(cxx_mpz_mat const& M, cxx_mpz_mat const& E, cxx_mpz_mat const& F)
276 {
277 /* Let O be an order in a degree n number field.
278 *
279 * Given the n*n^2 matrix M of the multiplication within that order,
280 * given matrices E and F of size k*n with coordinates of k elements
281 * of O, return the matrix of size k*n with i-th row giving
282 * coordinates of e_i times f_i in O.
283 */
284 unsigned int n = M->m;
285 unsigned int k = E->m;
286 ASSERT_ALWAYS(M->n == n * n);
287 ASSERT_ALWAYS(E->n == n);
288 ASSERT_ALWAYS(F->n == n);
289 ASSERT_ALWAYS(F->m == k);
290
291 cxx_mpz_mat EM;
292 mpz_mat_mul(EM, E, M);
293 cxx_mpz_mat R(k, n);
294 for(unsigned int ell = 0 ; ell < k ; ell++) {
295 for (unsigned int j = 0; j < n; j++) {
296 mpz_ptr Rlj = mpz_mat_entry(R, ell, j);
297 for(unsigned int i = 0 ; i < n ; i++) {
298 mpz_addmul(Rlj,
299 mpz_mat_entry_const(F, ell, i),
300 mpz_mat_entry_const(EM, ell, i*n+j));
301 }
302 }
303 }
304 return R;
305 }/*}}}*/
306
307 //{{{ frobenius_matrix
frobenius_matrix(cxx_mpz_mat const & M,cxx_mpz const & p)308 cxx_mpz_mat frobenius_matrix(cxx_mpz_mat const& M, cxx_mpz const& p)
309 {
310 // frobenius_matrix ; utility function for computing the p-radical
311 // Takes a matrix B containing the generators (w_0, ... w_{n-1}) of an
312 // order, expressed as polynomials in the root of the polynomial g,
313 // return the matrix U containing ((w_0)^p, ..., (w_{n-1})^p), and
314 // expressed as a linear transformation within the order B.
315
316 /* This version uses the multiplication table of the order (and does all
317 * arithmetic mod p), and binary powering. */
318 unsigned int n = M->m;
319 ASSERT_ALWAYS(M->n == n * n);
320 cxx_mpz_mat Mp = M;
321 mpz_mat_mod_mpz(Mp, Mp, p);
322
323 cxx_mpz_mat E(n, n, 1), F(E);
324
325 int k = mpz_sizeinbase(p, 2) - 1;
326 for( ; k-- ; ) {
327 F = multiply_elements_in_order(M, F, F);
328 if (mpz_tstbit(p, k))
329 F = multiply_elements_in_order(M, E, F);
330 mpz_mat_mod_mpz(F, F, p);
331 }
332 return F;
333 }//}}}
334
335 // {{{ cxx_mpz_mat p_radical_of_order
336 // Stores in I the p-radical of the order whose multiplication matrix is
337 // given by M. I is expressed with respect to the basis of the order.
p_radical_of_order(cxx_mpz_mat const & M,cxx_mpz const & p)338 cxx_mpz_mat p_radical_of_order(cxx_mpz_mat const& M, cxx_mpz const& p)
339 {
340 unsigned int n = M->m;
341 ASSERT_ALWAYS(M->n == n * n);
342
343 cxx_mpz_mat K;
344
345 // Now building the matrix U, containing all generators to the power
346 // of p Generators are polynomials, stored in the matrix B
347 cxx_mpz_mat T = frobenius_matrix(M, p);
348
349 /* Which power of p ? */
350 int k = 1;
351 for (cxx_mpz pk = p; mpz_cmp_ui(pk, n) < 0; k++)
352 mpz_mul(pk, pk, p);
353 mpz_mat_pow_ui_mod_mpz(T, T, k, p);
354
355 // Storing in K a basis of Ker((z -> (z^%p mod g mod p))^k)
356 mpz_mat_kernel_mod_mpz(K, T, p);
357
358 // Getting generators of the p radical from Ker(X) by computing HNF
359 // of the vertical block matrix (p*Id, K);
360 return join_HNF(K, p);
361 }// }}}
362
363 // {{{ cxx_mpq_mat p_maximal_order(cxx_mpz_poly const& f, cxx_mpz const& p)
364 //
365 // This functions runs the round-2 algorithm, starting from the trivial
366 // order defined by f (polynomial order generated by alpha_hat, which is
367 // leadingcoefficient(f)*alpha, where alpha is a root of f).
368 //
369 // The returned order is expressed with respect to the polynomial basis
370 // defined by f.
p_maximal_order(cxx_mpz_poly const & f,cxx_mpz const & p)371 cxx_mpq_mat p_maximal_order(cxx_mpz_poly const& f, cxx_mpz const& p)
372 {
373 unsigned int n = f->deg;
374
375 cxx_mpz_poly g;
376 mpz_poly_to_monic(g, f);
377 cxx_mpq_mat B(n, n, 1);
378
379 cxx_mpq_mat new_D = B;
380 cxx_mpq_mat D = B;
381
382 do {
383 D = new_D;
384
385 cxx_mpz_mat M = multiplication_table_of_order(D, g);
386
387 // Getting the p-radical of the order generated by D in I.
388 cxx_mpz_mat I = p_radical_of_order(M, p);
389
390 // Building the (n,n^2) matrix containing the integers mod p
391 // associated to all generators of O
392 cxx_mpz_mat M2 = multiplication_table_of_ideal(M, I);
393 mpz_mat_mod_mpz(M2, M2, p);
394 //printf("M is :\n"); mpz_mat_fprint(stdout, M); printf("\n");
395
396 // Computing Ker(M)
397 cxx_mpz_mat K; // Ker(M)
398 mpz_mat_kernel_mod_mpz(K, M2, p);
399 // printf("Ker(M) is :\n"); mpz_mat_fprint(stdout, K_M); printf("\n");
400
401 // Getting generators of p*O' by computing HNF of the vertical block matrix (p*Id, K_M);
402 cxx_mpz_mat J = join_HNF(K, p);
403
404 // Converting in the basis which is used to express elements of D
405 mpq_mat_mul(new_D, cxx_mpq_mat(J), D);
406
407 // Dividing by p
408 mpq_mat_div_mpz(new_D, new_D, p);
409
410 } while (D != new_D);
411
412 /* express D back in the basis of f */
413 cxx_mpz x;
414 mpz_set_ui(x, 1);
415 for(unsigned int j = 0 ; j < n ; j++) {
416 for(unsigned int i = 0; i < n; i++) {
417 mpq_ptr dij = mpq_mat_entry(D, i, j);
418 mpz_mul(mpq_numref(dij), mpq_numref(dij), x);
419 mpq_canonicalize(dij);
420 }
421 mpz_mul(x, x, f->coeff[f->deg]);
422 }
423 // Put D into HNF.
424 cxx_mpz_mat Dz;
425 cxx_mpz den;
426 mpq_mat_numden(Dz, den, D);
427 mpz_mat_hnf_backend_rev(Dz, NULL);
428 mpq_mat_set_mpz_mat_denom(D, Dz, den);
429 return D;
430 }
431 //}}}
432
433 //{{{ mpz_mat_minpoly_mod_ui
mpz_mat_minpoly_mod_mpz(cxx_mpz_mat M,cxx_mpz const & p)434 cxx_mpz_poly mpz_mat_minpoly_mod_mpz(cxx_mpz_mat M, cxx_mpz const& p)
435 {
436 ASSERT_ALWAYS(M->m == M->n);
437 unsigned long n = M->n;
438 cxx_mpz_mat B(n+1, n*n);
439 cxx_mpz_mat T(n,n,1);
440 for(unsigned int k = 0 ; k < n + 1 ; k++) {
441 for(unsigned int i = 0 ; i < n ; i++) {
442 mpz_mat_submat_set(B, n - k, i * n, T, i, 0, 1, n);
443 }
444 mpz_mat_mul_mod_mpz(T, T, M, p);
445 }
446 cxx_mpz_mat K;
447 mpz_mat_kernel_mod_mpz(K, B, p);
448 mpz_mat_gauss_backend_mod_mpz(K, NULL, p);
449
450 /* TODO: write down exactly what we need in terms of ordering from
451 * mpz_mat_gauss_backend_mod_ui. I take it that we're happy if it
452 * computes a RREF, so that the last row is the one with the largest
453 * number of column coordinates killed, corresponding to highest
454 * degree coefficients absent ?
455 */
456 cxx_mpz_poly f;
457 mpz_mat_row_to_poly_rev(f,K,K->m-1);
458
459 return f;
460 }
461 /*}}}*/
462
463 /*{{{ matrix_of_multmap */
matrix_of_multmap(cxx_mpz_mat const & M,cxx_mpz_mat const & J,cxx_mpz_mat const & c,cxx_mpz const & p)464 cxx_mpz_mat matrix_of_multmap(
465 cxx_mpz_mat const& M,
466 cxx_mpz_mat const& J,
467 cxx_mpz_mat const& c,
468 cxx_mpz const& p)
469 {
470 /* Let O be an order, represented here simply by its multiplication
471 * matrix. Let J be an m-dimensional subalgebra of O/pO. Let c be an
472 * element of J (a linear combination of the rows of J).
473 * This returns the m times m matrix which expresses multiplication by c
474 * within J.
475 */
476 /* We write the coordinates in O/pO of all products of the for c*j_k,
477 * j_k being one of the generators of J, and denote this matrix by
478 * CJ. Now we look for a matrix M such that M*J = CJ. To do this, we
479 * need to do gaussian reduction on transpose(J): Find a matrix K
480 * such that K*transpose(J) is verticaljoin(identity_matrix(m),
481 * zero_matrix(n-m,m)). Then, taking K' as the first m rows of K, we
482 * will have M = CJ * transpose(K').
483 */
484 unsigned int m = J->m;
485 unsigned int n = M->m;
486 ASSERT_ALWAYS(M->n == n*n);
487 ASSERT_ALWAYS(J->n == n);
488 ASSERT_ALWAYS(c->n == n);
489 ASSERT_ALWAYS(c->m == 1);
490 cxx_mpz_mat Mc;
491 cxx_mpz_mat Jt;
492 cxx_mpz_mat K;
493 mpz_mat_transpose(Jt, J);
494 cxx_mpz_mat T;
495 mpz_mat_gauss_backend_mod_mpz(Jt, T, p);
496 cxx_mpz_mat Kt(m, n);
497 mpz_mat_submat_swap(Kt, 0, 0, T, 0, 0, m, n);
498 mpz_mat_transpose(Kt, Kt); /* Kt is now n times m */
499 cxx_mpz_mat CJ(m, n);
500 for(unsigned int k = 0 ; k < m ; k++) {
501 cxx_mpz_mat jk(1,n);
502 mpz_mat_submat_set(jk,0,0,J,k,0,1,n);
503 cxx_mpz_mat w = multiply_elements_in_order(M, c, jk);
504 mpz_mat_submat_swap(CJ,k,0,w,0,0,1,n);
505 }
506 mpz_mat_mul_mod_mpz(Mc, CJ, Kt, p);
507 return Mc;
508 }
509 /*}}}*/
510
511 /*{{{ factorization_of_polynomial_mod_mpz */
factorization_of_polynomial_mod_mpz(cxx_mpz_poly const & f,cxx_mpz const & p,gmp_randstate_t state)512 vector<pair<cxx_mpz_poly, int> > factorization_of_polynomial_mod_mpz(cxx_mpz_poly const& f, cxx_mpz const& p, gmp_randstate_t state)
513 {
514 mpz_poly_factor_list lf;
515 mpz_poly_factor_list_init(lf);
516 mpz_poly_factor(lf,f,p,state);
517 vector<pair<cxx_mpz_poly, int> > res(lf->size);
518 for(int i = 0 ; i < lf->size ; i++) {
519 mpz_poly_swap(res[i].first, lf->factors[i]->f);
520 res[i].second = lf->factors[i]->m;
521 }
522 mpz_poly_factor_list_clear(lf);
523 return res;
524 }
525 /*}}}*/
526
527 /*{{{ template <typename T> void append_move(vector<T> &a, vector<T> &b) */
append_move(vector<T> & a,vector<T> & b)528 template <typename T> void append_move(vector<T> &a, vector<T> &b)
529 {
530 a.reserve(a.size() + b.size());
531 size_t na = a.size();
532 size_t nb = b.size();
533 if (&a == &b) {
534 for(size_t i = 0 ; i < nb ; i++) {
535 a.push_back(b[i]);
536 }
537 } else {
538 a.insert(a.end(), nb, T());
539 for(size_t i = 0 ; i < nb ; i++) {
540 std::swap(a[na + i], b[i]);
541 }
542 }
543 }
544 /*}}}*/
545
546 // {{{ factorization_of_prime
factorization_of_prime_inner(cxx_mpq_mat const & B,cxx_mpz_mat const & M,cxx_mpz const & p,cxx_mpz_mat const & Ip,cxx_mpz_mat const & I,cxx_mpz_mat const & J,gmp_randstate_t state)547 vector<pair<cxx_mpz_mat, int> > factorization_of_prime_inner(
548 cxx_mpq_mat const & B,
549 cxx_mpz_mat const & M,
550 cxx_mpz const& p,
551 cxx_mpz_mat const& Ip,
552 cxx_mpz_mat const& I,
553 cxx_mpz_mat const& J,
554 gmp_randstate_t state)
555 {
556 unsigned int m = J->m;
557 unsigned int n = J->n;
558 ASSERT_ALWAYS(I->m == n);
559 ASSERT_ALWAYS(I->n == n);
560 /* J represents an m-dimensional subalgebra of O/pO */
561
562 /* Pick a random element of J */
563 cxx_mpz_mat c(1, m);
564 for(unsigned int i = 0 ; i < m ; i++) {
565 mpz_urandomm(mpz_mat_entry(c,0,i), state, p);
566 }
567 mpz_mat_mul_mod_mpz(c, c, J, p);
568
569 cxx_mpz_mat Mc = matrix_of_multmap(M, J, c, p);
570
571 /* Now we would like to find the minimal polynomial of Mc, which is
572 * simply the matrix of multiplication by c. For this, we write an
573 * (m+1) times m^2 matrix, and compute a left kernel.
574 */
575 cxx_mpz_poly Pc = mpz_mat_minpoly_mod_mpz(Mc, p);
576
577 vector<pair<cxx_mpz_poly, int> > facP = factorization_of_polynomial_mod_mpz(Pc, p, state);
578
579 vector<pair<cxx_mpz_mat, int> > ideals;
580
581 vector<cxx_mpz_mat> characteristic_subspaces;
582
583 for(unsigned int i = 0 ; i < facP.size() ; i++) {
584 cxx_mpz_poly const& f(facP[i].first);
585 /* We need the basis of the kernel of f(Mc) */
586 cxx_mpz_mat E;
587 mpz_poly_eval_mpz_mat_mod_mpz(E, Mc, f, p);
588 mpz_mat_pow_ui_mod_mpz(E, E, facP[i].second, p);
589 mpz_mat_kernel_mod_mpz(E, E, p);
590 /* This line is just to be exactly in line with what magma says
591 */
592 mpz_mat_gauss_backend_mod_mpz(E, NULL, p);
593 mpz_mat_mul_mod_mpz(E, E, J, p);
594 characteristic_subspaces.push_back(E);
595 }
596
597 for(unsigned int i = 0 ; i < facP.size() ; i++) {
598 unsigned int e = facP[i].second;
599 cxx_mpz_poly const& f(facP[i].first);
600 cxx_mpz_mat const& Ci(characteristic_subspaces[i]);
601 /* We need to find an ideal which is smaller than Ip, and whose
602 * generators are p*the generators of Ci, as well as the
603 * unmodified generators of I and of the other characteristic
604 * subspaces.
605 */
606 cxx_mpz_mat Ix(n + J->m - Ci->m, n);
607 mpz_mat_submat_set(Ix,0,0,I,0,0,n,n);
608 for(unsigned int r = n, k = 0; k < facP.size() ; k++) {
609 if (k == i) continue;
610 cxx_mpz_mat const& Ck(characteristic_subspaces[k]);
611 unsigned int mk = Ck->m;
612 mpz_mat_submat_set(Ix, r, 0, Ck, 0, 0, mk, n);
613 r += mk;
614 }
615 cxx_mpz_mat Ihead(n, n);
616 if (Ci->m == e * f->deg) {
617 mpz_mat_vertical_join(Ix, Ix, Ip);
618 mpz_mat_hnf_backend_rev(Ix, NULL);
619 mpz_mat_submat_swap(Ihead,0,0,Ix,0,0,n,n);
620 ideals.push_back(make_pair(Ihead, e));
621 } else {
622 mpz_mat_hnf_backend_rev(Ix, NULL);
623 mpz_mat_submat_swap(Ihead,0,0,Ix,0,0,n,n);
624 vector<pair<cxx_mpz_mat, int> > more_ideals;
625 more_ideals = factorization_of_prime_inner(B,M,p,Ip,Ihead,Ci,state);
626 append_move(ideals, more_ideals);
627 }
628 }
629 sort(ideals.begin(), ideals.end(), ideal_comparator());
630 return ideals;
631 }
632
factorization_of_prime(cxx_mpq_mat & B,cxx_mpz_poly const & g,cxx_mpz const & p,gmp_randstate_t state)633 vector<pair<cxx_mpz_mat, int> > factorization_of_prime(
634 cxx_mpq_mat & B, cxx_mpz_poly const& g,
635 cxx_mpz const& p,
636 gmp_randstate_t state)
637 {
638 int n = g->deg;
639 cxx_mpz_mat M = multiplication_table_of_order(B, g);
640 cxx_mpz_mat Ip = p_radical_of_order(M, p);
641 return factorization_of_prime_inner(B, M, p, Ip,
642 cxx_mpz_mat(n, n, p),
643 cxx_mpz_mat(n, n, 1), state);
644 }
645 //}}}
646
647 // {{{ valuation_helper_for_ideal
648 //
649 // compute an uniformizing element for the prime ideal I of the order O,
650 // where I is above the rational prime p, and O is p-maximal and given by
651 // means of its multiplication matrix M.
652 //
653 // The uniformizing element is such that (a/p)*I is in O, yet a is not in
654 // p*O.
valuation_helper_for_ideal(cxx_mpz_mat const & M,cxx_mpz_mat const & I,cxx_mpz const & p)655 cxx_mpz_mat valuation_helper_for_ideal(cxx_mpz_mat const& M, cxx_mpz_mat const& I, cxx_mpz const& p)
656 {
657 unsigned int n = M->m;
658 ASSERT_ALWAYS(M->n == n * n);
659 ASSERT_ALWAYS(I->m == n);
660 ASSERT_ALWAYS(I->n == n);
661
662 /* We begin by something which is very similar to
663 * multiplication_table_of_ideal. It's a bit like computing I*M,
664 * except that we want it in a different order.
665 */
666 cxx_mpz_mat MI(n, n * n);
667 for(unsigned int i = 0 ; i < n ; i++) {
668 cxx_mpz_mat Tz(n, n);
669 for (unsigned int j = 0; j < n; j++) {
670 mpz_mat_submat_set(Tz, j, 0, M, i, j*n, 1, n);
671 }
672 mpz_mat_mul_mod_mpz(Tz, I, Tz, p);
673 for (unsigned int j = 0; j < n; j++) {
674 mpz_mat_submat_swap(MI, i, j*n, Tz, j, 0, 1, n);
675 }
676 }
677
678 cxx_mpz_mat ker;
679 mpz_mat_kernel_mod_mpz(ker, MI, p);
680 mpz_mat_hnf_backend(ker, NULL);
681
682 cxx_mpz_mat res(1, n);
683 mpz_mat_submat_swap(res, 0, 0, ker, 0, 0, 1, n);
684 return res;
685 }// }}}
686
687 // {{{ generate_ideal -- create an ideal from a set of generators
688 // generators (gens) are here given as elements of the field, not
689 // elements of the order. In case the ideal is fractional, its
690 // denominator is also returned. For an integral ideal, the denominator
691 // is always 1.
generate_ideal(cxx_mpq_mat const & O,cxx_mpz_mat const & M,cxx_mpq_mat const & gens)692 pair<cxx_mpz_mat, cxx_mpz> generate_ideal(cxx_mpq_mat const& O, cxx_mpz_mat const& M, cxx_mpq_mat const& gens)
693 {
694 unsigned int n = M->m;
695 ASSERT_ALWAYS(M->n == n * n);
696 ASSERT_ALWAYS(O->m == n);
697 ASSERT_ALWAYS(O->n == n);
698 ASSERT_ALWAYS(gens->n == n);
699 cxx_mpq_mat R;
700 mpq_mat_inv(R, O);
701 cxx_mpq_mat I0q;
702 mpq_mat_mul(I0q, gens, R);
703 cxx_mpz_mat I0;
704 cxx_mpz denom;
705 mpq_mat_numden(I0, denom, I0q);
706 /* Now create the full list of elements of O which are generated by
707 * the products I_i * O_j */
708 cxx_mpz_mat products(gens->m * O->m, n);
709 for(unsigned int i = 0 ; i < I0->m ; i++) {
710 cxx_mpz_mat a(1,n);
711 mpz_mat_submat_set(a,0,0,I0,i,0,1,n);
712 mpz_mat_mul(a,a,M);
713 for(unsigned int j = 0 ; j < n ; j++){
714 mpz_mat_submat_swap(a,0,j*n,products,i*n+j,0,1,n);
715 }
716 }
717 /* And put this in HNF */
718 mpz_mat_hnf_backend_rev(products, NULL);
719 cxx_mpz_mat I(n,n);
720 mpz_mat_submat_swap(I,0,0,products,0,0,n,n);
721 return make_pair(I, denom);
722 }//}}}
723
prime_ideal_inertia_degree(cxx_mpz_mat const & I)724 int prime_ideal_inertia_degree(cxx_mpz_mat const& I)/*{{{*/
725 {
726 unsigned int n = I->m;
727 ASSERT_ALWAYS(I->n == n);
728 int f = 0;
729 for(unsigned int i = 0 ; i < n ; i++) {
730 f += mpz_cmp_ui(mpz_mat_entry_const(I, i, i), 1) != 0;
731 }
732 return f;
733 }
734 /*}}}*/
valuation_of_ideal_at_prime_ideal(cxx_mpz_mat const & M,cxx_mpz_mat const & I,cxx_mpz_mat const & a,cxx_mpz const & p)735 int valuation_of_ideal_at_prime_ideal(cxx_mpz_mat const& M, cxx_mpz_mat const& I, cxx_mpz_mat const& a, cxx_mpz const& p)/*{{{*/
736 {
737 unsigned int n = M->m;
738 ASSERT_ALWAYS(M->n == n * n);
739 ASSERT_ALWAYS(I->m == n);
740 ASSERT_ALWAYS(I->n == n);
741 ASSERT_ALWAYS(a->m == 1);
742 ASSERT_ALWAYS(a->n == n);
743
744 cxx_mpz_mat Ia(I);
745 int val = 0;
746 for(;;val++) {
747 for(unsigned int i = 0 ; i < n ; i++) {
748 cxx_mpz_mat b(1,n);
749 mpz_mat_submat_set(b,0,0,Ia,i,0,1,n);
750 b = multiply_elements_in_order(M, a, b);
751 mpz_mat_submat_swap(b,0,0,Ia,i,0,1,n);
752 }
753 if (mpz_mat_p_valuation(Ia, p) < 1)
754 return val;
755 mpz_mat_divexact_mpz(Ia, Ia, p);
756 }
757 }
758 /*}}}*/
valuation_of_ideal_at_prime_ideal(cxx_mpz_mat const & M,pair<cxx_mpz_mat,cxx_mpz> const & Id,cxx_mpz_mat const & a,int e,cxx_mpz const & p)759 int valuation_of_ideal_at_prime_ideal(cxx_mpz_mat const& M, pair<cxx_mpz_mat,cxx_mpz> const& Id, cxx_mpz_mat const& a, int e, cxx_mpz const& p)/*{{{*/
760 {
761 int w = mpz_p_valuation(Id.second, p);
762 int v = valuation_of_ideal_at_prime_ideal(M, Id.first, a, p);
763 return v-w*e;
764 }
765 /*}}}*/
766
767 struct hypercube_walk {/*{{{*/
768 struct iterator {
769 int B;
770 vector<int> v;
771 vector<int> speed;
772 pair<int, int> last;
operator ++hypercube_walk::iterator773 iterator& operator++() {
774 for(unsigned int j = 0 ; j < v.size() ; j++) {
775 int s = v[j] + speed[j];
776 if (s >= 0 && s <= B) {
777 v[j] = s;
778 last = make_pair(j, speed[j]);
779 return *this;
780 }
781 speed[j]=-speed[j];
782 }
783 last = make_pair(-1, -1);
784 return *this;
785 }
operator !=hypercube_walk::iterator786 inline bool operator!=(iterator const& x) const { return last != x.last; }
operator *hypercube_walk::iterator787 vector<int> & operator*() { return v; }
788 };
789 int n,B;
hypercube_walkhypercube_walk790 hypercube_walk(int n, int B) : n(n), B(B) {}
beginhypercube_walk791 iterator begin() const {
792 iterator z;
793 z.B = B;
794 z.v.assign(n, 0);
795 z.speed.assign(n, 1);
796 return z;
797 }
middlehypercube_walk798 iterator middle(vector<int> const& x) const {
799 ASSERT_ALWAYS(x.size() == (unsigned int) n);
800 iterator z;
801 z.B = B;
802 z.v = x;
803 z.speed.assign(n, 1);
804 return z;
805 }
endhypercube_walk806 iterator end() const {
807 iterator z;
808 z.B = B;
809 z.v.assign(n,0);
810 z.speed.assign(n, 1);
811 z.last = make_pair(-1,-1);
812 return z;
813 }
814 };/*}}}*/
815
816 /* {{{ prime_ideal_two_element */
817 /* I must be in HNF */
prime_ideal_two_element(cxx_mpq_mat const & O,cxx_mpz_poly const & f,cxx_mpz_mat const & M,cxx_mpz_mat const & I)818 pair<cxx_mpz, cxx_mpz_mat> prime_ideal_two_element(cxx_mpq_mat const& O, cxx_mpz_poly const& f, cxx_mpz_mat const& M, cxx_mpz_mat const& I)
819 {
820 cxx_mpz p;
821 unsigned int n = M->m;
822 ASSERT_ALWAYS(M->n == n * n);
823 ASSERT_ALWAYS(I->m == n);
824 ASSERT_ALWAYS(I->n == n);
825 int inertia = 0;
826 for(unsigned int i = 0 ; i < n ; i++) {
827 mpz_srcptr di = mpz_mat_entry_const(I, i, i);
828 if (mpz_cmp_ui(di, 1) != 0) {
829 mpz_set(p, di);
830 inertia++;
831 }
832 }
833 ASSERT_ALWAYS(mpz_size(p) != 0);
834
835 /* We use a homebrew algorithm, loosely inspired on Cohen's. We
836 * prefer an enumeration which favors small generators. Based on
837 * [Cohen93, lemma 4.7.9], we're simply going to search for an
838 * element whose norm has p-valuation exactly f.
839 *
840 * Note though that we're only achieving small element size in terms
841 * of coordinates on the order.
842 */
843 /* There's a question about which generators we should combine in
844 * order to find our winning generator. Cohen takes "the generators
845 * of p*O, and the generators of I". In fact, we can simply take the
846 * generators of I as a Z-basis: that generates the right set. Then,
847 * we'd better work on some ordering for these. Elements in the HNF
848 * form of I with a p on the diagonal have in fact all other elements
849 * on that row equal to zero. That is because since I contains p*O,
850 * the subtraction of this generating row and the p-th multiple of
851 * the corresponding generator has to be zero, or I would not be in
852 * HNF. Therefore, these elements are of secondary importance in
853 * generating I (e.g.: alone, they can't).
854 */
855 unsigned int m = n; /* number of generators we take */
856 cxx_mpz_mat OI(m,n);
857 for(unsigned int i = 0, r = 0, s = 0 ; i < n ; i++) {
858 mpz_srcptr di = mpz_mat_entry_const(I, i, i);
859 if (mpz_cmp_ui(di, 1) != 0) {
860 mpz_mat_submat_set(OI,n - inertia + s++,0,I,i,0,1,n);
861 } else {
862 mpz_mat_submat_set(OI,r++,0,I,i,0,1,n);
863 }
864 }
865
866 cxx_mpq_mat OIOq;
867 mpq_mat_mul(OIOq, cxx_mpq_mat(OI), O);
868 for(int B = 1 ; ; B++) {
869 hypercube_walk H(m, B);
870 vector<int> H0(m,0);
871 cxx_mpq_mat gen(1, n);
872 // H0[n]=1;
873 // mpq_mat_submat_set(gen, 0, 0, OIOq, n, 0, 1, n);
874 cxx_mpq_mat temp(1, n);
875 for(hypercube_walk::iterator it = H.middle(H0) ; it != H.end() ; ++it) {
876 int j = it.last.first;
877 int s = it.last.second;
878 /* add (s=1) or subtract (s=-1) the j-th generator to gen */
879 mpq_mat_submat_swap(temp,0,0,OIOq,j,0,1,n);
880 if (s==1)
881 mpq_mat_add(gen, gen, temp);
882 else if (s==-1)
883 mpq_mat_sub(gen, gen, temp);
884 mpq_mat_submat_swap(temp,0,0,OIOq,j,0,1,n);
885
886 cxx_mpz_poly pgen;
887 cxx_mpz dgen;
888 cxx_mpz res;
889 mpq_mat_row_to_poly(pgen, dgen, gen, 0);
890 if (pgen->deg < 0) continue;
891 mpz_poly_resultant(res, pgen, f);
892 /* We want the absolute norm to have p-valuation equal to the
893 * inertia degree.
894 * absolute norm is galois norm. galois norm is product of
895 * all conjugate of pgen(alpha). Resultant(f,pgen) is
896 * lc(f)^deg(pgen) times the galois norm.
897 */
898 int v = mpz_p_valuation(res, p) - pgen->deg * mpz_p_valuation(f->coeff[f->deg], p) - f->deg * mpz_p_valuation(dgen, p);
899 ASSERT_ALWAYS(v >= inertia);
900 if (v == inertia) {
901 cxx_mpz_mat lambda(1, m);
902 vector<int> & v(*it);
903 for(unsigned int i = 0 ; i < m ; i++) {
904 mpz_set_si(mpz_mat_entry(lambda,0,i),v[i]);
905 }
906 mpz_mat_mul(lambda, lambda, OI);
907 return make_pair(p, lambda);
908 }
909 }
910 }
911 }
912 // }}}
913
write_element_as_polynomial(cxx_mpq_mat const & theta_q,string const & var)914 string write_element_as_polynomial(cxx_mpq_mat const& theta_q, string const& var)
915 {
916 ASSERT_ALWAYS(theta_q->m == 1);
917 cxx_mpz theta_denom;
918 cxx_mpz_poly theta;
919 mpq_mat_row_to_poly(theta, theta_denom, theta_q, 0);
920
921 string num = theta.print_poly(var);
922 /* first write the numerator as a string */
923 if (mpz_cmp_ui(theta_denom, 1) == 0) {
924 return num;
925 } else {
926 ostringstream os2;
927 os2 << "(" << num << ")/" << theta_denom;
928 return os2.str();
929 }
930 }
931