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 <http://www.gnu.org/licenses/>.
10 */
11 
12 #include <iostream>
13 
14 #include "nmod_poly_matxx.h"
15 
16 #include "flintxx/test/helpers.h"
17 
18 using namespace flint;
19 
20 void
test_init()21 test_init()
22 {
23     mp_limb_t M = 1039;
24     nmod_poly_matxx A(3, 4, M);
25     nmodxx_ctx_srcref ctx = A.estimate_ctx();
26     tassert(ctx.n() == M);
27     tassert((A + A).modulus() == M);
28     tassert(A.rows() == 3 && A.cols() == 4);
29     tassert(A.at(0, 0) == nmod_polyxx::from_ground(0, ctx));
30     A.at(0, 0) = nmod_polyxx::from_ground(1, ctx);
31 
32     nmod_poly_matxx B(A);
33     tassert(A == B);
34     tassert(B.rows() == 3 && B.cols() == 4);
35     tassert(B.at(0, 0) == nmod_polyxx::from_ground(1, ctx));
36     B.at(0, 0) = nmod_polyxx::from_ground(0, ctx);
37     tassert(A.at(0, 0) == nmod_polyxx::from_ground(1, ctx));
38     tassert(A != B);
39 
40     B = A;
41     tassert(A == B);
42 
43     frandxx state;
44     nmod_matxx C(A.rows(), A.cols(), A.modulus());
45     C.set_randtest(state);
46     A = nmod_poly_matxx::from_ground(C);
47     for(slong i = 0;i < A.rows();++i)
48         for(slong j = 0;j < A.cols();++j)
49             tassert(A.at(i, j) == nmod_polyxx::from_ground(C.at(i, j)));
50 
51     tassert(nmod_poly_matxx::zero(2, 2, M).is_zero()
52             && nmod_poly_matxx::one(2, 2, M).is_one());
53 }
54 
55 template<class Expr>
has_explicit_temporaries(const Expr &)56 bool has_explicit_temporaries(const Expr&)
57 {
58     return Expr::ev_traits_t::rule_t::temporaries_t::len != 0;
59 }
60 void
test_arithmetic()61 test_arithmetic()
62 {
63     mp_limb_t M = 1039;
64     nmod_poly_matxx A(10, 10, M);
65     nmod_poly_matxx v(10, 1, M);
66     nmodxx_ctx_srcref ctx = A.estimate_ctx();
67     for(unsigned i = 0;i < 10;++i)
68         v.at(i, 0) = nmod_polyxx::from_ground(i, ctx);
69     nmod_polyxx two = nmod_polyxx::from_ground(2, ctx);
70 
71     tassert(transpose(v).rows() == 1);
72     tassert(v.transpose().cols() == 10);
73     tassert((two*v).rows() == 10);
74     tassert((v*two).rows() == 10);
75     tassert((v*transpose(v)).rows() == 10 && (v*transpose(v)).cols() == 10);
76 
77     tassert(!has_explicit_temporaries(trace(transpose(v))));
78     tassert(!has_explicit_temporaries(trace(A + v*transpose(v))));
79     tassert(!has_explicit_temporaries(A + v*transpose(v)));
80     tassert(!has_explicit_temporaries(trace((v*transpose(v) + A))));
81     tassert(!has_explicit_temporaries(trace(v*transpose(v) + v*transpose(v))));
82     tassert(!has_explicit_temporaries(v*transpose(v) + v*transpose(v)));
83 
84     tassert(trace(transpose(v)) == nmod_polyxx::from_ground(0, ctx));
85     tassert(trace(A + v*transpose(v)) == nmod_polyxx::from_ground(285, ctx));
86     tassert(trace(v*transpose(v) + A) == nmod_polyxx::from_ground(285, ctx));
87     tassert(trace(v*transpose(v) + v*transpose(v))
88             == nmod_polyxx::from_ground(2*285, ctx));
89     tassert(trace((A+A)*(nmod_polyxx::from_ground(1, ctx)
90                     + nmod_polyxx::from_ground(1, ctx)))
91             == nmod_polyxx::from_ground(0, ctx));
92 
93     for(unsigned i = 0;i < 10; ++i)
94         for(unsigned j = 0; j < 10; ++j)
95             A.at(i, j) = nmod_polyxx::from_ground(i*j, ctx);
96     tassert(A == v*transpose(v));
97     tassert(A != transpose(v)*v);
98     A.at(0, 0) = nmod_polyxx::from_ground(15, ctx);
99     tassert(A != v*transpose(v));
100 
101     A.at(0, 0) = nmod_polyxx::from_ground(0, ctx);
102     for(unsigned i = 0;i < 10; ++i)
103         for(unsigned j = 0; j < 10; ++j)
104             A.at(i, j) *= two;
105     tassert(A == v*transpose(v) + v*transpose(v));
106     tassert(A - v*transpose(v) == v*transpose(v));
107     tassert(((-A) + A).is_zero());
108     tassert((A + A).at(0, 0) == A.at(0, 0) + A.at(0, 0));
109 
110     tassert(A * nmodxx::red(17, ctx) == A * nmod_polyxx::from_ground(17, ctx));
111 
112     frandxx rand;
113     nmodxx x = nmodxx::red(17, ctx);
114     A.set_randtest(rand, 5);
115     nmod_matxx B(A.rows(), A.cols(), M);
116     B = A(x);
117     for(slong i = 0;i < A.rows();++i)
118         for(slong j = 0;j < A.cols();++j)
119             tassert(B.at(i, j) == A.at(i, j)(x));
120     tassert(A(x) == evaluate(A, x));
121 }
122 
123 void
test_functions()124 test_functions()
125 {
126     mp_limb_t M = 1031;
127     nmod_poly_matxx A(2, 3, M), B(2, 2, M), empty(0, 15, M);
128     nmodxx_ctx_srcref ctx = A.estimate_ctx();
129     B.at(0, 0) = nmod_polyxx::from_ground(1, ctx);
130     tassert(A.is_zero() && !A.is_empty() && !A.is_square() && !A.is_one());
131     tassert(!B.is_zero() == B.is_square());
132     tassert(empty.is_zero() && empty.is_empty());
133     B.at(1, 1) = B.at(0, 0);
134     tassert(B.is_one());
135 
136     // transpose tested in arithmetic
137     // mul tested in arithmetic
138     // trace tested in arithmetic
139 
140     A.at(0, 0).set_coeff(35, 1);
141     tassert(A.max_length() == 36);
142 
143     frandxx rand;
144     A.set_randtest(rand, 5);
145     B.set_randtest(rand, 5);
146     tassert(B*A == B.mul_classical(A));
147     tassert(B*A == B.mul_KS(A));
148     tassert(B*A == B.mul_interpolate(A));
149 
150     tassert(B.sqr() == B*B);
151     tassert(B.sqr_classical() == B*B);
152     tassert(B.sqr_KS() == B*B);
153     tassert(B.sqr_interpolate() == B*B);
154 
155     tassert(B.pow(5u) == B*B.sqr().sqr());
156 
157     nmod_matxx Bp(B.rows(), B.cols(), B.modulus());
158     Bp.set_randrank(rand, 1);
159     tassert(nmod_poly_matxx::from_ground(Bp).det().is_zero());
160     Bp.set_randrank(rand, 2);
161     tassert(nmod_poly_matxx::from_ground(Bp).det()
162             == nmod_polyxx::from_ground(Bp.det()));
163 
164     Bp.set_randrank(rand, 1);
165     tassert(inv(nmod_poly_matxx::from_ground(Bp)).get<0>() == false);
166 
167     Bp.set_randrank(rand, 2);
168     bool worked;nmod_polyxx den(B.modulus());
169     ltupleref(worked, B, den) = inv(nmod_poly_matxx::from_ground(Bp));
170     tassert(worked && B*nmod_poly_matxx::from_ground(Bp)*A == A*den);
171 
172     tassert(rank(B) == 2);
173 
174     Bp.set_randrank(rand, 1);
175     tassert(nmod_poly_matxx::from_ground(Bp).solve(A).get<0>() == false);
176     Bp.set_randrank(rand, 2);
177     nmod_poly_matxx P(A.rows(), A.cols(), A.modulus());
178     ltupleref(worked, P, den) = nmod_poly_matxx::from_ground(Bp).solve(A);
179     tassert(worked && nmod_poly_matxx::from_ground(Bp)*P == A*den);
180     B = nmod_poly_matxx::from_ground(Bp);
181     tassert(B.solve(A) == B.solve_fflu(A));
182 
183     permxx perm(B.rows());
184     tassert(solve_fflu_precomp(perm, B.fflu(&perm, false).get<1>().evaluate(), A)
185             == B.solve_fflu(A).get<1>());
186 
187     Bp.set_randtest(rand);
188     B = nmod_poly_matxx::from_ground(Bp);
189     slong nullity;nmod_poly_matxx C(2, 2, M);
190     tassert(nullspace(B).get<1>().rows() == 2);
191     tassert(nullspace(B).get<1>().cols() == 2);
192     ltupleref(nullity, C) = nullspace(B);
193     tassert(nullity == 2 - rank(B));
194     tassert(C.rank() == nullity);
195     tassert((B*C).is_zero());
196 
197     B.set_zero();tassert(B.is_zero());
198 }
199 
200 void
test_randomisation()201 test_randomisation()
202 {
203     frandxx rand, rand2;
204     mp_limb_t M = 1031;
205     nmod_poly_matxx A(2, 2, M);
206 
207     A.set_randtest(rand, 17);
208     tassert(A.at(0, 0).length() <= 17);
209     tassert(A == nmod_poly_matxx::randtest(2, 2, M, rand2, 17));
210     A.set_randtest_sparse(rand, 17, 0.5);
211     tassert(A.at(0, 0).length() <= 17);
212     tassert(A == nmod_poly_matxx::randtest_sparse(2, 2, M, rand2, 17, 0.5));
213 }
214 
215 void
test_row_reduction()216 test_row_reduction()
217 {
218     frandxx state;
219     nmod_poly_matxx A = nmod_poly_matxx::randtest(5, 5, 1031, state, 7);
220     slong rank1, rank2;
221     nmod_polyxx den1(A.modulus()), den2(A.modulus());
222     nmod_poly_matxx res1(A.rows(), A.cols(), A.modulus());
223     nmod_poly_matxx res2(res1);
224 
225     tassert(find_pivot_any(A, 2, 4, 1)
226             == nmod_poly_mat_find_pivot_any(A._mat(), 2, 4, 1));
227     tassert(find_pivot_partial(A, 2, 4, 1)
228             == nmod_poly_mat_find_pivot_partial(A._mat(), 2, 4, 1));
229     tassert(A.fflu(0, false).get<1>().rows() == A.rows());
230     permxx p1(5), p2(5);
231     ltupleref(rank1, res1, den1) = fflu(A, &p1);
232     rank2 = nmod_poly_mat_fflu(res2._mat(), den2._poly(), p2._data(),
233             A._mat(), false);
234     tassert(rank1 == rank2 && res1 == res2 && p1 == p2 && den1 == den2);
235     tassert(rank1 == A.fflu(0, false).get<0>());
236 
237     ltupleref(rank1, res1, den1) = rref(A);
238     rank2 = nmod_poly_mat_rref(res2._mat(), den2._poly(), A._mat());
239     tassert(rank1 == rank2 && res1 == res2 && p1 == p2 && den1 == den2);
240 }
241 
242 
243 void
test_printing()244 test_printing()
245 {
246     if(0)
247         print_pretty(nmod_poly_matxx::zero(2, 2, 7), "x"); // make sure this compiles
248 }
249 
250 int
main()251 main()
252 {
253     std::cout << "nmod_poly_matxx....";
254 
255     test_init();
256     test_arithmetic();
257     test_functions();
258     test_randomisation();
259     test_row_reduction();
260     test_printing();
261 
262     std::cout << "PASS" << std::endl;
263     return 0;
264 }
265 
266