There are a couple of ways to construct matrices, with or without preset
elements:
+@cindex @code{lst_to_matrix()}
+@cindex @code{diag_matrix()}
+@cindex @code{unit_matrix()}
+@cindex @code{symbolic_matrix()}
@example
matrix::matrix(unsigned r, unsigned c);
matrix::matrix(unsigned r, unsigned c, const lst & l);
ex lst_to_matrix(const lst & l);
ex diag_matrix(const lst & l);
+ex unit_matrix(unsigned x);
+ex unit_matrix(unsigned r, unsigned c);
+ex symbolic_matrix(unsigned r, unsigned c, const string & base_name);
+ex symbolic_matrix(unsigned r, unsigned c, const string & base_name, const string & tex_base_name);
@end example
The first two functions are @code{matrix} constructors which create a matrix
with @samp{r} rows and @samp{c} columns. The matrix elements can be
initialized from a (flat) list of expressions @samp{l}. Otherwise they are
all set to zero. The @code{lst_to_matrix()} function constructs a matrix
-from a list of lists, each list representing a matrix row. Finally,
-@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
-elements. Note that the last two functions return expressions, not matrix
-objects.
+from a list of lists, each list representing a matrix row. @code{diag_matrix()}
+constructs a diagonal matrix given the list of diagonal elements.
+@code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r} by @samp{c})
+unit matrix. And finally, @code{symbolic_matrix} constructs a matrix filled
+with newly generated symbols made of the specified base name and the
+position of each element in the matrix.
Matrix elements can be accessed and set using the parenthesis (function call)
operator:
the @code{op()} method. But C++-style subscripting with square brackets
@samp{[]} is not available.
-Here are a couple of examples that all construct the same 2x2 diagonal
-matrix:
+Here are a couple of examples of constructing matrices:
@example
@{
symbol a("a"), b("b");
- ex e;
matrix M(2, 2);
M(0, 0) = a;
M(1, 1) = b;
- e = M;
-
- e = matrix(2, 2, lst(a, 0, 0, b));
+ cout << M << endl;
+ // -> [[a,0],[0,b]]
- e = lst_to_matrix(lst(lst(a, 0), lst(0, b)));
+ cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
+ // -> [[a,0],[0,b]]
- e = diag_matrix(lst(a, b));
+ cout << lst_to_matrix(lst(lst(a, 0), lst(0, b))) << endl;
+ // -> [[a,0],[0,b]]
- cout << e << endl;
+ cout << diag_matrix(lst(a, b)) << endl;
// -> [[a,0],[0,b]]
+
+ cout << unit_matrix(3) << endl;
+ // -> [[1,0,0],[0,1,0],[0,0,1]]
+
+ cout << symbolic_matrix(2, 3, "x") << endl;
+ // -> [[x00,x01,x02],[x10,x11,x12]]
@}
@end example
The @code{matrix} class provides a couple of additional methods for
computing determinants, traces, and characteristic polynomials:
+@cindex @code{determinant()}
+@cindex @code{trace()}
+@cindex @code{charpoly()}
@example
ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
ex matrix::trace(void) const;
@subsection Expanding and collecting
@cindex @code{expand()}
@cindex @code{collect()}
+@cindex @code{collect_common_factors()}
A polynomial in one or more variables has many equivalent
representations. Some useful ones serve a specific purpose. Consider
(1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
@end example
+Polynomials can often be brought into a more compact form by collecting
+common factors from the terms of sums. This is accomplished by the function
+
+@example
+ex collect_common_factors(const ex & e);
+@end example
+
+This function doesn't perform a full factorization but only looks for
+factors which are already explicitly present:
+
+@example
+> collect_common_factors(a*x+a*y);
+(x+y)*a
+> collect_common_factors(a*x^2+2*a*x*y+a*y^2);
+a*(2*x*y+y^2+x^2)
+> collect_common_factors(a*(b*(a+c)*x+b*((a+c)*x+(a+c)*y)*y));
+(c+a)*a*(x*y+y^2+x)*b
+@end example
+
@subsection Degree and coefficients
@cindex @code{degree()}
@cindex @code{ldegree()}
symminfo(const ex & symmterm_, const ex & orig_)
{
- if (is_a<mul>(orig_)) {
+ if (is_a<mul>(orig_) && is_a<numeric>(orig_.op(orig_.nops()-1))) {
ex tmp = orig_.op(orig_.nops()-1);
orig = orig_ / tmp;
} else
orig = orig_;
- if (is_a<mul>(symmterm_)) {
+ if (is_a<mul>(symmterm_) && is_a<numeric>(symmterm_.op(symmterm_.nops()-1))) {
coeff = symmterm_.op(symmterm_.nops()-1);
symmterm = symmterm_ / coeff;
} else {
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
+#include <string>
#include <iostream>
+#include <sstream>
#include <algorithm>
#include <map>
#include <stdexcept>
m2[r*col+c] = m[r*col+c].eval(level);
return (new matrix(row, col, m2))->setflag(status_flags::dynallocated |
- status_flags::evaluated );
+ status_flags::evaluated);
}
ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) const
cols = l.op(i).nops();
// Allocate and fill matrix
- matrix &m = *new matrix(rows, cols);
- m.setflag(status_flags::dynallocated);
+ matrix &M = *new matrix(rows, cols);
+ M.setflag(status_flags::dynallocated);
for (i=0; i<rows; i++)
for (j=0; j<cols; j++)
if (l.op(i).nops() > j)
- m(i, j) = l.op(i).op(j);
+ M(i, j) = l.op(i).op(j);
else
- m(i, j) = _ex0;
- return m;
+ M(i, j) = _ex0;
+ return M;
}
ex diag_matrix(const lst & l)
return Id;
}
+ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name)
+{
+ matrix &M = *new matrix(r, c);
+ M.setflag(status_flags::dynallocated | status_flags::evaluated);
+
+ bool long_format = (r > 10 || c > 10);
+ bool single_row = (r == 1 || c == 1);
+
+ for (unsigned i=0; i<r; i++) {
+ for (unsigned j=0; j<c; j++) {
+ std::ostringstream s1, s2;
+ s1 << base_name;
+ s2 << tex_base_name << "_{";
+ if (single_row) {
+ if (c == 1) {
+ s1 << i;
+ s2 << i << '}';
+ } else {
+ s1 << j;
+ s2 << j << '}';
+ }
+ } else {
+ if (long_format) {
+ s1 << '_' << i << '_' << j;
+ s2 << i << ';' << j << "}";
+ } else {
+ s1 << i << j;
+ s2 << i << j << '}';
+ }
+ }
+ M(i, j) = symbol(s1.str(), s2.str());
+ }
+ }
+
+ return M;
+}
+
} // namespace GiNaC
#define __GINAC_MATRIX_H__
#include <vector>
+#include <string>
#include "basic.h"
#include "ex.h"
/** Convert list of diagonal elements to matrix. */
extern ex diag_matrix(const lst & l);
-/** Create a r times c unit matrix. */
+/** Create an r times c unit matrix. */
extern ex unit_matrix(unsigned r, unsigned c);
/** Create a x times x unit matrix. */
inline ex unit_matrix(unsigned x)
{ return unit_matrix(x, x); }
+/** Create an r times c matrix of newly generated symbols consisting of the
+ * given base name plus the numeric row/column position of each element.
+ * The base name for LaTeX output is specified separately. */
+extern ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const std::string & tex_base_name);
+
+/** Create an r times c matrix of newly generated symbols consisting of the
+ * given base name plus the numeric row/column position of each element. */
+inline ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name)
+{ return symbolic_matrix(r, c, base_name, base_name); }
+
} // namespace GiNaC
#endif // ndef __GINAC_MATRIX_H__
c.s << "(";
}
- bool first = true;
-
// First print the overall numeric coefficient
const numeric &coeff = ex_to<numeric>(overall_coeff);
if (coeff.csgn() == -1)
// Then proceed with the remaining factors
epvector::const_iterator it = seq.begin(), itend = seq.end();
- while (it != itend) {
- if (!first) {
- if (is_a<print_latex>(c))
- c.s << ' ';
+ if (is_a<print_latex>(c)) {
+
+ // Separate factors into those with negative numeric exponent
+ // and all others
+ exvector neg_powers, others;
+ while (it != itend) {
+ GINAC_ASSERT(is_a<numeric>(it->coeff));
+ if (ex_to<numeric>(it->coeff).is_negative())
+ neg_powers.push_back(recombine_pair_to_ex(expair(it->rest, -(it->coeff))));
else
- c.s << '*';
+ others.push_back(recombine_pair_to_ex(*it));
+ ++it;
+ }
+
+ if (!neg_powers.empty()) {
+
+ // Factors with negative exponent are printed as a fraction
+ c.s << "\\frac{";
+ mul(others).eval().print(c);
+ c.s << "}{";
+ mul(neg_powers).eval().print(c);
+ c.s << "}";
+
} else {
- first = false;
+
+ // All other factors are printed in the ordinary way
+ exvector::const_iterator vit = others.begin(), vitend = others.end();
+ while (vit != vitend) {
+ c.s << ' ';
+ vit->print(c, precedence());
+ ++vit;
+ }
+ }
+
+ } else {
+
+ bool first = true;
+ while (it != itend) {
+ if (!first)
+ c.s << '*';
+ else
+ first = false;
+ recombine_pair_to_ex(*it).print(c, precedence());
+ ++it;
}
- recombine_pair_to_ex(*it).print(c, precedence());
- ++it;
}
if (precedence() <= level) {
}
+/** 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)
+{
+ if (is_a<add>(e)) {
+
+ unsigned num = e.nops();
+ exvector terms; terms.reserve(num);
+ lst repl;
+ ex gc;
+
+ // Find the common GCD
+ for (unsigned i=0; i<num; i++) {
+ ex x = e.op(i).to_rational(repl);
+ if (is_a<add>(x) || is_a<mul>(x)) {
+ ex f = 1;
+ x = find_common_factor(x, f);
+ x *= f;
+ }
+
+ if (i == 0)
+ gc = x;
+ else
+ gc = gcd(gc, x);
+
+ terms.push_back(x);
+ }
+
+ if (gc.is_equal(_ex1))
+ return e;
+
+ // The GCD is the factor we pull out
+ factor *= gc.subs(repl);
+
+ // Now divide all terms by the GCD
+ for (unsigned i=0; i<num; i++) {
+ ex x;
+
+ // Try to avoid divide() because it expands the polynomial
+ ex &t = terms[i];
+ if (is_a<mul>(t)) {
+ for (unsigned j=0; j<t.nops(); j++) {
+ if (t.op(j).is_equal(gc)) {
+ exvector v; v.reserve(t.nops());
+ for (unsigned k=0; k<t.nops(); k++) {
+ if (k == j)
+ v.push_back(_ex1);
+ else
+ v.push_back(t.op(k));
+ }
+ t = (new mul(v))->setflag(status_flags::dynallocated).subs(repl);
+ goto term_done;
+ }
+ }
+ }
+
+ divide(t, gc, x);
+ t = x.subs(repl);
+term_done: ;
+ }
+ return (new add(terms))->setflag(status_flags::dynallocated);
+
+ } else if (is_a<mul>(e)) {
+
+ exvector v;
+ for (unsigned i=0; i<e.nops(); i++)
+ v.push_back(find_common_factor(e.op(i), factor));
+ return (new mul(v))->setflag(status_flags::dynallocated);
+
+ } else
+ return e;
+}
+
+
+/** Collect common factors in sums. This converts expressions like
+ * 'a*(b*x+b*y)' to 'a*b*(x+y)'. */
+ex collect_common_factors(const ex & e)
+{
+ if (is_a<add>(e) || is_a<mul>(e)) {
+
+ ex factor = 1;
+ ex r = find_common_factor(e, factor);
+ return factor * r;
+
+ } else
+ return e;
+}
+
+
/*
* Normal form of rational functions
*/
// Square-free partial fraction decomposition of a rational function a(x)
extern ex sqrfree_parfrac(const ex & a, const symbol & x);
+// Collect common factors in sums.
+extern ex collect_common_factors(const ex & e);
+
} // namespace GiNaC
#endif // ndef __GINAC_NORMAL_H__
bool is_tex = is_a<print_latex>(c);
- if (exponent.is_equal(_ex1_2)) {
+ if (is_tex && is_a<numeric>(exponent) && ex_to<numeric>(exponent).is_negative()) {
+
+ // Powers with negative numeric exponents are printed as fractions in TeX
+ c.s << "\\frac{1}{";
+ power(basis, -exponent).eval().print(c);
+ c.s << "}";
+
+ } else if (exponent.is_equal(_ex1_2)) {
+
+ // Square roots are printed in a special way
c.s << (is_tex ? "\\sqrt{" : "sqrt(");
basis.print(c);
c.s << (is_tex ? '}' : ')');
+
} else {
+
+ // Ordinary output of powers using '^' or '**'
if (precedence() <= level)
c.s << (is_tex ? "{(" : "(");
basis.print(c, precedence());
.BI collect_distributed( expression ", " list )
\- collects coefficients of like powers (result in distributed form)
.br
+.BI collect_common_factors( expression )
+\- collects common factors from the terms of sums
+.br
.BI content( expression ", " symbol )
\- content part of a polynomial
.br
// Table of names and descriptors of functions to be added
static const fcn_init extended_fcns[] = {
- {NULL, fcn_desc(f_dummy, 0)} // End marker
+ {NULL, f_dummy, 0} // End marker
};
// Table of help strings for functions
typedef ex (*fcnp2)(const exprseq &e, int serial);
struct fcn_desc {
- fcn_desc() : p(NULL), num_params(0) {}
- fcn_desc(fcnp func, int num) : p(func), num_params(num), is_ginac(false) {}
+ fcn_desc() : p(NULL), num_params(0), is_ginac(false), serial(0) {}
+ fcn_desc(fcnp func, int num) : p(func), num_params(num), is_ginac(false), serial(0) {}
fcn_desc(fcnp2 func, int num, int ser) : p((fcnp)func), num_params(num), is_ginac(true), serial(ser) {}
fcnp p; // Pointer to function
static ex f_collect(const exprseq &e) {return e[0].collect(e[1]);}
static ex f_collect_distributed(const exprseq &e) {return e[0].collect(e[1], true);}
+static ex f_collect_common_factors(const exprseq &e) {return collect_common_factors(e[0]);}
static ex f_degree(const exprseq &e) {return e[0].degree(e[1]);}
static ex f_denom(const exprseq &e) {return e[0].denom();}
static ex f_eval1(const exprseq &e) {return e[0].eval();}
// Tables for initializing the "fcns" map and the function help topics
struct fcn_init {
const char *name;
- const fcn_desc desc;
+ fcnp p;
+ int num_params;
};
static const fcn_init builtin_fcns[] = {
- {"charpoly", fcn_desc(f_charpoly, 2)},
- {"coeff", fcn_desc(f_coeff, 3)},
- {"collect", fcn_desc(f_collect, 2)},
- {"collect_distributed", fcn_desc(f_collect_distributed, 2)},
- {"content", fcn_desc(f_content, 2)},
- {"decomp_rational", fcn_desc(f_decomp_rational, 2)},
- {"degree", fcn_desc(f_degree, 2)},
- {"denom", fcn_desc(f_denom, 1)},
- {"determinant", fcn_desc(f_determinant, 1)},
- {"diag", fcn_desc(f_diag, 0)},
- {"diff", fcn_desc(f_diff2, 2)},
- {"diff", fcn_desc(f_diff3, 3)},
- {"divide", fcn_desc(f_divide, 2)},
- {"eval", fcn_desc(f_eval1, 1)},
- {"eval", fcn_desc(f_eval2, 2)},
- {"evalf", fcn_desc(f_evalf1, 1)},
- {"evalf", fcn_desc(f_evalf2, 2)},
- {"evalm", fcn_desc(f_evalm, 1)},
- {"expand", fcn_desc(f_expand, 1)},
- {"find", fcn_desc(f_find, 2)},
- {"gcd", fcn_desc(f_gcd, 2)},
- {"has", fcn_desc(f_has, 2)},
- {"inverse", fcn_desc(f_inverse, 1)},
- {"is", fcn_desc(f_is, 1)},
- {"lcm", fcn_desc(f_lcm, 2)},
- {"lcoeff", fcn_desc(f_lcoeff, 2)},
- {"ldegree", fcn_desc(f_ldegree, 2)},
- {"lsolve", fcn_desc(f_lsolve, 2)},
- {"map", fcn_desc(f_map, 2)},
- {"match", fcn_desc(f_match, 2)},
- {"nops", fcn_desc(f_nops, 1)},
- {"normal", fcn_desc(f_normal1, 1)},
- {"normal", fcn_desc(f_normal2, 2)},
- {"numer", fcn_desc(f_numer, 1)},
- {"numer_denom", fcn_desc(f_numer_denom, 1)},
- {"op", fcn_desc(f_op, 2)},
- {"pow", fcn_desc(f_pow, 2)},
- {"prem", fcn_desc(f_prem, 3)},
- {"primpart", fcn_desc(f_primpart, 2)},
- {"quo", fcn_desc(f_quo, 3)},
- {"rem", fcn_desc(f_rem, 3)},
- {"series", fcn_desc(f_series, 3)},
- {"sprem", fcn_desc(f_sprem, 3)},
- {"sqrfree", fcn_desc(f_sqrfree1, 1)},
- {"sqrfree", fcn_desc(f_sqrfree2, 2)},
- {"sqrt", fcn_desc(f_sqrt, 1)},
- {"subs", fcn_desc(f_subs2, 2)},
- {"subs", fcn_desc(f_subs3, 3)},
- {"tcoeff", fcn_desc(f_tcoeff, 2)},
- {"time", fcn_desc(f_dummy, 0)},
- {"trace", fcn_desc(f_trace, 1)},
- {"transpose", fcn_desc(f_transpose, 1)},
- {"unassign", fcn_desc(f_unassign, 1)},
- {"unit", fcn_desc(f_unit, 2)},
- {NULL, fcn_desc(f_dummy, 0)} // End marker
+ {"charpoly", f_charpoly, 2},
+ {"coeff", f_coeff, 3},
+ {"collect", f_collect, 2},
+ {"collect_common_factors", f_collect_common_factors, 1},
+ {"collect_distributed", f_collect_distributed, 2},
+ {"content", f_content, 2},
+ {"decomp_rational", f_decomp_rational, 2},
+ {"degree", f_degree, 2},
+ {"denom", f_denom, 1},
+ {"determinant", f_determinant, 1},
+ {"diag", f_diag, 0},
+ {"diff", f_diff2, 2},
+ {"diff", f_diff3, 3},
+ {"divide", f_divide, 2},
+ {"eval", f_eval1, 1},
+ {"eval", f_eval2, 2},
+ {"evalf", f_evalf1, 1},
+ {"evalf", f_evalf2, 2},
+ {"evalm", f_evalm, 1},
+ {"expand", f_expand, 1},
+ {"find", f_find, 2},
+ {"gcd", f_gcd, 2},
+ {"has", f_has, 2},
+ {"inverse", f_inverse, 1},
+ {"is", f_is, 1},
+ {"lcm", f_lcm, 2},
+ {"lcoeff", f_lcoeff, 2},
+ {"ldegree", f_ldegree, 2},
+ {"lsolve", f_lsolve, 2},
+ {"map", f_map, 2},
+ {"match", f_match, 2},
+ {"nops", f_nops, 1},
+ {"normal", f_normal1, 1},
+ {"normal", f_normal2, 2},
+ {"numer", f_numer, 1},
+ {"numer_denom", f_numer_denom, 1},
+ {"op", f_op, 2},
+ {"pow", f_pow, 2},
+ {"prem", f_prem, 3},
+ {"primpart", f_primpart, 2},
+ {"quo", f_quo, 3},
+ {"rem", f_rem, 3},
+ {"series", f_series, 3},
+ {"sprem", f_sprem, 3},
+ {"sqrfree", f_sqrfree1, 1},
+ {"sqrfree", f_sqrfree2, 2},
+ {"sqrt", f_sqrt, 1},
+ {"subs", f_subs2, 2},
+ {"subs", f_subs3, 3},
+ {"tcoeff", f_tcoeff, 2},
+ {"time", f_dummy, 0},
+ {"trace", f_trace, 1},
+ {"transpose", f_transpose, 1},
+ {"unassign", f_unassign, 1},
+ {"unit", f_unit, 2},
+ {NULL, f_dummy, 0} // End marker
};
struct fcn_help_init {
static void insert_fcns(const fcn_init *p)
{
while (p->name) {
- fcns.insert(make_pair(string(p->name), p->desc));
+ fcns.insert(make_pair(string(p->name), fcn_desc(p->p, p->num_params)));
p++;
}
}
git://www.ginac.de / ginac.git / commitdiff