check/pgcd_relatively_prime_bug.cpp

   1 /** @file pgcd_relatively_prime_bug.cpp
   2  *
   3  * A program exposing historical bug in the pgcd() function.
   4  */
   5 #include <string>
   6 #include <iostream>
   7 #include <utility>
   8 #include "ginac.h"
   9 using namespace std;
  10 using namespace GiNaC;
  11
  12 int main(int argc, char** argv)
  13 {
  14         cout << "Checking for pgcd() bug regarding relatively prime polynomials: " << flush;
  15         const symbol q("q");
  16         parser reader;
  17         reader.get_syms().insert(make_pair(string("q"), q));
  18
  19         ex t = reader("-E20^16*E4^8*E5^8*E1^16*q^4"
  20                       "-(E10^24-E20^8*E5^16)*E4^16*E1^8"
  21                       "+E2^24*E20^16*E5^8*q^4");
  22         ex g = gcd(t.expand(), t.diff(q).expand()) - 1;
  23         if (!g.is_zero()) {
  24                 cout << " oops!" << endl <<
  25                         "** Error: should be 0, got " << g << endl << flush;
  26                 throw std::logic_error("gcd was miscalculated");
  27         }
  28         cout << "not found" << endl << flush;
  29         return 0;
  30 }
  31