From b42c582b86eb2c39b5c3d3c1a80dae44152de050 Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Thu, 21 Nov 2002 23:05:36 +0000 Subject: [PATCH] synced to 1.0 --- doc/tutorial/ginac.texi | 64 ++++++++++++++++----- ginac/indexed.cpp | 4 +- ginac/matrix.cpp | 51 +++++++++++++++-- ginac/matrix.h | 13 ++++- ginac/mul.cpp | 52 +++++++++++++---- ginac/normal.cpp | 90 +++++++++++++++++++++++++++++ ginac/normal.h | 3 + ginac/power.cpp | 14 ++++- ginsh/ginsh.1.in | 3 + ginsh/ginsh_extensions.h | 2 +- ginsh/ginsh_parser.yy | 121 ++++++++++++++++++++------------------- 11 files changed, 324 insertions(+), 93 deletions(-) diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index a642d656..113cf491 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -1363,21 +1363,31 @@ second one in the range 0@dots{}@math{n-1}. 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: @@ -1391,27 +1401,32 @@ It is also possible to access the matrix elements in a linear fashion with 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 @@ -1499,6 +1514,9 @@ general. 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; @@ -3415,6 +3433,7 @@ argument. You can not use functions like @samp{diff()}, @samp{op()}, @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 @@ -3473,6 +3492,25 @@ d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q (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()} diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 9cc01b47..7978327a 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -838,13 +838,13 @@ public: symminfo(const ex & symmterm_, const ex & orig_) { - if (is_a(orig_)) { + if (is_a(orig_) && is_a(orig_.op(orig_.nops()-1))) { ex tmp = orig_.op(orig_.nops()-1); orig = orig_ / tmp; } else orig = orig_; - if (is_a(symmterm_)) { + if (is_a(symmterm_) && is_a(symmterm_.op(symmterm_.nops()-1))) { coeff = symmterm_.op(symmterm_.nops()-1); symmterm = symmterm_ / coeff; } else { diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 1c316a35..33e38a63 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -20,7 +20,9 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ +#include #include +#include #include #include #include @@ -230,7 +232,7 @@ ex matrix::eval(int level) const 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 @@ -1455,15 +1457,15 @@ ex lst_to_matrix(const lst & l) 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 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) @@ -1486,4 +1488,41 @@ ex unit_matrix(unsigned r, unsigned c) 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 +#include #include "basic.h" #include "ex.h" @@ -144,13 +145,23 @@ extern ex lst_to_matrix(const lst & l); /** 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__ diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 34f2021a..aeb8c30e 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -193,8 +193,6 @@ void mul::print(const print_context & c, unsigned level) const c.s << "("; } - bool first = true; - // First print the overall numeric coefficient const numeric &coeff = ex_to(overall_coeff); if (coeff.csgn() == -1) @@ -220,17 +218,51 @@ void mul::print(const print_context & c, unsigned level) const // Then proceed with the remaining factors epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { - if (!first) { - if (is_a(c)) - c.s << ' '; + if (is_a(c)) { + + // Separate factors into those with negative numeric exponent + // and all others + exvector neg_powers, others; + while (it != itend) { + GINAC_ASSERT(is_a(it->coeff)); + if (ex_to(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) { diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 063ad69a..0b062a01 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1913,6 +1913,96 @@ ex sqrfree_parfrac(const ex & a, const symbol & x) } +/** 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(e)) { + + unsigned num = e.nops(); + exvector terms; terms.reserve(num); + lst repl; + ex gc; + + // Find the common GCD + for (unsigned i=0; i(x) || is_a(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(t)) { + for (unsigned j=0; jsetflag(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(e)) { + + exvector v; + for (unsigned i=0; isetflag(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(e) || is_a(e)) { + + ex factor = 1; + ex r = find_common_factor(e, factor); + return factor * r; + + } else + return e; +} + + /* * Normal form of rational functions */ diff --git a/ginac/normal.h b/ginac/normal.h index be67d8eb..9196f965 100644 --- a/ginac/normal.h +++ b/ginac/normal.h @@ -63,6 +63,9 @@ extern ex sqrfree(const ex &a, const lst &l = lst()); // 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__ diff --git a/ginac/power.cpp b/ginac/power.cpp index a2c8c1af..86c684c3 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -173,11 +173,23 @@ void power::print(const print_context & c, unsigned level) const bool is_tex = is_a(c); - if (exponent.is_equal(_ex1_2)) { + if (is_tex && is_a(exponent) && ex_to(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()); diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in index e5202138..f1615768 100644 --- a/ginsh/ginsh.1.in +++ b/ginsh/ginsh.1.in @@ -240,6 +240,9 @@ detail here. Please refer to the GiNaC documentation. .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 diff --git a/ginsh/ginsh_extensions.h b/ginsh/ginsh_extensions.h index 55c17ac5..21cd2d1f 100644 --- a/ginsh/ginsh_extensions.h +++ b/ginsh/ginsh_extensions.h @@ -24,7 +24,7 @@ // 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 diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index fb644351..8685f906 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -63,8 +63,8 @@ typedef ex (*fcnp)(const exprseq &e); 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 @@ -294,6 +294,7 @@ static void push(const ex &e) 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();} @@ -531,65 +532,67 @@ static ex f_dummy(const exprseq &e) // 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 { @@ -638,7 +641,7 @@ static const fcn_help_init builtin_help[] = { 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++; } } -- 2.46.1