Initial revision
[ginac.git] / check / poly_gcd.cpp
1 // check/poly_gcd.cpp
2
3 /* Some test with polynomial GCD calculations. See also the checks for
4  * rational function normalization in normalization.cpp. */
5
6 #include "ginac.h"
7
8 const int MAX_VARIABLES = 5;
9
10 static symbol x("x"), z("z");
11 static symbol y[MAX_VARIABLES];
12
13 // GCD = 1
14 static unsigned poly_gcd1(void)
15 {
16         for (int v=1; v<=MAX_VARIABLES; v++) {
17                 ex e1 = x;
18                 ex e2 = pow(x, 2);
19                 for (int i=0; i<v; i++) {
20                         e1 += y[i];
21                         e2 += pow(y[i], 2);
22                 }
23
24                 ex f = (e1 + 1) * (e1 + 2);
25                 ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1);
26                 ex r = gcd(f, g);
27                 if (r != exONE()) {
28                         clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl;
29                         return 1;
30                 }
31         }
32         return 0;
33 }
34
35 // Linearly dense quartic inputs with quadratic GCDs
36 static unsigned poly_gcd2(void)
37 {
38         for (int v=1; v<=MAX_VARIABLES; v++) {
39                 ex e1 = x;
40                 ex e2 = x;
41                 for (int i=0; i<v; i++) {
42                         e1 += y[i];
43                         e2 -= y[i];
44                 }
45
46                 ex d = pow(e1 + 1, 2);
47                 ex f = d * pow(e2 - 2, 2);
48                 ex g = d * pow(e1 + 2, 2);
49                 ex r = gcd(f, g);
50         ex re=r.expand();
51         ex df1=r-d;
52         ex df=(r-d).expand();
53                 if ((r - d).expand().compare(exZERO()) != 0) {
54                         clog << "case 2, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
55                         return 1;
56                 }
57         }
58         return 0;
59 }
60
61 // Sparse GCD and inputs where degrees are proportional to the number of variables
62 static unsigned poly_gcd3(void)
63 {
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);
68
69                 ex d = e1 + 1;
70                 ex f = d * (e1 - 2);
71                 ex g = d * (e1 + 2);
72                 ex r = gcd(f, g);
73                 if ((r - d).expand().compare(exZERO()) != 0) {
74                         clog << "case 3, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
75                         return 1;
76                 }
77         }
78         return 0;
79 }
80
81 // Variation of case 3; major performance degradation with PRS
82 static unsigned poly_gcd3p(void)
83 {
84         for (int v=1; v<=MAX_VARIABLES; v++) {
85                 ex e1 = pow(x, v + 1);
86                 ex e2 = pow(x, v);
87                 for (int i=0; i<v; i++) {
88                         e1 += pow(y[i], v + 1);
89                         e2 += pow(y[i], v);
90                 }
91
92                 ex d = e1 + 1;
93                 ex f = d * (e1 - 2);
94                 ex g = d * (e2 + 2);
95                 ex r = gcd(f, g);
96                 if ((r - d).expand().compare(exZERO()) != 0) {
97                         clog << "case 3p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
98                         return 1;
99                 }
100         }
101         return 0;
102 }
103
104 // Quadratic non-monic GCD; f and g have other quadratic factors
105 static unsigned poly_gcd4(void)
106 {
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);
110                 ex e3 = x * y[0];
111                 for (int i=1; i<v; i++) {
112                         e1 += pow(y[i], 2);
113                         e2 += pow(y[i], 2);
114                         e3 += y[i];
115                 }
116
117                 ex d = e1 + 1;
118                 ex f = d * (e2 - 1);
119                 ex g = d * pow(e3 + 2, 2);
120                 ex r = gcd(f, g);
121                 if ((r - d).expand().compare(exZERO()) != 0) {
122                         clog << "case 4, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
123                         return 1;
124                 }
125         }
126         return 0;
127 }
128
129 // Completely dense non-monic quadratic inputs with dense non-monic linear GCDs
130 static unsigned poly_gcd5(void)
131 {
132         for (int v=1; v<=MAX_VARIABLES; v++) {
133                 ex e1 = x + 1;
134                 ex e2 = x - 2;
135                 ex e3 = x + 2;
136                 for (int i=0; i<v; i++) {
137                         e1 *= y[i] + 1;
138                         e2 *= y[i] - 2;
139                         e3 *= y[i] + 2;
140                 }
141
142                 ex d = e1 - 3;
143                 ex f = d * (e2 + 3);
144                 ex g = d * (e3 - 3);
145                 ex r = gcd(f, g);
146                 if ((r - d).expand().compare(exZERO()) != 0) {
147                         clog << "case 5, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
148                         return 1;
149                 }
150         }
151         return 0;
152 }
153
154 // Sparse non-monic quadratic inputs with linear GCDs
155 static unsigned poly_gcd5p(void)
156 {
157         for (int v=1; v<=MAX_VARIABLES; v++) {
158                 ex e1 = x;
159                 for (int i=0; i<v; i++)
160                         e1 *= y[i];
161
162                 ex d = e1 - 1;
163                 ex f = d * (e1 + 3);
164                 ex g = d * (e1 - 3);
165                 ex r = gcd(f, g);
166                 if ((r - d).expand().compare(exZERO()) != 0) {
167                         clog << "case 5p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
168                         return 1;
169                 }
170         }
171         return 0;
172 }
173
174 // Trivariate inputs with increasing degrees
175 static unsigned poly_gcd6(void)
176 {
177         symbol y("y");
178
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);
183                 ex r = gcd(f, g);
184                 if ((r - d).expand().compare(exZERO()) != 0) {
185                         clog << "case 6, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
186                         return 1;
187                 }
188         }
189         return 0;
190 }
191
192 // Trivariate polynomials whose GCD has common factors with its cofactors
193 static unsigned poly_gcd7(void)
194 {
195         symbol y("y");
196         ex p = x - y * z + 1;
197         ex q = x - y + z * 3;
198
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); 
204                         ex r = gcd(f, g);
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;
207                                 return 1;
208                         }
209                 }
210         }
211         return 0;
212 }
213
214 unsigned poly_gcd(void)
215 {
216     unsigned result = 0;
217
218         cout << "checking polynomial GCD computation..." << flush;
219         clog << "---------polynomial GCD computation:" << endl;
220
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();
230
231         if (!result) {
232                 cout << " passed ";
233         clog << "(no output)" << endl;
234     } else {
235                 cout << " failed ";
236     }
237         return result;
238 }