From 3cbd851611c82629116046ecbe36fd85f373496c Mon Sep 17 00:00:00 2001
From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Mon, 28 Jun 2004 23:09:11 +0000
Subject: [PATCH] * Added resultant() function (by Ralf Stephan
 <ralf@ark.in-berlin.de>).

---
 doc/tutorial/ginac.texi | 34 +++++++++++++++++++++++++++++++++-
 ginac/normal.cpp        | 29 +++++++++++++++++++++++++++++
 ginac/normal.h          |  3 +++
 ginsh/ginsh.1.in        |  3 +++
 ginsh/ginsh_parser.yy   |  7 +++++++
 5 files changed, 75 insertions(+), 1 deletion(-)

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index b58d631b..4e789a51 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -4723,7 +4723,7 @@ content parts). The product of unit, content, and primitive part is the
 original polynomial.
 
 
-@subsection GCD and LCM
+@subsection GCD, LCM and resultant
 @cindex GCD
 @cindex LCM
 @cindex @code{gcd()}
@@ -4760,6 +4760,38 @@ int main()
 @}
 @end example
 
+@cindex resultant
+@cindex @code{resultant()}
+
+The resultant of two expressions only makes sense with polynomials.
+It is always computed with respect to a specific symbol within the
+expressions. The function has the interface
+
+@example
+ex resultant(const ex & a, const ex & b, const symbol & s);
+@end example
+
+Resultants are symmetric in @code{a} and @code{b}. The following example
+computes the resultant of two expressions with respect to @code{x} and
+@code{y}, respectively:
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+    symbol x("x"), y("y");
+
+    ex e1 = x+pow(y,2), e2 = 2*pow(x,3)-1; // x+y^2, 2*x^3-1
+    ex r;
+    
+    r = resultant (e1, e2, x); 
+    // -> 1+2*y^6
+    r = resultant (e1, e2, y); 
+    // -> 1-4*x^3+4*x^6
+@}
+@end example
 
 @subsection Square-free decomposition
 @cindex square-free decomposition
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 98efe985..8ba5df0d 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -2375,4 +2375,33 @@ ex collect_common_factors(const ex & e)
 }
 
 
+/** Resultant of two expressions e1,e2 with respect to symbol s.
+ *  Method: Compute determinant of Sylvester matrix of e1,e2,s.  */
+ex resultant(const ex & e1, const ex & e2, const ex & s)
+{
+	const ex ee1 = e1.expand();
+	const ex ee2 = e2.expand();
+	const int h1 = ee1.degree(s);
+	const int l1 = ee1.ldegree(s);
+	const int h2 = ee2.degree(s);
+	const int l2 = ee2.ldegree(s);
+
+	const int msize = h1 + h2;
+	matrix m(msize, msize);
+
+	for (int l = h1; l >= l1; --l) {
+		const ex e = ee1.coeff(s, l);
+		for (int k = 0; k < h2; ++k)
+			m(k, k+h1-l) = e;
+	}
+	for (int l = h2; l >= l2; --l) {
+		const ex e = ee2.coeff(s, l);
+		for (int k = 0; k < h1; ++k)
+			m(k+h2, k+h2-l) = e;
+	}
+
+	return m.determinant();
+}
+
+
 } // namespace GiNaC
diff --git a/ginac/normal.h b/ginac/normal.h
index 40f291ea..ff418bbe 100644
--- a/ginac/normal.h
+++ b/ginac/normal.h
@@ -66,6 +66,9 @@ extern ex sqrfree_parfrac(const ex & a, const symbol & x);
 // Collect common factors in sums.
 extern ex collect_common_factors(const ex & e);
 
+// Resultant of two polynomials e1,e2 with respect to symbol s.
+extern ex resultant(const ex & e1, const ex & e2, const ex & s);
+
 } // namespace GiNaC
 
 #endif // ndef __GINAC_NORMAL_H__
diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in
index 1662ce29..21db17e6 100644
--- a/ginsh/ginsh.1.in
+++ b/ginsh/ginsh.1.in
@@ -356,6 +356,9 @@ detail here. Please refer to the GiNaC documentation.
 .BI rem( expression ", " expression ", " symbol )
 \- remainder of polynomials
 .br
+.BI resultant( expression ", " expression ", " symbol )
+\- resultant of two polynomials with respect to symbol s
+.br
 .BI series( expression ", " relation-or-symbol ", " order )
 \- series expansion
 .br
diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy
index ed2fcda6..be65cf8e 100644
--- a/ginsh/ginsh_parser.yy
+++ b/ginsh/ginsh_parser.yy
@@ -480,6 +480,12 @@ static ex f_rem(const exprseq &e)
 	return rem(e[0], e[1], e[2]);
 }
 
+static ex f_resultant(const exprseq &e)
+{
+	CHECK_ARG(2, symbol, resultant);
+	return resultant(e[0], e[1], ex_to<symbol>(e[2]));
+}
+
 static ex f_series(const exprseq &e)
 {
 	CHECK_ARG(2, numeric, series);
@@ -589,6 +595,7 @@ static const fcn_init builtin_fcns[] = {
 	{"quo", f_quo, 3},
 	{"rank", f_rank, 1},
 	{"rem", f_rem, 3},
+	{"resultant", f_resultant, 3},
 	{"series", f_series, 3},
 	{"sprem", f_sprem, 3},
 	{"sqrfree", f_sqrfree1, 1},
-- 
2.44.0