pgcd(): detect relatively prime polynomials properly...

... so pgcd() doesn't loop forever any more. Division test does not handle
relatively prime polynomials (because C = 1 divides any polynomial). So we
should stop interpolation (and return gcd of contents) if we got relatively
prime images (we should do that for efficiency reasons too).

Thanks to Jörg Arndt for a bug report.

diff --git a/check/Makefile.am b/check/Makefile.am
index 5b9dd588a903c931d164d64e4f0224798d474d5a..375161336806c7f89653cb5b63916220104a4d32 100644 (file)
--- a/check/Makefile.am
+++ b/check/Makefile.am
@@ -30,6 +30,7 @@ EXAMS = exam_paranoia \
        exam_misc \
        exam_mod_gcd \
        bugme_chinrem_gcd \
        exam_misc \
        exam_mod_gcd \
        bugme_chinrem_gcd \

+       pgcd_relatively_prime_bug \

        exam_cra
 
 TIMES = time_dennyfliegner \
        exam_cra
 
 TIMES = time_dennyfliegner \


@@ -255,6 +256,8 @@ time_parser_LDADD = ../ginac/libginac.la
 bugme_chinrem_gcd_SOURCES = bugme_chinrem_gcd.cpp
 bugme_chinrem_gcd_LDADD = ../ginac/libginac.la
 
 bugme_chinrem_gcd_SOURCES = bugme_chinrem_gcd.cpp
 bugme_chinrem_gcd_LDADD = ../ginac/libginac.la
 

+pgcd_relatively_prime_bug_SOURCES = pgcd_relatively_prime_bug.cpp
+pgcd_relatively_prime_bug_LDADD = ../ginac/libginac.la

 
 AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac -DIN_GINAC
 
 
 AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac -DIN_GINAC
 

diff --git a/check/pgcd_relatively_prime_bug.cpp b/check/pgcd_relatively_prime_bug.cpp
new file mode 100644 (file)
index 0000000..1d4c287
--- /dev/null
+++ b/check/pgcd_relatively_prime_bug.cpp
@@ -0,0 +1,31 @@
+/** @file pgcd_relatively_prime_bug.cpp
+ *
+ * A program exposing historical bug in the pgcd() function. 
+ */
+#include <string>
+#include <iostream>
+#include <utility>
+#include "ginac.h"
+using namespace std;
+using namespace GiNaC;
+
+int main(int argc, char** argv)
+{
+       cout << "Checking for pgcd() bug regarding relatively prime polynomials: " << flush;
+       const symbol q("q");
+       parser reader;
+       reader.get_syms().insert(make_pair(string("q"), q));
+
+       ex t = reader("-E20^16*E4^8*E5^8*E1^16*q^4"
+                     "-(E10^24-E20^8*E5^16)*E4^16*E1^8"
+                     "+E2^24*E20^16*E5^8*q^4");
+       ex g = gcd(t.expand(), t.diff(q).expand()) - 1;
+       if (!g.is_zero()) {
+               cout << " oops!" << endl <<
+                       "** Error: should be 0, got " << g << endl << flush;
+               throw std::logic_error("gcd was miscalculated");
+       }
+       cout << "not found" << endl << flush;
+       return 0;
+}
+

diff --git a/ginac/polynomial/pgcd.cpp b/ginac/polynomial/pgcd.cpp
index 6e5e6b2ebc4bbaadc6adc1bcc20401f672d5af9d..ef3748dcd49139a0d384c1a4f87ae4169a5bb647 100644 (file)
--- a/ginac/polynomial/pgcd.cpp
+++ b/ginac/polynomial/pgcd.cpp
@@ -121,6 +121,8 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
                        // evaluation point is bad. Skip it.
                        continue;
                }
                        // evaluation point is bad. Skip it.
                        continue;
                }

+               if (img_gcd_deg == 0)
+                       return cont_gcd;

 
                // Image has the same degree as the previous one
                // (or at least not higher than the limit)
 
                // Image has the same degree as the previous one
                // (or at least not higher than the limit)


@@ -145,8 +147,6 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
                        if (divide_in_z_p(Aprim, C, dummy1, vars, p) &&
                                        divide_in_z_p(Bprim, C, dummy2, vars, p))
                                return (cont_gcd*C).expand().smod(pn);
                        if (divide_in_z_p(Aprim, C, dummy1, vars, p) &&
                                        divide_in_z_p(Bprim, C, dummy2, vars, p))
                                return (cont_gcd*C).expand().smod(pn);

-                       else if (img_gcd_deg == 0)
-                               return cont_gcd;

                        // else continue building the candidate
                } 
        } while(true);
                        // else continue building the candidate
                } 
        } while(true);