3 /* Some test with polynomial GCD calculations. See also the checks for
4 * rational function normalization in normalization.cpp. */
8 const int MAX_VARIABLES = 5;
10 static symbol x("x"), z("z");
11 static symbol y[MAX_VARIABLES];
14 static unsigned poly_gcd1(void)
16 for (int v=1; v<=MAX_VARIABLES; v++) {
19 for (int i=0; i<v; i++) {
24 ex f = (e1 + 1) * (e1 + 2);
25 ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1);
28 clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl;
35 // Linearly dense quartic inputs with quadratic GCDs
36 static unsigned poly_gcd2(void)
38 for (int v=1; v<=MAX_VARIABLES; v++) {
41 for (int i=0; i<v; i++) {
46 ex d = pow(e1 + 1, 2);
47 ex f = d * pow(e2 - 2, 2);
48 ex g = d * pow(e1 + 2, 2);
53 if ((r - d).expand().compare(exZERO()) != 0) {
54 clog << "case 2, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
61 // Sparse GCD and inputs where degrees are proportional to the number of variables
62 static unsigned poly_gcd3(void)
64 for (int v=1; v<=MAX_VARIABLES; v++) {
65 ex e1 = pow(x, v + 1);
66 for (int i=0; i<v; i++)
67 e1 += pow(y[i], v + 1);
73 if ((r - d).expand().compare(exZERO()) != 0) {
74 clog << "case 3, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
81 // Variation of case 3; major performance degradation with PRS
82 static unsigned poly_gcd3p(void)
84 for (int v=1; v<=MAX_VARIABLES; v++) {
85 ex e1 = pow(x, v + 1);
87 for (int i=0; i<v; i++) {
88 e1 += pow(y[i], v + 1);
96 if ((r - d).expand().compare(exZERO()) != 0) {
97 clog << "case 3p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
104 // Quadratic non-monic GCD; f and g have other quadratic factors
105 static unsigned poly_gcd4(void)
107 for (int v=1; v<=MAX_VARIABLES; v++) {
108 ex e1 = pow(x, 2) * pow(y[0], 2);
109 ex e2 = pow(x, 2) - pow(y[0], 2);
111 for (int i=1; i<v; i++) {
119 ex g = d * pow(e3 + 2, 2);
121 if ((r - d).expand().compare(exZERO()) != 0) {
122 clog << "case 4, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
129 // Completely dense non-monic quadratic inputs with dense non-monic linear GCDs
130 static unsigned poly_gcd5(void)
132 for (int v=1; v<=MAX_VARIABLES; v++) {
136 for (int i=0; i<v; i++) {
146 if ((r - d).expand().compare(exZERO()) != 0) {
147 clog << "case 5, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
154 // Sparse non-monic quadratic inputs with linear GCDs
155 static unsigned poly_gcd5p(void)
157 for (int v=1; v<=MAX_VARIABLES; v++) {
159 for (int i=0; i<v; i++)
166 if ((r - d).expand().compare(exZERO()) != 0) {
167 clog << "case 5p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
174 // Trivariate inputs with increasing degrees
175 static unsigned poly_gcd6(void)
179 for (int j=1; j<=MAX_VARIABLES; j++) {
180 ex d = pow(x, j) * y * (z - 1);
181 ex f = d * (pow(x, j) + pow(y, j + 1) * pow(z, j) + 1);
182 ex g = d * (pow(x, j + 1) + pow(y, j) * pow(z, j + 1) - 7);
184 if ((r - d).expand().compare(exZERO()) != 0) {
185 clog << "case 6, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
192 // Trivariate polynomials whose GCD has common factors with its cofactors
193 static unsigned poly_gcd7(void)
196 ex p = x - y * z + 1;
197 ex q = x - y + z * 3;
199 for (int j=1; j<=3; j++) {
200 for (int k=j+1; k<=4; k++) {
201 ex d = pow(p, j) * pow(q, j);
202 ex f = pow(p, j) * pow(q, k);
203 ex g = pow(p, k) * pow(q, j);
205 if ((r - d).expand().compare(exZERO()) != 0 && (r + d).expand().compare(exZERO()) != 0) {
206 clog << "case 7, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
214 unsigned poly_gcd(void)
218 cout << "checking polynomial GCD computation..." << flush;
219 clog << "---------polynomial GCD computation:" << endl;
221 result += poly_gcd1();
222 result += poly_gcd2();
223 result += poly_gcd3();
224 // result += poly_gcd3p(); // takes extremely long (PRS "worst" case)
225 result += poly_gcd4();
226 result += poly_gcd5();
227 result += poly_gcd5p();
228 result += poly_gcd6();
229 result += poly_gcd7();
233 clog << "(no output)" << endl;