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