* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* function returns for a given expression.
*
* @param e expression to search
- * @param x pointer to first symbol found (returned)
+ * @param x first symbol found (returned)
* @return "false" if no symbol was found, "true" otherwise */
-static bool get_first_symbol(const ex &e, const symbol *&x)
+static bool get_first_symbol(const ex &e, ex &x)
{
if (is_a<symbol>(e)) {
- x = &ex_to<symbol>(e);
+ x = e;
return true;
} else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
for (size_t i=0; i<e.nops(); i++)
*
* @see get_symbol_stats */
struct sym_desc {
- /** Pointer to symbol */
- const symbol *sym;
+ /** Reference to symbol */
+ ex sym;
/** Highest degree of symbol in polynomial "a" */
int deg_a;
typedef std::vector<sym_desc> sym_desc_vec;
// Add symbol the sym_desc_vec (used internally by get_symbol_stats())
-static void add_symbol(const symbol *s, sym_desc_vec &v)
+static void add_symbol(const ex &s, sym_desc_vec &v)
{
sym_desc_vec::const_iterator it = v.begin(), itend = v.end();
while (it != itend) {
- if (it->sym->compare(*s) == 0) // If it's already in there, don't add it a second time
+ if (it->sym.is_equal(s)) // If it's already in there, don't add it a second time
return;
++it;
}
static void collect_symbols(const ex &e, sym_desc_vec &v)
{
if (is_a<symbol>(e)) {
- add_symbol(&ex_to<symbol>(e), v);
+ add_symbol(e, v);
} else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
for (size_t i=0; i<e.nops(); i++)
collect_symbols(e.op(i), v);
collect_symbols(b.eval(), v);
sym_desc_vec::iterator it = v.begin(), itend = v.end();
while (it != itend) {
- int deg_a = a.degree(*(it->sym));
- int deg_b = b.degree(*(it->sym));
+ int deg_a = a.degree(it->sym);
+ int deg_b = b.degree(it->sym);
it->deg_a = deg_a;
it->deg_b = deg_b;
it->max_deg = std::max(deg_a, deg_b);
- it->max_lcnops = std::max(a.lcoeff(*(it->sym)).nops(), b.lcoeff(*(it->sym)).nops());
- it->ldeg_a = a.ldegree(*(it->sym));
- it->ldeg_b = b.ldegree(*(it->sym));
+ it->max_lcnops = std::max(a.lcoeff(it->sym).nops(), b.lcoeff(it->sym).nops());
+ it->ldeg_a = a.ldegree(it->sym);
+ it->ldeg_b = b.ldegree(it->sym);
++it;
}
std::sort(v.begin(), v.end());
std::clog << "Symbols:\n";
it = v.begin(); itend = v.end();
while (it != itend) {
- std::clog << " " << *it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
- std::clog << " lcoeff_a=" << a.lcoeff(*(it->sym)) << ", lcoeff_b=" << b.lcoeff(*(it->sym)) << endl;
+ std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
+ std::clog << " lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << endl;
++it;
}
#endif
/** Compute the integer content (= GCD of all numeric coefficients) of an
* expanded polynomial.
*
- * @param e expanded polynomial
* @return integer content */
numeric ex::integer_content() const
{
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
* @return quotient of a and b in Q[x] */
-ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
+ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
{
if (b.is_zero())
throw(std::overflow_error("quo: division by zero"));
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
* @return remainder of a(x) and b(x) in Q[x] */
-ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
+ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
{
if (b.is_zero())
throw(std::overflow_error("rem: division by zero"));
* @param a rational function in x
* @param x a is a function of x
* @return decomposed function. */
-ex decomp_rational(const ex &a, const symbol &x)
+ex decomp_rational(const ex &a, const ex &x)
{
ex nd = numer_denom(a);
ex numer = nd.op(0), denom = nd.op(1);
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
* @return pseudo-remainder of a(x) and b(x) in Q[x] */
-ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
+ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
{
if (b.is_zero())
throw(std::overflow_error("prem: division by zero"));
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
* @return sparse pseudo-remainder of a(x) and b(x) in Q[x] */
-ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args)
+ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
{
if (b.is_zero())
throw(std::overflow_error("prem: division by zero"));
throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
// Find first symbol
- const symbol *x;
+ ex x;
if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
throw(std::invalid_argument("invalid expression in divide()"));
q = _ex0;
return true;
}
- int bdeg = b.degree(*x);
- int rdeg = r.degree(*x);
- ex blcoeff = b.expand().coeff(*x, bdeg);
+ int bdeg = b.degree(x);
+ int rdeg = r.degree(x);
+ ex blcoeff = b.expand().coeff(x, bdeg);
bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
while (rdeg >= bdeg) {
- ex term, rcoeff = r.coeff(*x, rdeg);
+ ex term, rcoeff = r.coeff(x, rdeg);
if (blcoeff_is_numeric)
term = rcoeff / blcoeff;
else
if (!divide(rcoeff, blcoeff, term, false))
return false;
- term *= power(*x, rdeg - bdeg);
+ term *= power(x, rdeg - bdeg);
v.push_back(term);
r -= (term * b).expand();
if (r.is_zero()) {
q = (new add(v))->setflag(status_flags::dynallocated);
return true;
}
- rdeg = r.degree(*x);
+ rdeg = r.degree(x);
}
return false;
}
/** Exact polynomial division of a(X) by b(X) in Z[X].
* This functions works like divide() but the input and output polynomials are
* in Z[X] instead of Q[X] (i.e. they have integer coefficients). Unlike
- * divide(), it doesnยดt check whether the input polynomials really are integer
+ * divide(), it doesn't check whether the input polynomials really are integer
* polynomials, so be careful of what you pass in. Also, you have to run
* get_symbol_stats() over the input polynomials before calling this function
* and pass an iterator to the first element of the sym_desc vector. This
#endif
// Main symbol
- const symbol *x = var->sym;
+ const ex &x = var->sym;
// Compare degrees
- int adeg = a.degree(*x), bdeg = b.degree(*x);
+ int adeg = a.degree(x), bdeg = b.degree(x);
if (bdeg > adeg)
return false;
numeric point = _num0;
ex c;
for (i=0; i<=adeg; i++) {
- ex bs = b.subs(*x == point, subs_options::no_pattern);
+ ex bs = b.subs(x == point, subs_options::no_pattern);
while (bs.is_zero()) {
point += _num1;
- bs = b.subs(*x == point, subs_options::no_pattern);
+ bs = b.subs(x == point, subs_options::no_pattern);
}
- if (!divide_in_z(a.subs(*x == point, subs_options::no_pattern), bs, c, var+1))
+ if (!divide_in_z(a.subs(x == point, subs_options::no_pattern), bs, c, var+1))
return false;
alpha.push_back(point);
u.push_back(c);
// Convert from Newton form to standard form
c = v[adeg];
for (k=adeg-1; k>=0; k--)
- c = c * (*x - alpha[k]) + v[k];
+ c = c * (x - alpha[k]) + v[k];
- if (c.degree(*x) == (adeg - bdeg)) {
+ if (c.degree(x) == (adeg - bdeg)) {
q = c.expand();
return true;
} else
return true;
int rdeg = adeg;
ex eb = b.expand();
- ex blcoeff = eb.coeff(*x, bdeg);
+ ex blcoeff = eb.coeff(x, bdeg);
exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
while (rdeg >= bdeg) {
- ex term, rcoeff = r.coeff(*x, rdeg);
+ ex term, rcoeff = r.coeff(x, rdeg);
if (!divide_in_z(rcoeff, blcoeff, term, var+1))
break;
- term = (term * power(*x, rdeg - bdeg)).expand();
+ term = (term * power(x, rdeg - bdeg)).expand();
v.push_back(term);
r -= (term * eb).expand();
if (r.is_zero()) {
#endif
return true;
}
- rdeg = r.degree(*x);
+ rdeg = r.degree(x);
}
#if USE_REMEMBER
dr_remember[ex2(a, b)] = exbool(q, false);
* @param x variable in which to compute the unit part
* @return unit part
* @see ex::content, ex::primpart */
-ex ex::unit(const symbol &x) const
+ex ex::unit(const ex &x) const
{
ex c = expand().lcoeff(x);
if (is_exactly_a<numeric>(c))
return c < _ex0 ? _ex_1 : _ex1;
else {
- const symbol *y;
+ ex y;
if (get_first_symbol(c, y))
- return c.unit(*y);
+ return c.unit(y);
else
throw(std::invalid_argument("invalid expression in unit()"));
}
* @param x variable in which to compute the content part
* @return content part
* @see ex::unit, ex::primpart */
-ex ex::content(const symbol &x) const
+ex ex::content(const ex &x) const
{
if (is_zero())
return _ex0;
* @param x variable in which to compute the primitive part
* @return primitive part
* @see ex::unit, ex::content */
-ex ex::primpart(const symbol &x) const
+ex ex::primpart(const ex &x) const
{
if (is_zero())
return _ex0;
* @param x variable in which to compute the primitive part
* @param c previously computed content part
* @return primitive part */
-ex ex::primpart(const symbol &x, const ex &c) const
+ex ex::primpart(const ex &x, const ex &c) const
{
if (is_zero())
return _ex0;
#endif
// The first symbol is our main variable
- const symbol &x = *(var->sym);
+ const ex &x = var->sym;
// Sort c and d so that c has higher degree
ex c, d;
/** Return maximum (absolute value) coefficient of a polynomial.
* This function is used internally by heur_gcd().
*
- * @param e expanded multivariate polynomial
* @return maximum coefficient
* @see heur_gcd */
numeric ex::max_coefficient() const
/** xi-adic polynomial interpolation */
-static ex interpolate(const ex &gamma, const numeric &xi, const symbol &x, int degree_hint = 1)
+static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint = 1)
{
exvector g; g.reserve(degree_hint);
ex e = gamma;
}
// The first symbol is our main variable
- const symbol &x = *(var->sym);
+ const ex &x = var->sym;
// Remove integer content
numeric gc = gcd(a.integer_content(), b.integer_content());
*
* @param a first multivariate polynomial
* @param b second multivariate polynomial
+ * @param ca pointer to expression that will receive the cofactor of a, or NULL
+ * @param cb pointer to expression that will receive the cofactor of b, or NULL
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
* @return the GCD as a new expression */
// The symbol with least degree is our main variable
sym_desc_vec::const_iterator var = sym_stats.begin();
- const symbol &x = *(var->sym);
+ const ex &x = var->sym;
// Cancel trivial common factor
int ldeg_a = var->ldeg_a;
/** Compute a square-free factorization of a multivariate polynomial in Q[X].
*
* @param a multivariate polynomial over Q[X]
- * @param x lst of variables to factor in, may be left empty for autodetection
+ * @param l lst of variables to factor in, may be left empty for autodetection
* @return a square-free factorization of \p a.
*
* \note
get_symbol_stats(a, _ex0, sdv);
sym_desc_vec::const_iterator it = sdv.begin(), itend = sdv.end();
while (it != itend) {
- args.append(*it->sym);
+ args.append(it->sym);
++it;
}
} else {
const ex tmp = multiply_lcm(a,lcm);
// find the factors
- exvector factors = sqrfree_yun(tmp,x);
+ exvector factors = sqrfree_yun(tmp, x);
// construct the next list of symbols with the first element popped
lst newargs = args;
// because subs() is not recursive
ex es = (new symbol)->setflag(status_flags::dynallocated);
ex e_replaced = e.subs(repl, subs_options::no_pattern);
- repl[es] = e_replaced;
- rev_lookup[e_replaced] = es;
+ repl.insert(std::make_pair(es, e_replaced));
+ rev_lookup.insert(std::make_pair(e_replaced, es));
return es;
}
/** Create a symbol for replacing the expression "e" (or return a previously
- * assigned symbol). An expression of the form "symbol == expression" is added
- * to repl_lst and the symbol is returned.
+ * assigned symbol). The symbol and expression are appended to repl, and the
+ * symbol is returned.
* @see basic::to_rational
* @see basic::to_polynomial */
-static ex replace_with_symbol(const ex & e, lst & repl_lst)
+static ex replace_with_symbol(const ex & e, exmap & repl)
{
- // Expression already in repl_lst? Then return the assigned symbol
- for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
- if (it->op(1).is_equal(e))
- return it->op(0);
+ // Expression already replaced? Then return the assigned symbol
+ for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
+ if (it->second.is_equal(e))
+ return it->first;
// Otherwise create new symbol and add to list, taking care that the
- // replacement expression doesn't itself contain symbols from the repl_lst,
+ // replacement expression doesn't itself contain symbols from repl,
// because subs() is not recursive
ex es = (new symbol)->setflag(status_flags::dynallocated);
- ex e_replaced = e.subs(repl_lst, subs_options::no_pattern);
- repl_lst.append(es == e_replaced);
+ ex e_replaced = e.subs(repl, subs_options::no_pattern);
+ repl.insert(std::make_pair(es, e_replaced));
return es;
}
den *= _ex_1;
}
} else {
- const symbol *x;
+ ex x;
if (get_first_symbol(den, x)) {
- GINAC_ASSERT(is_exactly_a<numeric>(den.unit(*x)));
- if (ex_to<numeric>(den.unit(*x)).is_negative()) {
+ GINAC_ASSERT(is_exactly_a<numeric>(den.unit(x)));
+ if (ex_to<numeric>(den.unit(x)).is_negative()) {
num *= _ex_1;
den *= _ex_1;
}
* on non-rational functions by applying to_rational() on the arguments,
* calling the desired function and re-substituting the temporary symbols
* in the result. To make the last step possible, all temporary symbols and
- * their associated expressions are collected in the list specified by the
- * repl_lst parameter in the form {symbol == expression}, ready to be passed
- * as an argument to ex::subs().
+ * their associated expressions are collected in the map specified by the
+ * repl parameter, ready to be passed as an argument to ex::subs().
*
- * @param repl_lst collects a list of all temporary symbols and their replacements
+ * @param repl collects all temporary symbols and their replacements
* @return rationalized expression */
-ex ex::to_rational(lst &repl_lst) const
+ex ex::to_rational(exmap & repl) const
{
- return bp->to_rational(repl_lst);
+ return bp->to_rational(repl);
}
-ex ex::to_polynomial(lst &repl_lst) const
+// GiNaC 1.1 compatibility function
+ex ex::to_rational(lst & repl_lst) const
{
- return bp->to_polynomial(repl_lst);
+ // Convert lst to exmap
+ exmap m;
+ for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
+ m.insert(std::make_pair(it->op(0), it->op(1)));
+
+ ex ret = bp->to_rational(m);
+
+ // Convert exmap back to lst
+ repl_lst.remove_all();
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
+ repl_lst.append(it->first == it->second);
+
+ return ret;
}
+ex ex::to_polynomial(exmap & repl) const
+{
+ return bp->to_polynomial(repl);
+}
+
+// GiNaC 1.1 compatibility function
+ex ex::to_polynomial(lst & repl_lst) const
+{
+ // Convert lst to exmap
+ exmap m;
+ for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
+ m.insert(std::make_pair(it->op(0), it->op(1)));
+
+ ex ret = bp->to_polynomial(m);
+
+ // Convert exmap back to lst
+ repl_lst.remove_all();
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
+ repl_lst.append(it->first == it->second);
+
+ return ret;
+}
/** Default implementation of ex::to_rational(). This replaces the object with
* a temporary symbol. */
-ex basic::to_rational(lst &repl_lst) const
+ex basic::to_rational(exmap & repl) const
{
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
-ex basic::to_polynomial(lst &repl_lst) const
+ex basic::to_polynomial(exmap & repl) const
{
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
/** Implementation of ex::to_rational() for symbols. This returns the
* unmodified symbol. */
-ex symbol::to_rational(lst &repl_lst) const
+ex symbol::to_rational(exmap & repl) const
{
return *this;
}
/** Implementation of ex::to_polynomial() for symbols. This returns the
* unmodified symbol. */
-ex symbol::to_polynomial(lst &repl_lst) const
+ex symbol::to_polynomial(exmap & repl) const
{
return *this;
}
/** Implementation of ex::to_rational() for a numeric. It splits complex
* numbers into re+I*im and replaces I and non-rational real numbers with a
* temporary symbol. */
-ex numeric::to_rational(lst &repl_lst) const
+ex numeric::to_rational(exmap & repl) const
{
if (is_real()) {
if (!is_rational())
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
} else { // complex
numeric re = real();
numeric im = imag();
- ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
- ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
- return re_ex + im_ex * replace_with_symbol(I, repl_lst);
+ ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
+ ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
+ return re_ex + im_ex * replace_with_symbol(I, repl);
}
return *this;
}
/** Implementation of ex::to_polynomial() for a numeric. It splits complex
* numbers into re+I*im and replaces I and non-integer real numbers with a
* temporary symbol. */
-ex numeric::to_polynomial(lst &repl_lst) const
+ex numeric::to_polynomial(exmap & repl) const
{
if (is_real()) {
if (!is_integer())
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
} else { // complex
numeric re = real();
numeric im = imag();
- ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl_lst);
- ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl_lst);
- return re_ex + im_ex * replace_with_symbol(I, repl_lst);
+ ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl);
+ ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl);
+ return re_ex + im_ex * replace_with_symbol(I, repl);
}
return *this;
}
/** Implementation of ex::to_rational() for powers. It replaces non-integer
* powers by temporary symbols. */
-ex power::to_rational(lst &repl_lst) const
+ex power::to_rational(exmap & repl) const
{
if (exponent.info(info_flags::integer))
- return power(basis.to_rational(repl_lst), exponent);
+ return power(basis.to_rational(repl), exponent);
else
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
/** Implementation of ex::to_polynomial() for powers. It replaces non-posint
* powers by temporary symbols. */
-ex power::to_polynomial(lst &repl_lst) const
+ex power::to_polynomial(exmap & repl) const
{
if (exponent.info(info_flags::posint))
- return power(basis.to_rational(repl_lst), exponent);
+ return power(basis.to_rational(repl), exponent);
else
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
/** Implementation of ex::to_rational() for expairseqs. */
-ex expairseq::to_rational(lst &repl_lst) const
+ex expairseq::to_rational(exmap & repl) const
{
epvector s;
s.reserve(seq.size());
epvector::const_iterator i = seq.begin(), end = seq.end();
while (i != end) {
- s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_rational(repl_lst)));
+ s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_rational(repl)));
++i;
}
- ex oc = overall_coeff.to_rational(repl_lst);
+ ex oc = overall_coeff.to_rational(repl);
if (oc.info(info_flags::numeric))
return thisexpairseq(s, overall_coeff);
else
}
/** Implementation of ex::to_polynomial() for expairseqs. */
-ex expairseq::to_polynomial(lst &repl_lst) const
+ex expairseq::to_polynomial(exmap & repl) const
{
epvector s;
s.reserve(seq.size());
epvector::const_iterator i = seq.begin(), end = seq.end();
while (i != end) {
- s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_polynomial(repl_lst)));
+ s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_polynomial(repl)));
++i;
}
- ex oc = overall_coeff.to_polynomial(repl_lst);
+ ex oc = overall_coeff.to_polynomial(repl);
if (oc.info(info_flags::numeric))
return thisexpairseq(s, overall_coeff);
else
/** Remove the common factor in the terms of a sum 'e' by calculating the GCD,
* and multiply it into the expression 'factor' (which needs to be initialized
* to 1, unless you're accumulating factors). */
-static ex find_common_factor(const ex & e, ex & factor, lst & repl)
+static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
{
if (is_exactly_a<add>(e)) {
{
if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
- lst repl;
+ exmap repl;
ex factor = 1;
ex r = find_common_factor(e, factor, repl);
return factor.subs(repl, subs_options::no_pattern) * r.subs(repl, subs_options::no_pattern);
}
+/** 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
git://www.ginac.de / ginac.git / blobdiff