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