1 /*++
2 Copyright (c) 2011 Microsoft Corporation
3
4 Module Name:
5
6 polynomial.cpp
7
8 Abstract:
9
10 Goodies for creating and handling polynomials.
11
12 Author:
13
14 Leonardo (leonardo) 2011-11-15
15
16 Notes:
17
18 --*/
19 #if !defined(__clang__)
20 #include "math/polynomial/polynomial.h"
21 #include "math/polynomial/polynomial_var2value.h"
22 #include "math/polynomial/polynomial_cache.h"
23 #include "math/polynomial/linear_eq_solver.h"
24 #include "util/rlimit.h"
25
tst1()26 static void tst1() {
27 std::cout << "\n----- Basic testing -------\n";
28 reslimit rl;
29 polynomial::numeral_manager nm;
30 polynomial::manager m(rl, nm);
31 polynomial_ref x0(m);
32 polynomial_ref x1(m);
33 polynomial_ref x2(m);
34 x0 = m.mk_polynomial(m.mk_var());
35 x1 = m.mk_polynomial(m.mk_var());
36 x2 = m.mk_polynomial(m.mk_var());
37 std::cout << x0 << " " << x1 << " " << x2 << "\n";
38 polynomial_ref p(m);
39 p = (x0^3) + x1*x0 + 2;
40 std::cout << p << "\n";
41 std::cout << "max_var(p): " << max_var(p) << "\n";
42 ENSURE(max_var(p) == 1);
43 std::cout << (2*x2 - x1*x0) << "\n";
44 std::cout << (p + (2*x2 - x1*x0)) << "\n";
45 std::cout << (p*p + 2*x2) << "\n";
46 std::cout << derivative(p*p + 2*x2, 0) << "\n";
47 polynomial_ref q(m);
48 q = (x0^4) + x0 + 1;
49 std::cout << "q(x): " << q << "\n";
50 std::cout << "q(y): " << compose_y(q, 2) << "\n";
51 std::cout << "q(x-y): " << compose_x_minus_y(q, 2) << "\n";
52 q = (x0 - 1)*(x0 - 2)*(x0 - 1)*(x0 + 2);
53 std::cout << "q: " << q << "\n";
54 polynomial_ref s(m);
55 s = (x0 - 1)*((x0 + 3)^2);
56 std::cout << "s: " << s << "\n";
57 }
58
tst_pseudo_div(polynomial_ref const & A,polynomial_ref const & B,polynomial::var x)59 static void tst_pseudo_div(polynomial_ref const & A, polynomial_ref const & B, polynomial::var x) {
60 reslimit rl;
61 polynomial::manager & m = A.m();
62 std::cout << "---- Pseudo-division test ----\n";
63 std::cout << "A: " << A << "\n";
64 std::cout << "B: " << B << "\n";
65 std::cout << "x: " << x << "\n";
66 polynomial_ref Q(m);
67 polynomial_ref R(m);
68 unsigned d;
69 Q = pseudo_division(A, B, x, d, R);
70 std::cout << "d: " << d << "\n";
71 std::cout << "Q: " << Q << "\n";
72 std::cout << "R: " << R << "\n";
73 polynomial_ref l_B(m);
74 l_B = coeff(B, x, degree(B, x));
75 std::cout << "l_B: " << l_B << "\n";
76 polynomial_ref l_B_d(m);
77 l_B_d = l_B^d;
78 polynomial_ref t(m);
79 std::cout << "l_B^d: " << l_B_d << "\n";
80 std::cout << "Q * B + R: " << Q * B + R << "\n";
81 std::cout << "l_B_d * A: " << l_B_d * A << "\n";
82 ENSURE(eq((Q * B + R), (l_B_d * A)));
83 }
84
tst2()85 static void tst2() {
86 reslimit rl;
87 polynomial::numeral_manager nm;
88 polynomial::manager m(rl, nm);
89 polynomial_ref x0(m);
90 polynomial_ref x1(m);
91 polynomial_ref x2(m);
92 x0 = m.mk_polynomial(m.mk_var());
93 x1 = m.mk_polynomial(m.mk_var());
94 x2 = m.mk_polynomial(m.mk_var());
95 polynomial_ref p(m);
96 polynomial_ref q(m);
97 p = ((x0 - 1)^2)*x2 + (x1^2)*((x2 - 2)^2) + 1;
98 q = (x0 - 1)*x2 + (x1^3)*(x2 - 2) + (x0 - 2)*(x1 - 2) + 10;
99 tst_pseudo_div(p, q, 0);
100 }
101
102
tst3()103 static void tst3() {
104 reslimit rl;
105 polynomial::numeral_manager nm;
106 polynomial::manager m(rl, nm);
107 polynomial_ref x0(m);
108 polynomial_ref x1(m);
109 x0 = m.mk_polynomial(m.mk_var());
110 x1 = m.mk_polynomial(m.mk_var());
111 polynomial_ref p(m);
112 polynomial_ref q(m);
113 p = (x1^2) + (x0^2) - 1;
114 q = (x1*x0) - 1;
115 tst_pseudo_div(p, q, 1);
116 }
117
tst4()118 static void tst4() {
119 std::cout << "---- Testing renaming/reordering ----\n";
120 reslimit rl;
121 polynomial::numeral_manager nm;
122 polynomial::manager m(rl, nm);
123 polynomial_ref x0(m);
124 polynomial_ref x1(m);
125 polynomial_ref x2(m);
126 x0 = m.mk_polynomial(m.mk_var());
127 x1 = m.mk_polynomial(m.mk_var());
128 x2 = m.mk_polynomial(m.mk_var());
129 polynomial_ref p1 = x2 + ((x0 - 1)^2)*x1 + (x2^3) + 10;
130 polynomial_ref p2 = x0*x1*x2 + x1*(x2^3) + ((x0 - 2)^2);
131 std::cout << "p1: " << p1 << "\n";
132 std::cout << "p2: " << p2 << "\n";
133 polynomial::var new_order[3] = { 2, 0, 1 };
134 m.rename(3, new_order);
135 std::cout << "----- x0 -> x2, x1 -> x0, x2 -> x1 \n";
136 std::cout << "p1: " << p1 << "\n";
137 std::cout << "p2: " << p2 << "\n";
138 }
139
tst_quasi_resultant(polynomial_ref const & p,polynomial_ref const & q,polynomial::var x)140 static void tst_quasi_resultant(polynomial_ref const & p, polynomial_ref const & q, polynomial::var x) {
141 std::cout << "---- Testing quasi-resultants ---- \n";
142 std::cout << "p : " << p << "\n";
143 std::cout << "q : " << q << "\n";
144 std::cout << "x : " << x << "\n--->\n";
145 std::cout << quasi_resultant(p, q, x) << "\n";
146 }
147
tst5()148 static void tst5() {
149 reslimit rl;
150 polynomial::numeral_manager nm;
151 polynomial::manager m(rl, nm);
152 polynomial_ref x0(m);
153 polynomial_ref x1(m);
154 polynomial_ref x2(m);
155 x0 = m.mk_polynomial(m.mk_var());
156 x1 = m.mk_polynomial(m.mk_var());
157 x2 = m.mk_polynomial(m.mk_var());
158 polynomial_ref p(m);
159 polynomial_ref q(m);
160 p = ((x0 - x1)^2) - 2;
161 q = (x1^2) - 3;
162 // sqrt(2) + sqrt(3) must be a root of the quasi-resultant
163 tst_quasi_resultant(p, q, 1);
164 }
165
tst6()166 static void tst6() {
167 reslimit rl;
168 polynomial::numeral_manager nm;
169 polynomial::manager m(rl, nm);
170 polynomial_ref x0(m);
171 polynomial_ref x1(m);
172 polynomial_ref x2(m);
173 x0 = m.mk_polynomial(m.mk_var());
174 x1 = m.mk_polynomial(m.mk_var());
175 polynomial_ref p(m);
176 polynomial_ref q(m);
177 p = (x0 - 2)*(x0 - 3)*(x0 + 2);
178 std::cout << "p(x0): " << p << "\n";
179 std::cout << "p(-x0): " << compose_minus_x(p) << "\n";
180 std::cout << "x^3*p(1/x0): " << compose_1_div_x(p) << "\n";
181 std::cout << "p(x0 - x1): " << compose_x_minus_y(p, 1) << "\n";
182 std::cout << "x1^3*p(x0/x1): " << compose_x_div_y(p, 1) << "\n";
183 }
184
tst7()185 static void tst7() {
186 reslimit rl;
187 polynomial::numeral_manager nm;
188 polynomial::manager m(rl, nm);
189 polynomial_ref x0(m);
190 polynomial_ref x1(m);
191 polynomial_ref x2(m);
192 x0 = m.mk_polynomial(m.mk_var());
193 x1 = m.mk_polynomial(m.mk_var());
194 x2 = m.mk_polynomial(m.mk_var());
195 polynomial_ref p(m);
196 polynomial_ref q1(m);
197 polynomial_ref q2(m);
198 p = (x0 - x1)*(x2 - 1);
199 q1 = (x0^2) - 2;
200 q2 = (x1^2) - 2;
201 polynomial_ref r(m);
202 r = quasi_resultant(p, q1, 0);
203 std::cout << "1) r: " << r << "\n";
204 r = quasi_resultant(r, q2, 1);
205 std::cout << "2) r: " << r << "\n";
206 }
207
tst8()208 static void tst8() {
209 reslimit rl;
210 polynomial::numeral_manager nm;
211 polynomial::manager m(rl, nm);
212 polynomial_ref x0(m);
213 polynomial_ref x1(m);
214 polynomial_ref x2(m);
215 x0 = m.mk_polynomial(m.mk_var());
216 x1 = m.mk_polynomial(m.mk_var());
217 x2 = m.mk_polynomial(m.mk_var());
218 polynomial_ref p(m);
219 polynomial_ref sqrt2(m);
220 polynomial_ref sqrt3(m);
221 p = x2 - (x0*x1 - (x0^2) + 1);
222 sqrt3 = (x0^2) - 3;
223 sqrt2 = (x1^2) - 2;
224 polynomial_ref r(m);
225 r = quasi_resultant(p, sqrt2, 1);
226 r = quasi_resultant(r, sqrt3, 0);
227 std::cout << "r: " << r << "\n";
228 }
229
230
tst9()231 static void tst9() {
232 reslimit rl;
233 polynomial::numeral_manager nm;
234 polynomial::manager m(rl, nm);
235 polynomial_ref x0(m);
236 polynomial_ref x1(m);
237 polynomial_ref x2(m);
238 x0 = m.mk_polynomial(m.mk_var());
239 x1 = m.mk_polynomial(m.mk_var());
240 x2 = m.mk_polynomial(m.mk_var());
241 polynomial_ref p(m);
242 polynomial_ref sqrt2(m);
243 p = ((x0^3) - 1)*(x1^2) - 1;
244 sqrt2 = ((x0^2) - 2)*(x0 - 2); // added garbage to polynomial
245 // sqrt2 = (x0^2) - 2; // added garbage to polynomial
246 // p(sqrt(2), x1) has the roots -0.7395 and 0.7395
247 polynomial_ref r(m);
248 r = quasi_resultant(p, sqrt2, 0);
249 std::cout << "p: " << p << "\n";
250 std::cout << "sqrt2: " << sqrt2 << "\n";
251 std::cout << "r: " << r << "\n";
252 // r contains the roots -0.7395 and 0.7395, plus garbage roots: 0, -0.3779, 0.3779
253 polynomial_ref q(m);
254 q = x2 - (((x0^3) - 1)*(x1^2) - 1);
255 std::cout << "q: " << q << "\n";
256 polynomial_ref r2(m);
257 TRACE("polynomial", tout << "QUASI_RESULTANT: q, sqrt2.....\n";);
258 r2 = quasi_resultant(q, sqrt2, 0);
259 // TRACE("polynomial", tout << "QUASI_RESULTANT: sqrt2, q.....\n";);
260 // std::cout << "r2: " << r2 << "\n";
261 // r2 = quasi_resultant(sqrt2, q, 0);
262 // std::cout << "r2: " << r2 << "\n";
263 // return;
264 std::cout << "r2: " << r2 << "\n";
265 r2 = normalize(quasi_resultant(r2, r, 1));
266 std::cout << "r2: " << r2 << "\n";
267 polynomial_ref_vector seq(m);
268 r2 = normalize(quasi_resultant(sqrt2, q, 0));
269 // sturm_seq(r2, seq);
270 std::cout << "r2:\n" << r2 << "\n";
271 }
272
tst10()273 static void tst10() {
274 reslimit rl;
275 polynomial::numeral_manager nm;
276 polynomial::manager m(rl, nm);
277 polynomial_ref x0(m);
278 polynomial_ref x1(m);
279 polynomial_ref x2(m);
280 x0 = m.mk_polynomial(m.mk_var());
281 x1 = m.mk_polynomial(m.mk_var());
282 x2 = m.mk_polynomial(m.mk_var());
283 polynomial_ref A(m);
284 polynomial_ref B(m);
285 polynomial_ref g(m);
286 polynomial_ref h(m);
287 polynomial_ref R(m);
288 A = x2 - (((x0^3) - 1)*(x1^2) - 1);
289 B = ((x0^2) - 2)*(x0 - 2);
290 std::cout << "A: " << A << "\nB: " << B << "\n";
291 unsigned d;
292 R = pseudo_remainder(A, B, 0, d);
293 std::cout << "R: " << R << "\n";
294 // second iteration
295 std::cout << "second iteration...\n";
296 A = B;
297 B = R;
298 g = coeff(A, 0, degree(A, 0));
299 std::cout << "A: " << A << "\nB: " << B << "\ng: " << g << "\n";
300 R = pseudo_remainder(A, B, 0, d);
301 std::cout << "R: " << R << "\n";
302 // third iteration
303 std::cout << "third iteration...\n";
304 A = B;
305 B = R;
306 g = coeff(A, 0, degree(A, 0));
307 std::cout << "A: " << A << "\nB: " << B << "\ng: " << g << "\n";
308 R = pseudo_remainder(A, B, 0, d);
309 std::cout << "R: " << R << "\n";
310
311 }
312
tst11()313 static void tst11() {
314 reslimit rl;
315 polynomial::numeral_manager nm;
316 polynomial::manager m(rl, nm);
317 polynomial_ref x0(m);
318 polynomial_ref x1(m);
319 polynomial_ref x2(m);
320 polynomial_ref x3(m);
321 x0 = m.mk_polynomial(m.mk_var());
322 x1 = m.mk_polynomial(m.mk_var());
323 x2 = m.mk_polynomial(m.mk_var());
324 x3 = m.mk_polynomial(m.mk_var());
325 polynomial_ref p(m);
326 polynomial_ref q(m);
327 q = ((x1^2) + 1)*(x2 + 1);
328 p = (x3 + 1)*q;
329 polynomial_ref d(m);
330 d = exact_div(p, q);
331 std::cout << "p: " << p << "\nq: " << q << "\nd: " << d << "\n";
332 ENSURE(eq(q * d, p));
333
334 q = ((x1^3) + x1 + 1)*((x2^2) + x2 + x2 + 1)*((x3^2) + 2);
335 p = (x1 + (x3^2) + x3 + x2 + (x2^2) + 1)*((x1^3) + x1 + 1)*((x2^2) + x2 + x2 + 1)*((x3^2) + 2);
336 d = exact_div(p, q);
337 std::cout << "p: " << p << "\nq: " << q << "\nd: " << d << "\n";
338 ENSURE(eq(q * d, p));
339 }
340
tst_discriminant(polynomial_ref const & p,polynomial::var x,polynomial_ref const & expected)341 static void tst_discriminant(polynomial_ref const & p, polynomial::var x, polynomial_ref const & expected) {
342 polynomial::manager & m = p.m();
343 polynomial_ref r(m);
344 r = discriminant(p, x);
345 std::cout << "r: " << r << "\n";
346 std::cout << "expected: " << expected << "\n";
347 ENSURE(eq(r, expected));
348 m.lex_sort(r);
349 std::cout << "r (sorted): " << r << "\n";
350 }
351
tst_discriminant(polynomial_ref const & p,polynomial_ref const & expected)352 static void tst_discriminant(polynomial_ref const & p, polynomial_ref const & expected) {
353 tst_discriminant(p, max_var(p), expected);
354 }
355
tst_discriminant()356 static void tst_discriminant() {
357 reslimit rl;
358 polynomial::numeral_manager nm;
359 polynomial::manager m(rl, nm);
360 polynomial_ref a(m);
361 polynomial_ref b(m);
362 polynomial_ref c(m);
363 polynomial_ref d(m);
364 polynomial_ref x(m);
365 a = m.mk_polynomial(m.mk_var());
366 b = m.mk_polynomial(m.mk_var());
367 c = m.mk_polynomial(m.mk_var());
368 d = m.mk_polynomial(m.mk_var());
369 x = m.mk_polynomial(m.mk_var());
370 tst_discriminant(a*(x^2) + b*x + c,
371 (b^2) - 4*a*c);
372 tst_discriminant(a*(x^3) + b*(x^2) + c*x + d,
373 (b^2)*(c^2) - 4*a*(c^3) - 4*(b^3)*d + 18*a*b*c*d - 27*(a^2)*(d^2));
374 tst_discriminant(a*(x^3) + b*(x^2) + c*(x^2) + d,
375 -4*(b^3)*d - 12*(b^2)*c*d - 12*b*(c^2)*d - 4*(c^3)*d - 27*(a^2)*(d^2));
376 tst_discriminant(a*(x^3) + b*(x^2) + c*(x^2) + d,
377 -4*(b^3)*d - 12*(b^2)*c*d - 12*b*(c^2)*d - 4*(c^3)*d - 27*(a^2)*(d^2));
378 tst_discriminant(a*(x^3) + (b^2)*d*(x^2) + c*(x^2) + d,
379 -4*(b^6)*(d^4) - 12*(b^4)*c*(d^3) - 12*(b^2)*(c^2)*(d^2) - 4*(c^3)*d - 27*(a^2)*(d^2));
380 tst_discriminant(a*(x^4) + b*(x^2) + c,
381 16*a*(b^4)*c - 128*(a^2)*(b^2)*(c^2) + 256*(a^3)*(c^3));
382 polynomial_ref one(m);
383 one = m.mk_const(rational(1));
384 tst_discriminant(x, one);
385 tst_discriminant(3*x, one);
386 polynomial_ref zero(m);
387 zero = m.mk_zero();
388 tst_discriminant(one, zero);
389 tst_discriminant(a*(x^7) + b,
390 -823543*(a^6)*(b^6));
391 tst_discriminant(((a^2)+(b^2)+c)*(x^4) + (d + a*b)*x + a,
392 -27*(a^8)*(b^4) - 54*(a^6)*(b^6) - 27*(a^4)*(b^8) - 54*(a^6)*(b^4)*c - 54*(a^4)*(b^6)*c -
393 108*(a^7)*(b^3)*d - 216*(a^5)*(b^5)*d - 108*(a^3)*(b^7)*d - 27*(a^4)*(b^4)*(c^2) -
394 216*(a^5)*(b^3)*c*d - 216*(a^3)*(b^5)*c*d - 162*(a^6)*(b^2)*(d^2) - 324*(a^4)*(b^4)*(d^2) -
395 162*(a^2)*(b^6)*(d^2) + 256*(a^9) + 768*(a^7)*(b^2) + 768*(a^5)*(b^4) + 256*(a^3)*(b^6) -
396 108*(a^3)*(b^3)*(c^2)*d - 324*(a^4)*(b^2)*c*(d^2) - 324*(a^2)*(b^4)*c*(d^2) -
397 108*(a^5)*b*(d^3) - 216*(a^3)*(b^3)*(d^3) - 108*a*(b^5)*(d^3) + 768*(a^7)*c +
398 1536*(a^5)*(b^2)*c + 768*(a^3)*(b^4)*c - 162*(a^2)*(b^2)*(c^2)*(d^2) - 216*(a^3)*b*c*(d^3) -
399 216*a*(b^3)*c*(d^3) - 27*(a^4)*(d^4) - 54*(a^2)*(b^2)*(d^4) - 27*(b^4)*(d^4) + 768*(a^5)*(c^2)
400 + 768*(a^3)*(b^2)*(c^2) - 108*a*b*(c^2)*(d^3) - 54*(a^2)*c*(d^4) - 54*(b^2)*c*(d^4) +
401 256*(a^3)*(c^3) - 27*(c^2)*(d^4));
402 tst_discriminant((x^5) + a*(x^2) + a,
403 108*(a^6) + 3125*(a^4));
404 tst_discriminant((x^5) + (a*b)*(x^2) + a,
405 108*(a^6)*(b^5) + 3125*(a^4));
406 tst_discriminant((x^5) + (a*b*c)*(x^2) + a,
407 108*(a^6)*(b^5)*(c^5) + 3125*(a^4));
408 tst_discriminant((x^5) + (a*b*c + d)*(x^2) + a,
409 108*(a^6)*(b^5)*(c^5) + 540*(a^5)*(b^4)*(c^4)*d + 1080*(a^4)*(b^3)*(c^3)*(d^2) +
410 1080*(a^3)*(b^2)*(c^2)*(d^3) + 540*(a^2)*b*c*(d^4) + 108*a*(d^5) + 3125*(a^4));
411 tst_discriminant((x^4) + a*(x^2) + (a + c)*x + (c^2),
412 16*(a^4)*(c^2) - 128*(a^2)*(c^4) + 256*(c^6) - 4*(a^5) - 8*(a^4)*c + 140*(a^3)*(c^2) +
413 288*(a^2)*(c^3) + 144*a*(c^4) - 27*(a^4) - 108*(a^3)*c - 162*(a^2)*(c^2) - 108*a*(c^3) -
414 27*(c^4));
415 tst_discriminant((x^4) + (a + b)*(x^2) + (a + c)*x,
416 -4*(a^5) - 12*(a^4)*b - 12*(a^3)*(b^2) - 4*(a^2)*(b^3) - 8*(a^4)*c - 24*(a^3)*b*c -
417 24*(a^2)*(b^2)*c - 8*a*(b^3)*c - 4*(a^3)*(c^2) - 12*(a^2)*b*(c^2) - 12*a*(b^2)*(c^2) -
418 4*(b^3)*(c^2) - 27*(a^4) - 108*(a^3)*c - 162*(a^2)*(c^2) - 108*a*(c^3) - 27*(c^4));
419 tst_discriminant((x^4) + (a + c)*x + (c^2),
420 256*(c^6) - 27*(a^4) - 108*(a^3)*c - 162*(a^2)*(c^2) - 108*a*(c^3) - 27*(c^4)
421 );
422 tst_discriminant((x^4) + (a + b)*(x^2) + (a + c)*x + (c^2),
423 16*(a^4)*(c^2) + 64*(a^3)*b*(c^2) + 96*(a^2)*(b^2)*(c^2) + 64*a*(b^3)*(c^2) + 16*(b^4)*(c^2) -
424 128*(a^2)*(c^4) - 256*a*b*(c^4) - 128*(b^2)*(c^4) + 256*(c^6) - 4*(a^5) - 12*(a^4)*b -
425 12*(a^3)*(b^2) - 4*(a^2)*(b^3) - 8*(a^4)*c - 24*(a^3)*b*c - 24*(a^2)*(b^2)*c - 8*a*(b^3)*c
426 + 140*(a^3)*(c^2) + 132*(a^2)*b*(c^2) - 12*a*(b^2)*(c^2) - 4*(b^3)*(c^2) + 288*(a^2)*(c^3) +
427 288*a*b*(c^3) + 144*a*(c^4) + 144*b*(c^4) - 27*(a^4) - 108*(a^3)*c - 162*(a^2)*(c^2) -
428 108*a*(c^3) - 27*(c^4));
429 tst_discriminant((a + c)*(x^3) + (a + b)*(x^2) + (a + c)*x + (c^2),
430 -27*(a^2)*(c^4) - 54*a*(c^5) - 27*(c^6) + 14*(a^3)*(c^2) + 6*(a^2)*b*(c^2) -
431 12*a*(b^2)*(c^2) - 4*(b^3)*(c^2) + 36*(a^2)*(c^3) + 36*a*b*(c^3) + 18*a*(c^4) + 18*b*(c^4)
432 - 3*(a^4) + 2*(a^3)*b + (a^2)*(b^2) - 14*(a^3)*c + 4*(a^2)*b*c + 2*a*(b^2)*c -
433 23*(a^2)*(c^2) + 2*a*b*(c^2) + (b^2)*(c^2) - 16*a*(c^3) - 4*(c^4));
434 tst_discriminant((a^4) - 2*(a^3) + (a^2) - 3*(b^2)*a + 2*(b^4),
435 max_var(a),
436 2048*(b^12) - 4608*(b^10) + 37*(b^8) + 12*(b^6));
437 tst_discriminant((a^4) - 2*(a^3) + (a^2) - 3*(b^2)*a + 2*(b^4),
438 max_var(b),
439 2048*(a^12) - 12288*(a^11) + 26112*(a^10) - 22528*(a^9) + 5664*(a^8) + 960*(a^7) +
440 32*(a^6));
441 tst_discriminant((x^4) + a*(x^2) + b*x + c,
442 -4*(a^3)*(b^2) + 16*(a^4)*c - 27*(b^4) + 144*a*(b^2)*c - 128*(a^2)*(c^2) + 256*(c^3));
443 tst_discriminant((((a-1)^2) + a*b + ((b-1)^2) - 1)*(x^3) + (a*b)*(x^2) + ((a^2) - (b^2))*x + c*a,
444 -4*(a^8) - 4*(a^7)*b + 9*(a^6)*(b^2) + 12*(a^5)*(b^3) - 2*(a^4)*(b^4) - 12*(a^3)*(b^5) -
445 7*(a^2)*(b^6) + 4*a*(b^7) + 4*(b^8) + 18*(a^6)*b*c + 18*(a^5)*(b^2)*c - 4*(a^4)*(b^3)*c -
446 18*(a^3)*(b^4)*c - 18*(a^2)*(b^5)*c - 27*(a^6)*(c^2) - 54*(a^5)*b*(c^2) - 81*(a^4)*(b^2)*(c^2)
447 - 54*(a^3)*(b^3)*(c^2) - 27*(a^2)*(b^4)*(c^2) + 8*(a^7) + 8*(a^6)*b - 24*(a^5)*(b^2) -
448 24*(a^4)*(b^3) + 24*(a^3)*(b^4) + 24*(a^2)*(b^5) - 8*a*(b^6) - 8*(b^7) - 36*(a^5)*b*c -
449 36*(a^4)*(b^2)*c + 36*(a^3)*(b^3)*c + 36*(a^2)*(b^4)*c + 108*(a^5)*(c^2) + 216*(a^4)*b*(c^2)
450 + 216*(a^3)*(b^2)*(c^2) + 108*(a^2)*(b^3)*(c^2) - 4*(a^6) + 12*(a^4)*(b^2) - 12*(a^2)*(b^4) +
451 4*(b^6) + 18*(a^4)*b*c - 18*(a^2)*(b^3)*c - 162*(a^4)*(c^2) - 270*(a^3)*b*(c^2) -
452 162*(a^2)*(b^2)*(c^2) + 108*(a^3)*(c^2) + 108*(a^2)*b*(c^2) - 27*(a^2)*(c^2));
453 }
454
tst_resultant(polynomial_ref const & p,polynomial_ref const & q,polynomial::var x,polynomial_ref const & expected)455 static void tst_resultant(polynomial_ref const & p, polynomial_ref const & q, polynomial::var x, polynomial_ref const & expected) {
456 polynomial::manager & m = p.m();
457 polynomial_ref r(m);
458 std::cout << "----------------\n";
459 std::cout << "p: " << p << "\n";
460 std::cout << "q: " << q << std::endl;
461 r = resultant(p, q, x);
462 std::cout << "r: " << r << "\n";
463 std::cout << "expected: " << expected << "\n";
464 if (degree(p, x) > 0 && degree(q, x) > 0)
465 std::cout << "quasi-resultant: " << quasi_resultant(p, q, x) << "\n";
466 ENSURE(eq(r, expected));
467 m.lex_sort(r);
468 std::cout << "r (sorted): " << r << "\n";
469 }
470
tst_resultant(polynomial_ref const & p,polynomial_ref const & q,polynomial_ref const & expected)471 static void tst_resultant(polynomial_ref const & p, polynomial_ref const & q, polynomial_ref const & expected) {
472 tst_resultant(p, q, max_var(p), expected);
473 }
474
tst_resultant()475 static void tst_resultant() {
476 reslimit rl;
477 polynomial::numeral_manager nm;
478 polynomial::manager m(rl, nm);
479 polynomial_ref a(m);
480 polynomial_ref b(m);
481 polynomial_ref c(m);
482 polynomial_ref d(m);
483 polynomial_ref x(m);
484 a = m.mk_polynomial(m.mk_var());
485 b = m.mk_polynomial(m.mk_var());
486 c = m.mk_polynomial(m.mk_var());
487 d = m.mk_polynomial(m.mk_var());
488 x = m.mk_polynomial(m.mk_var());
489
490 tst_resultant((((a-1)^2) + a*b + ((b-1)^2) - 1)*(x^3) + (a*b)*(x^2) + ((a^2) - (b^2))*x + c*a,
491 a*b*(x^2) - (a^2) - (b^2),
492 -4*(a^9)*b - (a^10) - 9*(a^8)*(b^2) - 11*(a^7)*(b^3) - 14*(a^6)*(b^4) - 10*(a^5)*(b^5) -
493 10*(a^4)*(b^6) - 3*(a^3)*(b^7) - 5*(a^2)*(b^8) - (b^10) + 2*(a^6)*(b^3)*c + 2*(a^4)*(b^5)*c +
494 (a^5)*(b^3)*(c^2) + 4*(a^9) + 12*(a^8)*b + 24*(a^7)*(b^2) + 32*(a^6)*(b^3) + 40*(a^5)*(b^4) +
495 32*(a^4)*(b^5) + 24*(a^3)*(b^6) + 16*(a^2)*(b^7) + 4*a*(b^8) + 4*(b^9) - 6*(a^8) -
496 12*(a^7)*b - 24*(a^6)*(b^2) - 32*(a^5)*(b^3) - 36*(a^4)*(b^4) - 28*(a^3)*(b^5) -
497 24*(a^2)*(b^6) - 8*a*(b^7) - 6*(b^8) + 4*(a^7) + 4*(a^6)*b + 12*(a^5)*(b^2) + 12*(a^4)*(b^3)
498 + 12*(a^3)*(b^4) + 12*(a^2)*(b^5) + 4*a*(b^6) + 4*(b^7) - (a^6) - 3*(a^4)*(b^2) -
499 3*(a^2)*(b^4) - (b^6));
500
501
502 tst_resultant(a*(x^5) + b,
503 c*x + d,
504 a*(d^5) - b*(c^5));
505 tst_resultant(a*(x^5) + 3*(c + d)*(x^2) + 2*b,
506 c*x + d,
507 -2*b*(c^5) - 3*(c^4)*(d^2) - 3*(c^3)*(d^3) + a*(d^5));
508 tst_resultant(c*x + d,
509 a*(x^5) + 3*(c + d)*(x^2) + 2*b,
510 2*b*(c^5) + 3*(c^4)*(d^2) + 3*(c^3)*(d^3) - a*(d^5));
511 tst_resultant((x^2) - (a^3)*(x^2) + b + 1,
512 -49*(x^10) + 21*(x^8) + 5*(x^6) - (x^4),
513 (a^18)*(b^4) + 4*(a^18)*(b^3) + 6*(a^18)*(b^2) - 10*(a^15)*(b^5) + 4*(a^18)*b -
514 56*(a^15)*(b^4) + (a^18) - 124*(a^15)*(b^3) - 17*(a^12)*(b^6) - 136*(a^15)*(b^2) -
515 52*(a^12)*(b^5) - 74*(a^15)*b + 10*(a^12)*(b^4) + 308*(a^9)*(b^7) - 16*(a^15) +
516 220*(a^12)*(b^3) + 2224*(a^9)*(b^6) + 335*(a^12)*(b^2) + 6776*(a^9)*(b^5) - 49*(a^6)*(b^8) +
517 208*(a^12)*b + 11280*(a^9)*(b^4) - 1316*(a^6)*(b^7) + 48*(a^12) + 11060*(a^9)*(b^3) -
518 7942*(a^6)*(b^6) - 2058*(a^3)*(b^9) + 6368*(a^9)*(b^2) - 22660*(a^6)*(b^5) -
519 18424*(a^3)*(b^8) + 1984*(a^9)*b - 36785*(a^6)*(b^4) - 72380*(a^3)*(b^7) + 2401*(b^10) +
520 256*(a^9) - 36064*(a^6)*(b^3) - 163592*(a^3)*(b^6) + 26068*(b^9) - 21216*(a^6)*(b^2) -
521 234058*(a^3)*(b^5) + 126518*(b^8) - 6912*(a^6)*b - 219344*(a^3)*(b^4) + 361508*(b^7) -
522 960*(a^6) - 134208*(a^3)*(b^3) + 673537*(b^6) - 51456*(a^3)*(b^2) + 855056*(b^5) -
523 11136*(a^3)*b + 749104*(b^4) - 1024*(a^3) + 447232*(b^3) + 174144*(b^2) + 39936*b + 4096);
524 tst_resultant(((a - x)^2) + 2,
525 (x^5) - x - 1,
526 (a^10) + 10*(a^8) + 38*(a^6) - 2*(a^5) + 100*(a^4) + 40*(a^3) + 121*(a^2) - 38*a + 19);
527 tst_resultant(c - (((a^3) - 1)*(b^2) - 1),
528 ((a^2) - 2)*(a - 2),
529 max_var(a),
530 -49*(b^6) + 21*(b^4)*c + 21*(b^4) + 5*(b^2)*(c^2) + 10*(b^2)*c - (c^3) + 5*(b^2) - 3*(c^2) - 3*c - 1);
531 tst_resultant(-49*(b^6) + 21*(b^4)*c + 21*(b^4) + 5*(b^2)*(c^2) + 10*(b^2)*c - (c^3) + 5*(b^2) - 3*(c^2) - 3*c - 1,
532 (7*(b^4) - 2*(b^2) - 1),
533 max_var(b),
534 117649*(c^12) + 1075648*(c^11) + 1651888*(c^10) - 12293120*(c^9) - 46560192*(c^8)
535 - 9834496*(c^7) + 186855424*(c^6) + 314703872*(c^5) + 157351936*(c^4));
536 tst_resultant(144*(b^2) + 96*(a^2)*b + 9*(a^4) + 105*(a^2) + 70*a - 98,
537 a*(b^2) + 6*a*b + (a^3) + 9*a,
538 max_var(b),
539 81*(a^10) + 3330*(a^8) + 1260*(a^7) - 37395*(a^6) - 45780*(a^5) - 32096*(a^4) +
540 167720*(a^3) + 1435204*(a^2));
541 tst_resultant(144*(b^2) + 96*(a^2)*b + 9*(a^4) + 105*(a^2) + 70*a - 98,
542 a*(b^2) + 6*a*b + (a^3) + 9*a,
543 max_var(a),
544 11664*(b^10) + 31104*(b^9) - 119394*(b^8) - 1550448*(b^7) - 2167524*(b^6) +
545 7622712*(b^5) + 46082070*(b^4) + 46959720*(b^3) - 9774152*(b^2) - 35007168*b -
546 13984208);
547 polynomial_ref n1(m);
548 polynomial_ref n2(m);
549 polynomial_ref one(m);
550 n1 = m.mk_const(rational(10));
551 n2 = m.mk_const(rational(100));
552 one = m.mk_const(rational(1));
553 tst_resultant(n1, (x^2) + 2*x + 1, max_var(x), n2);
554 tst_resultant(n1, 2*x + 1, max_var(x), n1);
555 tst_resultant(n1, n2, 0, one);
556 tst_resultant((x^2) + 2*x + 1, n1, max_var(x), n2);
557 tst_resultant(2*x + 1, n1, max_var(x), n1);
558 tst_resultant((x^2) + 8*x + 1, n1, max_var(x), n2);
559 }
560
tst_compose()561 static void tst_compose() {
562 reslimit rl;
563 polynomial::numeral_manager nm;
564 polynomial::manager m(rl, nm);
565 polynomial_ref x0(m);
566 polynomial_ref x1(m);
567 x0 = m.mk_polynomial(m.mk_var());
568 x1 = m.mk_polynomial(m.mk_var());
569 polynomial_ref p(m);
570 p = (x0^3) - x0 + 3;
571 std::cout << "p: " << p << "\np(x - y): " << compose_x_minus_y(p, 1)
572 << "\np(x + y): " << compose_x_plus_y(p, 1) << "\np(x - x): " << compose_x_minus_y(p, 0) << "\np(x + x): " << compose_x_plus_y(p, 0) << "\n";
573 ENSURE(eq(compose_x_minus_y(p, 1), (x0^3) - 3*(x0^2)*x1 + 3*x0*(x1^2) - (x1^3) - x0 + x1 + 3));
574 ENSURE(eq(compose_x_plus_y(p, 1), (x0^3) + 3*(x0^2)*x1 + 3*x0*(x1^2) + (x1^3) - x0 - x1 + 3));
575 }
576
tst_prem()577 void tst_prem() {
578 reslimit rl;
579 polynomial::numeral_manager nm;
580 polynomial::manager m(rl, nm);
581 polynomial_ref x(m);
582 polynomial_ref y(m);
583 x = m.mk_polynomial(m.mk_var());
584 y = m.mk_polynomial(m.mk_var());
585 polynomial_ref p(m);
586 polynomial_ref q(m);
587 p = (x^2) - 2;
588 q = y*(x^3);
589 std::cout << "p: " << p << "\n";
590 std::cout << "q: " << q << "\n";
591 // unsigned d;
592 std::cout << "srem: " << exact_pseudo_remainder(q, p, 0) << "\n";
593 }
594
tst_sqrt()595 void tst_sqrt() {
596 reslimit rl;
597 polynomial::numeral_manager nm;
598 polynomial::manager m(rl, nm);
599 polynomial_ref x(m);
600 polynomial_ref y(m);
601 x = m.mk_polynomial(m.mk_var());
602 y = m.mk_polynomial(m.mk_var());
603 polynomial_ref p(m);
604 p = (4*x*y + 3*(x^2)*y + (y^2) + 3)^4;
605 polynomial_ref q(m);
606 VERIFY(sqrt(p, q));
607 ENSURE(eq(p, q*q));
608 std::cout << "p: " << p << "\n";
609 std::cout << "q: " << q << "\n";
610 p = p - 1;
611 ENSURE(!sqrt(p, q));
612 }
613
tst_content(polynomial_ref const & p,polynomial::var x,polynomial_ref const & expected)614 static void tst_content(polynomial_ref const & p, polynomial::var x, polynomial_ref const & expected) {
615 std::cout << "---------------\n";
616 std::cout << "p: " << p << std::endl;
617 std::cout << "content(p): " << content(p, x) << std::endl;
618 std::cout << "expected: " << expected << std::endl;
619 ENSURE(eq(content(p, x), expected));
620 }
621
tst_primitive(polynomial_ref const & p,polynomial::var x,polynomial_ref const & expected)622 static void tst_primitive(polynomial_ref const & p, polynomial::var x, polynomial_ref const & expected) {
623 std::cout << "---------------\n";
624 std::cout << "p: " << p << std::endl;
625 std::cout << "primitive(p): " << primitive(p, x) << std::endl;
626 std::cout << "expected: " << expected << std::endl;
627 ENSURE(eq(primitive(p, x), expected));
628 }
629
tst_gcd(polynomial_ref const & p,polynomial_ref const & q,polynomial_ref const & expected)630 static void tst_gcd(polynomial_ref const & p, polynomial_ref const & q, polynomial_ref const & expected) {
631 std::cout << "---------------\n";
632 std::cout << "p: " << p << std::endl;
633 std::cout << "q: " << q << std::endl;
634 polynomial_ref r(p.m());
635 r = gcd(p, q);
636 std::cout << "gcd(p, q): " << r << std::endl;
637 std::cout << "expected: " << expected << std::endl;
638 ENSURE(eq(r, expected));
639 }
640
tst_gcd()641 static void tst_gcd() {
642 reslimit rl;
643 polynomial::numeral_manager nm;
644 polynomial::manager m(rl, nm);
645 polynomial_ref x0(m);
646 polynomial_ref x1(m);
647 polynomial_ref x2(m);
648 polynomial_ref x3(m);
649 x0 = m.mk_polynomial(m.mk_var());
650 x1 = m.mk_polynomial(m.mk_var());
651 x2 = m.mk_polynomial(m.mk_var());
652 x3 = m.mk_polynomial(m.mk_var());
653 polynomial_ref three(m);
654 three = m.mk_const(mpz(3));
655
656 std::cout << "tst_gcd\n======================\n";
657
658 tst_gcd(((x0^2) + x0*x1 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*x1*x1 + x1*x2 + 1),
659 ((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1)*(x3*x1 + x1*x2 + 17),
660 ((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1));
661
662 tst_gcd((-1)*((x0^2) + x0*x1 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*x1*x1 + x1*x2 + 1),
663 ((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1)*(x3*x1 + x1*x2 + 17),
664 ((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1));
665
666 tst_gcd(((x0^2) + x0*x1 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*x1*x1 + x1*x2 + 1),
667 (-1)*((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1)*(x3*x1 + x1*x2 + 17),
668 ((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1));
669
670 tst_gcd((-1)*((x0^2) + x0*x1 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*x1*x1 + x1*x2 + 1),
671 (-1)*((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1)*(x3*x1 + x1*x2 + 17),
672 ((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1));
673
674 tst_gcd(21*(x0 + 1), 6*x0^2, three);
675 tst_content(x0*x1 + x0, 1, x0);
676 tst_primitive(x0*x1 + x0, 1, x1 + 1);
677 tst_primitive((x1^2) + x0*x1 + x0, 1, (x1^2) + x0*x1 + x0);
678 tst_primitive((x0 + 1)*(2*x1) + 1, 1, (x0 + 1)*(2*x1) + 1);
679 tst_primitive((x0 + 1)*(2*x1) + (x0^2)*(x0 + 1), 1, 2*x1 + (x0^2));
680 tst_primitive((x0 + 1)*(x2 + 1)*(x2^2)*(x0 + 1)*(x1^2) + (x0 + 1)*(x2^2)*x1 + (x0+1)*(x0+1), 1,
681 (x2 + 1)*(x2^2)*(x0 + 1)*(x1^2) + (x2^2)*x1 + (x0+1));
682 tst_primitive((x0 + (x3^2))*(x2 + x3 + 1)*(x2^2)*(x1^2) +
683 (x0 + (x3^2))*(x2 + x3 + 1)*x1 +
684 (x0 + (x3^2))*(x2 + x3 + 1)*(x3^2),
685 1,
686 (x2^2)*(x1^2) + x1 + (x3^2));
687 tst_content((x0 + (x3^2))*(x2 + x3 + 1)*(x2^2)*(x1^2) +
688 (x0 + (x3^2))*(x2 + x3 + 1)*x1 +
689 (x0 + (x3^2))*(x2 + x3 + 1)*(x3^2),
690 1,
691 (x0 + (x3^2))*(x2 + x3 + 1));
692 tst_primitive(4*(x0 + (x3^2))*(x2 + x3 + 1)*(x2^2)*(x1^2) +
693 2*(x0 + (x3^2))*(x2 + x3 + 1)*x1 +
694 4*(x0 + (x3^2))*(x2 + x3 + 1)*(x3^2),
695 1,
696 2*(x2^2)*(x1^2) + x1 + 2*(x3^2));
697 tst_gcd(63*((x0^2) + x0*x1 + 1)*(x2*x2 + x3 + 2)*(x3*x1 + 2)*(x3*x1*x1 + x1*x2 + 1),
698 14*((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1)*(x3*x1 + x1*x2 + 17),
699 7*((x0^2) + x0*x1 + 1)*(x3*x1*x1 + x1*x2 + 1));
700 }
701
tst_psc(polynomial_ref const & p,polynomial_ref const & q,polynomial::var x,polynomial_ref const & first,polynomial_ref const & second)702 static void tst_psc(polynomial_ref const & p, polynomial_ref const & q, polynomial::var x, polynomial_ref const & first, polynomial_ref const & second) {
703 polynomial::manager & m = p.m();
704 polynomial_ref_vector S(m);
705 std::cout << "---------" << std::endl;
706 std::cout << "p: " << p << std::endl;
707 std::cout << "q: " << q << std::endl;
708 m.psc_chain(p, q, x, S);
709 unsigned sz = S.size();
710 for (unsigned i = 0; i < sz; i++) {
711 std::cout << "S_" << i << ": " << polynomial_ref(S.get(i), m) << std::endl;
712 }
713 if (sz > 0) {
714 ENSURE(m.eq(S.get(0), first) || m.eq(S.get(0), neg(first)));
715 }
716 if (sz > 1) {
717 ENSURE(m.eq(S.get(1), second) || m.eq(S.get(1), neg(second)));
718 }
719 if (sz > 0) {
720 polynomial_ref Res(m);
721 Res = resultant(p, q, x);
722 ENSURE(m.eq(Res, S.get(0)) || m.eq(S.get(0), neg(Res)));
723 }
724 }
725
726 #if 0
727 static void tst_psc_perf(polynomial_ref const & p, polynomial_ref const & q, polynomial::var x) {
728 polynomial::manager & m = p.m();
729 polynomial_ref_vector S(m);
730 std::cout << "---------" << std::endl;
731 std::cout << "p: " << p << std::endl;
732 std::cout << "q: " << q << std::endl;
733 m.psc_chain(p, q, x, S);
734 unsigned sz = S.size();
735 for (unsigned i = 0; i < sz; i++) {
736 std::cout << "S_" << i << ": " << m.size(S.get(i)) << std::endl; // polynomial_ref(S.get(i), m) << std::endl;
737 }
738 }
739 #endif
740
tst_psc()741 static void tst_psc() {
742 reslimit rl;
743 polynomial::numeral_manager nm;
744 polynomial::manager m(rl, nm);
745 polynomial_ref x0(m);
746 polynomial_ref x1(m);
747 polynomial_ref x2(m);
748 polynomial_ref x3(m);
749 polynomial_ref x4(m);
750 polynomial_ref x5(m), x6(m), x7(m), x8(m), x9(m), x10(m);
751 x0 = m.mk_polynomial(m.mk_var());
752 x1 = m.mk_polynomial(m.mk_var());
753 x2 = m.mk_polynomial(m.mk_var());
754 x3 = m.mk_polynomial(m.mk_var());
755 x4 = m.mk_polynomial(m.mk_var());
756 x5 = m.mk_polynomial(m.mk_var());
757 x6 = m.mk_polynomial(m.mk_var());
758 x7 = m.mk_polynomial(m.mk_var());
759 x8 = m.mk_polynomial(m.mk_var());
760 x9 = m.mk_polynomial(m.mk_var());
761 x10 = m.mk_polynomial(m.mk_var());
762 tst_psc(x0*(x1^2) + (x0 + 1)*x1 + 2, x0*x1 + 3, 1,
763 6*x0 - (x0^2), x0);
764 tst_psc(x0*(x1^4) + (x0 + 1)*(x1^3) + 2, x0*(x1^3) + 3, 1,
765 72*(x0^3) - (x0^4) - 27*(x0^2) - 27*(x0), 9*(x0^3));
766 polynomial_ref & a = x0;
767 polynomial_ref & b = x1;
768 polynomial_ref & c = x2;
769 polynomial_ref & d = x3;
770 polynomial_ref & e = x4;
771 polynomial_ref & f = x5;
772 polynomial_ref & x = x9;
773 tst_psc((x^4) + a*(x^2) + b*x + c, 4*(x^3) + 2*a*x + b, 9,
774 16*(a^4)*c - 4*(a^3)*(b^2) - 128*(a^2)*(c^2) + 144*a*(b^2)*c - 27*(b^4) + 256*(c^3), 8*(a^3) - 32*a*c + 36*(b^2));
775 polynomial_ref & y = x10;
776
777 tst_psc(((y^2) + 6)*(x - 1) - y*((x^2) + 1), ((x^2) + 6)*(y - 1) - x*((y^2) + 1), 10,
778 2*(x^6) - 22*(x^5) + 102*(x^4) - 274*(x^3) + 488*(x^2) - 552*x + 288,
779 5*x - (x^2) - 6
780 );
781
782 tst_psc(((y^3) + 6)*(x - 1) - y*((x^3) + 1), ((x^3) + 6)*(y - 1) - x*((y^3) + 1), 10,
783 3*(x^11) - 3*(x^10) - 37*(x^9) + 99*(x^8) + 51*(x^7) - 621*(x^6) + 1089*(x^5) - 39*(x^4) - 3106*(x^3) + 5868*(x^2) - 4968*x + 1728,
784 (x^6) - 10*(x^4) + 12*(x^3) + 25*(x^2) - 60*x + 36);
785
786 polynomial_ref p = (x^6) + a * (x^3) + b;
787 polynomial_ref q = (x^6) + c * (x^3) + d;
788
789 tst_psc(p, q, 9,
790 (b^6) - 3*a*(b^5)*c + 3*(a^2)*(b^4)*(c^2) + 3*(b^5)*(c^2) - (a^3)*(b^3)*(c^3) -
791 6*a*(b^4)*(c^3) + 3*(a^2)*(b^3)*(c^4) + 3*(b^4)*(c^4) - 3*a*(b^3)*(c^5) + (b^3)*(c^6) +
792 3*(a^2)*(b^4)*d - 6*(b^5)*d - 6*(a^3)*(b^3)*c*d + 9*a*(b^4)*c*d +
793 3*(a^4)*(b^2)*(c^2)*d + 6*(a^2)*(b^3)*(c^2)*d - 12*(b^4)*(c^2)*d - 9*(a^3)*(b^2)*(c^3)*d +
794 6*a*(b^3)*(c^3)*d + 9*(a^2)*(b^2)*(c^4)*d - 6*(b^3)*(c^4)*d - 3*a*(b^2)*(c^5)*d +
795 3*(a^4)*(b^2)*(d^2) - 12*(a^2)*(b^3)*(d^2) + 15*(b^4)*(d^2) - 3*(a^5)*b*c*(d^2) +
796 6*(a^3)*(b^2)*c*(d^2) - 6*a*(b^3)*c*(d^2) + 9*(a^4)*b*(c^2)*(d^2) -
797 18*(a^2)*(b^2)*(c^2)*(d^2) + 18*(b^3)*(c^2)*(d^2) - 9*(a^3)*b*(c^3)*(d^2) +
798 6*a*(b^2)*(c^3)*(d^2) + 3*(a^2)*b*(c^4)*(d^2) + 3*(b^2)*(c^4)*(d^2) + (a^6)*(d^3) -
799 6*(a^4)*b*(d^3) + 18*(a^2)*(b^2)*(d^3) - 20*(b^3)*(d^3) - 3*(a^5)*c*(d^3) +
800 6*(a^3)*b*c*(d^3) - 6*a*(b^2)*c*(d^3) + 3*(a^4)*(c^2)*(d^3) + 6*(a^2)*b*(c^2)*(d^3) -
801 12*(b^2)*(c^2)*(d^3) - (a^3)*(c^3)*(d^3) - 6*a*b*(c^3)*(d^3) + 3*(a^4)*(d^4) -
802 12*(a^2)*b*(d^4) + 15*(b^2)*(d^4) - 6*(a^3)*c*(d^4) + 9*a*b*c*(d^4) +
803 3*(a^2)*(c^2)*(d^4) + 3*b*(c^2)*(d^4) + 3*(a^2)*(d^5) - 6*b*(d^5) -
804 3*a*c*(d^5) + (d^6),
805 3*(a^2)*c - (a^3) - 3*a*(c^2) + (c^3)
806 );
807
808
809 tst_psc(x,
810 a * x + b * c + d - e,
811 9,
812 b*c + d - e, a);
813
814 polynomial_ref zero(m);
815 zero = m.mk_zero();
816
817 tst_psc( a*d*x + a*c*f + a*e - b*a,
818 d*x + c*f + e - b,
819 9, zero, zero);
820
821
822 #if 0
823 tst_psc_perf((x^7) + a*(x^3) + b*(x^2) + c*x + d,
824 (x^7) + e*(x^3) + f*(x^2) + g*x + h,
825 9);
826
827 tst_psc_perf((x^15) + a * (x^10) + b,
828 (x^15) + c * (x^10) + d,
829 9);
830
831 tst_psc_perf((y^5) + a * (y^4) + b * (y^3) + c * (y^2) + d * y + e,
832 (y^5) + f * (y^4) + g * (y^3) + h * (y^2) + i * y + x,
833 10);
834 #endif
835 }
836
tst_vars(polynomial_ref const & p,unsigned sz,polynomial::var * xs)837 static void tst_vars(polynomial_ref const & p, unsigned sz, polynomial::var * xs) {
838 polynomial::var_vector r;
839 p.m().vars(p, r);
840 std::cout << "---------------\n";
841 std::cout << "p: " << p << "\nvars: ";
842 for (unsigned i = 0; i < r.size(); i++) {
843 std::cout << r[i] << " ";
844 }
845 std::cout << std::endl;
846 ENSURE(r.size() == sz);
847 std::sort(r.begin(), r.end());
848 std::sort(xs, xs + sz);
849 for (unsigned i = 0; i < r.size(); i++) {
850 ENSURE(r[i] == xs[i]);
851 }
852 }
853
tst_vars()854 static void tst_vars() {
855 polynomial::numeral_manager nm;
856 reslimit rl; polynomial::manager m(rl, nm);
857 polynomial_ref x0(m);
858 polynomial_ref x1(m);
859 polynomial_ref x2(m);
860 polynomial_ref x3(m);
861 polynomial_ref x4(m);
862 x0 = m.mk_polynomial(m.mk_var());
863 x1 = m.mk_polynomial(m.mk_var());
864 x2 = m.mk_polynomial(m.mk_var());
865 x3 = m.mk_polynomial(m.mk_var());
866 x4 = m.mk_polynomial(m.mk_var());
867 polynomial::var s023[3] = {0, 2, 3};
868 polynomial::var s14[2] = {1, 4};
869 polynomial::var s012[3] = {0, 1, 2};
870 polynomial::var s3[1] = {3};
871 polynomial::var s01234[5] = {0, 1, 2, 3, 4};
872
873 tst_vars((x0 + 1)*((x0^2) + (x3^2))*(x2*x3), 3, s023);
874 tst_vars((x0 + x2)*((x0^2) + (x3^2))*(x2*x3), 3, s023);
875 tst_vars(((x1 + x4 + 1)^5), 2, s14);
876 tst_vars(((x1 + x4*x2 + 1)^4) + x0 + (x3^2), 5, s01234);
877 tst_vars((x3 + 1)^5, 1, s3);
878 tst_vars(x0*x1*x2, 3, s012);
879 tst_vars(x0*x1*x2 + 1, 3, s012);
880 }
881
tst_sqf(polynomial_ref const & p,polynomial_ref const & expected)882 static void tst_sqf(polynomial_ref const & p, polynomial_ref const & expected) {
883 polynomial_ref r(p.m());
884 std::cout << "---------------\n";
885 std::cout << "p: " << p << std::endl;
886 r = square_free(p);
887 std::cout << "sqf(p): " << r << std::endl;
888 std::cout << "expected: " << expected << std::endl;
889 ENSURE(is_square_free(r));
890 ENSURE(!eq(r, p) || is_square_free(p));
891 ENSURE(eq(expected, r));
892 }
893
tst_sqf()894 static void tst_sqf() {
895 polynomial::numeral_manager nm;
896 reslimit rl; polynomial::manager m(rl, nm);
897 polynomial_ref x0(m);
898 polynomial_ref x1(m);
899 polynomial_ref x2(m);
900 polynomial_ref x3(m);
901 x0 = m.mk_polynomial(m.mk_var());
902 x1 = m.mk_polynomial(m.mk_var());
903 x2 = m.mk_polynomial(m.mk_var());
904 x3 = m.mk_polynomial(m.mk_var());
905 tst_sqf(((x0 + x1)^2)*((x0^2) - 3)*((x2*x2 + x3 + 1)^3),
906 (x0 + x1)*((x0^2) - 3)*(x2*x2 + x3 + 1));
907 tst_sqf((x0 + x1)*(x0 - x1), (x0 + x1)*(x0 - x1));
908 tst_sqf(((x0 + x1)^3)*(x0 - x1), (x0 + x1)*(x0 - x1));
909 polynomial_ref c1(m);
910 c1 = m.mk_const(rational(3));
911 tst_sqf(c1, c1);
912 polynomial_ref z(m);
913 z = m.mk_zero();
914 tst_sqf(z, z);
915 tst_sqf((x0 + x1 + x2 + x3)^5, (x0 + x1 + x2 + x3));
916 tst_sqf(((x0 + x1 + x2 + x3)^5) + 1, ((x0 + x1 + x2 + x3)^5) + 1);
917 }
918
tst_substitute(polynomial_ref const & p,polynomial::var x1,mpz const & v1,polynomial::var x2,mpz const & v2,polynomial_ref const & expected)919 static void tst_substitute(polynomial_ref const & p,
920 polynomial::var x1, mpz const & v1,
921 polynomial::var x2, mpz const & v2,
922 polynomial_ref const & expected) {
923 polynomial::numeral_manager & nm = p.m().m();
924 polynomial::var xs[2] = { x1, x2 };
925 scoped_mpz_vector vs(nm);
926 vs.push_back(v1);
927 vs.push_back(v2);
928 std::cout << "---------------\n";
929 std::cout << "p: " << p << std::endl;
930 polynomial_ref r(p.m());
931 r = p.m().substitute(p, 2, xs, vs.data());
932 std::cout << "r: " << r << std::endl;
933 std::cout << "expected: " << expected << std::endl;
934 ENSURE(eq(r, expected));
935 p.m().lex_sort(r);
936 std::cout << "r (sorted): " << r << std::endl;
937 }
938
tst_substitute()939 static void tst_substitute() {
940 polynomial::numeral_manager nm;
941 reslimit rl; polynomial::manager m(rl, nm);
942 polynomial_ref x0(m);
943 polynomial_ref x1(m);
944 polynomial_ref x2(m);
945 polynomial_ref x3(m);
946 x0 = m.mk_polynomial(m.mk_var());
947 x1 = m.mk_polynomial(m.mk_var());
948 x2 = m.mk_polynomial(m.mk_var());
949 x3 = m.mk_polynomial(m.mk_var());
950 tst_substitute(x0 + x1*x0 + x3, 1, mpz(1), 3, mpz(2), 2*x0 + 2);
951 tst_substitute((x0^2) + x1*x0 + x3, 0, mpz(2), 3, mpz(2), 2*x1 + 6);
952 tst_substitute((x0 + x1 + x2)^3, 0, mpz(2), 2, mpz(3), (x1 + 5)^3);
953 tst_substitute(((x0 + x1 + x2)^3) + ((x0*x1 + x3)^2), 0, mpz(2), 2, mpz(3), ((x1 + 5)^3) + ((2*x1 + x3)^2));
954 tst_substitute((x0 + x1 + 1)^5, 2, mpz(2), 3, mpz(3), (x0 + x1 + 1)^5);
955 polynomial_ref zero(m);
956 zero = m.mk_zero();
957 tst_substitute(zero, 2, mpz(2), 3, mpz(3), zero);
958 }
959
tst_qsubstitute(polynomial_ref const & p,unsynch_mpq_manager & qm,polynomial::var x1,rational const & v1,polynomial::var x2,rational const & v2,polynomial_ref const & expected)960 static void tst_qsubstitute(polynomial_ref const & p,
961 unsynch_mpq_manager & qm,
962 polynomial::var x1, rational const & v1,
963 polynomial::var x2, rational const & v2,
964 polynomial_ref const & expected) {
965 polynomial::var xs[2] = { x1, x2 };
966 scoped_mpq_vector vs(qm);
967 vs.push_back(v1.to_mpq());
968 vs.push_back(v2.to_mpq());
969 std::cout << "---------------\n";
970 std::cout << "p: " << p << std::endl;
971 polynomial_ref r(p.m());
972 r = p.m().substitute(p, 2, xs, vs.data());
973 std::cout << "r: " << r << std::endl;
974 std::cout << "expected (modulo a constant): " << expected << std::endl;
975 ENSURE(eq(r, normalize(expected)));
976 p.m().lex_sort(r);
977 std::cout << "r (sorted): " << r << std::endl;
978 }
979
tst_qsubstitute()980 static void tst_qsubstitute() {
981 unsynch_mpq_manager qm;
982 reslimit rl; polynomial::manager m(rl, qm);
983 polynomial_ref x0(m);
984 polynomial_ref x1(m);
985 polynomial_ref x2(m);
986 polynomial_ref x3(m);
987 x0 = m.mk_polynomial(m.mk_var());
988 x1 = m.mk_polynomial(m.mk_var());
989 x2 = m.mk_polynomial(m.mk_var());
990 x3 = m.mk_polynomial(m.mk_var());
991 tst_qsubstitute(x0 + x1 + x2, qm, 0, rational(1)/rational(2), 1, rational(3), 2*x2 + 2*3 + 1);
992 tst_qsubstitute((x0^2)*x2 + x1 + x2*x0, qm, 0, rational(3)/rational(2), 1, rational(3), 9*x2 + (4*3) + (2*3)*x2);
993 tst_qsubstitute(x0*x1*x2*(x3^2) + (x0^3)*x3 + (x1^2)*x0, qm,
994 0, rational(5)/rational(2),
995 1, rational(7)/rational(3),
996 (2*2*3*7*5)*x2*(x3^2) + (5*5*5*3*3)*x3 + (7*7*5*2*2));
997 tst_qsubstitute((x2 + x3)^3, qm,
998 0, rational(5)/rational(2),
999 1, rational(7)/rational(3),
1000 (x2 + x3)^3);
1001 tst_qsubstitute((x0 + x1 + x2 + x3)^3, qm,
1002 0, rational(5)/rational(2),
1003 2, rational(7)/rational(3),
1004 (6*x1 + 6*x3 + 15 + 14)^3);
1005 tst_qsubstitute((x0 + 2*x1 + 11*(x2^2)*x3 + x2 + (x3^2))^3, qm,
1006 0, rational(5)/rational(2),
1007 3, rational(7)/rational(3),
1008 ((2*3*3*2)*x1 + (11*2*3*7)*(x2^2) + (2*3*3)*x2 + (5*9 + 7*7*2))^3);
1009 polynomial_ref zero(m);
1010 zero = m.mk_zero();
1011 tst_qsubstitute(zero, qm,
1012 0, rational(5)/rational(2),
1013 3, rational(7)/rational(3),
1014 zero);
1015 }
1016
tst_mfact(polynomial_ref const & p,unsigned num_distinct_factors)1017 void tst_mfact(polynomial_ref const & p, unsigned num_distinct_factors) {
1018 std::cout << "---------------\n";
1019 std::cout << "p: " << p << std::endl;
1020 polynomial::factors fs(p.m());
1021 factor(p, fs);
1022 std::cout << "factors:\n";
1023 std::cout << p.m().m().to_string(fs.get_constant()) << "\n";
1024 for (unsigned i = 0; i < fs.distinct_factors(); i++) {
1025 std::cout << "*(" << fs[i] << ")^" << fs.get_degree(i) << std::endl;
1026 }
1027 ENSURE(fs.distinct_factors() == num_distinct_factors);
1028 polynomial_ref p2(p.m());
1029 fs.multiply(p2);
1030 ENSURE(eq(p, p2));
1031 }
1032
tst_mfact()1033 static void tst_mfact() {
1034 polynomial::numeral_manager nm;
1035 reslimit rl; polynomial::manager m(rl, nm);
1036 polynomial_ref x0(m);
1037 polynomial_ref x1(m);
1038 polynomial_ref x2(m);
1039 polynomial_ref x3(m);
1040 polynomial_ref x4(m);
1041 polynomial_ref x5(m);
1042 x0 = m.mk_polynomial(m.mk_var());
1043 x1 = m.mk_polynomial(m.mk_var());
1044 x2 = m.mk_polynomial(m.mk_var());
1045 x3 = m.mk_polynomial(m.mk_var());
1046 x4 = m.mk_polynomial(m.mk_var());
1047 x5 = m.mk_polynomial(m.mk_var());
1048 polynomial_ref & x = x0;
1049
1050 tst_mfact((x0 - (x1^3))*(x0 - ((x2^3) - 2)), 2);
1051 tst_mfact((((x3^3) + 1)*x0 - (x1^3))*(x0 - ((x2^3) - 2)), 2);
1052 tst_mfact((((x3^3) + 1)*x0 - (x1^3))*((x3 - 1)*x0 - ((x2^3) - 2)), 2);
1053 tst_mfact((((x3^3) + 1)*x0 - (x1^3))*((-1)*x0 - ((x2^3) - 2)), 2);
1054 tst_mfact((-1)*(x0 - (x1^3))*(x0 - ((x2^3) - 2)), 2);
1055 tst_mfact((-1)*(x0 - (x1^3) + (x1^2) + 2)*(x0 - ((x2^3) - 2)), 2);
1056
1057 tst_mfact(((x0 - (x1^3))*(x0 - ((x2^3) - 2)))^2, 2);
1058 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*(x0 - ((x2^3) - 2)))^2, 2);
1059 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*((x3 - 1)*x0 - ((x2^3) - 2)))^2, 2);
1060 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*((-1)*x0 - ((x2^3) - 2)))^2, 2);
1061 tst_mfact(((-1)*(x0 - (x1^3))*(x0 - ((x2^3) - 2)))^2, 2);
1062 tst_mfact(((-1)*(x0 - (x1^3) + (x1^2) + 2)*(x0 - ((x2^3) - 2)))^2, 2);
1063
1064 tst_mfact(((x0 - (x1^3))*(x0 - ((x2^3) - 2)))^3, 2);
1065 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*(x0 - ((x2^3) - 2)))^3, 2);
1066 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*((x3 - 1)*x0 - ((x2^3) - 2)))^3, 2);
1067 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*((-1)*x0 - ((x2^3) - 2)))^3, 2);
1068 tst_mfact(((-1)*(x0 - (x1^3))*(x0 - ((x2^3) - 2)))^3, 2);
1069 tst_mfact(((-1)*(x0 - (x1^3) + (x1^2) + 2)*(x0 - ((x2^3) - 2)))^3, 2);
1070
1071 tst_mfact(((x0 - (x1^3))*(x0 - ((x2^3) - 2)))^4, 2);
1072 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*(x0 - ((x2^3) - 2)))^4, 2);
1073 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*((x3 - 1)*x0 - ((x2^3) - 2)))^4, 2);
1074 tst_mfact(((((x3^3) + 1)*x0 - (x1^3))*((-1)*x0 - ((x2^3) - 2)))^4, 2);
1075 tst_mfact(((-1)*(x0 - (x1^3))*(x0 - ((x2^3) - 2)))^4, 2);
1076 tst_mfact(((-1)*(x0 - (x1^3) + (x1^2) + 2)*(x0 - ((x2^3) - 2)))^4, 2);
1077
1078 tst_mfact(((x^5) - (x^2) + 1)*((-1)*x + 1)*((x^2) - 2*x + 3), 3);
1079 tst_mfact(11*((x^5) - (x^2) + 1)*((-1)*x + 1)*((x^2) - 2*x + 3), 3);
1080 tst_mfact(11*(7*(x^5) - (x^2) + 1)*((-1)*(x^2) + 1)*((x^2) - 2*x + 3), 4);
1081 tst_mfact(11*(7*(x^5) - (x^2) + 1)*((-1)*(x^2) + 1)*((x^2) - 2*x + 3)*((x^7) - x +2), 5);
1082 tst_mfact((7*(x^5) - (x^2) + 1)*((-1)*(x^2) + 1)*((x^2) - 2*x + 3)*((x^7) - x +2)*((x^3) - x + 1), 6);
1083 tst_mfact((7*(x^5) - (x^2) + 1)*((-1)*(x^2) + 1)*((x^2) - 2*x + 3)*((x^7) - x +2)*((x^3) - x + 1)*((x^7) - (x^5) + (x^3) + (x^2) + x + 3), 7);
1084 tst_mfact((7*(x^5) - (x^2) + 1)*((-1)*(x^2) + 1)*((x^2) - 2*x + 3)*((x^7) - x +2)*
1085 ((x^3) - x + 1)*((x^7) - (x^5) + (x^3) + (x^2) + x + 3)*(x - (x^3) + 11)*
1086 (x - 10)*(x - 9)*(33*x + 12)*((x^5) - x + 1),
1087 12);
1088 tst_mfact((x^4) + (x^2) - 20, 3);
1089 tst_mfact((-11)*((x^5) - (x^2) + 1)*((-1)*x + 1)*((x^2) - 2*x + 3), 3);
1090 tst_mfact(x0 - 2*(x0^2) + 1, 2);
1091 tst_mfact((x0 + 1)*(x0 - 1)*(x0 + 2)*(((x1^5) - x1 - 1)^2), 4);
1092 tst_mfact((x0 + 1)*((x1 + 2)^2), 2);
1093 tst_mfact(7*(x0 + 1)*((x1 + 2)^2), 2);
1094 tst_mfact(11*(x0 + 1)*((x1 + 2)^2)*((x1 - x3)^4), 3);
1095 tst_mfact(11*(x0 + 1)*((x1 + 2)^2)*((x1 - x3)^4)*(((x0*(x2^2) + (x0 + x1)*x2 + x1))^3), 5);
1096 tst_mfact((11*(x0 + 1)*((x1 + 2)^2))^3, 2);
1097 tst_mfact((3*(x0 + 1)*(x1 + 2))^3, 2);
1098 tst_mfact((3*(x0 + 1)*(2*x1 + 2))^3, 2);
1099 tst_mfact(3*(2*(x0^2) + 4*(x1^2))*x2, 2);
1100 tst_mfact(13*((x0 - x2)^6)*((x1 - x2)^5)*((x0 - x3)^7), 3);
1101 tst_mfact((x0+1)^100, 1);
1102 tst_mfact((x0^70) - 6*(x0^65) - (x0^60) + 60*(x0^55) - 54*(x0^50) - 230*(x0^45) + 274*(x0^40) + 542*(x0^35) - 615*(x0^30) - 1120*(x0^25) + 1500*(x0^20) - 160*(x0^15) - 395*(x0^10) + 76*(x0^5) + 34, 3);
1103 tst_mfact(((x0^4) - 8*(x0^2)), 2);
1104 tst_mfact((x0^5) - 2*(x0^3) + x0 - 1, 1);
1105 tst_mfact( (x0^25) - 4*(x0^21) - 5*(x0^20) + 6*(x0^17) + 11*(x0^16) + 10*(x0^15) - 4*(x0^13) - 7*(x0^12) - 9*(x0^11) - 10*(x0^10) +
1106 (x0^9) + (x0^8) + (x0^7) + (x0^6) + 3*(x0^5) + x0 - 1, 2);
1107 tst_mfact( (x0^25) - 10*(x0^21) - 10*(x0^20) - 95*(x0^17) - 470*(x0^16) - 585*(x0^15) - 40*(x0^13) - 1280*(x0^12) - 4190*(x0^11) - 3830*(x0^10) + 400*(x0^9)+ 1760*(x0^8) + 760*(x0^7) - 2280*(x0^6) + 449*(x0^5) + 640*(x0^3) - 640*(x0^2) + 240*x0 - 32, 2);
1108 tst_mfact( x0^10, 1);
1109 polynomial_ref c(m);
1110 c = m.mk_zero();
1111 tst_mfact(c, 0);
1112 c = m.mk_const(mpz(3));
1113 tst_mfact(c, 0);
1114 tst_mfact(x0, 1);
1115 tst_mfact(x0 + x1, 1);
1116 tst_mfact(x0 - x1, 1);
1117 tst_mfact( (x0^10) - 10*(x0^8) + 38*(x0^6) - 2*(x0^5) - 100*(x0^4) - 40*(x0^3) + 121*(x0^2) - 38*x0 - 17, 1);
1118 tst_mfact( (x0^50) - 10*(x0^40) + 38*(x0^30) - 2*(x0^25) - 100*(x0^20) - 40*(x0^15) + 121*(x0^10) - 38*(x0^5) - 17, 1);
1119 polynomial_ref & y = x0;
1120 tst_mfact( (((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^10)
1121 + 10*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^9)
1122 + 35*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^8)
1123 + 40*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^7)
1124 - 32*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^6)
1125 - 82*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^5)
1126 - 30*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^4)
1127 - 140*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^3)
1128 - 284*(((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)^2)
1129 - 168*((y^5) + 5*(y^4) + 10*(y^3) + 10*(y^2) + 5*y)
1130 - 47, 1);
1131 tst_mfact( (y^4) - 404*(y^2) + 39204, 2);
1132 tst_mfact( ((y^5) - 15552)*
1133 ((y^20)- 15708*(y^15) + rational("138771724")*(y^10)- rational("432104148432")*(y^5) + rational("614198284585616")),
1134 2);
1135 tst_mfact( (y^25) -
1136 rational("3125")*(y^21) -
1137 rational("15630")*(y^20) +
1138 rational("3888750")*(y^17) +
1139 rational("38684375")*(y^16) +
1140 rational("95765635")*(y^15) -
1141 rational("2489846500")*(y^13) -
1142 rational("37650481875")*(y^12) -
1143 rational("190548065625")*(y^11) -
1144 rational("323785250010")*(y^10) +
1145 rational("750249453025")*(y^9) +
1146 rational("14962295699875")*(y^8) +
1147 rational("111775113235000")*(y^7) +
1148 rational("370399286731250")*(y^6) +
1149 rational("362903064503129")*(y^5) -
1150 rational("2387239013984400")*(y^4) -
1151 rational("23872390139844000")*(y^3) -
1152 rational("119361950699220000")*(y^2) -
1153 rational("298404876748050000")*y -
1154 rational("298500366308609376"), 2);
1155
1156 tst_mfact( rational("54")*(y^24) - (y^27) - 324*(y^21) + rational("17496")*(y^18) - 34992*(y^15)+ rational("1889568")*(y^12)- 1259712*(y^9) + rational("68024448")*(y^6), 3);
1157
1158 tst_mfact( ((y^3)- 432)*(((y^3)+54)^2)*((y^6)+108)*((y^6)+6912)*((y^6)- 324*(y^3)+37044),
1159 5);
1160
1161 tst_mfact( ((y^6)- 6*(y^4) - 864*(y^3) + 12*(y^2) - 5184*y + 186616)*
1162 (((y^6) - 6*(y^4) + 108*(y^3) + 12*(y^2) + 648*y + 2908)^2)*
1163 ((y^12) - 12*(y^10) + 60*(y^8) + 56*(y^6) + 6720*(y^4) + 12768*(y^2) + 13456)*
1164 ((y^12) - 12*(y^10) + 60*(y^8) + 13664*(y^6) + 414960*(y^4) + 829248*(y^2) + 47886400)*
1165 ((y^12) - 12*(y^10) - 648*(y^9)+ 60*(y^8) + 178904*(y^6) + 15552*(y^5) + 1593024*(y^4) - 24045984*(y^3) + 5704800*(y^2) - 143995968*y + 1372010896),
1166 5);
1167
1168 {
1169 polynomial_ref q1(m);
1170 polynomial_ref q2(m);
1171 polynomial_ref q3(m);
1172 polynomial_ref q4(m);
1173 polynomial_ref q5(m);
1174 polynomial_ref p(m);
1175 q1 = (x0^3) - 2;
1176 q2 = (x1^3) - 2;
1177 q3 = (x2^3) - 2;
1178 q4 = (x3^2) - 2;
1179 q5 = (x4^7) - x4 + 3;
1180 p = x5 - x0 - 2*x1 /* - 3*x2 - x3 */ + x4;
1181 p = resultant(p, q1, 0);
1182 std::cout << "finished resultant 1... size: " << size(p) << std::endl;
1183 p = resultant(p, q2, 1);
1184 std::cout << "finished resultant 2... size: " << size(p) << std::endl;
1185 // p = resultant(p, q3, 2);
1186 std::cout << "finished resultant 3... size: " << size(p) << std::endl;
1187 // p = resultant(p, q4, 3);
1188 std::cout << "finished resultant 4... size: " << size(p) << std::endl;
1189 p = resultant(p, q5, 4);
1190 tst_mfact(p, 2);
1191 }
1192 }
1193
1194
tst_zp()1195 static void tst_zp() {
1196 unsynch_mpz_manager z;
1197 reslimit rl; polynomial::manager pm(rl, z);
1198
1199 polynomial_ref x(pm);
1200 polynomial_ref y(pm);
1201 x = pm.mk_polynomial(pm.mk_var());
1202 y = pm.mk_polynomial(pm.mk_var());
1203
1204 polynomial_ref p(pm);
1205 polynomial_ref q(pm);
1206 p = (x^4) + 2*(x^3) + 2*(x^2) + x;
1207 q = (x^3) + x + 1;
1208 std::cout << "p: " << p << "\n";
1209 std::cout << "q: " << q << "\n";
1210 std::cout << "gcd: " << gcd(p, q) << "\n";
1211
1212 {
1213 polynomial::scoped_set_zp setZ3(pm, 3);
1214 polynomial_ref p3(pm);
1215 polynomial_ref q3(pm);
1216 p3 = normalize(p);
1217 q3 = normalize(q);
1218 std::cout << "p[Z_3]: " << p3 << "\n";
1219 std::cout << "q[Z_3]: " << q3 << "\n";
1220 std::cout << "gcd[Z_3]: " << gcd(p3, q3) << "\n";
1221 }
1222
1223 std::cout << "back into Z[x,y]\ngcd: " << gcd(p, q) << "\n";
1224
1225 p = 5*(x^2)*(y^2) + 3*(x^3) + 7*(y^3) + 3;
1226 {
1227 polynomial::scoped_set_zp setZ11(pm, 11);
1228 polynomial_ref p11(pm);
1229
1230 std::cout << "---------------\n";
1231 p11 = normalize(p);
1232 std::cout << "p[Z_11]: " << p11 << "\n";
1233 p11 = pm.mk_glex_monic(p11);
1234 std::cout << "monic p[Z_11]: " << p11 << "\n";
1235 }
1236 std::cout << "back into Z[x,y]\n";
1237 std::cout << "p: " << p << "\n";
1238 std::cout << "gcd: " << gcd(p, q) << "\n";
1239 }
1240
tst_translate(polynomial_ref const & p,polynomial::var x0,int v0,polynomial::var x1,int v1,polynomial::var x2,int v2,polynomial_ref const & expected)1241 static void tst_translate(polynomial_ref const & p, polynomial::var x0, int v0, polynomial::var x1, int v1, polynomial::var x2, int v2,
1242 polynomial_ref const & expected) {
1243 std::cout << "---------------\n";
1244 std::cout << "p: " << p << std::endl;
1245 polynomial::var xs[3] = { x0, x1, x2 };
1246 mpz vs[3] = { mpz(v0), mpz(v1), mpz(v2) };
1247 polynomial_ref r(p.m());
1248 p.m().translate(p, 3, xs, vs, r);
1249 std::cout << "r: " << r << std::endl;
1250 ENSURE(eq(expected, r));
1251 }
1252
tst_translate()1253 static void tst_translate() {
1254 unsynch_mpq_manager qm;
1255 reslimit rl; polynomial::manager m(rl, qm);
1256 polynomial_ref x0(m);
1257 polynomial_ref x1(m);
1258 polynomial_ref x2(m);
1259 polynomial_ref x3(m);
1260 x0 = m.mk_polynomial(m.mk_var());
1261 x1 = m.mk_polynomial(m.mk_var());
1262 x2 = m.mk_polynomial(m.mk_var());
1263 x3 = m.mk_polynomial(m.mk_var());
1264 tst_translate((x0^2) + x1 + 1, 0, 1, 1, 2, 3, 0,
1265 (x0^2) + 2*x0 + x1 + 4
1266 );
1267 tst_translate(x3 + 1, 0, 1, 1, 2, 2, 3,
1268 x3 + 1
1269 );
1270 tst_translate(x3 + 1, 0, 1, 1, 2, 3, 0,
1271 x3 + 1
1272 );
1273 tst_translate(x3 + 1, 0, 1, 1, 2, 3, 10,
1274 x3 + 11
1275 );
1276 tst_translate((x0^3)*(x1^2) + (x0^2)*(x1^3) + 10, 0, -3, 1, -2, 3, 0,
1277 (x0^3)*(x1^2) + (x0^2)*(x1^3) - 4*(x0^3)*x1 - 15*(x0^2)*(x1^2) - 6*x0*(x1^3) + 4*(x0^3) +
1278 48*(x0^2)*x1 + 63*x0*(x1^2) + 9*(x1^3) - 44*(x0^2) - 180*x0*x1 - 81*(x1^2) +
1279 156*x0 + 216*x1 - 170
1280 );
1281 }
1282
1283 #if 0
1284 static void tst_p25() {
1285 unsynch_mpq_manager qm;
1286 reslimit rl; polynomial::manager m(rl, qm);
1287 polynomial_ref x0(m);
1288 polynomial_ref x1(m);
1289 polynomial_ref x2(m);
1290 polynomial_ref x3(m);
1291 polynomial_ref x4(m);
1292 polynomial_ref x5(m);
1293 polynomial_ref x6(m);
1294 x0 = m.mk_polynomial(m.mk_var());
1295 x1 = m.mk_polynomial(m.mk_var());
1296 x2 = m.mk_polynomial(m.mk_var());
1297 x3 = m.mk_polynomial(m.mk_var());
1298 x4 = m.mk_polynomial(m.mk_var());
1299 x5 = m.mk_polynomial(m.mk_var());
1300 x6 = m.mk_polynomial(m.mk_var());
1301 polynomial_ref p(m);
1302 p = (x0 + x1 + x2 + x3 + x4 + x5 + x6)^25;
1303 std::cout << "size(p): " << size(p) << "\n";
1304 }
1305 #endif
1306
tst_mm()1307 static void tst_mm() {
1308 unsynch_mpq_manager qm;
1309 // pm1 and pm2 share the same monomial manager
1310 reslimit rl;
1311 polynomial::manager * pm1_ptr = alloc(polynomial::manager, rl, qm);
1312 polynomial::manager & pm1 = *pm1_ptr;
1313 polynomial::manager pm2(rl, qm, &pm1.mm());
1314 polynomial::manager pm3(rl, qm); // pm3 has its own manager
1315 polynomial_ref p2(pm2);
1316 {
1317 polynomial_ref x0(pm1);
1318 polynomial_ref x1(pm1);
1319 polynomial_ref x2(pm1);
1320 x0 = pm1.mk_polynomial(pm1.mk_var());
1321 x1 = pm1.mk_polynomial(pm1.mk_var());
1322 x2 = pm1.mk_polynomial(pm1.mk_var());
1323 polynomial_ref p1(pm1);
1324 p1 = (x0 + x1 + x2)^2;
1325
1326 std::cout << "p1: " << p1 << "\n";
1327 p2 = convert(pm1, p1, pm2);
1328 std::cout << "p2: " << p2 << "\n";
1329 ENSURE(pm1.get_monomial(p1, 0) == pm2.get_monomial(p2, 0));
1330
1331 polynomial_ref p3(pm3);
1332 p3 = convert(pm1, p1, pm3);
1333 ENSURE(pm1.get_monomial(p1, 0) != pm3.get_monomial(p3, 0));
1334 }
1335 dealloc(pm1_ptr);
1336 // p2 is still ok
1337 std::cout << "p2: " << p2 << "\n";
1338 }
1339
tst_eval(polynomial_ref const & p,polynomial::var x0,rational v0,polynomial::var x1,rational v1,polynomial::var x2,rational v2,rational expected)1340 static void tst_eval(polynomial_ref const & p, polynomial::var x0, rational v0, polynomial::var x1, rational v1, polynomial::var x2, rational v2,
1341 rational expected) {
1342 TRACE("eval_bug", tout << "tst_eval, " << p << "\n";);
1343 std::cout << "p: " << p << "\nx" << x0 << " -> " << v0 << "\nx" << x1 << " -> " << v1 << "\nx" << x2 << " -> " << v2 << "\n";
1344 unsynch_mpq_manager qm;
1345 polynomial::simple_var2value<unsynch_mpq_manager> x2v(qm);
1346 x2v.push_back(x0, v0.to_mpq());
1347 x2v.push_back(x1, v1.to_mpq());
1348 x2v.push_back(x2, v2.to_mpq());
1349 scoped_mpq r(qm);
1350 p.m().eval(p, x2v, r);
1351 std::cout << "r: " << r << "\nexpected: " << expected << "\n";
1352 scoped_mpq ex(qm);
1353 qm.set(ex, expected.to_mpq());
1354 ENSURE(qm.eq(r, ex));
1355 }
1356
tst_eval()1357 static void tst_eval() {
1358 polynomial::numeral_manager nm;
1359 reslimit rl; polynomial::manager m(rl, nm);
1360 polynomial_ref x0(m);
1361 polynomial_ref x1(m);
1362 polynomial_ref x2(m);
1363 x0 = m.mk_polynomial(m.mk_var());
1364 x1 = m.mk_polynomial(m.mk_var());
1365 x2 = m.mk_polynomial(m.mk_var());
1366 tst_eval(2000*x1 - x0, 0, rational(1), 1, rational(2), 2, rational(3), rational(3999));
1367 tst_eval(x0 + 1, 0, rational(1), 1, rational(2), 2, rational(3), rational(2));
1368 tst_eval((x0^3) + x0 + 1, 0, rational(2), 1, rational(2), 2, rational(3), rational(11));
1369 tst_eval((x0^3) - 2*x0 + 1, 0, rational(2), 1, rational(2), 2, rational(3), rational(5));
1370 tst_eval((x0^3) - 2*x0 + 1, 0, rational(-2), 1, rational(2), 2, rational(3), rational(-3));
1371 tst_eval((x0^4) - 2*x0 + x1 + 1, 0, rational(-2), 1, rational(10), 2, rational(3), rational(31));
1372 tst_eval((x0^4) - 2*x0 + ((x0^3) + 1)*x1 + 1, 0, rational(-2), 1, rational(10), 2, rational(3), rational(-49));
1373 tst_eval(((x0^4) - 2*x0)*(x1^2) + ((x0^3) + 1)*x1 + (x0^2) + 1, 0, rational(-2), 1, rational(10), 2, rational(3), rational(1935));
1374 tst_eval(((x0^4) - 2*x0)*(x1^2)*(x2^3) + ((x0^3) + 1)*x1 + (x0^2) + 1, 0, rational(-2), 1, rational(10), 2, rational(0), rational(-65));
1375 tst_eval(((x0^4) - 2*x0)*((x1^2) + 1)*(x2^3) + ((x0^3) + 1)*x1 + (x0^2) + 1, 0, rational(-2), 1, rational(10), 2, rational(1, 2), rational(375, 2));
1376 tst_eval(x0*x1*x2, 0, rational(2), 1, rational(3), 2, rational(1), rational(6));
1377 tst_eval(x0*x1*x2 + 1, 0, rational(2), 1, rational(3), 2, rational(1), rational(7));
1378 polynomial_ref one(m);
1379 one = x0 - x0 + 1;
1380 tst_eval(one, 0, rational(2), 1, rational(3), 2, rational(1), rational(1));
1381 tst_eval(x0*(x1^2)*x2 + 1, 0, rational(2), 1, rational(3), 2, rational(1), rational(19));
1382 tst_eval(x0*(x1^2)*x2 + x1 + 1, 0, rational(2), 1, rational(3), 2, rational(1), rational(22));
1383 tst_eval(x0*(x1^2)*x2 + x1 + 1 + (x2^2)*(2*x1 - 1), 0, rational(2), 1, rational(3), 2, rational(1), rational(27));
1384 tst_eval((x0^5) + 1, 0, rational(2), 1, rational(3), 2, rational(1), rational(33));
1385 tst_eval((x0^5) + x0*x1 + 1, 0, rational(2), 1, rational(1), 2, rational(5), rational(35));
1386 tst_eval((x1^5) + x0*x1 + 1, 0, rational(2), 1, rational(1), 2, rational(5), rational(4));
1387 tst_eval((x1^5) + x0*(x1^2) + 1, 0, rational(2), 1, rational(-2), 2, rational(5), rational(-23));
1388 }
1389
tst_mk_unique()1390 static void tst_mk_unique() {
1391 polynomial::numeral_manager nm;
1392 reslimit rl; polynomial::manager m(rl, nm);
1393 polynomial_ref x0(m);
1394 polynomial_ref x1(m);
1395 polynomial_ref x2(m);
1396 x0 = m.mk_polynomial(m.mk_var());
1397 x1 = m.mk_polynomial(m.mk_var());
1398 x2 = m.mk_polynomial(m.mk_var());
1399 polynomial::cache uniq(m);
1400 polynomial_ref p(m);
1401 polynomial_ref q(m);
1402 polynomial_ref r(m);
1403
1404 p = (x0^3) + (x2^5) + x0*x1 + x0*x1*x1 + 3*x0*x0 + 5;
1405 q = x0*x1*x1 + (x0^3) + 3*x0*x0 + (x2^5) + 5 + x0*x1;
1406 r = x0*x1*x1 + (x0^3) + 3*x0*x0 + (x2^5) + 6 + x0*x1;
1407 std::cout << "p: " << p << "\n";
1408 std::cout << "q: " << q << "\n";
1409 std::cout << "r: " << r << "\n";
1410 ENSURE(m.eq(p, q));
1411 ENSURE(!m.eq(p, r));
1412 ENSURE(p.get() != q.get());
1413 q = uniq.mk_unique(q);
1414 p = uniq.mk_unique(p);
1415 r = uniq.mk_unique(r);
1416 std::cout << "after mk_unique\np: " << p << "\n";
1417 std::cout << "q: " << q << "\n";
1418 std::cout << "r: " << r << "\n";
1419 ENSURE(m.eq(p, q));
1420 ENSURE(!m.eq(r, q));
1421 ENSURE(p.get() == q.get());
1422 }
1423
1424 struct dummy_del_eh : public polynomial::manager::del_eh {
1425 unsigned m_counter;
dummy_del_ehdummy_del_eh1426 dummy_del_eh():m_counter(0) {}
operator ()dummy_del_eh1427 virtual void operator()(polynomial::polynomial * p) {
1428 m_counter++;
1429 }
1430 };
1431
tst_del_eh()1432 static void tst_del_eh() {
1433 dummy_del_eh eh1;
1434 dummy_del_eh eh2;
1435
1436 polynomial::numeral_manager nm;
1437 reslimit rl; polynomial::manager m(rl, nm);
1438 polynomial_ref x0(m);
1439 polynomial_ref x1(m);
1440 x0 = m.mk_polynomial(m.mk_var());
1441 x1 = m.mk_polynomial(m.mk_var());
1442
1443 m.add_del_eh(&eh1);
1444 x1 = 0;
1445 ENSURE(eh1.m_counter == 1);
1446
1447 m.add_del_eh(&eh2);
1448 x1 = m.mk_polynomial(m.mk_var());
1449 x1 = 0;
1450 ENSURE(eh1.m_counter == 2);
1451 ENSURE(eh2.m_counter == 1);
1452 m.remove_del_eh(&eh1);
1453 x0 = 0;
1454 x1 = m.mk_polynomial(m.mk_var());
1455 x1 = 0;
1456 ENSURE(eh1.m_counter == 2);
1457 ENSURE(eh2.m_counter == 3);
1458 m.remove_del_eh(&eh2);
1459 x1 = m.mk_polynomial(m.mk_var());
1460 x1 = 0;
1461 ENSURE(eh1.m_counter == 2);
1462 ENSURE(eh2.m_counter == 3);
1463 }
1464
tst_const_coeff()1465 static void tst_const_coeff() {
1466 polynomial::numeral_manager nm;
1467 reslimit rl; polynomial::manager m(rl, nm);
1468 polynomial_ref x0(m);
1469 polynomial_ref x1(m);
1470 x0 = m.mk_polynomial(m.mk_var());
1471 x1 = m.mk_polynomial(m.mk_var());
1472
1473 scoped_mpz c(nm);
1474
1475 polynomial_ref p(m);
1476
1477 p = (x0^2)*x1 + 3*x0 + x1;
1478 ENSURE(!m.const_coeff(p, 0, 2, c));
1479 ENSURE(m.const_coeff(p, 0, 1, c) && c == 3);
1480 ENSURE(!m.const_coeff(p, 0, 0, c));
1481
1482 p = (x0^2)*x1 + 3*x0 + x1 + 1;
1483 ENSURE(!m.const_coeff(p, 0, 2, c));
1484 ENSURE(m.const_coeff(p, 0, 1, c) && c == 3);
1485 ENSURE(!m.const_coeff(p, 0, 0, c));
1486
1487 p = (x0^2)*x1 + 3*x0 + 1;
1488 ENSURE(!m.const_coeff(p, 0, 2, c));
1489 ENSURE(m.const_coeff(p, 0, 1, c) && c == 3);
1490 ENSURE(m.const_coeff(p, 0, 0, c) && c == 1);
1491
1492 p = x1 + 3*x0 + 1;
1493 ENSURE(m.const_coeff(p, 0, 2, c) && c == 0);
1494 ENSURE(m.const_coeff(p, 0, 1, c) && c == 3);
1495 ENSURE(!m.const_coeff(p, 0, 0, c));
1496
1497 p = 5*(x0^2) + 3*x0 + 7;
1498 ENSURE(m.const_coeff(p, 0, 5, c) && c == 0);
1499 ENSURE(m.const_coeff(p, 0, 2, c) && c == 5);
1500 ENSURE(m.const_coeff(p, 0, 1, c) && c == 3);
1501 ENSURE(m.const_coeff(p, 0, 0, c) && c == 7);
1502
1503 p = 5*(x0^2) + 3*x0;
1504 ENSURE(m.const_coeff(p, 0, 0, c) && c == 0);
1505
1506 p = - x0*x1 - x1 + 1;
1507 ENSURE(!m.const_coeff(p, 0, 1, c));
1508 }
1509
tst_gcd2()1510 static void tst_gcd2() {
1511 // enable_trace("mgcd");
1512 polynomial::numeral_manager nm;
1513 reslimit rl; polynomial::manager m(rl, nm);
1514 polynomial_ref x0(m);
1515 polynomial_ref x1(m);
1516 polynomial_ref x2(m);
1517 polynomial_ref x3(m);
1518 x0 = m.mk_polynomial(m.mk_var());
1519 x1 = m.mk_polynomial(m.mk_var());
1520 x2 = m.mk_polynomial(m.mk_var());
1521 x3 = m.mk_polynomial(m.mk_var());
1522 polynomial_ref p(m);
1523 polynomial_ref one(m);
1524 one = m.mk_const(mpz(1));
1525
1526 tst_gcd(one, one, one);
1527
1528 tst_gcd((5*x0 + 3*x1)*(x1 + 2)*(x2 + 3),
1529 (10*x0 + 6*x1)*(x1 + 2)*(x2 + 5),
1530 (5*x0 + 3*x1)*(x1 + 2));
1531
1532 tst_gcd((x0 + 1)*(x1 + 2)*(x2 + 3),
1533 (x0 + 1)*(x1 + 2)*(x2 + 5),
1534 (x0 + 1)*(x1 + 2));
1535
1536 tst_gcd(((x1^2) + 2*(x2^2) + 3*(x3^3))*(x1 + x2 + x3 + 1),
1537 ((x1^2) + 2*(x2^2) + 3*(x3^3))*(x1 + x2 + x3 + 2),
1538 ((x1^2) + 2*(x2^2) + 3*(x3^3)));
1539
1540 tst_gcd(5*(x1^3) + 11 + 7*(x0^2),
1541 5*(x1^3) + 13 + 7*(x0^2),
1542 one);
1543
1544 tst_gcd((5*3*(x1^2) + 5*6*(x2^2) + 5*21*(x3^3))*(5*(x1^3) + 7*(x0^2) + 11),
1545 (7*3*(x1^2) + 7*6*(x2^2) + 7*21*(x3^3))*(5*(x1^3) + 7*(x0^2) + 13),
1546 (3*(x1^2) + 6*(x2^2) + 21*(x3^3)));
1547
1548 tst_gcd((x2^6)*(x3^6) - 4*(x2^3)*(x3^6) + 2*(x2^6)*(x3^3) - 8*(x2^3)*(x3^3) + 4*(x1^3)*(x2^3)*(x3^3) - 8*(x1^3)*(x3^3) +
1549 4*(x3^6) + 8*(x3^3) + (x2^6) - 4*(x2^3) + 4*(x1^3)*(x2^3) - 8*(x1^3) + 4 + (x1^6),
1550 (-2)*(x2^3)*(x3^6) - 4*(x2^3)*(x3^3) + 4*(x3^6) + 8*(x3^3) - 2*(x1^3)*(x3^3) - 2*(x2^3) + 4 - 2*(x1^3),
1551 one);
1552
1553 tst_gcd((x1^2) - 2*x0 + 1 + (x0^2) + x0*x1 - 2*x1,
1554 x0*x1,
1555 one);
1556
1557 tst_gcd((5*3*(x1^2) + 5*6*(x2^2) + 5*21*(x3^3))*(x1 + x2 + x3 + 1)*(5*(x1^3) + 7*(x0^2) + 11),
1558 (7*3*(x1^2) + 7*6*(x2^2) + 7*21*(x3^3))*(x1 + x2 + x3 + 2)*(5*(x1^3) + 7*(x0^2) + 13),
1559 (3*(x1^2) + 6*(x2^2) + 21*(x3^3)));
1560
1561 p = 169*(x1^12)*(x2^16) - 468*x0*(x1^11)*(x2^16) + 428*(x0^2)*(x1^10)*(x2^16) - 92*(x0^3)*(x1^9)*(x2^16) - 82*(x0^4)*(x1^8)*(x2^16) + 52*(x0^5)*(x1^7)*(x2^16) - 4*(x0^6)*(x1^6)*(x2^16) - 4*(x0^7)*(x1^5)*(x2^16) + (x0^8)*(x1^4)*(x2^16) - 581*(x1^14)*(x2^14) + 1828*x0*(x1^13)*(x2^14) - 2452*(x0^2)*(x1^12)*(x2^14) + 548*(x0^3)*(x1^11)*(x2^14) + 1002*(x0^4)*(x1^10)*(x2^14) - 756*(x0^5)*(x1^9)*(x2^14) + 124*(x0^6)*(x1^8)*(x2^14) + 44*(x0^7)*(x1^7)*(x2^14) - 13*(x0^8)*(x1^6)*(x2^14) + 895*(x1^16)*(x2^12) - 1556*x0*(x1^15)*(x2^12) + 2864*(x0^2)*(x1^14)*(x2^12);
1562 tst_gcd(p, derivative(p, 2), (x1^4)*(x2^11));
1563
1564 tst_gcd((11*5*3)*((x0^2) + 1)*(x1 + 3),
1565 (11*5*7)*((x0^2) + 1)*(x1 + 5),
1566 (11*5)*((x0^2) + 1));
1567
1568 p = (x0^4)*(x3^8) - 2*(x0^4)*(x3^7) - 2*(x0^3)*(x2^3)*(x3^7) + 4*(x0^3)*(x3^7) + 2*(x0^4)*(x3^5) - 4*(x0^4)*(x3^4) - 4*(x0^3)*(x2^3)*(x3^4) + 8*(x0^3)*(x3^4) - 2*(x0^3)*(x1^3)*(x3^5) + 4*(x0^3)*(x1^3)*(x3^4) + 4*(x0^2)*(x1^3)*(x2^3)*(x3^4) - 8*(x0^2)*(x1^3)*(x3^4) + (x0^4)*(x3^6) + 2*(x0^3)*(x2^3)*(x3^6) - 4*(x0^3)*(x3^6) + 2*(x0^4)*(x3^3) + 4*(x0^3)*(x2^3)*(x3^3) - 8*(x0^3)*(x3^3) - 2*(x0^3)*(x1^3)*(x3^3) - 4*(x0^2)*(x1^3)*(x2^3)*(x3^3) + 8*(x0^2)*(x1^3)*(x3^3) + (x0^2)*(x2^6)*(x3^6) - 4*(x0^2)*(x2^3)*(x3^6) + 2*(x0^2)*(x2^6)*(x3^3) - 8*(x0^2)*(x2^3)*(x3^3) - 2*x0*(x1^3)*(x2^6)*(x3^3) + 8*x0*(x1^3)*(x2^3)*(x3^3) + 4*(x0^2)*(x3^6) + 8*(x0^2)*(x3^3) - 8*x0*(x1^3)*(x3^3) + (x0^4)*(x3^2) - 2*(x0^4)*x3 - 2*(x0^3)*(x2^3)*x3 + 4*(x0^3)*x3 - 2*(x0^3)*(x1^3)*(x3^2) + 4*(x0^3)*(x1^3)*x3 + 4*(x0^2)*(x1^3)*(x2^3)*x3 - 8*(x0^2)*(x1^3)*x3 + (x0^4) + 2*(x0^3)*(x2^3) - 4*(x0^3) - 2*(x0^3)*(x1^3) - 4*(x0^2)*(x1^3)*(x2^3) + 8*(x0^2)*(x1^3) + (x0^2)*(x2^6) - 4*(x0^2)*(x2^3) - 2*x0*(x1^3)*(x2^6) + 8*x0*(x1^3)*(x2^3) + 4*(x0^2) - 8*x0*(x1^3) + (x0^2)*(x1^6)*(x3^2) - 2*(x0^2)*(x1^6)*x3 - 2*x0*(x1^6)*(x2^3)*x3 + 4*x0*(x1^6)*x3 + (x0^2)*(x1^6) + 2*x0*(x1^6)*(x2^3) - 4*x0*(x1^6) + (x1^6)*(x2^6) - 4*(x1^6)*(x2^3) + 4*(x1^6);
1569
1570 // polynomial_ref p1(m);
1571 // p1 = derivative(p, 0);
1572 // polynomial_ref g(m);
1573 // for (unsigned i = 0; i < 50; i++)
1574 // g = gcd(p, p1);
1575 // return;
1576
1577 tst_gcd(p, derivative(p, 1),
1578 x0*(x2^6)*(x3^3) - (x1^3)*(x2^6) - 2*(x0^2)*(x2^3)*(x3^4) + 2*x0*(x1^3)*(x2^3)*x3 + 2*(x0^2)*(x2^3)*(x3^3) + (x0^3)*(x3^5) - 2*x0*(x1^3)*(x2^3) + x0*(x2^6) - (x0^2)*(x1^3)*(x3^2) - 4*x0*(x2^3)*(x3^3) - 2*(x0^3)*(x3^4) + 4*(x1^3)*(x2^3) + 2*(x0^2)*(x1^3)*x3 - 2*(x0^2)*(x2^3)*x3 + (x0^3)*(x3^3) + 4*(x0^2)*(x3^4) - (x0^2)*(x1^3) + 2*(x0^2)*(x2^3) - 4*x0*(x1^3)*x3 + (x0^3)*(x3^2) - 4*(x0^2)*(x3^3) + 4*x0*(x1^3) - 4*x0*(x2^3) - 2*(x0^3)*x3 + 4*x0*(x3^3) + (x0^3) - 4*(x1^3) + 4*(x0^2)*x3 - 4*(x0^2) + 4*x0
1579 );
1580
1581 tst_gcd(p, derivative(p, 0),
1582 neg((-1)*x0*(x2^3)*(x3^3) + (x1^3)*(x2^3) + (x0^2)*(x3^4) - x0*(x1^3)*x3 - (x0^2)*(x3^3) + x0*(x1^3) - x0*(x2^3) + 2*x0*(x3^3) - 2*(x1^3) + (x0^2)*x3 - (x0^2) + 2*x0));
1583
1584 tst_gcd(p, derivative(p, 2),
1585 neg((-1)*(x0^2)*(x2^3)*(x3^6) + 2*x0*(x1^3)*(x2^3)*(x3^3) + (x0^3)*(x3^7) - (x1^6)*(x2^3) - 2*(x0^2)*(x1^3)*(x3^4) - (x0^3)*(x3^6) + x0*(x1^6)*x3 + 2*(x0^2)*(x1^3)*(x3^3) - 2*(x0^2)*(x2^3)*(x3^3) + 2*(x0^2)*(x3^6) - x0*(x1^6) + 2*x0*(x1^3)*(x2^3) - 4*x0*(x1^3)*(x3^3) + 2*(x0^3)*(x3^4) + 2*(x1^6) - 2*(x0^2)*(x1^3)*x3 - 2*(x0^3)*(x3^3) + 2*(x0^2)*(x1^3) - (x0^2)*(x2^3) + 4*(x0^2)*(x3^3) - 4*x0*(x1^3) + (x0^3)*x3 - (x0^3) + 2*(x0^2))
1586 );
1587
1588 tst_gcd(((11*5*3)*(x0^2) + 1)*(x1 + 3),
1589 ((11*5*3)*(x0^2) + 1)*(x1 + 5),
1590 ((11*5*3)*(x0^2) + 1));
1591
1592 return;
1593 p = 169*(x1^12)*(x2^16) - 468*x0*(x1^11)*(x2^16) + 428*(x0^2)*(x1^10)*(x2^16) - 92*(x0^3)*(x1^9)*(x2^16) - 82*(x0^4)*(x1^8)*(x2^16) + 52*(x0^5)*(x1^7)*(x2^16) - 4*(x0^6)*(x1^6)*(x2^16) - 4*(x0^7)*(x1^5)*(x2^16) + (x0^8)*(x1^4)*(x2^16) - 581*(x1^14)*(x2^14) + 1828*x0*(x1^13)*(x2^14) - 2452*(x0^2)*(x1^12)*(x2^14) + 548*(x0^3)*(x1^11)*(x2^14) + 1002*(x0^4)*(x1^10)*(x2^14) - 756*(x0^5)*(x1^9)*(x2^14) + 124*(x0^6)*(x1^8)*(x2^14) + 44*(x0^7)*(x1^7)*(x2^14) - 13*(x0^8)*(x1^6)*(x2^14) + 895*(x1^16)*(x2^12) - 1556*x0*(x1^15)*(x2^12) + 2864*(x0^2)*(x1^14)*(x2^12) + 520*(x0^3)*(x1^13)*(x2^12) - 5402*(x0^4)*(x1^12)*(x2^12) + 3592*(x0^5)*(x1^11)*(x2^12) - 156*(x0^6)*(x1^10)*(x2^12) - 680*(x0^7)*(x1^9)*(x2^12) + 171*(x0^8)*(x1^8)*(x2^12) + 12*(x0^9)*(x1^7)*(x2^12) - 4*(x0^10)*(x1^6)*(x2^12) - 957*(x1^18)*(x2^10) - 1132*x0*(x1^17)*(x2^10) + 206*(x0^2)*(x1^16)*(x2^10) + 588*(x0^3)*(x1^15)*(x2^10) + 6861*(x0^4)*(x1^14)*(x2^10) - 5016*(x0^5)*(x1^13)*(x2^10) - 2756*(x0^6)*(x1^12)*(x2^10) + 3952*(x0^7)*(x1^11)*(x2^10) - 1143*(x0^8)*(x1^10)*(x2^10) - 124*(x0^9)*(x1^9)*(x2^10) + 30*(x0^10)*(x1^8)*(x2^10) + 4*(x0^11)*(x1^7)*(x2^10) - (x0^12)*(x1^6)*(x2^10) + 1404*(x1^20)*(x2^8) + 684*x0*(x1^19)*(x2^8) - 1224*(x0^2)*(x1^18)*(x2^8) - 4412*(x0^3)*(x1^17)*(x2^8) - 1442*(x0^4)*(x1^16)*(x2^8) + 4164*(x0^5)*(x1^15)*(x2^8) + 4116*(x0^6)*(x1^14)*(x2^8) - 5308*(x0^7)*(x1^13)*(x2^8) + 392*(x0^8)*(x1^12)*(x2^8) + 1600*(x0^9)*(x1^11)*(x2^8) - 468*(x0^10)*(x1^10)*(x2^8) - 24*(x0^11)*(x1^9)*(x2^8) + 6*(x0^12)*(x1^8)*(x2^8) - 594*(x1^22)*(x2^6) - 324*x0*(x1^21)*(x2^6) + 1980*(x0^2)*(x1^20)*(x2^6) + 1136*(x0^3)*(x1^19)*(x2^6) + 405*(x0^4)*(x1^18)*(x2^6) - 3916*(x0^5)*(x1^17)*(x2^6) - 396*(x0^6)*(x1^16)*(x2^6) + 1876*(x0^7)*(x1^15)*(x2^6) + 1108*(x0^8)*(x1^14)*(x2^6) - 2064*(x0^9)*(x1^13)*(x2^6) + 248*(x0^10)*(x1^12)*(x2^6) + 380*(x0^11)*(x1^11)*(x2^6) - 95*(x0^12)*(x1^10)*(x2^6) + 81*(x1^24)*(x2^4) + 108*x0*(x1^23)*(x2^4) - 432*(x0^2)*(x1^22)*(x2^4) - 276*(x0^3)*(x1^21)*(x2^4) + 481*(x0^4)*(x1^20)*(x2^4) + 144*(x0^5)*(x1^19)*(x2^4) + 788*(x0^6)*(x1^18)*(x2^4) - 1152*(x0^7)*(x1^17)*(x2^4) + 231*(x0^8)*(x1^16)*(x2^4) + 244*(x0^9)*(x1^15)*(x2^4) + 396*(x0^10)*(x1^14)*(x2^4) - 476*(x0^11)*(x1^13)*(x2^4) + 119*(x0^12)*(x1^12)*(x2^4) + 72*(x0^4)*(x1^22)*(x2^2) - 96*(x0^5)*(x1^21)*(x2^2) - 40*(x0^6)*(x1^20)*(x2^2) - 32*(x0^7)*(x1^19)*(x2^2) + 340*(x0^8)*(x1^18)*(x2^2) - 368*(x0^9)*(x1^17)*(x2^2) + 112*(x0^10)*(x1^16)*(x2^2) + 16*(x0^11)*(x1^15)*(x2^2) - 4*(x0^12)*(x1^14)*(x2^2) + 16*(x0^8)*(x1^20) - 64*(x0^9)*(x1^19) + 96*(x0^10)*(x1^18) - 64*(x0^11)*(x1^17) + 16*(x0^12)*(x1^16);
1594 polynomial_ref p_prime(m);
1595 p_prime = derivative(p, 2);
1596 tst_gcd(p, p_prime, x1^4);
1597 }
1598
1599 #if 0
1600 static void tst_gcd3() {
1601 enable_trace("polynomial_gcd");
1602 enable_trace("polynomial_gcd_detail");
1603 enable_trace("mpzzp");
1604 polynomial::numeral_manager nm;
1605 reslimit rl; polynomial::manager m(rl, nm);
1606 polynomial_ref x(m);
1607 x = m.mk_polynomial(m.mk_var());
1608 polynomial_ref p(m);
1609 polynomial_ref q(m);
1610 p = (x^8) + (x^6) + (x^4) + (x^3) + (x^2) + 1;
1611 q = (x^6) + (x^4) + x + 1;
1612 {
1613 polynomial::scoped_set_zp setZ2(m, 2);
1614 std::cout << "Z_p: " << nm.to_string(m.p()) << "\n";
1615 tst_gcd(normalize(p), normalize(q), x + 1);
1616 }
1617 {
1618 polynomial::scoped_set_zp setZ3(m, 3);
1619 std::cout << "Z_p: " << nm.to_string(m.p()) << "\n";
1620 polynomial_ref one(m);
1621 one = m.mk_const(mpz(1));
1622 tst_gcd(normalize(p), normalize(q), one);
1623 }
1624 }
1625
1626 static void tst_gcd4() {
1627 enable_trace("mgcd");
1628 // enable_trace("CRA");
1629 polynomial::numeral_manager nm;
1630 reslimit rl; polynomial::manager m(rl, nm);
1631 polynomial_ref x(m);
1632 x = m.mk_polynomial(m.mk_var());
1633 polynomial_ref p(m);
1634 polynomial_ref q(m);
1635 p = (15*x + 15)*((x + 2)^8)*(10000*x + 1)*(x + 3);
1636 q = (6*x + 6)*((x + 20)^8)*(10000*x + 3)*(x + 30);
1637 tst_gcd(p, q, 3*x + 3);
1638 p = (3*x + 2)*((x + 2)^8)*(10000*x + 1)*(x + 3);
1639 q = (3*x + 2)*((x + 20)^8)*(10000*x + 3)*(x + 30);
1640 tst_gcd(p, q, 3*x + 2);
1641 p = ((3*x + 2)*((x + 2)^8)*(10000*x + 1)*(x + 3))^3;
1642 q = ((3*x + 2)*((x + 20)^8)*(10000*x + 3)*(x + 30))^3;
1643 tst_gcd(p, q, (3*x + 2)^3);
1644 p = ((x + 3)^10)*((x^5) - x - 1)*(x + 1)*(x + 2)*(x + 4)*(10000*x + 33)*(x + 6)*(x + 11)*(x+33)*
1645 ((x^16) - 136*(x^14) + 6476*(x^12) - 141912*(x^10) + 1513334*(x^8) - 7453176*(x^6) + 13950764*(x^4) - 5596840*(x^2) + 46225)*
1646 (1000000*x + 1)*(333333333*x + 1)*(77777777*x + 1)*(11111111*x + 1)*(x + 128384747)*(x + 82837437)*(x + 22848481);
1647 tst_gcd(p, derivative(p, 0), (x + 3)^9);
1648 }
1649 #endif
1650
tst_newton_interpolation()1651 static void tst_newton_interpolation() {
1652 // enable_trace("newton");
1653 polynomial::numeral_manager nm;
1654 reslimit rl; polynomial::manager m(rl, nm);
1655 polynomial_ref x(m);
1656 polynomial_ref y(m);
1657 x = m.mk_polynomial(m.mk_var());
1658 y = m.mk_polynomial(m.mk_var());
1659 polynomial_ref p1(m), p2(m), p3(m);
1660 p1 = (-9)*y - 21;
1661 p2 = (-3)*y + 20;
1662 p3 = 5*y - 36;
1663 scoped_mpz_vector ins(nm);
1664 ins.push_back(mpz(0)); ins.push_back(mpz(1)); ins.push_back(mpz(2));
1665 polynomial::polynomial * outs[3] = { p1.get(), p2.get(), p3.get() };
1666 polynomial_ref r(m);
1667 {
1668 polynomial::scoped_set_zp setZ97(m, 97);
1669 m.newton_interpolation(0, 2, ins.data(), outs, r);
1670 }
1671 std::cout << "interpolation result: " << r << "\n";
1672 ENSURE(m.eq((x^2)*y + 5*x*y + 41*x - 9*y - 21, r));
1673 }
1674
tst_slow_mod_gcd()1675 static void tst_slow_mod_gcd() {
1676 polynomial::numeral_manager nm;
1677 reslimit rl; polynomial::manager m(rl, nm);
1678 polynomial_ref x0(m), x1(m), x2(m), x3(m), x4(m), x5(m);
1679 x0 = m.mk_polynomial(m.mk_var());
1680 x1 = m.mk_polynomial(m.mk_var());
1681 x2 = m.mk_polynomial(m.mk_var());
1682 x3 = m.mk_polynomial(m.mk_var());
1683 x4 = m.mk_polynomial(m.mk_var());
1684 x5 = m.mk_polynomial(m.mk_var());
1685 polynomial_ref p(m), q(m), b(m);
1686 polynomial_ref p_prime(m);
1687
1688 p = ((x0^3)*x1*x2*x3*x4*x5 + x0*x1*x2*x3*(x4^3)*x5 + x0*x1*x2*x3*x4*(x5^3) - x0*x1*x2*x3*x4*x5 - 2)^2;
1689 q = derivative(p, 0);
1690 tst_gcd(p, q, (x0^3)*x1*x2*x3*x4*x5 + x0*x1*x2*x3*(x4^3)*x5 + x0*x1*x2*x3*x4*(x5^3) - x0*x1*x2*x3*x4*x5 - 2);
1691
1692 b = (x0^10) + (x1^10) + (x2^10) + (x3^10);
1693 p = b*(x0 + 1);
1694 q = b*(x0 + 2);
1695 tst_gcd(p, q, b);
1696
1697 return;
1698 p = (x0^8) *
1699 (((x0^3)*x1*x2*x3*x4*x5 + x0*(x1^3)*x2*x3*x4*x5 + x0*x1*(x2^3)*x3*x4*x5 + x0*x1*x2*(x3^3)*x4*x5 +
1700 x0*x1*x2*x3*(x4^3)*x5 + x0*x1*x2*x3*x4*(x5^3) - x0*x1*x2*x3*x4*x5 - 2)^2) *
1701 (((x0^3)*x1*x2*x3*x4*x5 + x0*(x1^3)*x2*x3*x4*x5 + x0*x1*(x2^3)*x3*x4*x5 + x0*x1*x2*(x3^3)*x4*x5 +
1702 x0*x1*x2*x3*(x4^3)*x5 + x0*x1*x2*x3*x4*(x5^3) - x0*x1*x2*x3*x4*x5 + 2)^2);
1703 p_prime = derivative(p, 0);
1704 tst_gcd(p, p_prime,
1705 (x0^7) *
1706 ((x0^3)*x1*x2*x3*x4*x5 + x0*(x1^3)*x2*x3*x4*x5 + x0*x1*(x2^3)*x3*x4*x5 + x0*x1*x2*(x3^3)*x4*x5 +
1707 x0*x1*x2*x3*(x4^3)*x5 + x0*x1*x2*x3*x4*(x5^3) - x0*x1*x2*x3*x4*x5 - 2) *
1708 ((x0^3)*x1*x2*x3*x4*x5 + x0*(x1^3)*x2*x3*x4*x5 + x0*x1*(x2^3)*x3*x4*x5 + x0*x1*x2*(x3^3)*x4*x5 +
1709 x0*x1*x2*x3*(x4^3)*x5 + x0*x1*x2*x3*x4*(x5^3) - x0*x1*x2*x3*x4*x5 + 2));
1710 }
1711
tst_linear_solver()1712 void tst_linear_solver() {
1713 unsynch_mpq_manager qm;
1714 scoped_mpq_vector as(qm);
1715 scoped_mpq b(qm);
1716 scoped_mpq_vector xs(qm);
1717 linear_eq_solver<unsynch_mpq_manager> solver(qm);
1718
1719 solver.resize(3);
1720 xs.resize(3);
1721
1722 as.reset();
1723 as.push_back(mpq(2)); as.push_back(mpq(1)); as.push_back(mpq(-1)); qm.set(b, 8);
1724 solver.add(0, as.data(), b);
1725
1726 as.reset();
1727 as.push_back(mpq(-3)); as.push_back(mpq(-1)); as.push_back(mpq(2)); qm.set(b, -11);
1728 solver.add(1, as.data(), b);
1729
1730 as.reset();
1731 as.push_back(mpq(-2)); as.push_back(mpq(1)); as.push_back(mpq(2)); qm.set(b, -3);
1732 solver.add(2, as.data(), b);
1733
1734 VERIFY(solver.solve(xs.data()));
1735 ENSURE(qm.eq(xs[0], mpq(2)));
1736 ENSURE(qm.eq(xs[1], mpq(3)));
1737 ENSURE(qm.eq(xs[2], mpq(-1)));
1738 }
1739
tst_lex(polynomial_ref const & p1,polynomial_ref const & p2,int lex_expected,polynomial::var min,int lex2_expected)1740 static void tst_lex(polynomial_ref const & p1, polynomial_ref const & p2, int lex_expected, polynomial::var min, int lex2_expected) {
1741 polynomial::manager & m = p1.m();
1742 std::cout << "compare ";
1743 m.display(std::cout, m.get_monomial(p1, 0));
1744 std::cout << " ";
1745 m.display(std::cout, m.get_monomial(p2, 0));
1746 std::cout << " "; std::cout.flush();
1747 int r1 = lex_compare(m.get_monomial(p1, 0), m.get_monomial(p2, 0));
1748 int r2 = lex_compare2(m.get_monomial(p1, 0), m.get_monomial(p2, 0), min);
1749 ENSURE(r1 == lex_expected);
1750 ENSURE(r2 == lex2_expected);
1751 std::cout << r1 << " " << r2 << "\n";
1752 ENSURE(lex_compare(m.get_monomial(p2, 0), m.get_monomial(p1, 0)) == -r1);
1753 ENSURE(lex_compare2(m.get_monomial(p2, 0), m.get_monomial(p1, 0), min) == -r2);
1754 }
1755
tst_lex()1756 static void tst_lex() {
1757 polynomial::numeral_manager nm;
1758 reslimit rl; polynomial::manager m(rl, nm);
1759 polynomial_ref x0(m), x1(m), x2(m), x3(m), x4(m), x5(m);
1760 x0 = m.mk_polynomial(m.mk_var());
1761 x1 = m.mk_polynomial(m.mk_var());
1762 x2 = m.mk_polynomial(m.mk_var());
1763 x3 = m.mk_polynomial(m.mk_var());
1764 x4 = m.mk_polynomial(m.mk_var());
1765 x5 = m.mk_polynomial(m.mk_var());
1766
1767 polynomial_ref one(m);
1768 one = m.mk_const(mpz(1));
1769
1770 tst_lex(x0*x4*x1, (x0^10)*(x1^3), 1, 4, -1);
1771 tst_lex(x0*x3*(x1^2)*x4, x0*(x3^2)*(x1^2)*x4, -1, 3, -1);
1772 tst_lex((x0^2)*x3*(x1^2)*x4, x0*(x3^2)*(x1^2)*x4, -1, 3, 1);
1773 tst_lex(x0*x3*(x1^2)*x4, x0*x3*(x1^2)*x4, 0, 3, 0);
1774 tst_lex(x0*(x3^2)*(x1^2)*x4, x0*x3*(x1^2)*x4, 1, 3, 1);
1775 tst_lex((x1^2)*x4, x0*x2*x3*x4*x5, -1, 1, -1);
1776 tst_lex((x1^2)*x3*x4, x0*x1, 1, 1, 1);
1777 tst_lex(x1*x3*x4, x2*x3*x4, -1, 2, 1);
1778 tst_lex(x1*x3*x4, x2*x3*x4, -1, 1, -1);
1779 tst_lex(x1*x3*x4, x0*x2*x3*x4, -1, 1, -1);
1780 tst_lex(x3, x4, -1, 1, -1);
1781 tst_lex(x3, x4, -1, 4, 1);
1782 tst_lex(x2*x3, x4, -1, 1, -1);
1783 tst_lex(x2*x3, x4, -1, 4, 1);
1784 tst_lex(x3, x2*x4, -1, 1, -1);
1785 tst_lex(x3, x2*x4, -1, 4, 1);
1786 tst_lex(x3, x2*x4, -1, 4, 1);
1787 tst_lex(one, x3, -1, 1, -1);
1788 tst_lex(one, x3, -1, 3, -1);
1789 tst_lex(x3, one, 1, 3, 1);
1790 tst_lex(x4*x5, (x4^3)*x5, -1, 4, -1);
1791 }
1792
tst_divides()1793 static void tst_divides() {
1794 polynomial::numeral_manager nm;
1795 reslimit rl; polynomial::manager m(rl, nm);
1796 polynomial_ref x0(m);
1797 x0 = m.mk_polynomial(m.mk_var());
1798 polynomial_ref q(m);
1799 polynomial_ref p(m);
1800
1801 q = 16*(x0^27) - 1984*(x0^26) + 1762*(x0^25) + 17351*(x0^24) - 14165*(x0^23) + 16460*(x0^22) + 2919*(x0^21) - 16823*(x0^20) + 1530*(x0^19) +
1802 10646*(x0^18) + 19217*(x0^17);
1803 p = 16*(x0^39) - 3648*(x0^38) + 338136*(x0^37) - 16037936*(x0^36) + 392334357*(x0^35) - rational("3851617443")*(x0^34) -
1804 rational("14636221526")*(x0^33) + rational("377151717618")*(x0^32) + rational("677140776981")*(x0^31) - rational("4308280094419")*(x0^30) +
1805 rational("312708087606")*(x0^29) + rational("8205543533730")*(x0^28) + rational("3331586202704")*(x0^27) - rational("15291636627072")*(x0^26) +
1806 rational("433482645282")*(x0^25) + rational("7397104817486")*(x0^24) + rational("1021197979053")*(x0^23) - rational("1373737505247")*(x0^22) -
1807 rational("639394669026")*(x0^21) - rational("118513560618")*(x0^20) - rational("10405319535")*(x0^19) - rational("358722675")*(x0^18);
1808 std::cout << "----------------------\n";
1809 std::cout << "q: " << q << "\n";
1810 std::cout << "p: " << p << std::endl;
1811 std::cout << "divides(q, p): " << m.divides(q, p) << "\n";
1812 }
1813
tst_polynomial()1814 void tst_polynomial() {
1815 set_verbosity_level(1000);
1816 // enable_trace("factor");
1817 // enable_trace("poly_bug");
1818 // enable_trace("factor_bug");
1819 disable_trace("polynomial");
1820 enable_trace("psc_chain_classic");
1821 enable_trace("Lazard");
1822 // enable_trace("eval_bug");
1823 // enable_trace("mgcd");
1824 tst_psc();
1825 return;
1826 tst_eval();
1827 tst_divides();
1828 tst_gcd2();
1829 tst_slow_mod_gcd();
1830 tst_gcd();
1831
1832 tst_lex();
1833 tst_linear_solver();
1834 tst_newton_interpolation();
1835 tst_resultant();
1836 //
1837 // tst_gcd4();
1838 // tst_gcd3();
1839 tst_zp();
1840 tst_const_coeff();
1841 tst_psc();
1842 tst_del_eh();
1843 tst_mk_unique();
1844 tst_qsubstitute();
1845 tst_substitute();
1846 tst_discriminant();
1847 tst_mfact();
1848 tst_mm();
1849 // tst_p25();
1850 // return;
1851 tst_translate();
1852 // enable_trace("mpz_gcd");
1853 tst_vars();
1854 tst_sqf();
1855 enable_trace("resultant");
1856 enable_trace("psc");
1857 disable_trace("polynomial");
1858 enable_trace("pseudo_remainder");
1859 enable_trace("resultant_bug");
1860 tst_sqrt();
1861 tst_prem();
1862 tst_compose();
1863 tst11();
1864 tst10();
1865 tst9();
1866 tst8();
1867 tst7();
1868 tst6();
1869 tst5();
1870 tst3();
1871 tst2();
1872 tst1();
1873 tst4();
1874 }
1875 #else
tst_polynomial()1876 void tst_polynomial() {
1877 // it takes forever to compiler these regressions using clang++
1878 }
1879 #endif
1880