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 <stdio.h>
13 #include <stdlib.h>
14 #include "flint.h"
15 #include "nmod_poly.h"
16 #include "nmod_poly_mat.h"
17 #include "fmpz.h"
18
19 int
main(void)20 main(void)
21 {
22 slong i;
23
24 FLINT_TEST_INIT(state);
25
26 flint_printf("inv....");
27 fflush(stdout);
28
29 /* Test aliasing */
30 for (i = 0; i < 40 * flint_test_multiplier(); i++)
31 {
32 nmod_poly_mat_t A, Ainv;
33 nmod_poly_t den1, den2;
34 slong n, deg;
35 float density;
36 int ns1, ns2, result;
37 mp_limb_t mod;
38
39 mod = n_randtest_prime(state, 0);
40 n = n_randint(state, 8);
41 deg = 1 + n_randint(state, 5);
42 density = n_randint(state, 100) * 0.01;
43
44 nmod_poly_mat_init(A, n, n, mod);
45 nmod_poly_mat_init(Ainv, n, n, mod);
46 nmod_poly_init(den1, mod);
47 nmod_poly_init(den2, mod);
48
49 nmod_poly_mat_randtest_sparse(A, state, deg, density);
50
51 ns1 = nmod_poly_mat_inv(Ainv, den1, A);
52 ns2 = nmod_poly_mat_inv(A, den2, A);
53
54 result = ns1 == ns2;
55
56 if (result && ns1 != 0)
57 {
58 result = nmod_poly_equal(den1, den2) &&
59 nmod_poly_mat_equal(A, Ainv);
60 }
61
62 if (!result)
63 {
64 flint_printf("FAIL (aliasing)!\n");
65 nmod_poly_mat_print(A, "x"); flint_printf("\n");
66 nmod_poly_mat_print(Ainv, "x"); flint_printf("\n");
67 abort();
68 }
69
70 nmod_poly_mat_clear(A);
71 nmod_poly_mat_clear(Ainv);
72 nmod_poly_clear(den1);
73 nmod_poly_clear(den2);
74 }
75
76 /* Check A^(-1) = A = 1 */
77 for (i = 0; i < 100 * flint_test_multiplier(); i++)
78 {
79 nmod_poly_mat_t A, Ainv, B, Iden;
80 nmod_poly_t den, det;
81 slong n, deg;
82 float density;
83 int nonsingular;
84 mp_limb_t mod;
85
86 mod = n_randtest_prime(state, 0);
87 n = n_randint(state, 10);
88 deg = 1 + n_randint(state, 5);
89 density = n_randint(state, 100) * 0.01;
90
91 nmod_poly_mat_init(A, n, n, mod);
92 nmod_poly_mat_init(Ainv, n, n, mod);
93 nmod_poly_mat_init(B, n, n, mod);
94 nmod_poly_mat_init(Iden, n, n, mod);
95 nmod_poly_init(den, mod);
96 nmod_poly_init(det, mod);
97
98 nmod_poly_mat_randtest_sparse(A, state, deg, density);
99 nonsingular = nmod_poly_mat_inv(Ainv, den, A);
100 nmod_poly_mat_det_interpolate(det, A);
101
102 if (n == 0)
103 {
104 if (nonsingular == 0 || !nmod_poly_is_one(den))
105 {
106 flint_printf("FAIL: expected empty matrix to pass\n");
107 abort();
108 }
109 }
110 else
111 {
112 if (!nmod_poly_equal(den, det))
113 {
114 nmod_poly_neg(det, det);
115 flint_printf("FAIL: den != det(A)\n");
116 abort();
117 }
118
119 nmod_poly_mat_mul(B, Ainv, A);
120 nmod_poly_mat_one(Iden);
121 nmod_poly_mat_scalar_mul_nmod_poly(Iden, Iden, den);
122
123 if (!nmod_poly_mat_equal(B, Iden))
124 {
125 flint_printf("FAIL:\n");
126 flint_printf("A:\n");
127 nmod_poly_mat_print(A, "x");
128 flint_printf("Ainv:\n");
129 nmod_poly_mat_print(Ainv, "x");
130 flint_printf("B:\n");
131 nmod_poly_mat_print(B, "x");
132 flint_printf("den:\n");
133 nmod_poly_print(den);
134 abort();
135 }
136 }
137
138 nmod_poly_clear(den);
139 nmod_poly_clear(det);
140 nmod_poly_mat_clear(A);
141 nmod_poly_mat_clear(Ainv);
142 nmod_poly_mat_clear(B);
143 nmod_poly_mat_clear(Iden);
144 }
145
146 FLINT_TEST_CLEANUP(state);
147
148 flint_printf("PASS\n");
149 return 0;
150 }
151