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