1 /*
2 Copyright (C) 2011 Fredrik Johansson
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <stdlib.h>
13 #include <gmp.h>
14 #include "flint.h"
15 #include "nmod_mat.h"
16
17 slong
nmod_mat_nullspace(nmod_mat_t X,const nmod_mat_t A)18 nmod_mat_nullspace(nmod_mat_t X, const nmod_mat_t A)
19 {
20 slong i, j, k, m, n, rank, nullity;
21 slong * p;
22 slong * pivots;
23 slong * nonpivots;
24 nmod_mat_t tmp;
25
26 m = A->r;
27 n = A->c;
28
29 p = flint_malloc(sizeof(slong) * FLINT_MAX(m, n));
30
31 nmod_mat_init_set(tmp, A);
32 rank = nmod_mat_rref(tmp);
33 nullity = n - rank;
34
35 nmod_mat_zero(X);
36
37 if (rank == 0)
38 {
39 for (i = 0; i < nullity; i++)
40 nmod_mat_entry(X, i, i) = UWORD(1);
41 }
42 else if (nullity)
43 {
44 pivots = p; /* length = rank */
45 nonpivots = p + rank; /* length = nullity */
46
47 for (i = j = k = 0; i < rank; i++)
48 {
49 while (nmod_mat_entry(tmp, i, j) == UWORD(0))
50 {
51 nonpivots[k] = j;
52 k++;
53 j++;
54 }
55 pivots[i] = j;
56 j++;
57 }
58 while (k < nullity)
59 {
60 nonpivots[k] = j;
61 k++;
62 j++;
63 }
64
65 for (i = 0; i < nullity; i++)
66 {
67 for (j = 0; j < rank; j++)
68 {
69 mp_limb_t c = nmod_mat_entry(tmp, j, nonpivots[i]);
70 nmod_mat_entry(X, pivots[j], i) = nmod_neg(c, A->mod);
71 }
72
73 nmod_mat_entry(X, nonpivots[i], i) = UWORD(1);
74 }
75 }
76
77 flint_free(p);
78 nmod_mat_clear(tmp);
79
80 return nullity;
81 }
82