check/poly_gcd.cpp

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