1 /*
2 Copyright (C) 2009 William Hart
3 Copyright (C) 2010 Sebastian Pancratz
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <gmp.h>
16 #include "flint.h"
17 #include "fmpz.h"
18 #include "fmpz_poly.h"
19 #include "ulong_extras.h"
20
21 /*
22 Tests whether the polynomial is suitably normalised for the
23 result of a GCD operation, that is, whether it's leading
24 coefficient is non-negative.
25 */
26 static
_t_gcd_is_canonical(const fmpz_poly_t poly)27 int _t_gcd_is_canonical(const fmpz_poly_t poly)
28 {
29 return fmpz_poly_is_zero(poly) || (fmpz_sgn(fmpz_poly_lead(poly)) > 0);
30 }
31
32 int
main(void)33 main(void)
34 {
35 int i, result, d1, d2;
36 FLINT_TEST_INIT(state);
37
38 flint_printf("gcd_heuristic....");
39 fflush(stdout);
40
41 /* Check aliasing of a and b */
42 for (i = 0; i < 100 * flint_test_multiplier(); i++)
43 {
44 fmpz_poly_t a, b, c;
45
46 fmpz_poly_init(a);
47 fmpz_poly_init(b);
48 fmpz_poly_init(c);
49 fmpz_poly_randtest(b, state, n_randint(state, 40), 80);
50 fmpz_poly_randtest(c, state, n_randint(state, 40), 80);
51
52 d1 = fmpz_poly_gcd_heuristic(a, b, c);
53 d2 = fmpz_poly_gcd_heuristic(b, b, c);
54
55 result = ((d1 == 0 && d2 == 0) || (fmpz_poly_equal(a, b)
56 && _t_gcd_is_canonical(a)));
57 if (!result)
58 {
59 flint_printf("FAIL (aliasing a and b):\n");
60 flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n\n");
61 flint_printf("b = "), fmpz_poly_print(b), flint_printf("\n\n");
62 abort();
63 }
64
65 fmpz_poly_clear(a);
66 fmpz_poly_clear(b);
67 fmpz_poly_clear(c);
68 }
69
70 /* Check aliasing of a and c */
71 for (i = 0; i < 100 * flint_test_multiplier(); i++)
72 {
73 fmpz_poly_t a, b, c;
74
75 fmpz_poly_init(a);
76 fmpz_poly_init(b);
77 fmpz_poly_init(c);
78 fmpz_poly_randtest(b, state, n_randint(state, 40), 80);
79 fmpz_poly_randtest(c, state, n_randint(state, 40), 80);
80
81 d1 = fmpz_poly_gcd_heuristic(a, b, c);
82 d2 = fmpz_poly_gcd_heuristic(c, b, c);
83
84 result = ((d1 == 0 && d2 == 0) || (fmpz_poly_equal(a, c)
85 && _t_gcd_is_canonical(a)));
86 if (!result)
87 {
88 flint_printf("FAIL (aliasing a and c):\n");
89 flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n\n");
90 flint_printf("c = "), fmpz_poly_print(c), flint_printf("\n\n");
91 abort();
92 }
93
94 fmpz_poly_clear(a);
95 fmpz_poly_clear(b);
96 fmpz_poly_clear(c);
97 }
98
99 /* Check that a divides GCD(af, ag) */
100 for (i = 0; i < 300 * flint_test_multiplier(); i++)
101 {
102 fmpz_poly_t a, d, f, g, q, r;
103
104 fmpz_poly_init(a);
105 fmpz_poly_init(d);
106 fmpz_poly_init(f);
107 fmpz_poly_init(g);
108 fmpz_poly_init(q);
109 fmpz_poly_init(r);
110 fmpz_poly_randtest_not_zero(a, state, n_randint(state, 100) + 1, 40);
111 fmpz_poly_randtest(f, state, n_randint(state, 100), 40);
112 fmpz_poly_randtest(g, state, n_randint(state, 100), 40);
113
114 fmpz_poly_mul(f, a, f);
115 fmpz_poly_mul(g, a, g);
116 d1 = fmpz_poly_gcd_heuristic(d, f, g);
117
118 if (d1)
119 {
120 fmpz_poly_divrem_divconquer(q, r, d, a);
121
122 result = fmpz_poly_is_zero(r) && _t_gcd_is_canonical(d);
123 if (!result)
124 {
125 flint_printf("FAIL (check a | gcd(af, ag)):\n");
126 flint_printf("f = "), fmpz_poly_print(f), flint_printf("\n");
127 flint_printf("g = "), fmpz_poly_print(g), flint_printf("\n");
128 flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n");
129 flint_printf("d = "), fmpz_poly_print(d), flint_printf("\n");
130 abort();
131 }
132 }
133
134 fmpz_poly_clear(a);
135 fmpz_poly_clear(d);
136 fmpz_poly_clear(f);
137 fmpz_poly_clear(g);
138 fmpz_poly_clear(q);
139 fmpz_poly_clear(r);
140 }
141
142 /* Check that a == GCD(af, ag) when GCD(f, g) = 1 */
143 for (i = 0; i < 300 * flint_test_multiplier(); i++)
144 {
145 fmpz_poly_t a, d, f, g, q, r;
146
147 fmpz_poly_init(a);
148 fmpz_poly_init(d);
149 fmpz_poly_init(f);
150 fmpz_poly_init(g);
151 fmpz_poly_init(q);
152 fmpz_poly_init(r);
153 fmpz_poly_randtest_not_zero(a, state, n_randint(state, 100) + 1, 200);
154 do {
155 fmpz_poly_randtest(f, state, n_randint(state, 100), 200);
156 fmpz_poly_randtest(g, state, n_randint(state, 100), 200);
157 fmpz_poly_gcd_heuristic(d, f, g);
158 } while (!(d->length == 1 && fmpz_is_one(d->coeffs)));
159
160 fmpz_poly_mul(f, a, f);
161 fmpz_poly_mul(g, a, g);
162 d1 = fmpz_poly_gcd_heuristic(d, f, g);
163
164 if (d1)
165 {
166 if (!_t_gcd_is_canonical(a)) fmpz_poly_neg(a, a);
167
168 result = fmpz_poly_equal(d, a) && _t_gcd_is_canonical(d);
169 if (!result)
170 {
171 flint_printf("FAIL (check a == gcd(af, ag) when gcd(f, g) = 1):\n");
172 flint_printf("f = "), fmpz_poly_print(f), flint_printf("\n");
173 flint_printf("g = "), fmpz_poly_print(g), flint_printf("\n");
174 flint_printf("a = "), fmpz_poly_print(a), flint_printf("\n");
175 flint_printf("d = "), fmpz_poly_print(d), flint_printf("\n");
176 abort();
177 }
178 }
179
180 fmpz_poly_clear(a);
181 fmpz_poly_clear(d);
182 fmpz_poly_clear(f);
183 fmpz_poly_clear(g);
184 fmpz_poly_clear(q);
185 fmpz_poly_clear(r);
186 }
187
188 /*
189 Check that gcd(f, ga) divides f and ga for small generic f, g
190 and a small linear factor a. Exercises a bug found by Anton Mellit.
191 */
192 for (i = 0; i < 1000 * flint_test_multiplier(); i++)
193 {
194 fmpz_poly_t a, d, f, g, q, r;
195
196 fmpz_poly_init(d);
197 fmpz_poly_init(f);
198 fmpz_poly_init(g);
199 fmpz_poly_init(q);
200 fmpz_poly_init(r);
201 fmpz_poly_init(a);
202 fmpz_poly_randtest(f, state, n_randint(state, 10), 8);
203 fmpz_poly_randtest(g, state, n_randint(state, 10), 4);
204
205 /* multiply by small linear factor */
206 fmpz_poly_set_coeff_si(a, 0, n_randint(state, 2) ? 1 : -1);
207 fmpz_poly_set_coeff_si(a, 1, 1);
208 fmpz_poly_mul(g, g, a);
209
210 d1 = fmpz_poly_gcd_heuristic(d, f, g);
211
212 if (d1)
213 {
214 if (fmpz_poly_is_zero(d))
215 result = fmpz_poly_is_zero(f) && fmpz_poly_is_zero(g);
216 else
217 {
218 fmpz_poly_divrem_divconquer(q, r, f, d);
219 result = fmpz_poly_is_zero(r);
220 fmpz_poly_divrem_divconquer(q, r, g, d);
221 result &= fmpz_poly_is_zero(r);
222 }
223
224 if (!result)
225 {
226 flint_printf("FAIL (gcd(f, g) | f and g):\n");
227 flint_printf("f = "), fmpz_poly_print(f), flint_printf("\n");
228 flint_printf("g = "), fmpz_poly_print(g), flint_printf("\n");
229 flint_printf("d = "), fmpz_poly_print(d), flint_printf("\n");
230 abort();
231 }
232 }
233
234 fmpz_poly_clear(a);
235 fmpz_poly_clear(d);
236 fmpz_poly_clear(f);
237 fmpz_poly_clear(g);
238 fmpz_poly_clear(q);
239 fmpz_poly_clear(r);
240 }
241
242 /* Sebastian's test case */
243 {
244 fmpz_poly_t a, b, d;
245
246 fmpz_poly_init(a);
247 fmpz_poly_init(b);
248 fmpz_poly_init(d);
249
250 fmpz_poly_set_coeff_ui(b, 2, 1);
251 fmpz_poly_set_coeff_si(a, 0, -32);
252 fmpz_poly_set_coeff_si(a, 1, 24);
253
254 fmpz_poly_gcd_heuristic(d, a, b);
255
256 result = (d->length == 1 && fmpz_is_one(d->coeffs));
257 if (!result)
258 {
259 flint_printf("FAIL (check 1 == gcd(x^2, 24*x - 32):\n");
260 fmpz_poly_print(d); flint_printf("\n");
261 abort();
262 }
263
264 fmpz_poly_clear(a);
265 fmpz_poly_clear(b);
266 fmpz_poly_clear(d);
267 }
268
269 /* Anton Mellit's test case */
270 {
271 fmpz_poly_t a, b, d;
272 int heuristic;
273
274 fmpz_poly_init(a);
275 fmpz_poly_init(b);
276 fmpz_poly_init(d);
277
278 /*
279 b = 3*q^12 - 8*q^11 - 24*q^10 - 48*q^9 - 84*q^8 - 92*q^7 - 92*q^6 -
280 70*q^5 - 50*q^4 - 27*q^3 - 13*q^2 - 4*q - 1
281 a = q^13 - 2*q^12 + 2*q^10 - q^9
282 */
283 fmpz_poly_set_str(b, "13 -1 -4 -13 -27 -50 -70 -92 -92 -84 -48 -24 -8 3");
284 fmpz_poly_set_str(a, "14 0 0 0 0 0 0 0 0 0 -1 2 0 -2 1");
285
286 heuristic = fmpz_poly_gcd_heuristic(d, a, b);
287
288 result = (heuristic == 0 || (d->length == 1 && fmpz_is_one(d->coeffs)));
289 if (!result)
290 {
291 flint_printf("FAIL Mellit test case:\n");
292 fmpz_poly_print(d); flint_printf("\n");
293 abort();
294 }
295
296 fmpz_poly_clear(a);
297 fmpz_poly_clear(b);
298 fmpz_poly_clear(d);
299 }
300
301 /* Daniel's test case */
302 {
303 fmpz_poly_t a, b, g; fmpz_poly_init(a); fmpz_poly_init(b);
304 fmpz_poly_init(g);
305 fmpz_poly_set_str(a, "40 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -7609399 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 44");
306 fmpz_poly_set_str(b, "40 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 54909036 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -59769402");
307 fmpz_poly_gcd(g, a, b);
308 fmpz_poly_clear(a);
309 fmpz_poly_clear(b);
310 fmpz_poly_clear(g);
311 }
312
313 FLINT_TEST_CLEANUP(state);
314
315 flint_printf("PASS\n");
316 return 0;
317 }
318