1 /*
2 Copyright (C) 2012 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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "nmod_vec.h"
17 #include "nmod_mat.h"
18 #include "ulong_extras.h"
19 #include "perm.h"
20
check_rref_form(slong * perm,nmod_mat_t A,slong rank)21 int check_rref_form(slong * perm, nmod_mat_t A, slong rank)
22 {
23 slong i, j, k, prev_pivot;
24
25 /* bottom should be zero */
26 for (i = rank; i < A->r; i++)
27 for (j = 0; j < A->c; j++)
28 if (nmod_mat_entry(A, i, j) != 0)
29 return 0;
30
31 prev_pivot = -1;
32
33 for (i = 0; i < rank; i++)
34 {
35 for (j = 0; j < A->c; j++)
36 {
37 if (nmod_mat_entry(A, i, j) != 0)
38 {
39 /* pivot should have a higher column index than previous */
40 if (j <= prev_pivot)
41 return 0;
42
43 /* column should be 0 ... 0 1 0 ... 0 */
44 for (k = 0; k < rank; k++)
45 if (nmod_mat_entry(A, k, j) != (i == k))
46 return 0;
47
48 prev_pivot = j;
49 break;
50 }
51 }
52 }
53
54 return 1;
55 }
56
57 int
main(void)58 main(void)
59 {
60 slong i;
61
62 FLINT_TEST_INIT(state);
63
64 flint_printf("rref....");
65 fflush(stdout);
66
67 for (i = 0; i < 1000 * flint_test_multiplier(); i++)
68 {
69 nmod_mat_t A, B, C, D;
70 mp_limb_t mod;
71 slong j, k, m, n, rank1, rank2;
72 slong *perm;
73 int equal;
74 mp_limb_t c;
75
76 mod = n_randtest_prime(state, 0);
77
78 m = n_randint(state, 20);
79 n = n_randint(state, 20);
80 perm = _perm_init(2*m);
81
82 nmod_mat_init(A, m, n, mod);
83 nmod_mat_init(D, 2*m, n, mod);
84
85 nmod_mat_randtest(A, state);
86 nmod_mat_init_set(B, A);
87 nmod_mat_init_set(C, A);
88
89 rank1 = nmod_mat_rref(B);
90
91 if (!check_rref_form(perm, B, rank1))
92 {
93 flint_printf("FAIL (malformed rref)\n");
94 nmod_mat_print_pretty(A); flint_printf("\n\n");
95 nmod_mat_print_pretty(B); flint_printf("\n\n");
96 abort();
97 }
98
99 /* Concatenate the original matrix with the rref, scramble the rows,
100 and check that the rref is the same */
101 _perm_randtest(perm, 2 * m, state);
102
103 for (j = 0; j < m; j++)
104 {
105 do { c = n_randint(state, mod); } while (c == 0);
106 for (k = 0; k < n; k++)
107 nmod_mat_entry(D, perm[j], k) =
108 nmod_mul(nmod_mat_entry(A, j, k), c, A->mod);
109 }
110
111 for (j = 0; j < m; j++)
112 {
113 do { c = n_randint(state, mod); } while (c == 0);
114 for (k = 0; k < n; k++)
115 nmod_mat_entry(D, perm[m + j], k) =
116 nmod_mul(nmod_mat_entry(B, j, k), c, A->mod);
117 }
118
119 rank2 = nmod_mat_rref(D);
120 equal = (rank1 == rank2);
121
122 if (equal)
123 {
124 for (j = 0; j < rank2; j++)
125 for (k = 0; k < n; k++)
126 equal = equal && (nmod_mat_entry(B, j, k) ==
127 nmod_mat_entry(D, j, k));
128 for (j = rank2; j < 2 * rank2; j++)
129 for (k = 0; k < n; k++)
130 equal = equal && (nmod_mat_entry(D, j, k) == 0);
131 }
132
133 if (!equal)
134 {
135 flint_printf("FAIL (rank1 = %wd, rank2 = %wd)!\n", rank1, rank2);
136 nmod_mat_print_pretty(A); flint_printf("\n\n");
137 nmod_mat_print_pretty(B); flint_printf("\n\n");
138 nmod_mat_print_pretty(D); flint_printf("\n\n");
139 abort();
140 }
141
142 _perm_clear(perm);
143 nmod_mat_clear(A);
144 nmod_mat_clear(B);
145 nmod_mat_clear(C);
146 nmod_mat_clear(D);
147 }
148
149 FLINT_TEST_CLEANUP(state);
150
151 flint_printf("PASS\n");
152 return 0;
153 }
154