check/exam_pgcd.cpp

   1 /** @file exam_pgcd.cpp
   2  *
   3  * Exam GCD over prime fields computations.
   4  */
   5 #include <string>
   6 #include <iostream>
   7 #include <utility>
   8 #include "ginac.h"
   9 using namespace std;
  10 using namespace GiNaC;
  11
  12 // Check for an infite loop in PGCD, fixed 2010-02-23.
  13 static unsigned pgcd_relatively_prime_bug()
  14 {
  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                 clog << " oops!" << endl <<
  25                         "** Error: should be 0, got " << g << endl;
  26                 return 1;
  27         }
  28         return 0;
  29 }
  30
  31 // Check for an infinite loop in PGCD, fixed 2010-03-18.
  32 static unsigned pgcd_infinite_loop()
  33 {
  34         parser the_parser;
  35         ex e = the_parser(string(R"ex(792*z^8*w^4*x^3*y^4*u^7
  36 + 24*z^4*w^4*x^2*y^3*u^4 + 264*z^8*w^3*x^2*y^7*u^5 + 198*z^4*w^5*x^5*y*u^6
  37 + 110*z^2*w^3*x^5*y^4*u^6 - 120*z^8*w*x^4*u^6 - 480*z^5*w*x^4*y^6*u^8
  38 - 720*z^7*x^3*y^3*u^7 + 165*z^4*w^2*x^4*y*u^5 + 450*z^8*w^6*x^2*y*u^8
  39 + 40*z^2*w^3*x^3*y^3*u^6 - 288*z^7*w^2*x^3*y^6*u^6 + 250*z^6*w^4*x^2*y^4*u^8
  40 + 576*z^7*w^7*x^2*y^4*u^8 - 80*z^6*w^2*x^5*y^3*u^7 - 144*z^8*w^4*x^5*u^7
  41 + 120*z^4*w*x^2*y^6*u^6 + 320*z^5*w^5*x^2*y^7*u^8 + 192*z^7*w^6*x*y^7*u^6
  42 - 12*z^4*w^3*x^3*y^5*u^6 - 36*z^4*w^4*x^4*y^2*u^8 + 72*z^4*w^5*x^3*u^6
  43 - 20*z^2*w^2*x^4*y^5*u^8 + 660*z^8*w*x^2*y^4*u^6 + 66*z^4*w^4*x^4*y^4*u^4
  44 + 440*z^6*w^2*x^3*y^7*u^7 - 30*z^4*w*x^3*y^2*u^7 - 48*z^8*w^3*x^4*y^3*u^5
  45 + 72*z^6*w^2*x*y^6*u^4 - 864*z^7*w^3*x^4*y^3*u^8 + 480*z^7*w^4*x*y^4*u^7
  46 + 60*z^4*w^2*x^2*u^5 + 375*z^8*w^3*x*y*u^7 + 150*z^8*w^5*x*y^4*u^6
  47 + 180*z^6*x*y^3*u^5 + 216*z^6*w^3*x^2*y^3*u^6)ex"));
  48         const symbol x = ex_to<symbol>(the_parser.get_syms()["x"]);
  49         ex g = gcd(e, e.diff(x));
  50         ex should_be = the_parser(string("u^4*z^2"));
  51         if (!(g-should_be).expand().is_zero()) {
  52                 clog << "GCD was miscomputed. " << endl;
  53                 return 1;
  54         }
  55         return 0;
  56 }
  57
  58 int main()
  59 {
  60         unsigned result = 0;
  61
  62         cout << "Examining pgcd() bugs (infinite loop, miscalculation)" << flush;
  63
  64         result += pgcd_relatively_prime_bug();  cout << '.' << flush;
  65         result += pgcd_infinite_loop();  cout << '.' << flush;
  66
  67         return result;
  68 }