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