1 /*
2 Copyright (C) 2013 Tom Bachmann
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 <iostream>
13 #include <sstream>
14 #include <string>
15
16 #include "fmpq_matxx.h"
17 #include "flintxx/test/helpers.h"
18
19 using namespace flint;
20
21 void
test_init()22 test_init()
23 {
24 fmpq_matxx A(3, 4);
25 tassert(A.rows() == 3 && A.cols() == 4);
26 tassert(A.at(0, 0).is_zero());
27 A.at(0, 0) = 1;
28
29 fmpq_matxx B(A);
30 tassert(B.rows() == 3 && B.cols() == 4);
31 tassert(B.at(0, 0).is_one());
32 B.at(0, 0) = 0;
33 tassert(A.at(0, 0).is_one());
34
35 tassert(fmpq_matxx::zero(3, 4).is_zero());
36 fmpq_matxx eye = fmpq_matxx::one(4, 4);
37 for(slong i = 0;i < eye.rows();++i)
38 for(slong j = 0;j < eye.cols();++j)
39 tassert(eye.at(i, j) == fmpqxx::integer(int(i == j)));
40 }
41
42 void
test_assignment()43 test_assignment()
44 {
45 frandxx state;
46 fmpz_matxx A = fmpz_matxx::randtest(3, 4, state, 10);
47 fmpq_matxx Aq(A.rows(), A.cols());
48 Aq = A;
49 tassert(Aq == fmpq_matxx::integer_matrix(A));
50 }
51
52 void
test_conversion()53 test_conversion()
54 {
55 frandxx state;
56 fmpz_matxx A = fmpz_matxx::randtest(3, 4, state, 10);
57 fmpq_matxx Aq = fmpq_matxx::integer_matrix(A);
58 tassert(Aq.rows() == A.rows() && Aq.cols() == A.cols());
59 for(slong i = 0;i < A.rows();++i)
60 for(slong j = 0;j < A.cols();++j)
61 tassert(Aq.at(i, j) == fmpqxx::integer(A.at(i, j)));
62 tassert(A == fmpz_matxx::from_integral_fraction(Aq));
63 Aq.at(0, 0) = fmpqxx::frac(1, 2);
64 assert_exception(fmpz_matxx::from_integral_fraction(Aq));
65
66 tassert(Aq.numden_entrywise().get<0>().rows() == A.rows()
67 && Aq.numden_entrywise().get<0>().cols() == A.cols()
68 && Aq.numden_entrywise().get<1>().rows() == A.rows()
69 && Aq.numden_entrywise().get<1>().cols() == A.cols());
70
71 tassert(Aq.numden_matwise().get<0>().rows() == A.rows()
72 && Aq.numden_matwise().get<0>().cols() == A.cols());
73
74 tassert(Aq.numden_rowwise().get<0>().rows() == A.rows()
75 && Aq.numden_rowwise().get<0>().cols() == A.cols()
76 && Aq.numden_rowwise().get<1>().size() == A.rows());
77
78 tassert(Aq.numden_colwise().get<0>().rows() == A.rows()
79 && Aq.numden_colwise().get<0>().cols() == A.cols()
80 && Aq.numden_colwise().get<1>().size() == A.cols());
81
82 fmpz_matxx den(A.rows(), A.cols());
83 for(slong i = 0;i < A.rows();++i)
84 for(slong j = 0;j < A.cols();++j)
85 den.at(i, j) = 1u;
86 den.at(0, 0) = 2u;
87
88 A.at(0, 0) = 1;
89 tassert(Aq.numden_entrywise() == ltupleref(A, den));
90 tassert(Aq.num_rowwise().at(0, 0) == 1);
91 tassert(Aq.num_colwise().at(0, 0) == 1);
92 tassert(Aq.numden_matwise().get<1>() == 2);
93
94 fmpz_vecxx rowdens(A.rows());
95 rowdens[0] = 2;
96 for(slong i = 1;i < A.rows();++i)
97 rowdens[i] = 1;
98 for(slong i = 1;i < A.cols();++i)
99 A.at(0, i) *= 2;
100 tassert(Aq.numden_rowwise() == ltupleref(A, rowdens));
101 tassert(Aq.numden_colwise().get<1>()[0] == 2);
102 }
103
104 template<class Expr>
has_explicit_temporaries(const Expr &)105 bool has_explicit_temporaries(const Expr&)
106 {
107 return Expr::ev_traits_t::rule_t::temporaries_t::len != 0;
108 }
109 void
test_arithmetic()110 test_arithmetic()
111 {
112 fmpq_matxx A(10, 10);
113 fmpq_matxx v(10, 1);
114 for(int i = 0;i < 10;++i)
115 v.at(i, 0) = i;
116
117 fmpzxx two(2);
118 tassert(transpose(v).rows() == 1);
119 tassert(v.transpose().cols() == 10);
120 tassert((two*v).rows() == 10);
121 tassert((v*two).rows() == 10);
122 tassert((v*transpose(v)).rows() == 10 && (v*transpose(v)).cols() == 10);
123
124 tassert(!has_explicit_temporaries(trace(transpose(v))));
125 tassert(!has_explicit_temporaries(trace(A + v*transpose(v))));
126 tassert(!has_explicit_temporaries(A + v*transpose(v)));
127 tassert(!has_explicit_temporaries(trace((v*transpose(v) + A))));
128 tassert(!has_explicit_temporaries(trace(v*transpose(v) + v*transpose(v))));
129 tassert(!has_explicit_temporaries(v*transpose(v) + v*transpose(v)));
130
131 tassert(trace(transpose(v)).is_zero());
132 tassert(trace(A + v*transpose(v)) == fmpqxx(285, 1u));
133 tassert(trace(v*transpose(v) + A) == fmpqxx(285, 1u));
134 tassert(trace(v*transpose(v) + v*transpose(v)) == fmpqxx(2*285, 1u));
135 tassert(trace((A+A)*(two + two)).is_zero());
136
137 for(int i = 0;i < 10; ++i)
138 for(int j = 0; j < 10; ++j)
139 A.at(i, j) = i*j;
140 tassert(A == v*transpose(v));
141 tassert(A != transpose(v)*v);
142 A.at(0, 0) = 15;
143 tassert(A != v*transpose(v));
144
145 A.at(0, 0) = 0;
146 for(unsigned i = 0;i < 10; ++i)
147 for(unsigned j = 0; j < 10; ++j)
148 A.at(i, j) *= two;
149 tassert(A == v*transpose(v) + v*transpose(v));
150 tassert(A - v*transpose(v) == v*transpose(v));
151 tassert(((-A) + A).is_zero());
152 tassert((A + A).at(0, 0) == A.at(0, 0) + A.at(0, 0));
153
154 tassert((A + A) == A*two);
155 tassert((two*A) / two == A);
156
157 frandxx state;
158 A.set_randtest(state, 15);
159 fmpz_matxx B(A.rows(), A.cols());
160 B.set_randtest(state, 15);
161 fmpq_matxx C(A.rows(), A.cols());
162 for(slong i = 0;i < A.rows();++i)
163 for(slong j = 0;j < A.cols();++j)
164 C.at(i, j).num() = B.at(i, j);
165 tassert(C*A == B*A && A*C == A*B);
166 tassert(C.mul_direct(A) == C*A && C.mul_cleared(A) == C*A);
167 }
168
169 void
test_functions()170 test_functions()
171 {
172 fmpq_matxx A(2, 3), B(2, 2), empty(0, 15);
173 B.at(0, 0) = 1;
174 tassert(A.is_zero() && !A.is_empty() && !A.is_square());
175 tassert(!B.is_zero() == B.is_square());
176 tassert(empty.is_zero() && empty.is_empty());
177
178 // transpose tested in arithmetic
179 // mul tested in arithmetic
180 // trace tested in arithmetic
181
182 tassert(hilbert_matrix(4, 6).rows() == 4);
183 tassert(hilbert_matrix(4, 6).cols() == 6);
184 A.set_hilbert_matrix();
185 fmpq_matxx H(hilbert_matrix(2, 3));
186 tassert(A == H);
187 for(slong i = 0;i < A.rows();++i)
188 for(slong j = 0;j < A.rows();++j)
189 tassert(A.at(i, j).num() == 1 && A.at(i, j).den() == i + j + 1);
190
191 tassert(A.is_integral() == false);
192
193 frandxx rand;
194 fmpz_matxx Bp(B.rows(), B.cols());
195 Bp.set_randtest(rand, 10);
196 for(slong i = 0;i < B.rows();++i)
197 for(slong j = 0;j < B.rows();++j)
198 B.at(i, j) = fmpqxx(Bp.at(i, j), fmpzxx(1));
199 tassert(B.det().num() == Bp.det() && B.det().den() == 1);
200
201 B.at(0, 0) = fmpqxx::randtest_not_zero(rand, 10);
202 B.at(0, 1) = 0;
203 B.at(1, 0) = fmpqxx::randtest(rand, 10);
204 B.at(1, 1) = fmpqxx::randtest_not_zero(rand, 10);
205
206 tassert(B*B.solve_fraction_free(A) == A);
207 tassert(B*B.solve_dixon(A) == A);
208 fmpq_matxx eye(B.rows(), B.cols());
209 for(slong i = 0;i < B.rows();++i)
210 eye.at(i, i) = 1;
211 tassert(B*B.inv() == eye);
212
213 assert_exception(fmpq_matxx(2, 2).solve_fraction_free(A).evaluate());
214 assert_exception(fmpq_matxx(2, 2).solve_dixon(A).evaluate());
215
216 // make sure this compiles
217 if(0)
218 print(B);
219 }
220
221 void
test_extras()222 test_extras()
223 {
224 fmpq_matxx A(10, 10), B(10, 10);
225 frandxx rand;
226 A.set_randtest(rand, 15);
227 B.set_randtest(rand, 15);
228 A.at(0, 0) = B.at(0, 0) + fmpqxx(1, 1u);
229
230 fmpq_matxx_srcref Asr(A);
231 fmpq_matxx_ref Br(B);
232
233 tassert((A + A) + (B + B) == (Asr + Asr) + (Br + Br));
234
235 Br = Asr;
236 tassert(A == B);
237
238 fmpq_matxx C(Asr);
239 tassert(C == A);
240 C.at(0, 0) += fmpqxx(2, 1u);
241 tassert(C != A);
242 }
243
244 void
test_randomisation()245 test_randomisation()
246 {
247 frandxx rand;
248 fmpq_matxx A(2, 2);
249 A.set_randbits(rand, 5);
250 tassert(height(A.at(0, 0)) <= 31);
251 A.set_randtest(rand, 5);
252 tassert(height(A.at(0, 0)) <= 31);
253
254 frandxx rand2, rand3;
255 A.set_randbits(rand2, 5);
256 tassert(A == fmpq_matxx::randbits(2, 2, rand3, 5));
257 A.set_randtest(rand2, 5);
258 tassert(A == fmpq_matxx::randtest(2, 2, rand3, 5));
259 }
260
261 void
test_reduction_reconstruction()262 test_reduction_reconstruction()
263 {
264 fmpq_matxx A(4, 7);
265 frandxx state;
266 A.set_randtest(state, 4);
267 fmpzxx M(UWORD(123457891));
268 fmpz_matxx Ar = fmpz_matxx::reduce(A, M);
269 tassert(Ar.rows() == A.rows() && Ar.cols() == A.cols());
270 for(slong i = 0;i < A.rows();++i)
271 for(slong j = 0;j < A.cols();++j)
272 tassert(Ar.at(i, j) == A.at(i, j) % M);
273
274 tassert(A == fmpq_matxx::reconstruct(Ar, M));
275 // TODO test exception
276 }
277
278 void
test_row_reduction()279 test_row_reduction()
280 {
281 frandxx state;
282 fmpq_matxx A = fmpq_matxx::randtest(5, 5, state, 15);
283 slong rank1, rank2;
284 permxx p1(5), p2(5);
285 fmpq_matxx res1(A.rows(), A.cols()), res2(A.rows(), A.cols());
286
287 fmpq_matxx B(A);
288 tassert(B.pivot(1, 1, &p1) == fmpq_mat_pivot(p1._data(), A._mat(), 1, 1));
289 tassert(A == B);
290
291 ltupleref(rank1, res1) = rref(A);
292 rank2 = fmpq_mat_rref(res2._mat(), A._mat());
293 tassert(rank1 == rank2 && res1 == res2);
294
295 tassert(rref_classical(A) == rref(A) && rref_fraction_free(A) == rref(A));
296 }
297
298 void
test_unified_access()299 test_unified_access()
300 {
301 fmpq_matxx A(2, 2);
302 const fmpq_matxx& Ar = A;
303 const fmpq_matxx_ref Ar2(A);
304 Ar2.at(0, 0) = fmpqxx::one();
305 tassert(Ar.at(0, 0).is_one());
306 }
307
308 int
main()309 main()
310 {
311 std::cout << "fmpq_matxx....";
312
313 test_init();
314 test_assignment();
315 test_conversion();
316 test_arithmetic();
317 test_functions();
318 test_extras();
319 test_randomisation();
320 test_reduction_reconstruction();
321 test_row_reduction();
322 test_unified_access();
323
324 std::cout << "PASS" << std::endl;
325 return 0;
326 }
327