lexer: when switching to another output stream, clean last read character.
[ginac.git] / check / exam_polygcd.cpp
1 /** @file exam_polygcd.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-2008 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  */
23
24 #include <iostream>
25 #include "ginac.h"
26 using namespace std;
27 using namespace GiNaC;
28
29 const int MAX_VARIABLES = 3;
30
31 static symbol x("x"), z("z");
32 static symbol y[MAX_VARIABLES];
33
34 // GCD = 1
35 static unsigned poly_gcd1()
36 {
37         for (int v=1; v<=MAX_VARIABLES; v++) {
38                 ex e1 = x;
39                 ex e2 = pow(x, 2);
40                 for (int i=0; i<v; i++) {
41                         e1 += y[i];
42                         e2 += pow(y[i], 2);
43                 }
44
45                 ex f = (e1 + 1) * (e1 + 2);
46                 ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1);
47                 ex r = gcd(f, g);
48                 if (r != 1) {
49                         clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl;
50                         return 1;
51                 }
52         }
53         return 0;
54 }
55
56 // Linearly dense quartic inputs with quadratic GCDs
57 static unsigned poly_gcd2()
58 {
59         for (int v=1; v<=MAX_VARIABLES; v++) {
60                 ex e1 = x;
61                 ex e2 = x;
62                 for (int i=0; i<v; i++) {
63                         e1 += y[i];
64                         e2 -= y[i];
65                 }
66
67                 ex d = pow(e1 + 1, 2);
68                 ex f = d * pow(e2 - 2, 2);
69                 ex g = d * pow(e1 + 2, 2);
70                 ex r = gcd(f.expand(), g.expand());
71                 if (!(r - d).expand().is_zero()) {
72                         clog << "case 2, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
73                         return 1;
74                 }
75         }
76         return 0;
77 }
78
79 // Sparse GCD and inputs where degrees are proportional to the number of variables
80 static unsigned poly_gcd3()
81 {
82         for (int v=1; v<=MAX_VARIABLES; v++) {
83                 ex e1 = pow(x, v + 1);
84                 for (int i=0; i<v; i++)
85                         e1 += pow(y[i], v + 1);
86
87                 ex d = e1 + 1;
88                 ex f = d * (e1 - 2);
89                 ex g = d * (e1 + 2);
90                 ex r = gcd(f.expand(), g.expand());
91                 if (!(r - d).expand().is_zero()) {
92                         clog << "case 3, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
93                         return 1;
94                 }
95         }
96         return 0;
97 }
98
99 // Variation of case 3; major performance degradation with PRS
100 static unsigned poly_gcd3p()
101 {
102         for (int v=1; v<=MAX_VARIABLES; v++) {
103                 ex e1 = pow(x, v + 1);
104                 ex e2 = pow(x, v);
105                 for (int i=0; i<v; i++) {
106                         e1 += pow(y[i], v + 1);
107                         e2 += pow(y[i], v);
108                 }
109
110                 ex d = e1 + 1;
111                 ex f = d * (e1 - 2);
112                 ex g = d * (e2 + 2);
113                 ex r = gcd(f.expand(), g.expand());
114                 if (!(r - d).expand().is_zero()) {
115                         clog << "case 3p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
116                         return 1;
117                 }
118         }
119         return 0;
120 }
121
122 // Quadratic non-monic GCD; f and g have other quadratic factors
123 static unsigned poly_gcd4()
124 {
125         for (int v=1; v<=MAX_VARIABLES; v++) {
126                 ex e1 = pow(x, 2) * pow(y[0], 2);
127                 ex e2 = pow(x, 2) - pow(y[0], 2);
128                 ex e3 = x * y[0];
129                 for (int i=1; i<v; i++) {
130                         e1 += pow(y[i], 2);
131                         e2 += pow(y[i], 2);
132                         e3 += y[i];
133                 }
134
135                 ex d = e1 + 1;
136                 ex f = d * (e2 - 1);
137                 ex g = d * pow(e3 + 2, 2);
138                 ex r = gcd(f.expand(), g.expand());
139                 if (!(r - d).expand().is_zero()) {
140                         clog << "case 4, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
141                         return 1;
142                 }
143         }
144         return 0;
145 }
146
147 // Completely dense non-monic quadratic inputs with dense non-monic linear GCDs
148 static unsigned poly_gcd5()
149 {
150         for (int v=1; v<=MAX_VARIABLES; v++) {
151                 ex e1 = x + 1;
152                 ex e2 = x - 2;
153                 ex e3 = x + 2;
154                 for (int i=0; i<v; i++) {
155                         e1 *= y[i] + 1;
156                         e2 *= y[i] - 2;
157                         e3 *= y[i] + 2;
158                 }
159
160                 ex d = e1 - 3;
161                 ex f = d * (e2 + 3);
162                 ex g = d * (e3 - 3);
163                 ex r = gcd(f.expand(), g.expand());
164                 if (!(r - d).expand().is_zero()) {
165                         clog << "case 5, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
166                         return 1;
167                 }
168         }
169         return 0;
170 }
171
172 // Sparse non-monic quadratic inputs with linear GCDs
173 static unsigned poly_gcd5p()
174 {
175         for (int v=1; v<=MAX_VARIABLES; v++) {
176                 ex e1 = x;
177                 for (int i=0; i<v; i++)
178                         e1 *= y[i];
179
180                 ex d = e1 - 1;
181                 ex f = d * (e1 + 3);
182                 ex g = d * (e1 - 3);
183                 ex r = gcd(f.expand(), g.expand());
184                 if (!(r - d).expand().is_zero()) {
185                         clog << "case 5p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
186                         return 1;
187                 }
188         }
189         return 0;
190 }
191
192 // Trivariate inputs with increasing degrees
193 static unsigned poly_gcd6()
194 {
195         symbol y("y");
196
197         for (int j=1; j<=MAX_VARIABLES; j++) {
198                 ex d = pow(x, j) * y * (z - 1);
199                 ex f = d * (pow(x, j) + pow(y, j + 1) * pow(z, j) + 1);
200                 ex g = d * (pow(x, j + 1) + pow(y, j) * pow(z, j + 1) - 7);
201                 ex r = gcd(f.expand(), g.expand());
202                 if (!(r - d).expand().is_zero()) {
203                         clog << "case 6, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
204                         return 1;
205                 }
206         }
207         return 0;
208 }
209
210 // Trivariate polynomials whose GCD has common factors with its cofactors
211 static unsigned poly_gcd7()
212 {
213         symbol y("y");
214         ex p = x - y * z + 1;
215         ex q = x - y + z * 3;
216
217         for (int j=1; j<=MAX_VARIABLES; j++) {
218                 for (int k=j+1; k<=4; k++) {
219                         ex d = pow(p, j) * pow(q, j);
220                         ex f = pow(p, j) * pow(q, k);
221                         ex g = pow(p, k) * pow(q, j); 
222                         ex r = gcd(f, g);
223                         if (!(r - d).expand().is_zero() && !(r + d).expand().is_zero()) {
224                                 clog << "case 7, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
225                                 return 1;
226                         }
227                 }
228         }
229         return 0;
230 }
231
232 unsigned exam_polygcd()
233 {
234         unsigned result = 0;
235         
236         cout << "examining polynomial GCD computation" << flush;
237         
238         result += poly_gcd1();  cout << '.' << flush;
239         result += poly_gcd2();  cout << '.' << flush;
240         result += poly_gcd3();  cout << '.' << flush;
241         result += poly_gcd3p();  cout << '.' << flush; // PRS "worst" case
242         result += poly_gcd4();  cout << '.' << flush;
243         result += poly_gcd5();  cout << '.' << flush;
244         result += poly_gcd5p();  cout << '.' << flush;
245         result += poly_gcd6();  cout << '.' << flush;
246         result += poly_gcd7();  cout << '.' << flush;
247         
248         return result;
249 }
250
251 int main(int argc, char** argv)
252 {
253         return exam_polygcd();
254 }