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