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  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <ginac/ginac.h>
  24
  25 const int MAX_VARIABLES = 5;
  26
  27 static symbol x("x"), z("z");
  28 static symbol y[MAX_VARIABLES];
  29
  30 // GCD = 1
  31 static unsigned poly_gcd1(void)
  32 {
  33         for (int v=1; v<=MAX_VARIABLES; v++) {
  34                 ex e1 = x;
  35                 ex e2 = pow(x, 2);
  36                 for (int i=0; i<v; i++) {
  37                         e1 += y[i];
  38                         e2 += pow(y[i], 2);
  39                 }
  40
  41                 ex f = (e1 + 1) * (e1 + 2);
  42                 ex g = e2 * (-pow(x, 2) * y[0] * 3 + pow(y[0], 2) - 1);
  43                 ex r = gcd(f, g);
  44                 if (r != exONE()) {
  45                         clog << "case 1, gcd(" << f << "," << g << ") = " << r << " (should be 1)" << endl;
  46                         return 1;
  47                 }
  48         }
  49         return 0;
  50 }
  51
  52 // Linearly dense quartic inputs with quadratic GCDs
  53 static unsigned poly_gcd2(void)
  54 {
  55         for (int v=1; v<=MAX_VARIABLES; v++) {
  56                 ex e1 = x;
  57                 ex e2 = x;
  58                 for (int i=0; i<v; i++) {
  59                         e1 += y[i];
  60                         e2 -= y[i];
  61                 }
  62
  63                 ex d = pow(e1 + 1, 2);
  64                 ex f = d * pow(e2 - 2, 2);
  65                 ex g = d * pow(e1 + 2, 2);
  66                 ex r = gcd(f, g);
  67         ex re=r.expand();
  68         ex df1=r-d;
  69         ex df=(r-d).expand();
  70                 if ((r - d).expand().compare(exZERO()) != 0) {
  71                         clog << "case 2, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
  72                         return 1;
  73                 }
  74         }
  75         return 0;
  76 }
  77
  78 // Sparse GCD and inputs where degrees are proportional to the number of variables
  79 static unsigned poly_gcd3(void)
  80 {
  81         for (int v=1; v<=MAX_VARIABLES; v++) {
  82                 ex e1 = pow(x, v + 1);
  83                 for (int i=0; i<v; i++)
  84                         e1 += pow(y[i], v + 1);
  85
  86                 ex d = e1 + 1;
  87                 ex f = d * (e1 - 2);
  88                 ex g = d * (e1 + 2);
  89                 ex r = gcd(f, g);
  90                 if ((r - d).expand().compare(exZERO()) != 0) {
  91                         clog << "case 3, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
  92                         return 1;
  93                 }
  94         }
  95         return 0;
  96 }
  97
  98 // Variation of case 3; major performance degradation with PRS
  99 static unsigned poly_gcd3p(void)
 100 {
 101         for (int v=1; v<=MAX_VARIABLES; v++) {
 102                 ex e1 = pow(x, v + 1);
 103                 ex e2 = pow(x, v);
 104                 for (int i=0; i<v; i++) {
 105                         e1 += pow(y[i], v + 1);
 106                         e2 += pow(y[i], v);
 107                 }
 108
 109                 ex d = e1 + 1;
 110                 ex f = d * (e1 - 2);
 111                 ex g = d * (e2 + 2);
 112                 ex r = gcd(f, g);
 113                 if ((r - d).expand().compare(exZERO()) != 0) {
 114                         clog << "case 3p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
 115                         return 1;
 116                 }
 117         }
 118         return 0;
 119 }
 120
 121 // Quadratic non-monic GCD; f and g have other quadratic factors
 122 static unsigned poly_gcd4(void)
 123 {
 124         for (int v=1; v<=MAX_VARIABLES; v++) {
 125                 ex e1 = pow(x, 2) * pow(y[0], 2);
 126                 ex e2 = pow(x, 2) - pow(y[0], 2);
 127                 ex e3 = x * y[0];
 128                 for (int i=1; i<v; i++) {
 129                         e1 += pow(y[i], 2);
 130                         e2 += pow(y[i], 2);
 131                         e3 += y[i];
 132                 }
 133
 134                 ex d = e1 + 1;
 135                 ex f = d * (e2 - 1);
 136                 ex g = d * pow(e3 + 2, 2);
 137                 ex r = gcd(f, g);
 138                 if ((r - d).expand().compare(exZERO()) != 0) {
 139                         clog << "case 4, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
 140                         return 1;
 141                 }
 142         }
 143         return 0;
 144 }
 145
 146 // Completely dense non-monic quadratic inputs with dense non-monic linear GCDs
 147 static unsigned poly_gcd5(void)
 148 {
 149         for (int v=1; v<=MAX_VARIABLES; v++) {
 150                 ex e1 = x + 1;
 151                 ex e2 = x - 2;
 152                 ex e3 = x + 2;
 153                 for (int i=0; i<v; i++) {
 154                         e1 *= y[i] + 1;
 155                         e2 *= y[i] - 2;
 156                         e3 *= y[i] + 2;
 157                 }
 158
 159                 ex d = e1 - 3;
 160                 ex f = d * (e2 + 3);
 161                 ex g = d * (e3 - 3);
 162                 ex r = gcd(f, g);
 163                 if ((r - d).expand().compare(exZERO()) != 0) {
 164                         clog << "case 5, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
 165                         return 1;
 166                 }
 167         }
 168         return 0;
 169 }
 170
 171 // Sparse non-monic quadratic inputs with linear GCDs
 172 static unsigned poly_gcd5p(void)
 173 {
 174         for (int v=1; v<=MAX_VARIABLES; v++) {
 175                 ex e1 = x;
 176                 for (int i=0; i<v; i++)
 177                         e1 *= y[i];
 178
 179                 ex d = e1 - 1;
 180                 ex f = d * (e1 + 3);
 181                 ex g = d * (e1 - 3);
 182                 ex r = gcd(f, g);
 183                 if ((r - d).expand().compare(exZERO()) != 0) {
 184                         clog << "case 5p, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
 185                         return 1;
 186                 }
 187         }
 188         return 0;
 189 }
 190
 191 // Trivariate inputs with increasing degrees
 192 static unsigned poly_gcd6(void)
 193 {
 194         symbol y("y");
 195
 196         for (int j=1; j<=MAX_VARIABLES; j++) {
 197                 ex d = pow(x, j) * y * (z - 1);
 198                 ex f = d * (pow(x, j) + pow(y, j + 1) * pow(z, j) + 1);
 199                 ex g = d * (pow(x, j + 1) + pow(y, j) * pow(z, j + 1) - 7);
 200                 ex r = gcd(f, g);
 201                 if ((r - d).expand().compare(exZERO()) != 0) {
 202                         clog << "case 6, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
 203                         return 1;
 204                 }
 205         }
 206         return 0;
 207 }
 208
 209 // Trivariate polynomials whose GCD has common factors with its cofactors
 210 static unsigned poly_gcd7(void)
 211 {
 212         symbol y("y");
 213         ex p = x - y * z + 1;
 214         ex q = x - y + z * 3;
 215
 216         for (int j=1; j<=3; j++) {
 217                 for (int k=j+1; k<=4; k++) {
 218                         ex d = pow(p, j) * pow(q, j);
 219                         ex f = pow(p, j) * pow(q, k);
 220                         ex g = pow(p, k) * pow(q, j);
 221                         ex r = gcd(f, g);
 222                         if ((r - d).expand().compare(exZERO()) != 0 && (r + d).expand().compare(exZERO()) != 0) {
 223                                 clog << "case 7, gcd(" << f << "," << g << ") = " << r << " (should be " << d << ")" << endl;
 224                                 return 1;
 225                         }
 226                 }
 227         }
 228         return 0;
 229 }
 230
 231 unsigned poly_gcd(void)
 232 {
 233     unsigned result = 0;
 234
 235         cout << "checking polynomial GCD computation..." << flush;
 236         clog << "---------polynomial GCD computation:" << endl;
 237
 238         result += poly_gcd1();
 239         result += poly_gcd2();
 240         result += poly_gcd3();
 241 //      result += poly_gcd3p(); // takes extremely long (PRS "worst" case)
 242         result += poly_gcd4();
 243         result += poly_gcd5();
 244         result += poly_gcd5p();
 245         result += poly_gcd6();
 246         result += poly_gcd7();
 247
 248         if (!result) {
 249                 cout << " passed ";
 250         clog << "(no output)" << endl;
 251     } else {
 252                 cout << " failed ";
 253     }
 254         return result;
 255 }