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 }