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