... 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.
(cherry picked from commit
6afbb9793e359cf388462c471ea256a6662b0cd4)
exam_misc \
exam_mod_gcd \
bugme_chinrem_gcd \
+ pgcd_relatively_prime_bug \
exam_cra
TIMES = time_dennyfliegner \
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
--- /dev/null
+/** @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;
+}
+
// 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)
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);