1448ae3cc815dd0a9c06ac6f1353b7cc4d339c79
[ginac.git] / check / poly_gcd.cpp
1 /** @file 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 /*
7  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
8  *
9  *  This program is free software; you can redistribute it and/or modify
10  *  it under the terms of the GNU General Public License as published by
11  *  the Free Software Foundation; either version 2 of the License, or
12  *  (at your option) any later version.
13  *
14  *  This program is distributed in the hope that it will be useful,
15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  *  GNU General Public License for more details.
18  *
19  *  You should have received a copy of the GNU General Public License
20  *  along with this program; if not, write to the Free Software
21  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <ginac/ginac.h>
25
26 #ifndef NO_GINAC_NAMESPACE
27 using namespace GiNaC;
28 #endif // ndef NO_GINAC_NAMESPACE
29
30 const int MAX_VARIABLES = 5;
31
32 static symbol x("x"), z("z");
33 static symbol y[MAX_VARIABLES];
34
35 // GCD = 1
36 static unsigned poly_gcd1(void)
37 {
38         for (int v=1; v<=MAX_VARIABLES; v++) {
39                 ex e1 = x;
40                 ex e2 = pow(x, 2);
41                 for (int i=0; i<v; i++) {
42                         e1 += y[i];
43                         e2 += pow(y[i], 2);
44                 }
45
46                 ex f = (e1 + 1) * (e1 + 2);
47                 ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1);
48                 ex r = gcd(f, g);
49                 if (r != exONE()) {
50                         clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl;
51                         return 1;
52                 }
53         }
54         return 0;
55 }
56
57 // Linearly dense quartic inputs with quadratic GCDs
58 static unsigned poly_gcd2(void)
59 {
60         for (int v=1; v<=MAX_VARIABLES; v++) {
61                 ex e1 = x;
62                 ex e2 = x;
63                 for (int i=0; i<v; i++) {
64                         e1 += y[i];
65                         e2 -= y[i];
66                 }
67
68                 ex d = pow(e1 + 1, 2);
69                 ex f = d * pow(e2 - 2, 2);
70                 ex g = d * pow(e1 + 2, 2);
71                 ex r = gcd(f, g);
72         ex re=r.expand();
73         ex df1=r-d;
74         ex df=(r-d).expand();
75                 if ((r - d).expand().compare(exZERO()) != 0) {
76                         clog << "case 2, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
77                         return 1;
78                 }
79         }
80         return 0;
81 }
82
83 // Sparse GCD and inputs where degrees are proportional to the number of variables
84 static unsigned poly_gcd3(void)
85 {
86         for (int v=1; v<=MAX_VARIABLES; v++) {
87                 ex e1 = pow(x, v + 1);
88                 for (int i=0; i<v; i++)
89                         e1 += pow(y[i], v + 1);
90
91                 ex d = e1 + 1;
92                 ex f = d * (e1 - 2);
93                 ex g = d * (e1 + 2);
94                 ex r = gcd(f, g);
95                 if ((r - d).expand().compare(exZERO()) != 0) {
96                         clog << "case 3, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
97                         return 1;
98                 }
99         }
100         return 0;
101 }
102
103 // Variation of case 3; major performance degradation with PRS
104 static unsigned poly_gcd3p(void)
105 {
106         for (int v=1; v<=MAX_VARIABLES; v++) {
107                 ex e1 = pow(x, v + 1);
108                 ex e2 = pow(x, v);
109                 for (int i=0; i<v; i++) {
110                         e1 += pow(y[i], v + 1);
111                         e2 += pow(y[i], v);
112                 }
113
114                 ex d = e1 + 1;
115                 ex f = d * (e1 - 2);
116                 ex g = d * (e2 + 2);
117                 ex r = gcd(f, g);
118                 if ((r - d).expand().compare(exZERO()) != 0) {
119                         clog << "case 3p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
120                         return 1;
121                 }
122         }
123         return 0;
124 }
125
126 // Quadratic non-monic GCD; f and g have other quadratic factors
127 static unsigned poly_gcd4(void)
128 {
129         for (int v=1; v<=MAX_VARIABLES; v++) {
130                 ex e1 = pow(x, 2) * pow(y[0], 2);
131                 ex e2 = pow(x, 2) - pow(y[0], 2);
132                 ex e3 = x * y[0];
133                 for (int i=1; i<v; i++) {
134                         e1 += pow(y[i], 2);
135                         e2 += pow(y[i], 2);
136                         e3 += y[i];
137                 }
138
139                 ex d = e1 + 1;
140                 ex f = d * (e2 - 1);
141                 ex g = d * pow(e3 + 2, 2);
142                 ex r = gcd(f, g);
143                 if ((r - d).expand().compare(exZERO()) != 0) {
144                         clog << "case 4, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
145                         return 1;
146                 }
147         }
148         return 0;
149 }
150
151 // Completely dense non-monic quadratic inputs with dense non-monic linear GCDs
152 static unsigned poly_gcd5(void)
153 {
154         for (int v=1; v<=MAX_VARIABLES; v++) {
155                 ex e1 = x + 1;
156                 ex e2 = x - 2;
157                 ex e3 = x + 2;
158                 for (int i=0; i<v; i++) {
159                         e1 *= y[i] + 1;
160                         e2 *= y[i] - 2;
161                         e3 *= y[i] + 2;
162                 }
163
164                 ex d = e1 - 3;
165                 ex f = d * (e2 + 3);
166                 ex g = d * (e3 - 3);
167                 ex r = gcd(f, g);
168                 if ((r - d).expand().compare(exZERO()) != 0) {
169                         clog << "case 5, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
170                         return 1;
171                 }
172         }
173         return 0;
174 }
175
176 // Sparse non-monic quadratic inputs with linear GCDs
177 static unsigned poly_gcd5p(void)
178 {
179         for (int v=1; v<=MAX_VARIABLES; v++) {
180                 ex e1 = x;
181                 for (int i=0; i<v; i++)
182                         e1 *= y[i];
183
184                 ex d = e1 - 1;
185                 ex f = d * (e1 + 3);
186                 ex g = d * (e1 - 3);
187                 ex r = gcd(f, g);
188                 if ((r - d).expand().compare(exZERO()) != 0) {
189                         clog << "case 5p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
190                         return 1;
191                 }
192         }
193         return 0;
194 }
195
196 // Trivariate inputs with increasing degrees
197 static unsigned poly_gcd6(void)
198 {
199         symbol y("y");
200
201         for (int j=1; j<=MAX_VARIABLES; j++) {
202                 ex d = pow(x, j) * y * (z - 1);
203                 ex f = d * (pow(x, j) + pow(y, j + 1) * pow(z, j) + 1);
204                 ex g = d * (pow(x, j + 1) + pow(y, j) * pow(z, j + 1) - 7);
205                 ex r = gcd(f, g);
206                 if ((r - d).expand().compare(exZERO()) != 0) {
207                         clog << "case 6, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
208                         return 1;
209                 }
210         }
211         return 0;
212 }
213
214 // Trivariate polynomials whose GCD has common factors with its cofactors
215 static unsigned poly_gcd7(void)
216 {
217         symbol y("y");
218         ex p = x - y * z + 1;
219         ex q = x - y + z * 3;
220
221         for (int j=1; j<=3; j++) {
222                 for (int k=j+1; k<=4; k++) {
223                         ex d = pow(p, j) * pow(q, j);
224                         ex f = pow(p, j) * pow(q, k);
225                         ex g = pow(p, k) * pow(q, j); 
226                         ex r = gcd(f, g);
227                         if ((r - d).expand().compare(exZERO()) != 0 && (r + d).expand().compare(exZERO()) != 0) {
228                                 clog << "case 7, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
229                                 return 1;
230                         }
231                 }
232         }
233         return 0;
234 }
235
236 unsigned poly_gcd(void)
237 {
238     unsigned result = 0;
239
240         cout << "checking polynomial GCD computation..." << flush;
241         clog << "---------polynomial GCD computation:" << endl;
242
243         result += poly_gcd1();
244         result += poly_gcd2();
245         result += poly_gcd3();
246 //      result += poly_gcd3p(); // takes extremely long (PRS "worst" case)
247         result += poly_gcd4();
248         result += poly_gcd5();
249         result += poly_gcd5p();
250         result += poly_gcd6();
251         result += poly_gcd7();
252
253         if (!result) {
254                 cout << " passed ";
255         clog << "(no output)" << endl;
256     } else {
257                 cout << " failed ";
258     }
259         return result;
260 }