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-2011 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 "ginac.h"
  25 using namespace GiNaC;
  26
  27 #include <iostream>
  28 using namespace std;
  29
  30 const int MAX_VARIABLES = 3;
  31
  32 static symbol x("x"), z("z");
  33 static symbol y[MAX_VARIABLES];
  34
  35 // GCD = 1
  36 static unsigned poly_gcd1()
  37 {
  38         for (int v=1; v<=MAX_VARIABLES; v++) {
  39                 ex e1 = x;
  40                 ex e2 = pow(x, 2);
  41                 for (int i=0; i<v; i++) {
  42                         e1 += y[i];
  43                         e2 += pow(y[i], 2);
  44                 }
  45
  46                 ex f = (e1 + 1) * (e1 + 2);
  47                 ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1);
  48                 ex r = gcd(f, g);
  49                 if (r != 1) {
  50                         clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl;
  51                         return 1;
  52                 }
  53         }
  54         return 0;
  55 }
  56
  57 // Linearly dense quartic inputs with quadratic GCDs
  58 static unsigned poly_gcd2()
  59 {
  60         for (int v=1; v<=MAX_VARIABLES; v++) {
  61                 ex e1 = x;
  62                 ex e2 = x;
  63                 for (int i=0; i<v; i++) {
  64                         e1 += y[i];
  65                         e2 -= y[i];
  66                 }
  67
  68                 ex d = pow(e1 + 1, 2);
  69                 ex f = d * pow(e2 - 2, 2);
  70                 ex g = d * pow(e1 + 2, 2);
  71                 ex r = gcd(f.expand(), g.expand());
  72                 if (!(r - d).expand().is_zero()) {
  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()
  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.expand(), g.expand());
  92                 if (!(r - d).expand().is_zero()) {
  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()
 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.expand(), g.expand());
 115                 if (!(r - d).expand().is_zero()) {
 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()
 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.expand(), g.expand());
 140                 if (!(r - d).expand().is_zero()) {
 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()
 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.expand(), g.expand());
 165                 if (!(r - d).expand().is_zero()) {
 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()
 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.expand(), g.expand());
 185                 if (!(r - d).expand().is_zero()) {
 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()
 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.expand(), g.expand());
 203                 if (!(r - d).expand().is_zero()) {
 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()
 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<=MAX_VARIABLES; 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().is_zero() && !(r + d).expand().is_zero()) {
 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 exam_polygcd()
 234 {
 235         unsigned result = 0;
 236
 237         cout << "examining polynomial GCD computation" << flush;
 238
 239         result += poly_gcd1();  cout << '.' << flush;
 240         result += poly_gcd2();  cout << '.' << flush;
 241         result += poly_gcd3();  cout << '.' << flush;
 242         result += poly_gcd3p();  cout << '.' << flush; // PRS "worst" case
 243         result += poly_gcd4();  cout << '.' << flush;
 244         result += poly_gcd5();  cout << '.' << flush;
 245         result += poly_gcd5p();  cout << '.' << flush;
 246         result += poly_gcd6();  cout << '.' << flush;
 247         result += poly_gcd7();  cout << '.' << flush;
 248
 249         return result;
 250 }
 251
 252 int main(int argc, char** argv)
 253 {
 254         return exam_polygcd();
 255 }