From 9f0f6ede514e50a7722a1112c96b7a6e51668ddb Mon Sep 17 00:00:00 2001
From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Wed, 21 Jun 2000 19:06:40 +0000
Subject: [PATCH] - moved polynomial interpolation in heur_gcd() to its own
 function

---
 ginac/normal.cpp | 59 ++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 49 insertions(+), 10 deletions(-)

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 59323798..79d16461 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -59,7 +59,8 @@ namespace GiNaC {
 #define USE_REMEMBER 0
 
 // Set this if you want divide_in_z() to use trial division followed by
-// polynomial interpolation (usually slower except for very large problems)
+// polynomial interpolation (always slower except for completely dense
+// polynomials)
 #define USE_TRIAL_DIVISION 0
 
 // Set this to enable some statistical output for the GCD routines
@@ -1291,6 +1292,20 @@ ex mul::smod(const numeric &xi) const
 }
 
 
+/** xi-adic polynomial interpolation */
+static ex interpolate(const ex &gamma, const numeric &xi, const symbol &x)
+{
+	ex g = _ex0();
+	ex e = gamma;
+	numeric rxi = xi.inverse();
+	for (int i=0; !e.is_zero(); i++) {
+		ex gi = e.smod(xi);
+		g += gi * power(x, i);
+		e = (e - gi) * rxi;
+	}
+	return g;
+}
+
 /** Exception thrown by heur_gcd() to signal failure. */
 class gcdheu_failed {};
 
@@ -1355,21 +1370,17 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 		}
 
         // Apply evaluation homomorphism and calculate GCD
-        ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), NULL, NULL, var+1).expand();
+		ex cp, cq;
+        ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand();
         if (!is_ex_exactly_of_type(gamma, fail)) {
 
             // Reconstruct polynomial from GCD of mapped polynomials
-            ex g = _ex0();
-            numeric rxi = xi.inverse();
-            for (int i=0; !gamma.is_zero(); i++) {
-                ex gi = gamma.smod(xi);
-                g += gi * power(x, i);
-                gamma = (gamma - gi) * rxi;
-            }
+			ex g = interpolate(gamma, xi, x);
+
             // Remove integer content
             g /= g.integer_content();
 
-            // If the calculated polynomial divides both a and b, this is the GCD
+            // If the calculated polynomial divides both p and q, this is the GCD
             ex dummy;
             if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
                 g *= gc;
@@ -1379,6 +1390,34 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
                 else
                     return g;
             }
+#if 0
+			cp = interpolate(cp, xi, x);
+			if (divide_in_z(cp, p, g, var)) {
+				if (divide_in_z(g, q, cb ? *cb : dummy, var)) {
+					g *= gc;
+					if (ca)
+						*ca = cp;
+	                ex lc = g.lcoeff(x);
+	                if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative())
+	                    return -g;
+	                else
+	                    return g;
+				}
+			}
+			cq = interpolate(cq, xi, x);
+			if (divide_in_z(cq, q, g, var)) {
+				if (divide_in_z(g, p, ca ? *ca : dummy, var)) {
+					g *= gc;
+					if (cb)
+						*cb = cq;
+	                ex lc = g.lcoeff(x);
+	                if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative())
+	                    return -g;
+	                else
+	                    return g;
+				}
+			}
+#endif
         }
 
         // Next evaluation point
-- 
2.44.0