1 /* ---------------------------------------------------------------------
2 *
3 * -- Integer Matrix Library (IML)
4 * (C) Copyright 2004, 2006 All Rights Reserved
5 *
6 * -- IML routines -- Version 1.0.1 -- November, 2006
7 *
8 * Author : Zhuliang Chen
9 * Contributor(s) : Arne Storjohann
10 * University of Waterloo -- School of Computer Science
11 * Waterloo, Ontario, N2L3G1 Canada
12 *
13 * ---------------------------------------------------------------------
14 *
15 * -- Copyright notice and Licensing terms:
16 *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
19 * are met:
20 *
21 * 1. Redistributions of source code must retain the above copyright
22 * notice, this list of conditions and the following disclaimer.
23 * 2. Redistributions in binary form must reproduce the above copyright
24 * notice, this list of conditions, and the following disclaimer in
25 * the documentation and/or other materials provided with the distri-
26 * bution.
27 * 3. The name of the University, the IML group, or the names of its
28 * contributors may not be used to endorse or promote products deri-
29 * ved from this software without specific written permission.
30 *
31 * -- Disclaimer:
32 *
33 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
37 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
38 * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
39 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
40 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
41 * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
42 * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
43 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44 *
45 */
46
47
48
49 #include "mtrans.h"
50
51 /*
52 * Calling Sequence:
53 * RowEchelonTransform(p, A, n, m, frows, lrows, redflag, eterm, Q, rp, d)
54 *
55 * Summary:
56 * Compute a mod p row-echelon transform of a mod p input matrix
57 *
58 * Description:
59 * Given a n x m mod p matrix A, a row-echelon transform of A is a 4-tuple
60 * (U,P,rp,d) with rp the rank profile of A (the unique and strictly
61 * increasing list [j1,j2,...jr] of column indices of the row-echelon form
62 * which contain the pivots), P a permutation matrix such that all r leading
63 * submatrices of (PA)[0..r-1,rp] are nonsingular, U a nonsingular matrix
64 * such that UPA is in row-echelon form, and d the determinant of
65 * (PA)[0..r-1,rp].
66 *
67 * Generally, it is required that p be a prime, as inverses are needed, but
68 * in some cases it is possible to obtain an echelon transform when p is
69 * composite. For the cases where the echelon transform cannot be obtained
70 * for p composite, the function returns an error indicating that p is
71 * composite.
72 *
73 * The matrix U is structured, and has last n-r columns equal to the last n-r
74 * columns of the identity matrix, n the row dimension of A.
75 *
76 * The first r rows of UPA comprise a basis in echelon form for the row
77 * space of A, while the last n-r rows of U comprise a basis for the left
78 * nullspace of PA.
79 *
80 * For efficiency, this function does not output an echelon transform
81 * (U,P,rp,d) directly, but rather the expression sequence (Q,rp,d).
82 * Q, rp, d are the form of arrays and pointers in order to operate inplace,
83 * which require to preallocate spaces and initialize them. Initially,
84 * Q[i] = i (i=0..n), rp[i] = 0 (i=0..n), and *d = 1. Upon completion, rp[0]
85 * stores the rank r, rp[1..r] stores the rank profile. i<=Q[i]<=n for
86 * i=1..r. The input Matrix A is modified inplace and used to store U.
87 * Let A' denote the state of A on completion. Then U is obtained from the
88 * identity matrix by replacing the first r columns with those of A', and P
89 * is obtained from the identity matrix by swapping row i with row Q[i], for
90 * i=1..r in succession.
91 *
92 * Parameters flrows, lrows, redflag, eterm control the specific operations
93 * this function will perform. Let (U,P,rp,d) be as constructed above. If
94 * frows=0, the first r rows of U will not be correct. If lrows=0, the last
95 * n-r rows of U will not be correct. The computation can be up to four
96 * times faster if these flags are set to 0.
97 *
98 * If redflag=1, the row-echelon form is reduced, that is (UPA)[0..r-1,rp]
99 * will be the identity matrix. If redflag=0, the row-echelon form will not
100 * be reduced, that is (UPA)[1..r,rp] will be upper triangular and U is unit
101 * lower triangular. If frows=0 then redflag has no effect.
102 *
103 * If eterm=1, then early termination is triggered if a column of the
104 * input matrix is discovered that is linearly dependant on the previous
105 * columns. In case of early termination, the third return value d will be 0
106 * and the remaining components of the echelon transform will not be correct.
107 *
108 * Input:
109 * p: FiniteField, modulus
110 * A: 1-dim Double array length n*m, representation of a n x m input
111 * matrix
112 * n: long, row dimension of A
113 * m: long, column dimension of A
114 * frows: 1/0,
115 * - if frows = 1, the first r rows of U will be correct
116 * - if frows = 0, the first r rows of U will not be correct
117 * lrows: 1/0,
118 * - if lrows = 1, the last n-r rows of U will be correct
119 * - if lrows = 0, the last n-r rows of U will not be correct
120 * redflag: 1/0,
121 * - if redflag = 1, compute row-echelon form
122 * - if redflag = 0, not compute reow-echelon form
123 * eterm: 1/0,
124 * - if eterm = 1, terminate early if not in full rank
125 * - if eterm = 0, not terminate early
126 * Q: 1-dim long array length n+1, compact representation of
127 * permutation vector, initially Q[i] = i, 0 <= i <= n
128 * rp: 1-dim long array length n+1, representation of rank profile,
129 * initially rp[i] = 0, 0 <= i <= n
130 * d: pointer to FiniteField, storing determinant of the matrix,
131 * initially *d = 1
132 *
133 * Precondition:
134 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
135 *
136 */
137
138 void
RowEchelonTransform(const FiniteField p,Double * A,const long n,const long m,const long frows,const long lrows,const long redflag,const long eterm,long * Q,long * rp,FiniteField * d)139 RowEchelonTransform (const FiniteField p, Double *A, const long n, \
140 const long m, const long frows, const long lrows, \
141 const long redflag, const long eterm, long *Q, \
142 long *rp, FiniteField *d)
143 {
144 mpz_t mp_a, mp_p; /* preallocated temporary storage */
145
146 { mpz_init(mp_a); mpz_init_set_ui(mp_p, p); }
147
148 RowEchelonTransform_rec(p, A, n, m, 1, m, 0, 0, frows, lrows, redflag, \
149 eterm, Q, rp, d, mp_a, mp_p);
150
151 { mpz_clear(mp_a); mpz_clear(mp_p); }
152 }
153
154
155 long
RowEchelonTransform_rec(const FiniteField p,Double * A,const long n,const long m,long m1,long m2,long k,const long ks,long frows,long lrows,long redflag,long eterm,long * P,long * rp,FiniteField * d,mpz_t mp_a,mpz_t mp_p)156 RowEchelonTransform_rec (const FiniteField p, Double *A, const long n, \
157 const long m, long m1, long m2, long k, \
158 const long ks, long frows, long lrows, long redflag, \
159 long eterm, long *P, long *rp, FiniteField *d, \
160 mpz_t mp_a, mpz_t mp_p)
161 {
162 long i, j, r1, r2, r, ri, mm, inv;
163 FiniteField a;
164 Double *A1;
165 Double t, b;
166
167 if (m1 == m2)
168 {
169 for (i = k+1; i <= n; i++) { if (*(A+(i-1)*m+m1-1) != 0) { break; } }
170 if ((i > n) && (eterm == 0)) { return 0; }
171 else if ((i > n) && (eterm == 1))
172 {
173 *d = 0;
174 return 0;
175 }
176 if (i > k+1)
177 cblas_dswap(m-m1+1, A+k*m+m1-1, 1, A+(i-1)*m+m1-1, 1);
178 if (k-ks > 0)
179 cblas_dswap(k-ks, A+k*m, 1, A+(i-1)*m, 1);
180 P[k+1] = i;
181 mpz_set_d(mp_a, *(A+m*k+m1-1));
182 inv = mpz_invert(mp_a, mp_a, mp_p);
183
184 /* mp_a is not relatively prime with modulus */
185 if (!inv) { iml_fatal("in RowEchelonTransform: modulus is composite"); }
186 a = mpz_get_ui(mp_a);
187 b = fmod(*(A+k*m+m1-1), p);
188 if (b < 0) { b = p+b; }
189 if ((frows == 1) && (redflag == 1))
190 {
191 for (j = 1; j <= n; j++)
192 *(A+(j-1)*m+k-ks) = *(A+(j-1)*m+m1-1)*(p-a);
193 Dmod(p, A+k-ks, n, 1, m);
194 *(A+k*m+k-ks) = a;
195 }
196 else
197 {
198 if (k+2 <= n)
199 {
200 for (j = k+2; j <= n; j++)
201 *(A+(j-1)*m+k-ks) = *(A+(j-1)*m+m1-1)*(p-a);
202 Dmod(p, A+(k+1)*m+k-ks, n-k-1, 1, m);
203 }
204 if (k > 0)
205 {
206 for (j = 1; j <= k; j++)
207 *(A+(j-1)*m+k-ks) = 0;
208 }
209 *(A+k*m+k-ks) = 1;
210 }
211 ++rp[0];
212 *d = fmod( (*d)*b, p);
213 ri = rp[0];
214 rp[ri] = m1;
215 return 1;
216 }
217 /* Recursively solve the first subproblem */
218 mm = m1+(long)((m2-m1)/2);
219 r1 = RowEchelonTransform_rec(p, A, n, m, m1, mm, k, ks, frows, 1, redflag,\
220 eterm, P, rp, d, mp_a, mp_p);
221 if ((eterm == 1) && (r1 < mm-m1+1))
222 {
223 *d = 0;
224 return 0;
225 }
226 /* If r1=0 then don't need to construct second subproblem */
227 if (r1 > 0)
228 {
229 /* Compute U1.A2 by submatrix multiply */
230 if (k+r1 < n)
231 {
232 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n-k-r1, \
233 m2-mm, r1, 1.0, A+(k+r1)*m+k-ks, m, A+k*m+mm, m, \
234 1.0, A+(k+r1)*m+mm, m);
235 Dmod(p, A+(k+r1)*m+mm, n-k-r1, m2-mm, m);
236 }
237 if ((frows == 1) && (redflag == 1))
238 {
239 if (k > 0)
240 {
241 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, k, \
242 m2-mm, r1, 1.0, A+k-ks, m, A+k*m+mm, m, 1.0, \
243 A+mm, m);
244 Dmod(p, A+mm, k, m2-mm, m);
245 }
246 A1 = XMALLOC(Double, r1*(m2-mm));
247 DCopy(r1, m2-mm, A+k*m+mm, m, A1, m2-mm);
248 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, r1, m2-mm, \
249 r1, 1.0, A+k*m+k-ks, m, A1, m2-mm, 0.0, A+k*m+mm, m);
250 XFREE(A1);
251 Dmod(p, A+k*m+mm, r1, m2-mm, m);
252 }
253 }
254 /* Recursively solve the second subproblem */
255 r2 = RowEchelonTransform_rec(p, A, n, m, mm+1, m2, k+r1, ks, frows, lrows, \
256 redflag, eterm, P, rp, d, mp_a, mp_p);
257 r = r1+r2;
258 if ((eterm == 1) && (r < m2-m1+1))
259 {
260 *d = 0;
261 return 0;
262 }
263 /* Only need to combine when both subproblems nontrivial */
264 if ((r2 > 0) && (r1 > 0))
265 {
266 if ((k+r+1 <= n) && (lrows == 1))
267 {
268 /* Bottom block of U */
269 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n-k-r, \
270 r1, r-r1, 1.0, A+(k+r)*m+k-ks+r1, m, A+(k+r1)*m+k-ks, \
271 m, 1.0, A+(k+r)*m+k-ks, m);
272 Dmod(p, A+(k+r)*m+k-ks, n-k-r, r1, m);
273 }
274 if (frows == 1)
275 {
276 if (redflag == 1)
277 i = 1;
278 else
279 i = k+1;
280
281 /* Rows i..k of top block of U plus first r1 rows of middle block */
282 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, k+r1-i+1, \
283 r1, r-r1, 1.0, A+(i-1)*m+k-ks+r1, m, A+(k+r1)*m+k-ks, \
284 m, 1.0, A+(i-1)*m+k-ks, m);
285 Dmod(p, A+(i-1)*m+k-ks, k+r1-i+1, r1, m);
286
287 /* Last r2 rows of middle block */
288 A1 = XMALLOC(Double, r1*(r-r1));
289 DCopy(r-r1, r1, A+(k+r1)*m+k-ks, m, A1, r1);
290 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, r-r1, r1, \
291 r-r1, 1.0, A+(k+r1)*m+k-ks+r1, m, A1, r1, \
292 0.0, A+(k+r1)*m+k-ks, m);
293 XFREE(A1);
294 Dmod(p, A+(k+r1)*m+k-ks, r-r1, r1, m);
295 }
296 }
297 return r;
298 }
299
300
301
302 /*
303 * Calling Sequence:
304 * Adj <-- mAdjoint(p, A, n)
305 *
306 * Summary:
307 * Compute the adjoint of a mod p square matrix
308 *
309 * Description:
310 * Given a n x n mod p matrix A, the function computes adjoint of A. Input
311 * A is not modified upon completion.
312 *
313 * Input:
314 * p: FiniteField, prime modulus
315 * if p is a composite number, the routine will still work if no error
316 * message is returned
317 * A: 1-dim Double array length n*n, representation of a n x n mod p matrix.
318 * The entries of A are casted from integers
319 * n: long, dimension of A
320 *
321 * Return:
322 * 1-dim Double matrix length n*n, repesentation of a n x n mod p matrix,
323 * adjoint of A
324 *
325 * Precondition:
326 * n*(p-1)^2 <= 2^53-1 = 9007199254740991
327 *
328 */
329
330 Double *
mAdjoint(const FiniteField p,Double * A,const long n)331 mAdjoint (const FiniteField p, Double *A, const long n)
332 {
333 long i, j, k, r, count=0;
334 long *P, *rp;
335 FiniteField det, d[1];
336 Double p1, *B, *C;
337
338 p1 = (Double)p;
339 P = XMALLOC(long, n+1);
340 for (i = 0; i < n+1; i++) { P[i] = i; }
341 rp = XCALLOC(long, n+1);
342 d[0] = 1;
343 B = XMALLOC(Double, n*n);
344 DCopy(n, n, A, n, B, n);
345 RowEchelonTransform(p, B, n, n, 1, 1, 1, 0, P, rp, d);
346 det = d[0];
347 r = rp[0];
348 if (r < n-1)
349 {
350 for (i = 0; i < n*n; i++) { B[i] = 0; }
351 { XFREE(P); XFREE(rp); }
352 return B;
353 }
354 if (r == n)
355 {
356 for (i = r; i > 0; i--)
357 {
358 if (P[i] != i)
359 {
360 ++count;
361 cblas_dswap(n, B+i-1, n, B+P[i]-1, n);
362 }
363 }
364 if (count % 2 == 0)
365 for (i = 0; i < n*n; i++) { B[i] = fmod(det*B[i], p1); }
366 else
367 for (i = 0; i < n*n; i++) { B[i] = fmod((p-det)*B[i], p1); }
368 { XFREE(P); XFREE(rp); }
369 return B;
370 }
371 else
372 {
373 if (n == 1)
374 {
375 B[0] = 1;
376 { XFREE(P); XFREE(rp); }
377 return B;
378 }
379 for (i = 0; i < n; i++) { B[i*n+n-1] = 0; }
380 B[(n-1)*n+(n-1)] = 1;
381 for (i = r; i > 0; i--)
382 {
383 if (P[i] != i)
384 {
385 ++count;
386 cblas_dswap(n, B+i-1, n, B+P[i]-1, n);
387 }
388 }
389 for (j = 1; j < r+1; j++) { if (j != rp[j]) { break; } }
390 C = XMALLOC(Double, n);
391 cblas_dgemv(CblasRowMajor, CblasNoTrans, n, n, 1.0, B, n, \
392 A+j-1, n, 0.0, C, 1);
393 Dmod(p, C, 1, n, 1);
394 for (i = 0; i < n-1; i++) { C[i] = fmod((p-1)*C[i], p1); }
395 for (i = 0; i < n-1; i++) { for (k = 0; k < n; k++) { B[i*n+k] = 0; } }
396 cblas_dger(CblasRowMajor, n-1, n, 1.0, C, 1, B+(n-1)*n, 1, B, n);
397 Dmod(p, B, n-1, n, n);
398 for (i = n-2; i > j-2; i--)
399 {
400 ++count;
401 cblas_dswap(n, B+i*n, 1, B+(i+1)*n, 1);
402 }
403 if (count % 2 == 0)
404 for (i = 0; i < n*n; i++) { B[i] = fmod(det*B[i], p1); }
405 else
406 for (i = 0; i < n*n; i++) { B[i] = fmod((p-det)*B[i], p1); }
407 { XFREE(P); XFREE(rp); XFREE(C); }
408 return B;
409 }
410 }
411
412
413
414
415 /*
416 * Calling Sequence:
417 * r/-1 <-- mBasis(p, A, n, m, basis, nullsp, B, N)
418 *
419 * Summary:
420 * Compute a basis for the rowspace and/or a basis for the left nullspace
421 * of a mod p matrix
422 *
423 * Description:
424 * Given a n x m mod p matrix A, the function computes a basis for the
425 * rowspace B and/or a basis for the left nullspace N of A. Row vectors in
426 * the r x m matrix B consist of basis of A, where r is the rank of A in
427 * Z/pZ. If r is zero, then B will be NULL. Row vectors in the n-r x n
428 * matrix N consist of the left nullspace of A. N will be NULL if A is full
429 * rank.
430 *
431 * The pointers are passed into argument lists to store the computed basis
432 * and nullspace. Upon completion, the rank r will be returned. The
433 * parameters basis and nullsp control whether to compute basis and/or
434 * nullspace. If set basis and nullsp in the way that both basis and
435 * nullspace will not be computed, an error message will be printed and
436 * instead of rank r, -1 will be returned.
437 *
438 * Input:
439 * p: FiniteField, prime modulus
440 * if p is a composite number, the routine will still work if no
441 * error message is returned
442 * A: 1-dim Double array length n*m, representation of a n x m mod p
443 * matrix. The entries of A are casted from integers
444 * n: long, row dimension of A
445 * m: long, column dimension of A
446 * basis: 1/0, flag to indicate whether to compute basis for rowspace or
447 * not
448 * - basis = 1, compute the basis
449 * - basis = 0, not compute the basis
450 * nullsp: 1/0, flag to indicate whether to compute basis for left nullspace
451 * or not
452 * - nullsp = 1, compute the nullspace
453 * - nullsp = 0, not compute the nullspace
454 *
455 * Output:
456 * B: pointer to (Double *), if basis = 1, *B will be a 1-dim r*m Double
457 * array, representing the r x m basis matrix. If basis = 1 and r = 0,
458 * *B = NULL
459 *
460 * N: pointer to (Double *), if nullsp = 1, *N will be a 1-dim (n-r)*n Double
461 * array, representing the n-r x n nullspace matrix. If nullsp = 1 and
462 * r = n, *N = NULL.
463 *
464 * Return:
465 * - if basis and/or nullsp are set to be 1, then return the rank r of A
466 * - if both basis and nullsp are set to be 0, then return -1
467 *
468 * Precondition:
469 * n*(p-1)^2 <= 2^53-1 = 9007199254740991
470 *
471 * Note:
472 * - In case basis = 0, nullsp = 1, A will be destroyed inplace. Otherwise,
473 * A will not be changed.
474 * - Space of B and/or N will be allocated in the function
475 *
476 */
477
478 long
mBasis(const FiniteField p,Double * A,const long n,const long m,const long basis,const long nullsp,Double ** B,Double ** N)479 mBasis (const FiniteField p, Double *A, const long n, const long m, \
480 const long basis, const long nullsp, Double **B, Double **N)
481 {
482 long i, r;
483 long *P, *rp;
484 FiniteField d[1];
485 Double *A1, *B1, *N1, *U;
486
487 P = XMALLOC(long, n+1);
488 for (i = 0; i < n+1; i++) { P[i] = i; }
489 rp = XCALLOC(long, n+1);
490 d[0] = 1;
491 if ((basis == 1) && (nullsp == 1))
492 {
493 A1 = XMALLOC(Double, n*m);
494 DCopy(n, m, A, m, A1, m);
495 RowEchelonTransform(p, A1, n, m, 1, 1, 1, 0, P, rp, d);
496 r = rp[0];
497 U = XCALLOC(Double, n*n);
498 for (i = 0; i < n; i++) { U[i*n+i] = 1; }
499 if (r != 0) { DCopy(n, r, A1, m, U, n); }
500 for (i = r; i > 0; i--)
501 { if (P[i] != i) { cblas_dswap(n, U+i-1, n, U+P[i]-1, n); } }
502 if (r == 0) { *B = NULL; }
503 else
504 {
505 *B = XMALLOC(Double, r*m);
506 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, r, m, n, \
507 1.0, U, n, A, m, 0.0, *B, m);
508 Dmod(p, *B, r, m, m);
509 }
510 if (r ==n) { *N = NULL; }
511 else
512 {
513 *N = XMALLOC(Double, (n-r)*n);
514 DCopy(n-r, n, U+r*n, n, *N, n);;
515 }
516 { XFREE(A1); XFREE(U); XFREE(P); XFREE(rp); }
517 return r;
518 }
519 else if ((basis == 1) && (nullsp == 0))
520 {
521 A1 = XMALLOC(Double, n*m);
522 DCopy(n, m, A, m, A1, m);
523 RowEchelonTransform(p, A1, n, m, 1, 0, 1, 0, P, rp, d);
524 r = rp[0];
525 if (r == 0)
526 {
527 *B = NULL;
528 { XFREE(A1); XFREE(P); XFREE(rp); }
529 return r;
530 }
531 U = XCALLOC(Double, r*n);
532 DCopy(r, r, A1, m, U, n);
533 for (i = r; i > 0; i--)
534 { if (P[i] != i) { cblas_dswap(r, U+i-1, n, U+P[i]-1, n); } }
535 *B = XMALLOC(Double, r*m);
536 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, r, m, n, \
537 1.0, U, n, A, m, 0.0, *B, m);
538 Dmod(p, *B, r, m, m);
539 { XFREE(A1); XFREE(U); XFREE(P); XFREE(rp); }
540 return r;
541 }
542 else if ((basis == 0) && (nullsp == 1))
543 {
544 RowEchelonTransform(p, A, n, m, 0, 1, 1, 0, P, rp, d);
545 r = rp[0];
546 if (r == n)
547 {
548 *N = NULL;
549 { XFREE(P); XFREE(rp); }
550 return r;
551 }
552 *N = XCALLOC(Double, (n-r)*n);
553 if (r != 0) { DCopy(n-r, r, A+r*m, m, *N, n); }
554 for (i = 0; i < n-r; i++) { (*N)[i*n+r+i] = 1; }
555 for (i = r; i > 0; i--)
556 { if (P[i] != i) { cblas_dswap(n-r, *N+i-1, n, *N+P[i]-1, n); } }
557 { XFREE(P); XFREE(rp); }
558 return r;
559 }
560 else
561 {
562 iml_error("In mBasis, both basis and nullsp are zero.");
563 return -1;
564 }
565 }
566
567
568
569 /*
570 * Calling Sequence:
571 * det <-- mDeterminant(p, A, n)
572 *
573 * Summary:
574 * Compute the determinant of a square mod p matrix
575 *
576 * Input:
577 * p: FiniteField, prime modulus
578 * if p is a composite number, the routine will still work if no error
579 * message is returned
580 * A: 1-dim Double array length n*n, representation of a n x n mod p matrix.
581 * The entries of A are casted from integers
582 * n: long, dimension of A
583 *
584 * Output:
585 * det(A) mod p, the determinant of square matrix A
586 *
587 * Precondition:
588 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
589 *
590 * Note:
591 * A is destroyed inplace
592 *
593 */
594 long
mDeterminant(const FiniteField p,Double * A,const long n)595 mDeterminant (const FiniteField p, Double *A, const long n)
596 {
597 long i, count=0;
598 long *P, *rp;
599 FiniteField det, d[1];
600
601 P = XMALLOC(long, n+1);
602 for (i = 0; i < n+1; i++) { P[i] = i; }
603 rp = XCALLOC(long, n+1);
604 d[0] = 1;
605 RowEchelonTransform(p, A, n, n, 0, 0, 0, 1, P, rp, d);
606 det = d[0];
607 if (det != 0)
608 {
609 for (i = 1; i < n+1; i++) { if (P[i] != i) { ++count; } }
610 if (count % 2 == 0)
611 {
612 { XFREE(P); XFREE(rp); }
613 return det;
614 }
615 else
616 {
617 { XFREE(P); XFREE(rp); }
618 return p-det;
619 }
620 }
621 { XFREE(P); XFREE(rp); }
622 return det;
623 }
624
625
626 /*
627 * Calling Sequence:
628 * 1/0 <-- mInverse(p, A, n)
629 *
630 * Summary:
631 * Certified compute the inverse of a mod p matrix inplace
632 *
633 * Description:
634 * Given a n x n mod p matrix A, the function computes A^(-1) mod p
635 * inplace in case A is a nonsingular matrix in Z/Zp. If the inverse does
636 * not exist, the function returns 0.
637 *
638 * A will be destroyed at the end in both cases. If the inverse exists, A is
639 * inplaced by its inverse. Otherwise, the inplaced A is not the inverse.
640 *
641 * Input:
642 * p: FiniteField, prime modulus
643 * if p is a composite number, the routine will still work if no error
644 * message is returned
645 * A: 1-dim Double array length n*n, representation of a n x n mod p matrix.
646 * The entries of A are casted from integers
647 * n: long, dimension of A
648 *
649 * Return:
650 * - 1, if A^(-1) mod p exists
651 * - 0, if A^(-1) mod p does not exist
652 *
653 * Precondition:
654 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
655 *
656 * Note:
657 * A is destroyed inplace
658 *
659 */
660
661 long
mInverse(const FiniteField p,Double * A,const long n)662 mInverse (const FiniteField p, Double *A, const long n)
663 {
664 long i;
665 long *P, *rp;
666 FiniteField d[1];
667
668 P = XMALLOC(long, n+1);
669 for (i = 0; i < n+1; i++) { P[i] = i; }
670 rp = XCALLOC(long, n+1);
671 d[0] = 1;
672 RowEchelonTransform(p, A, n, n, 1, 1, 1, 0, P, rp, d);
673 if (rp[0] == n)
674 {
675 for (i = n; i > 0; i--)
676 if (P[i] != i)
677 cblas_dswap(n, A+i-1, n, A+P[i]-1, n);
678 { XFREE(P); XFREE(rp); }
679 return 1;
680 }
681 { XFREE(P); XFREE(rp); }
682 return 0;
683 }
684
685
686
687
688 /*
689 * Calling Sequence:
690 * r <-- mRank(p, A, n, m)
691 *
692 * Summary:
693 * Compute the rank of a mod p matrix
694 *
695 * Input:
696 * p: FiniteField, prime modulus
697 * if p is a composite number, the routine will still work if no
698 * error message is returned
699 * A: 1-dim Double array length n*m, representation of a n x m mod p
700 * matrix. The entries of A are casted from integers
701 * n: long, row dimension of A
702 * m: long, column dimension of A
703 *
704 * Return:
705 * r: long, rank of matrix A
706 *
707 * Precondition:
708 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
709 *
710 * Note:
711 * A is destroyed inplace
712 *
713 */
714
715 long
mRank(const FiniteField p,Double * A,const long n,const long m)716 mRank (const FiniteField p, Double *A, const long n, const long m)
717 {
718 long i, r;
719 long *P, *rp;
720 FiniteField d[1];
721
722 P = XMALLOC(long, n+1);
723 for (i = 0; i < n+1; i++) { P[i] = i; }
724 rp = XCALLOC(long, n+1);
725 d[0] = 1;
726 RowEchelonTransform(p, A, n, m, 0, 0, 0, 0, P, rp, d);
727 r = rp[0];
728 { XFREE(P); XFREE(rp); }
729 return r;
730 }
731
732
733
734
735 /*
736 * Calling Sequence:
737 * rp <-- mRankProfile(p, A, n, m)
738 *
739 * Summary:
740 * Compute the rank profile of a mod p matrix
741 *
742 * Input:
743 * p: FiniteField, prime modulus
744 * if p is a composite number, the routine will still work if no
745 * error message is returned
746 * A: 1-dim Double array length n*m, representation of a n x m mod p
747 * matrix. The entries of A are casted from integers
748 * n: long, row dimension of A
749 * m: long, column dimension of A
750 *
751 * Return:
752 * rp: 1-dim long array length n+1, where
753 * - rp[0] is the rank of matrix A
754 * - rp[1..r] is the rank profile of matrix A
755 *
756 * Precondition:
757 * ceil(n/2)*(p-1)^2+(p-1) <= 2^53-1 = 9007199254740991 (n >= 2)
758 *
759 * Note:
760 * A is destroyed inplace
761 *
762 */
763
764 long *
mRankProfile(const FiniteField p,Double * A,const long n,const long m)765 mRankProfile (const FiniteField p, Double *A, const long n, const long m)
766 {
767 long i;
768 long *P, *rp;
769 FiniteField d[1];
770
771 P = XMALLOC(long, n+1);
772 for (i = 0; i < n+1; i++) { P[i] = i; }
773 rp = XCALLOC(long, n+1);
774 d[0] = 1;
775 RowEchelonTransform(p, A, n, m, 0, 0, 0, 0, P, rp, d);
776 XFREE(P);
777 return rp;
778 }
779
780
781