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 }