1 #include "cado.h" // IWYU pragma: keep
2 #include <cstdio> // FILE // IWYU pragma: keep
3 #include <cstdlib>
4 #include <cstring>
5 #include <climits> /* for INT_MAX */
6 #include <ostream>
7 #include <gmp.h>
8 #include "mpz_mat.h"
9 #include "gmpxx.hpp"
10 #include "lll.h" // mat_Z, LLL
11 #include "macros.h"
12
13 /*{{{ entry access*/
mpz_mat_entry(mpz_mat_ptr M,unsigned int i,unsigned int j)14 mpz_ptr mpz_mat_entry(mpz_mat_ptr M, unsigned int i, unsigned int j)
15 {
16 return M->x[i * M->n + j];
17 }
18
mpz_mat_entry_const(mpz_mat_srcptr M,unsigned int i,unsigned int j)19 mpz_srcptr mpz_mat_entry_const(mpz_mat_srcptr M, unsigned int i, unsigned int j)
20 {
21 return M->x[i * M->n + j];
22 }
23
mpq_mat_entry(mpq_mat_ptr M,unsigned int i,unsigned int j)24 mpq_ptr mpq_mat_entry(mpq_mat_ptr M, unsigned int i, unsigned int j)
25 {
26 return M->x[i * M->n + j];
27 }
28
mpq_mat_entry_const(mpq_mat_srcptr M,unsigned int i,unsigned int j)29 mpq_srcptr mpq_mat_entry_const(mpq_mat_srcptr M, unsigned int i, unsigned int j)
30 {
31 return M->x[i * M->n + j];
32 }
33 /*}}}*/
34 /*{{{ init/clear/realloc*/
mpz_mat_init(mpz_mat_ptr M,unsigned int m,unsigned int n)35 void mpz_mat_init(mpz_mat_ptr M, unsigned int m, unsigned int n)
36 {
37 M->x = (mpz_t*) ((m && n) ? malloc(m * n * sizeof(mpz_t)) : NULL);
38
39 M->m = m;
40 M->n = n;
41 for(unsigned int i = 0 ; i < M->m ; i++)
42 for(unsigned int j = 0 ; j < M->n ; j++)
43 mpz_init(mpz_mat_entry(M, i, j));
44 }
45
mpz_mat_clear(mpz_mat_ptr M)46 void mpz_mat_clear(mpz_mat_ptr M)
47 {
48 for(unsigned int i = 0 ; i < M->m ; i++)
49 for(unsigned int j = 0 ; j < M->n ; j++)
50 mpz_clear(mpz_mat_entry(M, i, j));
51 free(M->x);
52 }
53
mpz_mat_realloc(mpz_mat_ptr M,unsigned int m,unsigned int n)54 void mpz_mat_realloc(mpz_mat_ptr M, unsigned int m, unsigned int n)
55 {
56 if (M->m == m && M->n == n) return;
57 mpz_mat_clear(M);
58 mpz_mat_init(M, m, n);
59 }
60
mpq_mat_init(mpq_mat_ptr M,unsigned int m,unsigned int n)61 void mpq_mat_init(mpq_mat_ptr M, unsigned int m, unsigned int n)
62 {
63 M->x = (mpq_t*) ((m && n) ? malloc(m * n * sizeof(mpq_t)) : NULL);
64 M->m = m;
65 M->n = n;
66 for(unsigned int i = 0 ; i < M->m ; i++)
67 for(unsigned int j = 0 ; j < M->n ; j++)
68 mpq_init(mpq_mat_entry(M, i, j));
69 }
70
mpq_mat_clear(mpq_mat_ptr M)71 void mpq_mat_clear(mpq_mat_ptr M)
72 {
73 for(unsigned int i = 0 ; i < M->m ; i++)
74 for(unsigned int j = 0 ; j < M->n ; j++)
75 mpq_clear(mpq_mat_entry(M, i, j));
76 free(M->x);
77 }
78
mpq_mat_realloc(mpq_mat_ptr M,unsigned int m,unsigned int n)79 void mpq_mat_realloc(mpq_mat_ptr M, unsigned int m, unsigned int n)
80 {
81 if (M->m == m && M->n == n) return;
82 mpq_mat_clear(M);
83 mpq_mat_init(M, m, n);
84 }
85
86 /*}}}*/
87 /*{{{ operations on submatrices, and swaps*/
mpz_mat_submat_swap(mpz_mat_ptr A0,unsigned int i0,unsigned int j0,mpz_mat_ptr A1,unsigned int i1,unsigned int j1,unsigned int dm,unsigned int dn)88 void mpz_mat_submat_swap(
89 mpz_mat_ptr A0, unsigned int i0, unsigned int j0,
90 mpz_mat_ptr A1, unsigned int i1, unsigned int j1,
91 unsigned int dm, unsigned int dn)
92 {
93 ASSERT_ALWAYS(i0 + dm <= A0->m);
94 ASSERT_ALWAYS(i1 + dm <= A1->m);
95 ASSERT_ALWAYS(j0 + dn <= A0->n);
96 ASSERT_ALWAYS(j1 + dn <= A1->n);
97 for(unsigned int i = 0 ; i < dm ; i++) {
98 for(unsigned int j = 0 ; j < dn ; j++) {
99 mpz_swap(
100 mpz_mat_entry(A0, i0 + i, j0 + j),
101 mpz_mat_entry(A1, i1 + i, j1 + j));
102 }
103 }
104 }
105
mpz_mat_submat_set(mpz_mat_ptr A0,unsigned int i0,unsigned int j0,mpz_mat_srcptr A1,unsigned int i1,unsigned int j1,unsigned int dm,unsigned int dn)106 void mpz_mat_submat_set(
107 mpz_mat_ptr A0, unsigned int i0, unsigned int j0,
108 mpz_mat_srcptr A1, unsigned int i1, unsigned int j1,
109 unsigned int dm, unsigned int dn)
110 {
111 ASSERT_ALWAYS(i0 + dm <= A0->m);
112 ASSERT_ALWAYS(i1 + dm <= A1->m);
113 ASSERT_ALWAYS(j0 + dn <= A0->n);
114 ASSERT_ALWAYS(j1 + dn <= A1->n);
115 for(unsigned int i = 0 ; i < dm ; i++) {
116 for(unsigned int j = 0 ; j < dn ; j++) {
117 mpz_set(
118 mpz_mat_entry(A0, i0 + i, j0 + j),
119 mpz_mat_entry_const(A1, i1 + i, j1 + j));
120 }
121 }
122 }
123
mpq_mat_swap(mpz_mat_ptr A,mpz_mat_ptr B)124 void mpq_mat_swap(mpz_mat_ptr A, mpz_mat_ptr B)
125 {
126 mpz_mat C;
127 memcpy(C, A, sizeof(mpz_mat));
128 memcpy(A, B, sizeof(mpz_mat));
129 memcpy(B, C, sizeof(mpz_mat));
130 }
131
mpq_mat_submat_swap(mpq_mat_ptr A0,unsigned int i0,unsigned int j0,mpq_mat_ptr A1,unsigned int i1,unsigned int j1,unsigned int dm,unsigned int dn)132 void mpq_mat_submat_swap(
133 mpq_mat_ptr A0, unsigned int i0, unsigned int j0,
134 mpq_mat_ptr A1, unsigned int i1, unsigned int j1,
135 unsigned int dm, unsigned int dn)
136 {
137 ASSERT_ALWAYS(i0 + dm <= A0->m);
138 ASSERT_ALWAYS(i1 + dm <= A1->m);
139 ASSERT_ALWAYS(j0 + dn <= A0->n);
140 ASSERT_ALWAYS(j1 + dn <= A1->n);
141 for(unsigned int i = 0 ; i < dm ; i++) {
142 for(unsigned int j = 0 ; j < dn ; j++) {
143 mpq_swap(
144 mpq_mat_entry(A0, i0 + i, j0 + j),
145 mpq_mat_entry(A1, i1 + i, j1 + j));
146 }
147 }
148 }
149
mpq_mat_submat_set(mpq_mat_ptr A0,unsigned int i0,unsigned int j0,mpq_mat_srcptr A1,unsigned int i1,unsigned int j1,unsigned int dm,unsigned int dn)150 void mpq_mat_submat_set(
151 mpq_mat_ptr A0, unsigned int i0, unsigned int j0,
152 mpq_mat_srcptr A1, unsigned int i1, unsigned int j1,
153 unsigned int dm, unsigned int dn)
154 {
155 ASSERT_ALWAYS(i0 + dm <= A0->m);
156 ASSERT_ALWAYS(i1 + dm <= A1->m);
157 ASSERT_ALWAYS(j0 + dn <= A0->n);
158 ASSERT_ALWAYS(j1 + dn <= A1->n);
159 for(unsigned int i = 0 ; i < dm ; i++) {
160 for(unsigned int j = 0 ; j < dn ; j++) {
161 mpq_set(
162 mpq_mat_entry(A0, i0 + i, j0 + j),
163 mpq_mat_entry_const(A1, i1 + i, j1 + j));
164 }
165 }
166 }
167
mpz_mat_swap(mpz_mat_ptr A,mpz_mat_ptr B)168 void mpz_mat_swap(mpz_mat_ptr A, mpz_mat_ptr B)
169 {
170 mpz_mat C;
171 memcpy(C, A, sizeof(mpz_mat));
172 memcpy(A, B, sizeof(mpz_mat));
173 memcpy(B, C, sizeof(mpz_mat));
174 }
175
mpq_mat_swap(mpq_mat_ptr A,mpq_mat_ptr B)176 void mpq_mat_swap(mpq_mat_ptr A, mpq_mat_ptr B)
177 {
178 mpq_mat C;
179 memcpy(C, A, sizeof(mpq_mat));
180 memcpy(A, B, sizeof(mpq_mat));
181 memcpy(B, C, sizeof(mpq_mat));
182 }
183
184 /*}}}*/
185 /*{{{ set/set_ui*/
mpz_mat_set(mpz_mat_ptr dst,mpz_mat_srcptr src)186 void mpz_mat_set(mpz_mat_ptr dst, mpz_mat_srcptr src)
187 {
188 mpz_mat_realloc(dst, src->m, src->n);
189 for(unsigned int i = 0 ; i < src->m ; i++)
190 for(unsigned int j = 0 ; j < src->n ; j++)
191 mpz_set(mpz_mat_entry(dst, i, j), mpz_mat_entry_const(src, i, j));
192 }
193
mpz_mat_set_ui(mpz_mat_ptr M,unsigned long a)194 void mpz_mat_set_ui(mpz_mat_ptr M, unsigned long a)
195 {
196 if (a) {
197 ASSERT_ALWAYS(M->m == M->n);
198 }
199 for(unsigned int i = 0 ; i < M->m ; i++)
200 for(unsigned int j = 0 ; j < M->n ; j++)
201 mpz_set_ui(mpz_mat_entry(M, i, j), i == j ? a : 0);
202 }
203
mpz_mat_set_mpz(mpz_mat_ptr M,mpz_srcptr a)204 void mpz_mat_set_mpz(mpz_mat_ptr M, mpz_srcptr a)
205 {
206 ASSERT_ALWAYS(M->m == M->n);
207 for(unsigned int i = 0 ; i < M->m ; i++) {
208 mpz_set(mpz_mat_entry(M, i, i), a);
209 }
210 }
211
mpz_mat_add_ui(mpz_mat_ptr M,unsigned long a)212 void mpz_mat_add_ui(mpz_mat_ptr M, unsigned long a)
213 {
214 ASSERT_ALWAYS(M->m == M->n);
215 for(unsigned int i = 0 ; i < M->m ; i++) {
216 mpz_ptr mii = mpz_mat_entry(M, i, i);
217 mpz_add_ui(mii, mii, a);
218 }
219 }
220
mpz_mat_add_mpz(mpz_mat_ptr M,mpz_srcptr a)221 void mpz_mat_add_mpz(mpz_mat_ptr M, mpz_srcptr a)
222 {
223 ASSERT_ALWAYS(M->m == M->n);
224 for(unsigned int i = 0 ; i < M->m ; i++) {
225 mpz_ptr mii = mpz_mat_entry(M, i, i);
226 mpz_add(mii, mii, a);
227 }
228 }
229
mpq_mat_set_ui(mpq_mat_ptr M,unsigned long a)230 void mpq_mat_set_ui(mpq_mat_ptr M, unsigned long a)
231 {
232 if (a) {
233 ASSERT_ALWAYS(M->m == M->n);
234 }
235 for(unsigned int i = 0 ; i < M->m ; i++)
236 for(unsigned int j = 0 ; j < M->n ; j++)
237 mpq_set_ui(mpq_mat_entry(M, i, j), i == j ? a : 0, 1);
238 }
239
mpq_mat_set(mpq_mat_ptr dst,mpq_mat_srcptr src)240 void mpq_mat_set(mpq_mat_ptr dst, mpq_mat_srcptr src)
241 {
242 mpq_mat_realloc(dst, src->m, src->n);
243 for(unsigned int i = 0 ; i < src->m ; i++)
244 for(unsigned int j = 0 ; j < src->n ; j++)
245 mpq_set(mpq_mat_entry(dst, i, j), mpq_mat_entry_const(src, i, j));
246 }
247
mpz_mat_urandomm(mpz_mat_ptr M,gmp_randstate_t state,mpz_srcptr p)248 void mpz_mat_urandomm(mpz_mat_ptr M, gmp_randstate_t state, mpz_srcptr p)
249 {
250 for(unsigned int i = 0 ; i < M->m ; i++)
251 for(unsigned int j = 0 ; j < M->n ; j++)
252 mpz_urandomm(mpz_mat_entry(M, i, j), state, p);
253 }
254
mpq_mat_urandomm(mpq_mat_ptr M,gmp_randstate_t state,mpz_srcptr p)255 void mpq_mat_urandomm(mpq_mat_ptr M, gmp_randstate_t state, mpz_srcptr p)
256 {
257 for(unsigned int i = 0 ; i < M->m ; i++)
258 for(unsigned int j = 0 ; j < M->n ; j++) {
259 mpz_urandomm(mpq_numref(mpq_mat_entry(M, i, j)), state, p);
260 mpz_set_ui(mpq_denref(mpq_mat_entry(M, i, j)), 1);
261 }
262 }
263
264
265
266 /*}}}*/
267 /*{{{ Joining matrices */
mpz_mat_vertical_join(mpz_mat_ptr N,mpz_mat_srcptr M1,mpz_mat_srcptr M2)268 void mpz_mat_vertical_join(mpz_mat_ptr N, mpz_mat_srcptr M1, mpz_mat_srcptr M2)
269 {
270 if (N == M1 || N == M2) {
271 mpz_mat Nx;
272 mpz_mat_init(Nx, 0, 0);
273 mpz_mat_vertical_join(Nx, M1, M2);
274 mpz_mat_swap(Nx, N);
275 mpz_mat_clear(Nx);
276 return;
277 }
278 ASSERT_ALWAYS(M1->n == M2->n);
279 mpz_mat_realloc(N,M1->m+M2->m,M1->n);
280 mpz_mat_submat_set(N,0,0,M1,0,0,M1->m,M1->n);
281 mpz_mat_submat_set(N,M1->m,0,M2,0,0,M2->m,M2->n);
282 }
mpq_mat_vertical_join(mpq_mat_ptr N,mpq_mat_srcptr M1,mpq_mat_srcptr M2)283 void mpq_mat_vertical_join(mpq_mat_ptr N, mpq_mat_srcptr M1, mpq_mat_srcptr M2)
284 {
285 if (N == M1 || N == M2) {
286 mpq_mat Nx;
287 mpq_mat_init(Nx, 0, 0);
288 mpq_mat_vertical_join(Nx, M1, M2);
289 mpq_mat_swap(Nx, N);
290 mpq_mat_clear(Nx);
291 return;
292 }
293 ASSERT_ALWAYS(M1->n == M2->n);
294 mpq_mat_realloc(N,M1->m+M2->m,M1->n);
295 mpq_mat_submat_set(N,0,0,M1,0,0,M1->m,M1->n);
296 mpq_mat_submat_set(N,M1->m,0,M2,0,0,M2->m,M2->n);
297 }
mpz_mat_horizontal_join(mpz_mat_ptr N,mpz_mat_srcptr M1,mpz_mat_srcptr M2)298 void mpz_mat_horizontal_join(mpz_mat_ptr N, mpz_mat_srcptr M1, mpz_mat_srcptr M2)
299 {
300 if (N == M1 || N == M2) {
301 mpz_mat Nx;
302 mpz_mat_init(Nx, 0, 0);
303 mpz_mat_horizontal_join(Nx, M1, M2);
304 mpz_mat_swap(Nx, N);
305 mpz_mat_clear(Nx);
306 return;
307 }
308 ASSERT_ALWAYS(M1->m == M2->m);
309 mpz_mat_realloc(N,M1->m,M1->n + M2->n);
310 mpz_mat_submat_set(N,0,0,M1,0,0,M1->m,M1->n);
311 mpz_mat_submat_set(N,0,M1->n,M2,0,0,M2->m,M2->n);
312 }
mpq_mat_horizontal_join(mpq_mat_ptr N,mpq_mat_srcptr M1,mpq_mat_srcptr M2)313 void mpq_mat_horizontal_join(mpq_mat_ptr N, mpq_mat_srcptr M1, mpq_mat_srcptr M2)
314 {
315 if (N == M1 || N == M2) {
316 mpq_mat Nx;
317 mpq_mat_init(Nx, 0, 0);
318 mpq_mat_horizontal_join(Nx, M1, M2);
319 mpq_mat_swap(Nx, N);
320 mpq_mat_clear(Nx);
321 return;
322 }
323 ASSERT_ALWAYS(M1->m == M2->m);
324 mpq_mat_realloc(N,M1->m,M1->n + M2->n);
325 mpq_mat_submat_set(N,0,0,M1,0,0,M1->m,M1->n);
326 mpq_mat_submat_set(N,0,M1->n,M2,0,0,M2->m,M2->n);
327 }
328 /*}}}*/
329 /*{{{ determinant, trace and transposition */
mpz_mat_trace(mpz_ptr t,mpz_mat_srcptr M)330 void mpz_mat_trace(mpz_ptr t, mpz_mat_srcptr M)/*{{{*/
331 {
332 ASSERT_ALWAYS(M->m == M->n);
333 mpz_set_ui(t, 0);
334 for(unsigned int i = 0 ; i < M->n ; i++)
335 mpz_add(t, t, mpz_mat_entry_const(M,i,i));
336 }
337 /*}}}*/
mpz_mat_determinant_triangular(mpz_ptr d,mpz_mat_srcptr M)338 void mpz_mat_determinant_triangular(mpz_ptr d, mpz_mat_srcptr M)/*{{{*/
339 {
340 // We assume that M is triangular
341 ASSERT_ALWAYS(M->m == M->n);
342 mpz_set_ui(d, 1);
343 for(unsigned int i = 0 ; i < M->n ; i++)
344 mpz_mul(d, d, mpz_mat_entry_const(M,i,i));
345 }
346 /*}}}*/
mpz_mat_transpose(mpz_mat_ptr D,mpz_mat_srcptr M)347 void mpz_mat_transpose(mpz_mat_ptr D, mpz_mat_srcptr M)/*{{{*/
348 {
349 if (D != M) {
350 mpz_mat_realloc(D,M->n,M->m);
351 for(unsigned int i = 0 ; i < M->m ; i++){
352 for(unsigned int j = 0 ; j < M->n ; j++){
353 mpz_srcptr mij = mpz_mat_entry_const(M,i,j);
354 mpz_ptr dji = mpz_mat_entry(D,j,i);
355 mpz_set(dji, mij);
356 }
357 }
358 } else if (M->m != M->n) {
359 /* transpose a rectangular matrix in place. Rather annoying to do
360 * with real swaps, right ? */
361 mpz_mat Mc;
362 mpz_mat_init(Mc, 0, 0);
363 mpz_mat_set(Mc, M);
364 mpz_mat_transpose(D, Mc);
365 mpz_mat_clear(Mc);
366 } else {
367 for(unsigned int i = 0 ; i < M->m ; i++){
368 for(unsigned int j = i + 1 ; j < M->n ; j++){
369 mpz_ptr dij = mpz_mat_entry(D,i,j);
370 mpz_ptr dji = mpz_mat_entry(D,j,i);
371 mpz_swap(dji, dij);
372 }
373 }
374 }
375 }/*}}}*/
376
mpz_mat_reverse_rows(mpz_mat_ptr B,mpz_mat_srcptr A)377 void mpz_mat_reverse_rows(mpz_mat_ptr B, mpz_mat_srcptr A)
378 {
379 if (A != B) {
380 mpz_mat_realloc(B, A->m,A->n);
381 for(unsigned int i = 0 ; i < A->m ; i++){
382 mpz_mat_submat_set(B, i, 0, A, A->m - 1 - i, 0, 1, A->n);
383 }
384 } else {
385 for(unsigned int i = 0 ; i < A->m - 1 - i ; i++){
386 mpz_mat_submat_swap(B, i, 0, B, A->m - 1 - i, 0, 1, A->n);
387 }
388 }
389 }
390
mpz_mat_reverse_columns(mpz_mat_ptr B,mpz_mat_srcptr A)391 void mpz_mat_reverse_columns(mpz_mat_ptr B, mpz_mat_srcptr A)
392 {
393 if (A != B) {
394 mpz_mat_realloc(B, A->m,A->n);
395 for(unsigned int j = 0 ; j < A->n ; j++){
396 mpz_mat_submat_set(B, 0, j, A, 0, A->n - 1 - j, A->m, 1);
397 }
398 } else {
399 for(unsigned int j = 0 ; j < A->n - 1 - j ; j++){
400 mpz_mat_submat_swap(B, 0, j, B, 0, A->n - 1 - j, A->m, 1);
401 }
402 }
403 }
404
mpq_mat_trace(mpq_ptr t,mpq_mat_srcptr M)405 void mpq_mat_trace(mpq_ptr t, mpq_mat_srcptr M)/*{{{*/
406 {
407 ASSERT_ALWAYS(M->m == M->n);
408 mpq_set_ui(t, 0, 1);
409 for(unsigned int i = 0 ; i < M->n ; i++)
410 mpq_add(t, t, mpq_mat_entry_const(M,i,i));
411 }
412 /*}}}*/
mpq_mat_determinant_triangular(mpq_ptr d,mpq_mat_srcptr M)413 void mpq_mat_determinant_triangular(mpq_ptr d, mpq_mat_srcptr M)/*{{{*/
414 {
415 // We assume that M is triangular
416 ASSERT_ALWAYS(M->m == M->n);
417 mpq_set_ui(d, 1, 1);
418 for(unsigned int i = 0 ; i < M->n ; i++)
419 mpq_mul(d, d, mpq_mat_entry_const(M,i,i));
420 }
421 /*}}}*/
mpq_mat_transpose(mpq_mat_ptr D,mpq_mat_srcptr M)422 void mpq_mat_transpose(mpq_mat_ptr D, mpq_mat_srcptr M)/*{{{*/
423 {
424 if (D != M) {
425 mpq_mat_realloc(D,M->n,M->m);
426 for(unsigned int i = 0 ; i < M->m ; i++){
427 for(unsigned int j = 0 ; j < M->n ; j++){
428 mpq_srcptr mij = mpq_mat_entry_const(M,i,j);
429 mpq_ptr dji = mpq_mat_entry(D,j,i);
430 mpq_set(dji, mij);
431 }
432 }
433 } else if (M->m != M->n) {
434 /* transpose a rectangular matrix in place. Rather annoying to do
435 * with real swaps, right ? */
436 mpq_mat Mc;
437 mpq_mat_init(Mc, 0, 0);
438 mpq_mat_set(Mc, M);
439 mpq_mat_transpose(D, Mc);
440 mpq_mat_clear(Mc);
441 } else {
442 for(unsigned int i = 0 ; i < M->m ; i++){
443 for(unsigned int j = i + 1 ; j < M->n ; j++){
444 mpq_ptr dij = mpq_mat_entry(D,i,j);
445 mpq_ptr dji = mpq_mat_entry(D,j,i);
446 mpq_swap(dji, dij);
447 }
448 }
449 }
450 }/*}}}*/
451 /*}}}*/
452 /*{{{ miscellaneous */
453
454 /* convert to integer matrix divided by lcm of denominator.
455 * return 1
456 * if den==NULL, return 0 if denominator happens to not be 1 (in which
457 * case the matrix returned is undefined).
458 */
459
mpq_mat_numden(mpz_mat_ptr num,mpz_ptr den,mpq_mat_srcptr M)460 int mpq_mat_numden(mpz_mat_ptr num, mpz_ptr den, mpq_mat_srcptr M)/*{{{*/
461 {
462 int ret = 1;
463 if (den) mpz_set_ui(den, 1);
464 for(unsigned int i = 0 ; i < M->m ; i++) {
465 for(unsigned int j = 0 ; j < M->n ; j++) {
466 mpz_srcptr Mij = mpq_denref(mpq_mat_entry_const(M, i, j));
467 if (den) {
468 mpz_lcm(den, den, Mij);
469 } else if (mpz_cmp_ui(Mij, 1) != 0) {
470 ret = 0;
471 }
472 }
473 }
474 if (!ret) return ret;
475 mpz_mat_realloc(num, M->m, M->n);
476 for(unsigned int i = 0 ; i < M->m ; i++) {
477 for(unsigned int j = 0 ; j < M->n ; j++) {
478 mpz_ptr dst = mpz_mat_entry(num, i, j);
479 mpq_srcptr src = mpq_mat_entry_const(M, i, j);
480 if (den) {
481 mpz_divexact(dst, den, mpq_denref(src));
482 mpz_mul(dst, dst, mpq_numref(src));
483 } else {
484 mpz_set(dst, mpq_numref(src));
485 }
486 }
487 }
488 return ret;
489 }
490 /*}}}*/
mpq_mat_set_mpz_mat(mpq_mat_ptr N,mpz_mat_srcptr M)491 void mpq_mat_set_mpz_mat(mpq_mat_ptr N, mpz_mat_srcptr M)/*{{{*/
492 {
493 mpq_mat_realloc(N,M->m,M->n);
494 for(unsigned int i = 0 ; i < M->m ; i++) {
495 for(unsigned int j = 0 ; j < M->n ; j++) {
496 mpq_set_z(mpq_mat_entry(N,i,j),mpz_mat_entry_const(M,i,j));
497 }
498 }
499 }
500 /*}}}*/
mpq_mat_set_mpz_mat_denom(mpq_mat_ptr N,mpz_mat_srcptr M,mpz_srcptr d)501 void mpq_mat_set_mpz_mat_denom(mpq_mat_ptr N, mpz_mat_srcptr M, mpz_srcptr d)/*{{{*/
502 {
503 mpq_mat_realloc(N,M->m,M->n);
504 for(unsigned int i = 0 ; i < M->m ; i++) {
505 for(unsigned int j = 0 ; j < M->n ; j++) {
506 mpq_ptr nij = mpq_mat_entry(N,i,j);
507 mpz_srcptr mij = mpz_mat_entry_const(M,i,j);
508 mpz_set(mpq_numref(nij), mij);
509 mpz_set(mpq_denref(nij), d);
510 mpq_canonicalize(nij);
511 }
512 }
513 }
514 /*}}}*/
mpz_mat_mod_ui(mpz_mat_ptr dst,mpz_mat_srcptr src,unsigned long p)515 void mpz_mat_mod_ui(mpz_mat_ptr dst, mpz_mat_srcptr src, unsigned long p)/*{{{*/
516 {
517 mpz_mat_realloc(dst, src->m, src->n); /* no-op if dst == src */
518 for (unsigned int i = 0 ; i < src->m ; i++){
519 for (unsigned int j = 0 ; j < src->n ; j++){
520 mpz_fdiv_r_ui(mpz_mat_entry(dst,i,j),mpz_mat_entry_const(src,i,j),p);
521 }
522 }
523 }/*}}}*/
mpz_mat_mod_mpz(mpz_mat_ptr dst,mpz_mat_srcptr src,mpz_srcptr p)524 void mpz_mat_mod_mpz(mpz_mat_ptr dst, mpz_mat_srcptr src, mpz_srcptr p)/*{{{*/
525 {
526 mpz_mat_realloc(dst, src->m, src->n); /* no-op if dst == src */
527 for (unsigned int i = 0 ; i < src->m ; i++){
528 for (unsigned int j = 0 ; j < src->n ; j++){
529 mpz_fdiv_r(mpz_mat_entry(dst,i,j),mpz_mat_entry_const(src,i,j),p);
530 }
531 }
532 }/*}}}*/
533
mpz_mat_p_valuation(mpz_mat_srcptr A,mpz_srcptr p)534 int mpz_mat_p_valuation(mpz_mat_srcptr A, mpz_srcptr p)
535 {
536 int val = INT_MAX;
537 mpz_t c;
538 mpz_init(c);
539 for (unsigned int i = 0 ; val && i < A->m ; i++){
540 for (unsigned int j = 0 ; val && j < A->n ; j++){
541 int v = 0;
542 mpz_srcptr aij = mpz_mat_entry_const(A, i, j);
543 if (mpz_size(aij) == 0) continue;
544 mpz_set(c, aij);
545 for( ; v < val && mpz_divisible_p(c, p) ; v++)
546 mpz_fdiv_q(c, c, p);
547 val = v;
548 }
549 }
550 mpz_clear(c);
551 return val;
552 }
553
mpz_mat_p_valuation_ui(mpz_mat_srcptr A,unsigned long p)554 int mpz_mat_p_valuation_ui(mpz_mat_srcptr A, unsigned long p)
555 {
556 int val = INT_MAX;
557 mpz_t c;
558 mpz_init(c);
559 for (unsigned int i = 0 ; val && i < A->m ; i++){
560 for (unsigned int j = 0 ; val && j < A->n ; j++){
561 int v = 0;
562 mpz_srcptr aij = mpz_mat_entry_const(A, i, j);
563 if (mpz_size(aij) == 0) continue;
564 mpz_set(c, aij);
565 for( ; v < val && mpz_divisible_ui_p(c, p) ; v++)
566 mpz_fdiv_q_ui(c, c, p);
567 val = v;
568 }
569 }
570 mpz_clear(c);
571 return val;
572 }
573
574 /*}}}*/
575
576 /* {{{ row-level operations (for Gauss and friends, mostly) */
577 // Return 1 if the k-th line of M is null, 0 else
mpz_mat_isnull_row(mpz_mat_srcptr M,unsigned int k)578 int mpz_mat_isnull_row(mpz_mat_srcptr M, unsigned int k){/*{{{*/
579 unsigned int j = 0;
580 ASSERT_ALWAYS(k < M->m);
581 while((j < M->n) && !mpz_cmp_si(mpz_mat_entry_const(M,k,j),0)){
582 j++;
583 }
584 if(j == M->n){
585 return 1;
586 }
587 return 0;
588 }
589 /*}}}*/
mpz_mat_swaprows(mpz_mat_ptr M,unsigned int i0,unsigned int i1)590 void mpz_mat_swaprows(mpz_mat_ptr M, unsigned int i0, unsigned int i1)/*{{{*/
591 {
592 ASSERT_ALWAYS(i0 < M->m);
593 ASSERT_ALWAYS(i1 < M->m);
594 if (i0 == i1) return;
595 for(unsigned int j = 0 ; j < M->n ; j++) {
596 mpz_swap(mpz_mat_entry(M, i0, j), mpz_mat_entry(M, i1, j));
597 }
598 }
599 /*}}}*/
mpq_mat_swaprows(mpq_mat_ptr M,unsigned int i0,unsigned int i1)600 void mpq_mat_swaprows(mpq_mat_ptr M, unsigned int i0, unsigned int i1)/*{{{*/
601 {
602 ASSERT_ALWAYS(i0 < M->m);
603 ASSERT_ALWAYS(i1 < M->m);
604 if (i0 == i1) return;
605 for(unsigned int j = 0 ; j < M->n ; j++) {
606 mpq_swap(mpq_mat_entry(M, i0, j), mpq_mat_entry(M, i1, j));
607 }
608 }
609 /*}}}*/
mpz_mat_rotatedownrows(mpz_mat_ptr M,unsigned int i0,unsigned int k)610 void mpz_mat_rotatedownrows(mpz_mat_ptr M, unsigned int i0, unsigned int k)/*{{{*/
611 {
612 /* apply a circular shift on rows [i0...i0+k[ */
613 ASSERT_ALWAYS(i0 < M->m);
614 ASSERT_ALWAYS(i0 + k <= M->m);
615 if (k <= 1) return;
616 for(unsigned int j = 0 ; j < M->n ; j++) {
617 for(unsigned int s = k - 1 ; s-- ; ) {
618 mpz_swap(mpz_mat_entry(M, i0 + s, j), mpz_mat_entry(M, i0 + s + 1, j));
619 }
620 }
621 }
622 /*}}}*/
mpq_mat_rotatedownrows(mpq_mat_ptr M,unsigned int i0,unsigned int k)623 void mpq_mat_rotatedownrows(mpq_mat_ptr M, unsigned int i0, unsigned int k)/*{{{*/
624 {
625 ASSERT_ALWAYS(i0 < M->m);
626 ASSERT_ALWAYS(i0 + k <= M->m);
627 if (k <= 1) return;
628 for(unsigned int j = 0 ; j < M->n ; j++) {
629 for(unsigned int s = k - 1 ; s-- ; ) {
630 mpq_swap(mpq_mat_entry(M, i0 + s, j), mpq_mat_entry(M, i0 + s + 1, j));
631 }
632 }
633 }
634 /*}}}*/
mpz_mat_permuterows(mpz_mat_ptr M,unsigned int * perm)635 void mpz_mat_permuterows(mpz_mat_ptr M, unsigned int * perm)/*{{{*/
636 {
637 /* put row perm[k] in row k */
638 mpz_mat Mx;
639 mpz_mat_init(Mx, M->m, M->n);
640 for(unsigned int i = 0 ; i < M->m ; i++) {
641 for(unsigned int j = 0 ; j < M->n ; j++) {
642 mpz_swap(mpz_mat_entry(Mx, i, j),
643 mpz_mat_entry(M, perm[i], j));
644 }
645 }
646 mpz_mat_swap(M, Mx);
647 mpz_mat_clear(Mx);
648 }
649 /*}}}*/
mpq_mat_permuterows(mpq_mat_ptr M,unsigned int * perm)650 void mpq_mat_permuterows(mpq_mat_ptr M, unsigned int * perm)/*{{{*/
651 {
652 mpq_mat Mx;
653 mpq_mat_init(Mx, M->m, M->n);
654 for(unsigned int i = 0 ; i < M->m ; i++) {
655 for(unsigned int j = 0 ; j < M->n ; j++) {
656 mpq_swap(mpq_mat_entry(Mx, i, j),
657 mpq_mat_entry(M, perm[i], j));
658 }
659 }
660 mpq_mat_swap(M, Mx);
661 mpq_mat_clear(Mx);
662 }
663 /*}}}*/
permutation_signature(const unsigned int * perm,unsigned int n)664 static int permutation_signature(const unsigned int * perm, unsigned int n)/*{{{*/
665 {
666 /* can this be done in O(n) only ?? */
667 int sign = 1;
668 for(unsigned int i = 0 ; i < n ; i++)
669 for(unsigned int j = i ; j < n ; j++)
670 if (perm[j] < perm[i]) sign*=-1;
671 return sign;
672 }
673 /*}}}*/
674
675 /* add lambda times row i1 to row i0 */
mpz_mat_addmulrow(mpz_mat_ptr M,unsigned int i0,unsigned int i1,mpz_srcptr lambda)676 void mpz_mat_addmulrow(mpz_mat_ptr M, unsigned int i0, unsigned int i1, mpz_srcptr lambda)/*{{{*/
677 {
678 ASSERT_ALWAYS(i0 < M->m);
679 ASSERT_ALWAYS(i1 < M->m);
680 for(unsigned int j = 0 ; j < M->n ; j++) {
681 mpz_addmul(mpz_mat_entry(M, i0, j),
682 lambda,
683 mpz_mat_entry(M, i1, j));
684 }
685 }
686 /*}}}*/
mpz_mat_addmulrow_mod(mpz_mat_ptr M,unsigned int i0,unsigned int i1,mpz_srcptr lambda,mpz_srcptr p)687 void mpz_mat_addmulrow_mod(mpz_mat_ptr M, unsigned int i0, unsigned int i1, mpz_srcptr lambda, mpz_srcptr p)/*{{{*/
688 {
689 ASSERT_ALWAYS(i0 < M->m);
690 ASSERT_ALWAYS(i1 < M->m);
691 for(unsigned int j = 0 ; j < M->n ; j++) {
692 mpz_addmul(mpz_mat_entry(M, i0, j),
693 lambda,
694 mpz_mat_entry(M, i1, j));
695 mpz_mod(mpz_mat_entry(M, i0, j), mpz_mat_entry(M, i0, j), p);
696 }
697 }
698 /*}}}*/
699 /* add lambda times row i1 to row i0 */
mpq_mat_addmulrow(mpq_mat_ptr M,unsigned int i0,unsigned int i1,mpq_srcptr lambda)700 void mpq_mat_addmulrow(mpq_mat_ptr M, unsigned int i0, unsigned int i1, mpq_srcptr lambda)/*{{{*/
701 {
702 ASSERT_ALWAYS(i0 < M->m);
703 ASSERT_ALWAYS(i1 < M->m);
704 mpq_t tmp;
705 mpq_init(tmp);
706 for(unsigned int j = 0 ; j < M->n ; j++) {
707 /* we have no mpq_addmul, unfortunately */
708 mpq_mul(tmp,
709 lambda,
710 mpq_mat_entry(M, i1, j));
711 mpq_add(mpq_mat_entry(M, i0, j),
712 mpq_mat_entry(M, i0, j),
713 tmp);
714 }
715 mpq_clear(tmp);
716 }
717 /*}}}*/
718 /* subtract lambda times row i1 to row i0 */
mpz_mat_submulrow(mpz_mat_ptr M,unsigned int i0,unsigned int i1,mpz_srcptr lambda)719 void mpz_mat_submulrow(mpz_mat_ptr M, unsigned int i0, unsigned int i1, mpz_srcptr lambda)/*{{{*/
720 {
721 ASSERT_ALWAYS(i0 < M->m);
722 ASSERT_ALWAYS(i1 < M->m);
723 for(unsigned int j = 0 ; j < M->n ; j++) {
724 mpz_submul(mpz_mat_entry(M, i0, j),
725 lambda,
726 mpz_mat_entry(M, i1, j));
727 }
728 }
729 /*}}}*/
mpz_mat_submulrow_mod(mpz_mat_ptr M,unsigned int i0,unsigned int i1,mpz_srcptr lambda,mpz_srcptr p)730 void mpz_mat_submulrow_mod(mpz_mat_ptr M, unsigned int i0, unsigned int i1, mpz_srcptr lambda, mpz_srcptr p)/*{{{*/
731 {
732 ASSERT_ALWAYS(i0 < M->m);
733 ASSERT_ALWAYS(i1 < M->m);
734 for(unsigned int j = 0 ; j < M->n ; j++) {
735 mpz_submul(mpz_mat_entry(M, i0, j),
736 lambda,
737 mpz_mat_entry(M, i1, j));
738 mpz_mod(mpz_mat_entry(M, i0, j), mpz_mat_entry(M, i0, j), p);
739 }
740 }
741 /*}}}*/
742 /* subtract lambda times row i1 to row i0 */
mpq_mat_submulrow(mpq_mat_ptr M,unsigned int i0,unsigned int i1,mpq_srcptr lambda)743 void mpq_mat_submulrow(mpq_mat_ptr M, unsigned int i0, unsigned int i1, mpq_srcptr lambda)/*{{{*/
744 {
745 ASSERT_ALWAYS(i0 < M->m);
746 ASSERT_ALWAYS(i1 < M->m);
747 mpq_t tmp;
748 mpq_init(tmp);
749 for(unsigned int j = 0 ; j < M->n ; j++) {
750 /* we have no mpq_submul, unfortunately */
751 mpq_mul(tmp,
752 lambda,
753 mpq_mat_entry(M, i1, j));
754 mpq_sub(mpq_mat_entry(M, i0, j),
755 mpq_mat_entry(M, i0, j),
756 tmp);
757 }
758 mpq_clear(tmp);
759 }
760 /*}}}*/
761
762 /* add row i1 to row i0 */
mpz_mat_addrow(mpz_mat_ptr M,unsigned int i0,unsigned int i1)763 void mpz_mat_addrow(mpz_mat_ptr M, unsigned int i0, unsigned int i1)/*{{{*/
764 {
765 ASSERT_ALWAYS(i0 < M->m);
766 ASSERT_ALWAYS(i1 < M->m);
767 for(unsigned int j = 0 ; j < M->n ; j++) {
768 mpz_add(mpz_mat_entry(M, i0, j),
769 mpz_mat_entry(M, i0, j),
770 mpz_mat_entry(M, i1, j));
771 }
772 }
773 /*}}}*/
774 /* add row i1 to row i0 */
mpq_mat_addrow(mpq_mat_ptr M,unsigned int i0,unsigned int i1)775 void mpq_mat_addrow(mpq_mat_ptr M, unsigned int i0, unsigned int i1)/*{{{*/
776 {
777 ASSERT_ALWAYS(i0 < M->m);
778 ASSERT_ALWAYS(i1 < M->m);
779 for(unsigned int j = 0 ; j < M->n ; j++) {
780 mpq_add(mpq_mat_entry(M, i0, j),
781 mpq_mat_entry(M, i0, j),
782 mpq_mat_entry(M, i1, j));
783 }
784 }
785 /*}}}*/
786
787 /* subtract row i1 to row i0 */
mpz_mat_subrow(mpz_mat_ptr M,unsigned int i0,unsigned int i1)788 void mpz_mat_subrow(mpz_mat_ptr M, unsigned int i0, unsigned int i1)/*{{{*/
789 {
790 ASSERT_ALWAYS(i0 < M->m);
791 ASSERT_ALWAYS(i1 < M->m);
792 for(unsigned int j = 0 ; j < M->n ; j++) {
793 mpz_sub(mpz_mat_entry(M, i0, j),
794 mpz_mat_entry(M, i0, j),
795 mpz_mat_entry(M, i1, j));
796 }
797 }
798 /*}}}*/
799 /* subtract row i1 to row i0 */
mpq_mat_subrow(mpq_mat_ptr M,unsigned int i0,unsigned int i1)800 void mpq_mat_subrow(mpq_mat_ptr M, unsigned int i0, unsigned int i1)/*{{{*/
801 {
802 ASSERT_ALWAYS(i0 < M->m);
803 ASSERT_ALWAYS(i1 < M->m);
804 for(unsigned int j = 0 ; j < M->n ; j++) {
805 mpq_sub(mpq_mat_entry(M, i0, j),
806 mpq_mat_entry(M, i0, j),
807 mpq_mat_entry(M, i1, j));
808 }
809 }
810 /*}}}*/
811 /* multiply row i0 by lambda */
mpz_mat_mulrow(mpz_mat_ptr M,unsigned int i0,mpz_srcptr lambda)812 void mpz_mat_mulrow(mpz_mat_ptr M, unsigned int i0, mpz_srcptr lambda)/*{{{*/
813 {
814 ASSERT_ALWAYS(i0 < M->m);
815 for(unsigned int j = 0 ; j < M->n ; j++) {
816 mpz_mul(mpz_mat_entry(M, i0, j), mpz_mat_entry(M, i0, j), lambda);
817 }
818 }
819 /*}}}*/
mpz_mat_mulrow_mod(mpz_mat_ptr M,unsigned int i0,mpz_srcptr lambda,mpz_srcptr p)820 void mpz_mat_mulrow_mod(mpz_mat_ptr M, unsigned int i0, mpz_srcptr lambda, mpz_srcptr p)/*{{{*/
821 {
822 ASSERT_ALWAYS(i0 < M->m);
823 for(unsigned int j = 0 ; j < M->n ; j++) {
824 mpz_mul(mpz_mat_entry(M, i0, j), mpz_mat_entry(M, i0, j), lambda);
825 mpz_mod(mpz_mat_entry(M, i0, j), mpz_mat_entry(M, i0, j), p);
826 }
827 }
828 /*}}}*/
829 /* multiply row i0 by lambda */
mpq_mat_mulrow(mpq_mat_ptr M,unsigned int i0,mpq_srcptr lambda)830 void mpq_mat_mulrow(mpq_mat_ptr M, unsigned int i0, mpq_srcptr lambda)/*{{{*/
831 {
832 ASSERT_ALWAYS(i0 < M->m);
833 for(unsigned int j = 0 ; j < M->n ; j++) {
834 mpq_mul(mpq_mat_entry(M, i0, j), mpq_mat_entry(M, i0, j), lambda);
835 }
836 }
837 /*}}}*/
838
839 #if 0
840 /* XXX never used, untested ! */
841 /* this computes an additive combination of n rows into row [didx] of the
842 * initial matrix. We require that this destination row be cleared
843 * initially.
844 */
845 void mpz_mat_combinerows(mpz_mat_ptr M, unsigned int didx, unsigned int sidx,/*{{{*/
846 mpz_srcptr * lambda, unsigned int n)
847 {
848 for(unsigned int j = 0 ; j < M->n ; j++) {
849 ASSERT_ALWAYS(mpz_cmp_ui(mpz_mat_entry(M, didx, j), 0) == 0);
850 }
851 for(unsigned int i = 0 ; i < n ; i++) {
852 mpz_mat_addmulrow(M, didx, sidx + i, lambda[i]);
853 }
854 }
855 /*}}}*/
856 #endif
857
858 /* }}} */
859 /*{{{ I/O*/
mpz_mat_fprint(FILE * stream,mpz_mat_srcptr M)860 void mpz_mat_fprint(FILE * stream, mpz_mat_srcptr M)
861 {
862 for(unsigned int i = 0 ; i < M->m ; i++) {
863 for(unsigned int j = 0 ; j < M->n ; j++) {
864 gmp_fprintf(stream, " %Zd", mpz_mat_entry_const(M, i, j));
865 }
866 fprintf(stream, "\n");
867 }
868 }
869
mpq_mat_fprint(FILE * stream,mpq_mat_srcptr M)870 void mpq_mat_fprint(FILE * stream, mpq_mat_srcptr M)
871 {
872 for(unsigned int i = 0 ; i < M->m ; i++) {
873 for(unsigned int j = 0 ; j < M->n ; j++) {
874 gmp_fprintf(stream, " %Qd", mpq_mat_entry_const(M, i, j));
875 }
876 fprintf(stream, "\n");
877 }
878 }
879
mpq_mat_fprint_as_mpz(FILE * f,mpq_mat_srcptr M)880 void mpq_mat_fprint_as_mpz(FILE* f, mpq_mat_srcptr M)
881 {
882 mpz_t denom;
883 mpz_mat N;
884
885 mpz_init(denom);
886 mpz_mat_init(N,0,0);
887
888 mpq_mat_numden(N, denom, M);
889
890 mpz_mat_fprint(f,N);
891 gmp_fprintf(f, "Denominator is : %Zd\n",denom);
892
893 mpz_mat_clear(N);
894 mpz_clear(denom);
895 }
896 /*}}}*/
897 /*{{{ comparison */
mpz_mat_cmp(mpz_mat_srcptr M,mpz_mat_srcptr N)898 int mpz_mat_cmp(mpz_mat_srcptr M, mpz_mat_srcptr N)/*{{{*/
899 {
900 /* shall we abort or return something for incompatible matrices ?? */
901 if (M->m != N->m) return M->m - N->m;
902 if (M->n != N->n) return M->n - N->n;
903 // ASSERT_ALWAYS((M->m == N->m) && (M->n == N->n));
904 unsigned int m = M->m;
905 unsigned int n = M->n;
906 unsigned int i,j;
907 for (i = 0 ; i < m ; i++){
908 for (j = 0 ; j < n ; j++){
909 int k = mpz_cmp(mpz_mat_entry_const(M,i,j),mpz_mat_entry_const(N,i,j));
910 if(k!=0) return k;
911 }
912 }
913 return 0;
914 }
915 /*}}}*/
mpq_mat_cmp(mpq_mat_srcptr M,mpq_mat_srcptr N)916 int mpq_mat_cmp(mpq_mat_srcptr M, mpq_mat_srcptr N)/*{{{*/
917 {
918 /* shall we abort or return something for incompatible matrices ?? */
919 if (M->m != N->m) return M->m - N->m;
920 if (M->n != N->n) return M->n - N->n;
921 // ASSERT_ALWAYS((M->m == N->m) && (M->n == N->n));
922 unsigned int m = M->m;
923 unsigned int n = M->n;
924 unsigned int i,j;
925 for (i = 0 ; i < m ; i++){
926 for (j = 0 ; j < n ; j++){
927 int k = mpq_cmp(mpq_mat_entry_const(M,i,j),mpq_mat_entry_const(N,i,j));
928 if(k!=0) return k;
929 }
930 }
931 return 0;
932 }
933 /*}}}*/
934
935 /* TODO (perhaps) :
936 * mp[qz]_mat_is_{upper,lower}_triangular
937 * mp[qz]_mat_is_diagonal
938 */
939 /*}}}*/
940 /*{{{ add and subtract */
mpz_mat_add(mpz_mat_ptr C,mpz_mat_srcptr A,mpz_mat_srcptr B)941 void mpz_mat_add(mpz_mat_ptr C, mpz_mat_srcptr A, mpz_mat_srcptr B)/*{{{*/
942 {
943 ASSERT_ALWAYS(A->m == B->m);
944 ASSERT_ALWAYS(A->n == B->n);
945 mpz_mat_realloc(C,A->m,A->n); /* no-op if C == A or C == B */
946
947 for(unsigned int i = 0 ; i < A->m ; i++){
948 for(unsigned int j = 0 ; j < A->n ; j++){
949 mpz_srcptr aij = mpz_mat_entry_const(A, i, j);
950 mpz_srcptr bij = mpz_mat_entry_const(B, i, j);
951 mpz_ptr cij = mpz_mat_entry(C, i, j);
952 mpz_add(cij,aij,bij);
953 }
954 }
955 }
956 /* }}} */
mpq_mat_add(mpq_mat_ptr C,mpq_mat_srcptr A,mpq_mat_srcptr B)957 void mpq_mat_add(mpq_mat_ptr C, mpq_mat_srcptr A, mpq_mat_srcptr B)/*{{{*/
958 {
959 ASSERT_ALWAYS(A->m == B->m);
960 ASSERT_ALWAYS(A->n == B->n);
961 mpq_mat_realloc(C,A->m,A->n); /* no-op if C == A or C == B */
962
963 for(unsigned int i = 0 ; i < A->m ; i++){
964 for(unsigned int j = 0 ; j < A->n ; j++){
965 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
966 mpq_srcptr bij = mpq_mat_entry_const(B, i, j);
967 mpq_ptr cij = mpq_mat_entry(C, i, j);
968 mpq_add(cij,aij,bij);
969 }
970 }
971 }
972 /*}}}*/
mpz_mat_sub(mpz_mat_ptr C,mpz_mat_srcptr A,mpz_mat_srcptr B)973 void mpz_mat_sub(mpz_mat_ptr C, mpz_mat_srcptr A, mpz_mat_srcptr B)/*{{{*/
974 {
975 ASSERT_ALWAYS(A->m == B->m);
976 ASSERT_ALWAYS(A->n == B->n);
977 mpz_mat_realloc(C,A->m,A->n); /* no-op if C == A or C == B */
978
979 for(unsigned int i = 0 ; i < A->m ; i++){
980 for(unsigned int j = 0 ; j < A->n ; j++){
981 mpz_srcptr aij = mpz_mat_entry_const(A, i, j);
982 mpz_srcptr bij = mpz_mat_entry_const(B, i, j);
983 mpz_ptr cij = mpz_mat_entry(C, i, j);
984 mpz_sub(cij,aij,bij);
985 }
986 }
987 }
988 /* }}} */
mpq_mat_sub(mpq_mat_ptr C,mpq_mat_srcptr A,mpq_mat_srcptr B)989 void mpq_mat_sub(mpq_mat_ptr C, mpq_mat_srcptr A, mpq_mat_srcptr B)/*{{{*/
990 {
991 ASSERT_ALWAYS(A->m == B->m);
992 ASSERT_ALWAYS(A->n == B->n);
993 mpq_mat_realloc(C,A->m,A->n); /* no-op if C == A or C == B */
994
995 for(unsigned int i = 0 ; i < A->m ; i++){
996 for(unsigned int j = 0 ; j < A->n ; j++){
997 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
998 mpq_srcptr bij = mpq_mat_entry_const(B, i, j);
999 mpq_ptr cij = mpq_mat_entry(C, i, j);
1000 mpq_sub(cij,aij,bij);
1001 }
1002 }
1003 }
1004 /*}}}*/
1005 /*}}}*/
1006 /*{{{ multiplication */
mpz_mat_mul(mpz_mat_ptr C,mpz_mat_srcptr A,mpz_mat_srcptr B)1007 void mpz_mat_mul(mpz_mat_ptr C, mpz_mat_srcptr A, mpz_mat_srcptr B)/*{{{*/
1008 {
1009 ASSERT_ALWAYS(A->n == B->m);
1010 if (C == A || C == B) {
1011 mpz_mat D;
1012 mpz_mat_init(D, A->m, B->n);
1013 mpz_mat_mul(D, A, B);
1014 mpz_mat_swap(D, C);
1015 mpz_mat_clear(D);
1016 return;
1017 }
1018 mpz_mat_realloc(C, A->m, B->n);
1019 for(unsigned int i = 0 ; i < A->m ; i++) {
1020 for(unsigned int j = 0 ; j < B->n ; j++) {
1021 mpz_set_ui(mpz_mat_entry(C, i, j), 0);
1022 for(unsigned int k = 0 ; k < B->m ; k++) {
1023 mpz_addmul(mpz_mat_entry(C, i, j),
1024 mpz_mat_entry_const(A, i, k),
1025 mpz_mat_entry_const(B, k, j));
1026 }
1027 }
1028 }
1029 }/*}}}*/
1030
mpz_mat_mul_mod_ui(mpz_mat_ptr C,mpz_mat_srcptr A,mpz_mat_srcptr B,unsigned long p)1031 void mpz_mat_mul_mod_ui(mpz_mat_ptr C, mpz_mat_srcptr A, mpz_mat_srcptr B, unsigned long p)/*{{{*/
1032 {
1033 ASSERT_ALWAYS(A->n == B->m);
1034 if (C == A || C == B) {
1035 mpz_mat D;
1036 mpz_mat_init(D, A->m, B->n);
1037 mpz_mat_mul_mod_ui(D, A, B, p);
1038 mpz_mat_swap(D, C);
1039 mpz_mat_clear(D);
1040 return;
1041 }
1042 mpz_mat_realloc(C, A->m, B->n);
1043 for(unsigned int i = 0 ; i < A->m ; i++) {
1044 for(unsigned int j = 0 ; j < B->n ; j++) {
1045 mpz_set_ui(mpz_mat_entry(C, i, j), 0);
1046 for(unsigned int k = 0 ; k < B->m ; k++) {
1047 mpz_addmul(mpz_mat_entry(C, i, j),
1048 mpz_mat_entry_const(A, i, k),
1049 mpz_mat_entry_const(B, k, j));
1050 }
1051 }
1052 }
1053 mpz_mat_mod_ui(C,C,p);
1054 }/*}}}*/
1055
mpz_mat_mul_mod_mpz(mpz_mat_ptr C,mpz_mat_srcptr A,mpz_mat_srcptr B,mpz_srcptr p)1056 void mpz_mat_mul_mod_mpz(mpz_mat_ptr C, mpz_mat_srcptr A, mpz_mat_srcptr B, mpz_srcptr p)/*{{{*/
1057 {
1058 ASSERT_ALWAYS(A->n == B->m);
1059 if (C == A || C == B) {
1060 mpz_mat D;
1061 mpz_mat_init(D, A->m, B->n);
1062 mpz_mat_mul_mod_mpz(D, A, B, p);
1063 mpz_mat_swap(D, C);
1064 mpz_mat_clear(D);
1065 return;
1066 }
1067 mpz_mat_realloc(C, A->m, B->n);
1068 for(unsigned int i = 0 ; i < A->m ; i++) {
1069 for(unsigned int j = 0 ; j < B->n ; j++) {
1070 mpz_set_ui(mpz_mat_entry(C, i, j), 0);
1071 for(unsigned int k = 0 ; k < B->m ; k++) {
1072 mpz_addmul(mpz_mat_entry(C, i, j),
1073 mpz_mat_entry_const(A, i, k),
1074 mpz_mat_entry_const(B, k, j));
1075 }
1076 }
1077 }
1078 mpz_mat_mod_mpz(C,C,p);
1079 }/*}}}*/
1080
mpz_mat_mul_mpz(mpz_mat_ptr B,mpz_mat_srcptr A,mpz_ptr k)1081 void mpz_mat_mul_mpz(mpz_mat_ptr B, mpz_mat_srcptr A, mpz_ptr k)/*{{{*/
1082 {
1083 mpz_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1084 for(unsigned int i = 0 ; i < B->m ; i++) {
1085 for(unsigned int j = 0 ; j < B->n ; j++) {
1086 mpz_mul(mpz_mat_entry(B, i, j),mpz_mat_entry_const(A,i,j),k);
1087 }
1088 }
1089 }/*}}}*/
1090
mpz_mat_divexact_mpz(mpz_mat_ptr B,mpz_mat_srcptr A,mpz_srcptr k)1091 void mpz_mat_divexact_mpz(mpz_mat_ptr B, mpz_mat_srcptr A, mpz_srcptr k)/*{{{*/
1092 {
1093 mpz_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1094 for(unsigned int i = 0 ; i < B->m ; i++) {
1095 for(unsigned int j = 0 ; j < B->n ; j++) {
1096 mpz_divexact(mpz_mat_entry(B, i, j),mpz_mat_entry_const(A,i,j),k);
1097 }
1098 }
1099 }/*}}}*/
1100
mpz_mat_divexact_ui(mpz_mat_ptr B,mpz_mat_srcptr A,unsigned long k)1101 void mpz_mat_divexact_ui(mpz_mat_ptr B, mpz_mat_srcptr A, unsigned long k)/*{{{*/
1102 {
1103 mpz_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1104 for(unsigned int i = 0 ; i < B->m ; i++) {
1105 for(unsigned int j = 0 ; j < B->n ; j++) {
1106 mpz_divexact_ui(mpz_mat_entry(B, i, j),mpz_mat_entry_const(A,i,j),k);
1107 }
1108 }
1109 }/*}}}*/
1110
mpz_mat_mul_si(mpz_mat_ptr B,mpz_mat_srcptr A,long k)1111 void mpz_mat_mul_si(mpz_mat_ptr B, mpz_mat_srcptr A, long k)/*{{{*/
1112 {
1113 mpz_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1114 for(unsigned int i = 0 ; i < B->m ; i++) {
1115 for(unsigned int j = 0 ; j < B->n ; j++) {
1116 mpz_mul_si(mpz_mat_entry(B, i, j),mpz_mat_entry_const(A,i,j),k);
1117 }
1118 }
1119 }/*}}}*/
1120
mpz_mat_mul_ui(mpz_mat_ptr B,mpz_mat_srcptr A,unsigned long k)1121 void mpz_mat_mul_ui(mpz_mat_ptr B, mpz_mat_srcptr A, unsigned long k)/*{{{*/
1122 {
1123 mpz_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1124 for(unsigned int i = 0 ; i < B->m ; i++) {
1125 for(unsigned int j = 0 ; j < B->n ; j++) {
1126 mpz_mul_ui(mpz_mat_entry(B, i, j),mpz_mat_entry_const(A,i,j),k);
1127 }
1128 }
1129 }/*}}}*/
1130
1131
mpq_mat_mul(mpq_mat_ptr C,mpq_mat_srcptr A,mpq_mat_srcptr B)1132 void mpq_mat_mul(mpq_mat_ptr C, mpq_mat_srcptr A, mpq_mat_srcptr B)/*{{{*/
1133 {
1134 ASSERT_ALWAYS(A->n == B->m);
1135 if (C == A || C == B) {
1136 mpq_mat D;
1137 mpq_mat_init(D, A->m, B->n);
1138 mpq_mat_mul(D, A, B);
1139 mpq_mat_swap(D, C);
1140 mpq_mat_clear(D);
1141 return;
1142 }
1143 mpq_t z;
1144 mpq_init(z);
1145 mpq_mat_realloc(C, A->m, B->n);
1146 for(unsigned int i = 0 ; i < A->m ; i++) {
1147 for(unsigned int j = 0 ; j < B->n ; j++) {
1148 mpq_ptr cij = mpq_mat_entry(C, i, j);
1149 mpq_set_ui(cij, 0, 1);
1150 for(unsigned int k = 0 ; k < B->m ; k++) {
1151 mpq_srcptr aik = mpq_mat_entry_const(A, i, k);
1152 mpq_srcptr bkj = mpq_mat_entry_const(B, k, j);
1153 mpq_mul(z, aik, bkj);
1154 mpq_add(cij, cij, z);
1155 }
1156 mpq_canonicalize(cij);
1157 }
1158 }
1159 mpq_clear(z);
1160 }/*}}}*/
1161
mpq_mat_mul_mpz(mpq_mat_ptr B,mpq_mat_srcptr A,mpz_srcptr k)1162 void mpq_mat_mul_mpz(mpq_mat_ptr B, mpq_mat_srcptr A, mpz_srcptr k)/*{{{*/
1163 {
1164 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1165 for(unsigned int i = 0 ; i < B->m ; i++) {
1166 for(unsigned int j = 0 ; j < B->n ; j++) {
1167 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1168 mpq_ptr bij = mpq_mat_entry(B, i, j);
1169 mpz_set(mpq_denref(bij), mpq_denref(aij));
1170 mpz_mul(mpq_numref(bij), mpq_numref(aij), k);
1171 mpq_canonicalize(bij);
1172 }
1173 }
1174 }/*}}}*/
mpq_mat_mul_mpq(mpq_mat_ptr B,mpq_mat_srcptr A,mpq_srcptr k)1175 void mpq_mat_mul_mpq(mpq_mat_ptr B, mpq_mat_srcptr A, mpq_srcptr k)/*{{{*/
1176 {
1177 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1178 for(unsigned int i = 0 ; i < B->m ; i++) {
1179 for(unsigned int j = 0 ; j < B->n ; j++) {
1180 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1181 mpq_ptr bij = mpq_mat_entry(B, i, j);
1182 mpq_mul(bij, aij, k);
1183 }
1184 }
1185 }
1186 /*}}}*/
mpq_mat_mul_si(mpq_mat_ptr B,mpq_mat_srcptr A,long k)1187 void mpq_mat_mul_si(mpq_mat_ptr B, mpq_mat_srcptr A, long k)/*{{{*/
1188 {
1189 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1190 for(unsigned int i = 0 ; i < B->m ; i++) {
1191 for(unsigned int j = 0 ; j < B->n ; j++) {
1192 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1193 mpq_ptr bij = mpq_mat_entry(B, i, j);
1194 mpz_set(mpq_denref(bij),mpq_denref(aij));
1195 mpz_mul_si(mpq_numref(bij),mpq_numref(aij),k);
1196 mpq_canonicalize(bij);
1197 }
1198 }
1199 }
1200 /*}}}*/
mpq_mat_mul_ui(mpq_mat_ptr B,mpq_mat_srcptr A,unsigned long k)1201 void mpq_mat_mul_ui(mpq_mat_ptr B, mpq_mat_srcptr A, unsigned long k)/*{{{*/
1202 {
1203 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1204 for(unsigned int i = 0 ; i < B->m ; i++) {
1205 for(unsigned int j = 0 ; j < B->n ; j++) {
1206 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1207 mpq_ptr bij = mpq_mat_entry(B, i, j);
1208 mpz_set(mpq_denref(bij),mpq_denref(aij));
1209 mpz_mul_ui(mpq_numref(bij),mpq_numref(aij),k);
1210 mpq_canonicalize(bij);
1211 }
1212 }
1213 }
1214 /*}}}*/
mpq_mat_div_mpz(mpq_mat_ptr B,mpq_mat_srcptr A,mpz_srcptr k)1215 void mpq_mat_div_mpz(mpq_mat_ptr B, mpq_mat_srcptr A, mpz_srcptr k)/*{{{*/
1216 {
1217 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1218 for(unsigned int i = 0 ; i < B->m ; i++) {
1219 for(unsigned int j = 0 ; j < B->n ; j++) {
1220 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1221 mpq_ptr bij = mpq_mat_entry(B, i, j);
1222 mpz_mul(mpq_denref(bij), mpq_denref(aij), k);
1223 mpz_set(mpq_numref(bij), mpq_numref(aij));
1224 mpq_canonicalize(bij);
1225 }
1226 }
1227 }/*}}}*/
mpq_mat_div_mpq(mpq_mat_ptr B,mpq_mat_srcptr A,mpq_srcptr k)1228 void mpq_mat_div_mpq(mpq_mat_ptr B, mpq_mat_srcptr A, mpq_srcptr k)/*{{{*/
1229 {
1230 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1231 for(unsigned int i = 0 ; i < B->m ; i++) {
1232 for(unsigned int j = 0 ; j < B->n ; j++) {
1233 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1234 mpq_ptr bij = mpq_mat_entry(B, i, j);
1235 mpz_mul(mpq_numref(bij), mpq_numref(aij), mpq_denref(k));
1236 mpz_mul(mpq_denref(bij), mpq_denref(aij), mpq_numref(k));
1237 mpq_canonicalize(bij);
1238 }
1239 }
1240 }
1241 /*}}}*/
mpq_mat_div_si(mpq_mat_ptr B,mpq_mat_srcptr A,long k)1242 void mpq_mat_div_si(mpq_mat_ptr B, mpq_mat_srcptr A, long k)/*{{{*/
1243 {
1244 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1245 for(unsigned int i = 0 ; i < B->m ; i++) {
1246 for(unsigned int j = 0 ; j < B->n ; j++) {
1247 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1248 mpq_ptr bij = mpq_mat_entry(B, i, j);
1249 mpz_mul_si(mpq_denref(bij),mpq_denref(aij), k);
1250 mpz_set(mpq_numref(bij),mpq_numref(aij));
1251 mpq_canonicalize(bij);
1252 }
1253 }
1254 }
1255 /*}}}*/
mpq_mat_div_ui(mpq_mat_ptr B,mpq_mat_srcptr A,unsigned long k)1256 void mpq_mat_div_ui(mpq_mat_ptr B, mpq_mat_srcptr A, unsigned long k)/*{{{*/
1257 {
1258 mpq_mat_realloc(B,A->m,A->n); /* no-op if A == B */
1259 for(unsigned int i = 0 ; i < B->m ; i++) {
1260 for(unsigned int j = 0 ; j < B->n ; j++) {
1261 mpq_srcptr aij = mpq_mat_entry_const(A, i, j);
1262 mpq_ptr bij = mpq_mat_entry(B, i, j);
1263 mpz_mul_ui(mpq_denref(bij),mpq_denref(aij), k);
1264 mpz_set(mpq_numref(bij),mpq_numref(aij));
1265 mpq_canonicalize(bij);
1266 }
1267 }
1268 }
1269 /*}}}*/
1270 /* }}} */
1271 /*{{{ powering */
1272 // Returns A^n for n >= 0, A being a square matrix
mpz_mat_pow_ui(mpz_mat_ptr B,mpz_mat_srcptr A,unsigned long n)1273 void mpz_mat_pow_ui(mpz_mat_ptr B, mpz_mat_srcptr A, unsigned long n)/*{{{*/
1274 {
1275 ASSERT_ALWAYS(A->n == A->m);
1276 if (n == 0) {
1277 mpz_mat_set_ui(B, 1);
1278 return;
1279 }
1280 if (B == A) {
1281 mpz_mat C;
1282 mpz_mat_init(C, A->m, A->n);
1283 mpz_mat_pow_ui(C, A, n);
1284 mpz_mat_swap(B, C);
1285 mpz_mat_clear(C);
1286 return;
1287 }
1288 unsigned long k = ((~0UL)>>1) + 1;
1289 for( ; k > n ; k >>= 1);
1290 mpz_mat_set(B, A);
1291 mpz_mat Q;
1292 mpz_mat_init(Q, A->m, A->n);
1293 for( ; k >>= 1 ; ) {
1294 mpz_mat_mul(Q, B, B);
1295 mpz_mat_swap(Q, B);
1296 if (n & k) {
1297 mpz_mat_mul(Q, B, A);
1298 mpz_mat_swap(Q, B);
1299 }
1300 }
1301 mpz_mat_clear(Q);
1302 }/*}}}*/
1303
mpz_mat_pow_ui_mod_ui(mpz_mat_ptr B,mpz_mat_srcptr A,unsigned long n,unsigned long p)1304 void mpz_mat_pow_ui_mod_ui(mpz_mat_ptr B, mpz_mat_srcptr A, unsigned long n, unsigned long p)/*{{{*/
1305 {
1306 ASSERT_ALWAYS(A->n == A->m);
1307 if (n == 0) {
1308 mpz_mat_set_ui(B, 1);
1309 return;
1310 }
1311 if (B == A) {
1312 mpz_mat C;
1313 mpz_mat_init(C, A->m, A->n);
1314 mpz_mat_pow_ui_mod_ui(C, A, n, p);
1315 mpz_mat_swap(B, C);
1316 mpz_mat_clear(C);
1317 return;
1318 }
1319 unsigned long k = ((~0UL)>>1) + 1;
1320 for( ; k > n ; k >>= 1);
1321 mpz_mat_set(B, A);
1322 mpz_mat Q;
1323 mpz_mat_init(Q, A->m, A->n);
1324 for( ; k >>= 1 ; ) {
1325 mpz_mat_mul_mod_ui(Q, B, B, p);
1326 mpz_mat_swap(Q, B);
1327 if (n & k) {
1328 mpz_mat_mul_mod_ui(Q, B, A, p);
1329 mpz_mat_swap(Q, B);
1330 }
1331 }
1332 mpz_mat_clear(Q);
1333 }/*}}}*/
1334
mpz_mat_pow_ui_mod_mpz(mpz_mat_ptr B,mpz_mat_srcptr A,unsigned long n,mpz_srcptr p)1335 void mpz_mat_pow_ui_mod_mpz(mpz_mat_ptr B, mpz_mat_srcptr A, unsigned long n, mpz_srcptr p)/*{{{*/
1336 {
1337 ASSERT_ALWAYS(A->n == A->m);
1338 if (n == 0) {
1339 mpz_mat_set_ui(B, 1);
1340 return;
1341 }
1342 if (B == A) {
1343 mpz_mat C;
1344 mpz_mat_init(C, A->m, A->n);
1345 mpz_mat_pow_ui_mod_mpz(C, A, n, p);
1346 mpz_mat_swap(B, C);
1347 mpz_mat_clear(C);
1348 return;
1349 }
1350 unsigned long k = ((~0UL)>>1) + 1;
1351 for( ; k > n ; k >>= 1);
1352 mpz_mat_set(B, A);
1353 mpz_mat Q;
1354 mpz_mat_init(Q, A->m, A->n);
1355 for( ; k >>= 1 ; ) {
1356 mpz_mat_mul_mod_mpz(Q, B, B, p);
1357 mpz_mat_swap(Q, B);
1358 if (n & k) {
1359 mpz_mat_mul_mod_mpz(Q, B, A, p);
1360 mpz_mat_swap(Q, B);
1361 }
1362 }
1363 mpz_mat_clear(Q);
1364 }/*}}}*/
1365 /*}}}*/
1366 /* {{{ polynomial evaluation */
1367
mpz_poly_eval_mpz_mat(mpz_mat_ptr D,mpz_mat_srcptr M,mpz_poly_srcptr f)1368 void mpz_poly_eval_mpz_mat(mpz_mat_ptr D, mpz_mat_srcptr M, mpz_poly_srcptr f)/*{{{*/
1369 {
1370 ASSERT_ALWAYS(M->m == M->n);
1371 unsigned int n = M->n;
1372 int d = f->deg;
1373 if (d < 0) {
1374 mpz_mat_realloc(D, n, n);
1375 mpz_mat_set_ui(D, 0);
1376 return;
1377 }
1378 if (D == M) {
1379 mpz_mat X;
1380 mpz_mat_init(X, 0, 0);
1381 mpz_poly_eval_mpz_mat(X, M, f);
1382 mpz_mat_swap(D, X);
1383 mpz_mat_clear(X);
1384 return;
1385 }
1386 mpz_mat_realloc(D, n, n);
1387 mpz_mat_set_mpz(D, f->coeff[f->deg]);
1388 for(int i = f->deg - 1 ; i >= 0 ; i--) {
1389 mpz_mat_mul(D, M, D);
1390 mpz_mat_add_mpz(D, f->coeff[i]);
1391 }
1392 }
1393 /*}}}*/
mpz_poly_eval_mpz_mat_mod_ui(mpz_mat_ptr D,mpz_mat_srcptr M,mpz_poly_srcptr f,unsigned long p)1394 void mpz_poly_eval_mpz_mat_mod_ui(mpz_mat_ptr D, mpz_mat_srcptr M, mpz_poly_srcptr f, unsigned long p)/*{{{*/
1395 {
1396 ASSERT_ALWAYS(M->m == M->n);
1397 unsigned int n = M->n;
1398 int d = f->deg;
1399 if (d < 0) {
1400 mpz_mat_realloc(D, n, n);
1401 mpz_mat_set_ui(D, 0);
1402 return;
1403 }
1404 if (D == M) {
1405 mpz_mat X;
1406 mpz_mat_init(X, 0, 0);
1407 mpz_poly_eval_mpz_mat_mod_ui(X, M, f, p);
1408 mpz_mat_swap(D, X);
1409 mpz_mat_clear(X);
1410 return;
1411 }
1412 mpz_mat_realloc(D, n, n);
1413 mpz_mat_set_ui(D, mpz_fdiv_ui(f->coeff[f->deg], p));
1414 for(int i = f->deg - 1 ; i >= 0 ; i--) {
1415 mpz_mat_mul_mod_ui(D, M, D, p);
1416 mpz_mat_add_ui(D, mpz_fdiv_ui(f->coeff[i], p));
1417 mpz_mat_mod_ui(D, D, p);
1418 }
1419 }
1420 /*}}}*/
mpz_poly_eval_mpz_mat_mod_mpz(mpz_mat_ptr D,mpz_mat_srcptr M,mpz_poly_srcptr f,mpz_srcptr p)1421 void mpz_poly_eval_mpz_mat_mod_mpz(mpz_mat_ptr D, mpz_mat_srcptr M, mpz_poly_srcptr f, mpz_srcptr p)/*{{{*/
1422 {
1423 ASSERT_ALWAYS(M->m == M->n);
1424 unsigned int n = M->n;
1425 int d = f->deg;
1426 if (d < 0) {
1427 mpz_mat_realloc(D, n, n);
1428 mpz_mat_set_ui(D, 0);
1429 return;
1430 }
1431 if (D == M) {
1432 mpz_mat X;
1433 mpz_mat_init(X, 0, 0);
1434 mpz_poly_eval_mpz_mat_mod_mpz(X, M, f, p);
1435 mpz_mat_swap(D, X);
1436 mpz_mat_clear(X);
1437 return;
1438 }
1439 mpz_mat_realloc(D, n, n);
1440 mpz_t tmp;
1441 mpz_init(tmp);
1442 mpz_fdiv_r(tmp, f->coeff[f->deg], p);
1443 mpz_mat_set_mpz(D, tmp);
1444 for(int i = f->deg - 1 ; i >= 0 ; i--) {
1445 mpz_mat_mul_mod_mpz(D, M, D, p);
1446 mpz_fdiv_r(tmp, f->coeff[i], p);
1447 mpz_mat_add_mpz(D, tmp);
1448 mpz_mat_mod_mpz(D, D, p);
1449 }
1450 mpz_clear(tmp);
1451 }
1452 /*}}}*/
1453 /*}}}*/
1454
1455 /* {{{ gaussian reduction over the rationals
1456 * this is a backend for row gaussian reduction. T receives the
1457 * transformation matrix. M is modified in place.
1458 * this includes reordering of the rows.
1459 *
1460 * this works reasonably well over Fp, but for rationals we suffer from
1461 * coefficient blowup.
1462 */
mpq_mat_gauss_backend(mpq_mat_ptr M,mpq_mat_ptr T)1463 void mpq_mat_gauss_backend(mpq_mat_ptr M, mpq_mat_ptr T)
1464 {
1465 unsigned int m = M->m;
1466 unsigned int n = M->n;
1467 ASSERT_ALWAYS(M != T);
1468 mpq_mat_realloc(T, m, m);
1469 mpq_mat_set_ui(T, 1);
1470 unsigned int rank = 0;
1471 mpq_t tmp1, tmp2;
1472 mpq_init(tmp1);
1473 mpq_init(tmp2);
1474 for(unsigned int j = 0 ; j < n ; j++) {
1475 unsigned int i1;
1476 for(i1 = rank ; i1 < M->m ; i1++) {
1477 if (mpq_cmp_ui(mpq_mat_entry(M, i1, j), 0, 1) != 0) break;
1478 }
1479 if (i1 == M->m) /* this column is zero */
1480 continue;
1481 mpq_mat_swaprows(M, rank, i1);
1482 mpq_mat_swaprows(T, rank, i1);
1483 i1 = rank;
1484 /* canonicalize this row */
1485 mpq_inv(tmp1, mpq_mat_entry(M, rank, j));
1486 mpq_mat_mulrow(M, rank, tmp1);
1487 mpq_mat_mulrow(T, rank, tmp1);
1488 for(unsigned int i0 = 0 ; i0 < m ; i0++) {
1489 if (i0 == rank) continue;
1490 /* we've normalized row rank, so all it takes is to multiply by
1491 * M[i0, j]
1492 */
1493 mpq_set(tmp2, mpq_mat_entry(M, i0, j));
1494 mpq_mat_submulrow(M, i0, rank, tmp2);
1495 mpq_mat_submulrow(T, i0, rank, tmp2);
1496 }
1497 rank++;
1498 }
1499 mpq_clear(tmp1);
1500 mpq_clear(tmp2);
1501 }
1502 /* }}} */
1503 /* {{{ Gaussian reduction over Z/pZ
1504 *
1505 * We sort of return something even in the case where p is not prime.
1506 *
1507 * M is the reduced matrix (modification in place of the input), T the
1508 * transformation matrix (so that T * input-M = output-M)
1509 *
1510 * Note that if M is rectangular, this the function we want to call in
1511 * order to get a row-reduced echelon form of M.
1512 *
1513 * T may be NULL in case we don't give a penny about the transformation
1514 * matrix.
1515 */
mpz_mat_gauss_backend_mod_mpz(mpz_mat_ptr M,mpz_mat_ptr T,mpz_srcptr p)1516 void mpz_mat_gauss_backend_mod_mpz(mpz_mat_ptr M, mpz_mat_ptr T, mpz_srcptr p)
1517 {
1518 unsigned int m = M->m;
1519 unsigned int n = M->n;
1520 ASSERT_ALWAYS(M != T);
1521 if (T) mpz_mat_realloc(T, m, m);
1522 if (T) mpz_mat_set_ui(T, 1);
1523 unsigned int rank = 0;
1524 mpz_t gcd, tmp2, tmp3;
1525 mpz_init(gcd);
1526 mpz_init(tmp2);
1527 mpz_init(tmp3);
1528 for(unsigned int j = 0 ; j < n ; j++) {
1529 unsigned int i1 = M->m;
1530 mpz_set(gcd, p);
1531 for(unsigned int i = rank ; i < M->m ; i++) {
1532 mpz_gcd(tmp2, mpz_mat_entry(M, i, j), p);
1533 if (mpz_cmp(tmp2, gcd) < 0) {
1534 i1 = i;
1535 mpz_set(gcd, tmp2);
1536 if (mpz_cmp_ui(gcd, 1) == 0) break;
1537 }
1538 }
1539 if (i1 == M->m) /* this column is zero */
1540 continue;
1541 mpz_mat_swaprows(M, rank, i1);
1542 if (T) mpz_mat_swaprows(T, rank, i1);
1543 i1 = rank;
1544 /* canonicalize this row */
1545 /* gcd is gcd(M[rank, j], p) */
1546 mpz_divexact(tmp2, mpz_mat_entry(M, rank, j), gcd);
1547 mpz_divexact(tmp3, p, gcd);
1548 mpz_invert(tmp2, tmp2, tmp3);
1549 mpz_mat_mulrow_mod(M, rank, tmp2, p);
1550 if (T) mpz_mat_mulrow_mod(T, rank, tmp2, p);
1551 ASSERT_ALWAYS(mpz_cmp(mpz_mat_entry(M, rank, j), gcd) == 0);
1552 for(unsigned int i0 = 0 ; i0 < m ; i0++) {
1553 if (i0 == rank) continue;
1554 /* we've normalized row rank, so all it takes is to multiply by
1555 * M[i0, j]
1556 */
1557 mpz_fdiv_q(tmp2, mpz_mat_entry(M, i0, j), gcd);
1558 /* for i0 > rank, we should have exact quotient */
1559 mpz_mat_submulrow_mod(M, i0, rank, tmp2, p);
1560 if (T) mpz_mat_submulrow_mod(T, i0, rank, tmp2, p);
1561 ASSERT_ALWAYS(mpz_cmp_ui(mpz_mat_entry(M, i0, j), 0) == 0);
1562 }
1563 rank++;
1564 }
1565 mpz_clear(gcd);
1566 mpz_clear(tmp2);
1567 mpz_clear(tmp3);
1568 }
1569 /* }}} */
1570
mpz_mat_gauss_backend_mod_ui(mpz_mat_ptr M,mpz_mat_ptr T,unsigned long p)1571 void mpz_mat_gauss_backend_mod_ui(mpz_mat_ptr M, mpz_mat_ptr T, unsigned long p)
1572 {
1573 mpz_t pz;
1574 mpz_init_set_ui(pz, p);
1575 mpz_mat_gauss_backend_mod_mpz(M, T, pz);
1576 mpz_clear(pz);
1577 }
1578
1579
1580 /* {{{ integer heap routines */
1581 /* input: heap[0..n-1[ a heap, and heap[n-1] unused.
1582 * output: heap[0..n[ a heap, includes x
1583 */
1584 typedef int (*cmp_t)(void *, unsigned int, unsigned int);
1585
push_heap(unsigned int * heap,unsigned int x,unsigned int n,cmp_t cmp,void * arg)1586 static void push_heap(unsigned int * heap, unsigned int x, unsigned int n, cmp_t cmp, void * arg)
1587 {
1588 unsigned int i = n - 1;
1589 for(unsigned int i0; i; i = i0) {
1590 i0 = (i-1) / 2;
1591 if (cmp(arg, x, heap[i0]) > 0)
1592 heap[i] = heap[i0];
1593 else
1594 break;
1595 }
1596 heap[i] = x;
1597 }
1598
1599 /* put heap[0] in x, and keep the heap property on the resulting
1600 * heap[0..n-1[ ; heap[n-1] is unused on output.
1601 */
pop_heap(unsigned int * x,unsigned int * heap,unsigned int n,cmp_t cmp,void * arg)1602 static void pop_heap(unsigned int * x, unsigned int * heap, unsigned int n, cmp_t cmp, void * arg)
1603 {
1604 *x = heap[0];
1605 unsigned int y = heap[n-1];
1606 /* make heap[0..n-1[ a heap, now */
1607 /* couldn't there be various strategies here ? */
1608 unsigned int i = 0;
1609 for(unsigned int i1 ; 2*i+2 < n - 1 ; i = i1) {
1610 unsigned int left = 2*i + 1;
1611 unsigned int right = 2*i + 2;
1612 if (cmp(arg, heap[left], heap[right]) > 0)
1613 i1 = left;
1614 else
1615 i1 = right;
1616 if (cmp(arg, y, heap[i1]) > 0)
1617 break;
1618 heap[i] = heap[i1];
1619 }
1620 /* there we know that i has no right child, maybe even no left child.
1621 */
1622 if (2*i + 1 < n - 1 && cmp(arg, y, heap[2*i+1]) <= 0) {
1623 heap[i] = heap[2*i+1];
1624 i=2*i+1;
1625 }
1626 heap[i] = y;
1627 }
1628
1629 /* make heap[0..n[ a heap. */
make_heap(unsigned int * heap,unsigned int n,cmp_t cmp,void * arg)1630 static void make_heap(unsigned int * heap, unsigned int n, cmp_t cmp, void * arg)
1631 {
1632 for(unsigned int i = 1 ; i <= n ; i++) {
1633 unsigned int t = heap[i-1];
1634 push_heap(heap, t, i, cmp, arg);
1635 }
1636 }
1637
1638 /* must be a heap already */
sort_heap(unsigned int * heap,unsigned int n,cmp_t cmp,void * arg)1639 static void sort_heap(unsigned int * heap, unsigned int n, cmp_t cmp, void * arg)
1640 {
1641 for(unsigned int i = n ; i ; i--) {
1642 unsigned int t;
1643 pop_heap(&t, heap, i, cmp, arg);
1644 heap[i-1] = t;
1645 }
1646 }
1647
sort_heap_inverse(unsigned int * heap,unsigned int n,cmp_t cmp,void * arg)1648 static void sort_heap_inverse(unsigned int * heap, unsigned int n, cmp_t cmp, void * arg)
1649 {
1650 sort_heap(heap, n, cmp, arg);
1651 for(unsigned int i = 0, j = n-1 ; i < j ; i++, j--) {
1652 int a = heap[i];
1653 int b = heap[j];
1654 heap[i] = b;
1655 heap[j] = a;
1656 }
1657 }
1658 /* }}} */
1659 /* {{{ Hermite normal form over the integers
1660 * The algorithm here is randomly cooked from a stupid idea I had. The
1661 * transformation matrices produced seem to be reasonably good, somewhat
1662 * better than what magma gives (for overdetermined matrices, of course
1663 * -- otherwise they're unique). However we do so in an embarrassingly
1664 * larger amount of time (more than 1000* for 100x100 matrices).
1665 * Fortunately, our intent is not to run this code on anything larger
1666 * than 30*30, in which case it's sufficient.
1667 */
1668 /*{{{ hnf helper algorithm (for column matrices) */
1669 struct mpz_mat_hnf_helper_heap_aux {
1670 double * dd;
1671 mpz_mat_ptr A;
1672 };
1673
mpz_mat_hnf_helper_heap_aux_cmp(struct mpz_mat_hnf_helper_heap_aux * data,unsigned int a,unsigned int b)1674 static int mpz_mat_hnf_helper_heap_aux_cmp(struct mpz_mat_hnf_helper_heap_aux * data, unsigned int a, unsigned int b)
1675 {
1676 /* The heap will have the largest elements on top according to this
1677 * measure. Thus the comparison on the ->a field must be as per
1678 * mpz_cmp
1679 */
1680 mpz_srcptr xa = mpz_mat_entry_const(data->A, a, 0);
1681 mpz_srcptr xb = mpz_mat_entry_const(data->A, b, 0);
1682
1683 ASSERT(mpz_sgn(xa) >= 0);
1684 ASSERT(mpz_sgn(xb) >= 0);
1685 int r = mpz_cmp(xa, xb);
1686 if (r) return r;
1687 /* among combinations giving the same result, we want the "best" to
1688 * be the one with *smallest* norm. So instead of the sign of da-db,
1689 * we want the sign of db-da, here.
1690 */
1691 double da = data->dd[a];
1692 double db = data->dd[b];
1693 return (db > da) - (da > db);
1694 }
1695
1696 /* this is quite the same as computing the hnf on a column matrix, with
1697 * the added functionality that entries in a[0..n0[ are reduced with
1698 * respect to the entries in a[n0..n[. This is equivalent to saying that
1699 * we apply a transformation matrix which is:
1700 * I_n0 R
1701 * T == 0 S
1702 * where S is unimodular. The resulting vector T*A has entries [n0..n[
1703 * equal to (g,0,...,0), while entries [0..n0[ are reduced modulo g.
1704 *
1705 * return +1 or -1 which is the determinant of T.
1706 */
mpz_mat_hnf_helper(mpz_mat_ptr T,mpz_mat_ptr dT,mpz_ptr * a,unsigned int n0,unsigned int n)1707 static int mpz_mat_hnf_helper(mpz_mat_ptr T, mpz_mat_ptr dT, mpz_ptr * a, unsigned int n0, unsigned int n)
1708 {
1709 int signdet=1;
1710 ASSERT_ALWAYS(dT != T);
1711 mpz_mat_realloc(dT, n, n);
1712 mpz_mat_set_ui(dT, 1);
1713
1714 if (n == n0) return signdet;
1715
1716 mpz_t q, r2;
1717 mpz_mat A;
1718 mpz_mat_init(A, n, 1);
1719 mpz_init(q);
1720 mpz_init(r2);
1721 unsigned int * heap = (unsigned int *) malloc(n * sizeof(unsigned int));
1722 struct mpz_mat_hnf_helper_heap_aux cmpdata[1];
1723 cmpdata->A = A;
1724 cmpdata->dd = (double*) malloc(n * sizeof(double));
1725 cmp_t cmp = (cmp_t) &mpz_mat_hnf_helper_heap_aux_cmp;
1726 for(unsigned int i = 0 ; i < n ; i++) {
1727 mpz_ptr Ai = mpz_mat_entry(A, i, 0);
1728 mpz_set(Ai, a[i]);
1729 cmpdata->dd[i] = 0;
1730 for(unsigned int j = 0 ; j < n ; j++) {
1731 double d = mpz_get_d(mpz_mat_entry(T, i, j));
1732 cmpdata->dd[i] += d * d;
1733 }
1734 heap[i] = i;
1735 if (i < n0)
1736 continue;
1737 if (mpz_sgn(Ai) == -1) {
1738 mpz_neg(Ai, Ai);
1739 signdet = -signdet;
1740 mpz_set_si(mpz_mat_entry(dT, i, i), -1);
1741 for(unsigned int j = 0 ; j < n ; j++) {
1742 mpz_neg(mpz_mat_entry(T, i, j), mpz_mat_entry(T, i, j));
1743 }
1744 }
1745 }
1746 make_heap(heap + n0, n - n0, cmp, cmpdata);
1747 unsigned int head;
1748 for( ; ; ) {
1749 /* extract the head */
1750 pop_heap(&head, heap + n0, n - n0, cmp, cmpdata);
1751 mpz_ptr r0 = mpz_mat_entry(A, head, 0);
1752 if (mpz_cmp_ui(r0, 0) == 0) {
1753 /* we're done ! */
1754 push_heap(heap + n0, head, n - n0, cmp, cmpdata);
1755 break;
1756 }
1757 /* reduce A[0..n0[ wrt r0 == A[head] */
1758 for(unsigned int i = 0 ; i < n0 ; i++) {
1759 mpz_fdiv_qr(q, r2, mpz_mat_entry(A, i, 0), r0);
1760 mpz_swap(r2, mpz_mat_entry(A, i, 0));
1761 mpz_mat_submulrow(T, i, head, q);
1762 mpz_mat_submulrow(dT, i, head, q);
1763 }
1764 if (n0 == n - 1) {
1765 /* we're actually computing the gcd of one single integer.
1766 * That makes no sense of course.
1767 */
1768 push_heap(heap + n0, head, n - n0, cmp, cmpdata);
1769 break;
1770 }
1771 mpz_ptr r1 = mpz_mat_entry(A, heap[n0], 0);
1772 if (mpz_cmp_ui(r1, 0) == 0) {
1773 /* we're done ! */
1774 push_heap(heap + n0, head, n - n0, cmp, cmpdata);
1775 break;
1776 }
1777 mpz_fdiv_qr(q, r2, r0, r1);
1778 mpz_swap(r0, r2);
1779 mpz_mat_submulrow(T, head, heap[n0], q);
1780 mpz_mat_submulrow(dT, head, heap[n0], q);
1781 /* adjust the norm of the row in T */
1782 cmpdata->dd[head] = 0;
1783 for(unsigned int j = 0 ; j < n ; j++) {
1784 double d = mpz_get_d(mpz_mat_entry(T, head, j));
1785 cmpdata->dd[head] += d * d;
1786 }
1787 push_heap(heap + n0, head, n - n0, cmp, cmpdata);
1788 }
1789 sort_heap_inverse(heap + n0, n - n0, cmp, cmpdata);
1790 mpz_mat_permuterows(T, heap);
1791 mpz_mat_permuterows(dT, heap);
1792 mpz_mat_permuterows(A, heap);
1793 signdet*= permutation_signature(heap, n);
1794 for(unsigned int i = 0 ; i < n ; i++) {
1795 mpz_ptr Ai = mpz_mat_entry(A, i, 0);
1796 mpz_swap(a[i], Ai);
1797 }
1798 free(cmpdata->dd);
1799 mpz_mat_clear(A);
1800 mpz_clear(q);
1801 mpz_clear(r2);
1802 free(heap);
1803
1804 return signdet;
1805 }
1806 /* {{{ this computes the row hnf on a column C.
1807 * The result is always of the form C'=(gcd(C),0,0,...,0). The transform
1808 * matrix T such that C'=T*C is most interesting. a is modified in
1809 * place.
1810 */
mpz_gcd_many(mpz_mat_ptr dT,mpz_ptr * a,unsigned int n)1811 void mpz_gcd_many(mpz_mat_ptr dT, mpz_ptr * a, unsigned int n)
1812 {
1813 mpz_mat T;
1814 mpz_mat_init(T, n, n);
1815 mpz_mat_set_ui(T, 1);
1816 mpz_mat_hnf_helper(T, dT, a, 0, n);
1817 mpz_mat_clear(T);
1818 }
1819 /*}}}*/
1820 /*}}}*/
1821
1822 /* return +1 or -1, which is the determinant of the transformation matrix
1823 * T */
mpz_mat_hnf_backend(mpz_mat_ptr M,mpz_mat_ptr T)1824 int mpz_mat_hnf_backend(mpz_mat_ptr M, mpz_mat_ptr T)
1825 {
1826 ASSERT_ALWAYS(M != T);
1827 if (T == NULL) {
1828 mpz_mat xT;
1829 mpz_mat_init(xT,0,0);
1830 int r = mpz_mat_hnf_backend(M, xT);
1831 mpz_mat_clear(xT);
1832 return r;
1833 }
1834 int signdet = 1;
1835 unsigned int m = M->m;
1836 unsigned int n = M->n;
1837 mpz_mat_realloc(T, m, m);
1838 mpz_mat_set_ui(T, 1);
1839 unsigned int rank = 0;
1840 mpz_t quo;
1841 mpz_init(quo);
1842 mpz_ptr * colm = (mpz_ptr *) malloc(m * sizeof(mpz_t));
1843 mpz_mat dT, Mx, My;
1844 mpz_mat_init(dT, 0, 0);
1845 mpz_mat_init(Mx, 0, 0);
1846 mpz_mat_init(My, 0, 0);
1847 for(unsigned int j = 0 ; j < n && rank < m; j++) {
1848 for(unsigned int i = 0 ; i < m ; i++) {
1849 colm[i] = mpz_mat_entry(M, i, j);
1850 }
1851 signdet *= mpz_mat_hnf_helper(T, dT, colm, rank, m);
1852
1853 /* apply dT to the submatrix of size m * (n - 1 - j) of M */
1854 mpz_mat_realloc(Mx, m, n - 1 - j);
1855 mpz_mat_submat_swap(Mx, 0, 0, M, 0, j + 1, m, n - 1 - j);
1856 mpz_mat_mul(My, dT, Mx);
1857 mpz_mat_submat_swap(M, 0, j + 1, My, 0, 0, m, n - 1 - j);
1858
1859 mpz_srcptr gcd = colm[rank];
1860 if (mpz_cmp_ui(gcd, 0) == 0) /* this column is zero -- nothing to do */
1861 continue;
1862 rank++;
1863 }
1864 free(colm);
1865 mpz_clear(quo);
1866 mpz_mat_clear(dT);
1867 mpz_mat_clear(Mx);
1868 mpz_mat_clear(My);
1869
1870 return signdet;
1871 }
1872 /* }}} */
1873 /*{{{ kernel*/
1874 // This is supposed to compute the Kernel of M mod p and to store it in
1875 // the matrix K. If r is the rank of M, and M is a square matrix n*n, K
1876 // is a n*(n-r) matrix
1877 //
1878 // self-assignment is ok.
mpz_mat_kernel_mod_mpz(mpz_mat_ptr K,mpz_mat_srcptr M,mpz_srcptr p)1879 void mpz_mat_kernel_mod_mpz(mpz_mat_ptr K, mpz_mat_srcptr M, mpz_srcptr p)
1880 {
1881 mpz_mat T, H;
1882 unsigned int r;
1883
1884 mpz_mat_init(T,M->m,M->m);
1885 mpz_mat_init(H,M->m,M->n);
1886 mpz_mat_set(H,M);
1887
1888 // Storing a reduced matrix of M in H with gauss
1889 mpz_mat_gauss_backend_mod_mpz(H,T,p);
1890
1891 // Finding the rank of M and H
1892 r = H->m;
1893 while((r>0) && mpz_mat_isnull_row(H,r-1)){
1894 r--;
1895 }
1896 // Kernel is of dimension n-r, and a basis of the kernel is given in the n-r last rows of T
1897 // We shall keep the convention of magma
1898
1899 // Reallocating K with n-r columns
1900 mpz_mat_realloc(K, H->m-r, H->m);
1901 mpz_mat_submat_swap( K, 0, 0,
1902 T, r, 0,
1903 H->m - r, H->m);
1904
1905 mpz_mat_clear(T);
1906 mpz_mat_clear(H);
1907 }
mpz_mat_kernel_mod_ui(mpz_mat_ptr K,mpz_mat_srcptr M,unsigned long p)1908 void mpz_mat_kernel_mod_ui(mpz_mat_ptr K, mpz_mat_srcptr M, unsigned long p)
1909 {
1910 mpz_t p2;
1911 mpz_init_set_ui(p2,p);
1912 mpz_mat_kernel_mod_mpz(K,M,p2);
1913 mpz_clear(p2);
1914 }
1915 /* }}} */
1916 /*{{{ inversion*/
mpq_mat_inv(mpq_mat_ptr dst,mpq_mat_srcptr src)1917 void mpq_mat_inv(mpq_mat_ptr dst, mpq_mat_srcptr src)
1918 {
1919 ASSERT_ALWAYS(src->m == src->n);
1920
1921 mpq_mat aux;
1922 mpq_mat_init(aux,src->m,src->n);
1923 mpq_mat_set(aux,src);
1924
1925 /* This is self-assignment compatible (done reading src before
1926 * dst is touched). */
1927 mpq_mat_gauss_backend(aux,dst);
1928
1929 mpq_mat_clear(aux);
1930 }
1931 /* }}} */
1932
1933 /*
1934 * For now, it is just a wrapper to use utils/lll.c.
1935 */
mpz_mat_LLL(mpz_ptr det,mpz_mat_ptr M,mpz_mat_ptr U,mpz_srcptr a,mpz_srcptr b)1936 void mpz_mat_LLL(mpz_ptr det, mpz_mat_ptr M, mpz_mat_ptr U, mpz_srcptr a,
1937 mpz_srcptr b)
1938 {
1939 mat_Z M_tmp;
1940 LLL_init(&M_tmp, M->m, M->n);
1941 for (unsigned int row = 1; row < M->m + 1; row++) {
1942 for (unsigned int col = 1; col < M->n + 1; col++) {
1943 mpz_set(M_tmp.coeff[row][col], mpz_mat_entry_const(M, row - 1, col - 1));
1944 }
1945 }
1946
1947 if (U) {
1948 mat_Z U_tmp;
1949 LLL_init(&U_tmp, U->m, U->m);
1950
1951 LLL(det, M_tmp, &U_tmp, a, b);
1952
1953 for (unsigned int row = 1; row < U->m + 1; row++) {
1954 for (unsigned int col = 1; col < U->n + 1; col++) {
1955 mpz_set(mpz_mat_entry(U, row - 1, col - 1), U_tmp.coeff[row][col]);
1956 }
1957 }
1958 LLL_clear(&U_tmp);
1959 } else {
1960 LLL(det, M_tmp, NULL, a, b);
1961 }
1962
1963 for (unsigned int row = 1; row < M->m + 1; row++) {
1964 for (unsigned int col = 1; col < M->n + 1; col++) {
1965 mpz_set(mpz_mat_entry(M, row - 1, col - 1), M_tmp.coeff[row][col]);
1966 }
1967 }
1968
1969 LLL_clear(&M_tmp);
1970 }
1971
1972 using namespace std;
1973
operator <<(ostream & os,cxx_mpz_mat const & M)1974 ostream& operator<<(ostream& os, cxx_mpz_mat const& M)/*{{{*/
1975 {
1976 for(unsigned int i = 0 ; i < M->m ; i++) {
1977 for(unsigned int j = 0 ; j < M->n ; j++) {
1978 if (j) os << " ";
1979 os << mpz_mat_entry_const(M, i, j);
1980 }
1981 if (i < M->m - 1) os << "\n";
1982 }
1983 return os;
1984 }
1985 /*}}}*/
operator <<(ostream & os,cxx_mpq_mat const & M)1986 ostream& operator<<(ostream& os, cxx_mpq_mat const& M)/*{{{*/
1987 {
1988 for(unsigned int i = 0 ; i < M->m ; i++) {
1989 for(unsigned int j = 0 ; j < M->n ; j++) {
1990 if (j) os << " ";
1991 os << mpq_mat_entry_const(M, i, j);
1992 }
1993 if (i < M->m - 1) os << "\n";
1994 }
1995 return os;
1996 }
1997 /*}}}*/
1998