- // Sort c and d so that c has higher degree
- ex c, d;
- int adeg = a.degree(*x), bdeg = b.degree(*x);
- int cdeg, ddeg;
- if (adeg >= bdeg) {
- c = a;
- d = b;
- cdeg = adeg;
- ddeg = bdeg;
- } else {
- c = b;
- d = a;
- cdeg = bdeg;
- ddeg = adeg;
- }
-
- // Remove content from c and d, to be attached to GCD later
- ex cont_c = c.content(*x);
- ex cont_d = d.content(*x);
- ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
- if (ddeg == 0)
- return gamma;
- c = c.primpart(*x, cont_c);
- d = d.primpart(*x, cont_d);
-
- // First element of subresultant sequence
- ex r = _ex0(), ri = _ex1(), psi = _ex1();
- int delta = cdeg - ddeg;
-
- for (;;) {
- // Calculate polynomial pseudo-remainder
- r = prem(c, d, *x, false);
- if (r.is_zero())
- return gamma * d.primpart(*x);
- c = d;
- cdeg = ddeg;
- if (!divide(r, ri * power(psi, delta), d, false))
- throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
- ddeg = d.degree(*x);
- if (ddeg == 0) {
- if (is_ex_exactly_of_type(r, numeric))
- return gamma;
- else
- return gamma * r.primpart(*x);
- }
-
- // Next element of subresultant sequence
- ri = c.expand().lcoeff(*x);
- if (delta == 1)
- psi = ri;
- else if (delta)
- divide(power(ri, delta), power(psi, delta-1), psi, false);
- delta = cdeg - ddeg;
- }
+//std::clog << "red_gcd(" << a << "," << b << ")\n";
+
+ // Sort c and d so that c has higher degree
+ ex c, d;
+ int adeg = a.degree(*x), bdeg = b.degree(*x);
+ int cdeg, ddeg;
+ if (adeg >= bdeg) {
+ c = a;
+ d = b;
+ cdeg = adeg;
+ ddeg = bdeg;
+ } else {
+ c = b;
+ d = a;
+ cdeg = bdeg;
+ ddeg = adeg;
+ }
+
+ // Remove content from c and d, to be attached to GCD later
+ ex cont_c = c.content(*x);
+ ex cont_d = d.content(*x);
+ ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
+ if (ddeg == 0)
+ return gamma;
+ c = c.primpart(*x, cont_c);
+ d = d.primpart(*x, cont_d);
+
+ // First element of divisor sequence
+ ex r, ri = _ex1();
+ int delta = cdeg - ddeg;
+
+ for (;;) {
+ // Calculate polynomial pseudo-remainder
+//std::clog << " d = " << d << endl;
+ r = prem(c, d, *x, false);
+ if (r.is_zero())
+ return gamma * d.primpart(*x);
+ c = d;
+ cdeg = ddeg;
+
+ if (!divide(r, pow(ri, delta), d, false))
+ throw(std::runtime_error("invalid expression in red_gcd(), division failed"));
+ ddeg = d.degree(*x);
+ if (ddeg == 0) {
+ if (is_ex_exactly_of_type(r, numeric))
+ return gamma;
+ else
+ return gamma * r.primpart(*x);
+ }
+
+ ri = c.expand().lcoeff(*x);
+ delta = cdeg - ddeg;
+ }
+}
+
+
+/** Compute GCD of multivariate polynomials using the subresultant PRS
+ * algorithm. This function is used internally by gcd().
+ *
+ * @param a first multivariate polynomial
+ * @param b second multivariate polynomial
+ * @param var iterator to first element of vector of sym_desc structs
+ * @return the GCD as a new expression
+ * @see gcd */
+
+static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
+{
+//std::clog << "sr_gcd(" << a << "," << b << ")\n";
+#if STATISTICS
+ sr_gcd_called++;
+#endif
+
+ // The first symbol is our main variable
+ const symbol &x = *(var->sym);
+
+ // Sort c and d so that c has higher degree
+ ex c, d;
+ int adeg = a.degree(x), bdeg = b.degree(x);
+ int cdeg, ddeg;
+ if (adeg >= bdeg) {
+ c = a;
+ d = b;
+ cdeg = adeg;
+ ddeg = bdeg;
+ } else {
+ c = b;
+ d = a;
+ cdeg = bdeg;
+ ddeg = adeg;
+ }
+
+ // Remove content from c and d, to be attached to GCD later
+ ex cont_c = c.content(x);
+ ex cont_d = d.content(x);
+ ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
+ if (ddeg == 0)
+ return gamma;
+ c = c.primpart(x, cont_c);
+ d = d.primpart(x, cont_d);
+//std::clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n";
+
+ // First element of subresultant sequence
+ ex r = _ex0(), ri = _ex1(), psi = _ex1();
+ int delta = cdeg - ddeg;
+
+ for (;;) {
+ // Calculate polynomial pseudo-remainder
+//std::clog << " start of loop, psi = " << psi << ", calculating pseudo-remainder...\n";
+//std::clog << " d = " << d << endl;
+ r = prem(c, d, x, false);
+ if (r.is_zero())
+ return gamma * d.primpart(x);
+ c = d;
+ cdeg = ddeg;
+//std::clog << " dividing...\n";
+ if (!divide_in_z(r, ri * pow(psi, delta), d, var))
+ throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
+ ddeg = d.degree(x);
+ if (ddeg == 0) {
+ if (is_ex_exactly_of_type(r, numeric))
+ return gamma;
+ else
+ return gamma * r.primpart(x);
+ }
+
+ // Next element of subresultant sequence
+//std::clog << " calculating next subresultant...\n";
+ ri = c.expand().lcoeff(x);
+ if (delta == 1)
+ psi = ri;
+ else if (delta)
+ divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1);
+ delta = cdeg - ddeg;
+ }