X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=bd276c56e454f035eb6043d1c9bf8b30ce769f72;hp=0d24a04f950f312c54934ae5b8d2a0e994e7e117;hb=fa13c6831b809abdae9eb2c17b33d2eff8e4878c;hpb=68fdf425abf14d016d5f95ee7b9d06a19a3c5926

diff --git a/ginac/power.cpp b/ginac/power.cpp
index 0d24a04f..bd276c56 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
 
 /*
- *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2005 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
@@ -17,7 +17,7 @@
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include <vector>
@@ -38,13 +38,17 @@
 #include "indexed.h"
 #include "symbol.h"
 #include "lst.h"
-#include "print.h"
 #include "archive.h"
 #include "utils.h"
 
 namespace GiNaC {
 
-GINAC_IMPLEMENT_REGISTERED_CLASS(power, basic)
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic,
+  print_func<print_dflt>(&power::do_print_dflt).
+  print_func<print_latex>(&power::do_print_latex).
+  print_func<print_csrc>(&power::do_print_csrc).
+  print_func<print_python>(&power::do_print_python).
+  print_func<print_python_repr>(&power::do_print_python_repr))
 
 typedef std::vector<int> intvector;
 
@@ -85,11 +89,58 @@ DEFAULT_UNARCHIVE(power)
 
 // public
 
+void power::print_power(const print_context & c, const char *powersymbol, const char *openbrace, const char *closebrace, unsigned level) const
+{
+	// Ordinary output of powers using '^' or '**'
+	if (precedence() <= level)
+		c.s << openbrace << '(';
+	basis.print(c, precedence());
+	c.s << powersymbol;
+	c.s << openbrace;
+	exponent.print(c, precedence());
+	c.s << closebrace;
+	if (precedence() <= level)
+		c.s << ')' << closebrace;
+}
+
+void power::do_print_dflt(const print_dflt & c, unsigned level) const
+{
+	if (exponent.is_equal(_ex1_2)) {
+
+		// Square roots are printed in a special way
+		c.s << "sqrt(";
+		basis.print(c);
+		c.s << ')';
+
+	} else
+		print_power(c, "^", "", "", level);
+}
+
+void power::do_print_latex(const print_latex & c, unsigned level) const
+{
+	if (is_exactly_a<numeric>(exponent) && ex_to<numeric>(exponent).is_negative()) {
+
+		// Powers with negative numeric exponents are printed as fractions
+		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 << "\\sqrt{";
+		basis.print(c);
+		c.s << '}';
+
+	} else
+		print_power(c, "^", "{", "}", level);
+}
+
 static void print_sym_pow(const print_context & c, const symbol &x, int exp)
 {
 	// Optimal output of integer powers of symbols to aid compiler CSE.
 	// C.f. ISO/IEC 14882:1998, section 1.9 [intro execution], paragraph 15
-	// to learn why such a parenthisation is really necessary.
+	// to learn why such a parenthesation is really necessary.
 	if (exp == 1) {
 		x.print(c);
 	} else if (exp == 2) {
@@ -109,96 +160,58 @@ static void print_sym_pow(const print_context & c, const symbol &x, int exp)
 	}
 }
 
-void power::print(const print_context & c, unsigned level) const
+void power::do_print_csrc(const print_csrc & c, unsigned level) const
 {
-	if (is_a<print_tree>(c)) {
-
-		inherited::print(c, level);
-
-	} else if (is_a<print_csrc>(c)) {
-
-		// Integer powers of symbols are printed in a special, optimized way
-		if (exponent.info(info_flags::integer)
-		 && (is_a<symbol>(basis) || is_a<constant>(basis))) {
-			int exp = ex_to<numeric>(exponent).to_int();
-			if (exp > 0)
-				c.s << '(';
-			else {
-				exp = -exp;
-				if (is_a<print_csrc_cl_N>(c))
-					c.s << "recip(";
-				else
-					c.s << "1.0/(";
-			}
-			print_sym_pow(c, ex_to<symbol>(basis), exp);
-			c.s << ')';
-
-		// <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
-		} else if (exponent.is_equal(_ex_1)) {
+	// Integer powers of symbols are printed in a special, optimized way
+	if (exponent.info(info_flags::integer)
+	 && (is_a<symbol>(basis) || is_a<constant>(basis))) {
+		int exp = ex_to<numeric>(exponent).to_int();
+		if (exp > 0)
+			c.s << '(';
+		else {
+			exp = -exp;
 			if (is_a<print_csrc_cl_N>(c))
 				c.s << "recip(";
 			else
 				c.s << "1.0/(";
-			basis.print(c);
-			c.s << ')';
-
-		// Otherwise, use the pow() or expt() (CLN) functions
-		} else {
-			if (is_a<print_csrc_cl_N>(c))
-				c.s << "expt(";
-			else
-				c.s << "pow(";
-			basis.print(c);
-			c.s << ',';
-			exponent.print(c);
-			c.s << ')';
 		}
+		print_sym_pow(c, ex_to<symbol>(basis), exp);
+		c.s << ')';
 
-	} else if (is_a<print_python_repr>(c)) {
+	// <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
+	} else if (exponent.is_equal(_ex_1)) {
+		if (is_a<print_csrc_cl_N>(c))
+			c.s << "recip(";
+		else
+			c.s << "1.0/(";
+		basis.print(c);
+		c.s << ')';
 
-		c.s << class_name() << '(';
+	// Otherwise, use the pow() or expt() (CLN) functions
+	} else {
+		if (is_a<print_csrc_cl_N>(c))
+			c.s << "expt(";
+		else
+			c.s << "pow(";
 		basis.print(c);
 		c.s << ',';
 		exponent.print(c);
 		c.s << ')';
+	}
+}
 
-	} else {
-
-		bool is_tex = is_a<print_latex>(c);
-
-		if (is_tex && is_exactly_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 {
+void power::do_print_python(const print_python & c, unsigned level) const
+{
+	print_power(c, "**", "", "", level);
+}
 
-			// Ordinary output of powers using '^' or '**'
-			if (precedence() <= level)
-				c.s << (is_tex ? "{(" : "(");
-			basis.print(c, precedence());
-			if (is_a<print_python>(c))
-				c.s << "**";
-			else
-				c.s << '^';
-			if (is_tex)
-				c.s << '{';
-			exponent.print(c, precedence());
-			if (is_tex)
-				c.s << '}';
-			if (precedence() <= level)
-				c.s << (is_tex ? ")}" : ")");
-		}
-	}
+void power::do_print_python_repr(const print_python_repr & c, unsigned level) const
+{
+	c.s << class_name() << '(';
+	basis.print(c);
+	c.s << ',';
+	exponent.print(c);
+	c.s << ')';
 }
 
 bool power::info(unsigned inf) const
@@ -209,12 +222,14 @@ bool power::info(unsigned inf) const
 		case info_flags::cinteger_polynomial:
 		case info_flags::rational_polynomial:
 		case info_flags::crational_polynomial:
-			return exponent.info(info_flags::nonnegint);
+			return exponent.info(info_flags::nonnegint) &&
+			       basis.info(inf);
 		case info_flags::rational_function:
-			return exponent.info(info_flags::integer);
+			return exponent.info(info_flags::integer) &&
+			       basis.info(inf);
 		case info_flags::algebraic:
-			return (!exponent.info(info_flags::integer) ||
-					basis.info(inf));
+			return !exponent.info(info_flags::integer) ||
+			       basis.info(inf);
 	}
 	return inherited::info(inf);
 }
@@ -233,7 +248,14 @@ ex power::op(size_t i) const
 
 ex power::map(map_function & f) const
 {
-	return (new power(f(basis), f(exponent)))->setflag(status_flags::dynallocated);
+	const ex &mapped_basis = f(basis);
+	const ex &mapped_exponent = f(exponent);
+
+	if (!are_ex_trivially_equal(basis, mapped_basis)
+	 || !are_ex_trivially_equal(exponent, mapped_exponent))
+		return (new power(mapped_basis, mapped_exponent))->setflag(status_flags::dynallocated);
+	else
+		return *this;
 }
 
 int power::degree(const ex & s) const
@@ -433,7 +455,7 @@ ex power::eval(int level) const
 	
 		// ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
 		if (num_exponent->is_integer() && is_exactly_a<mul>(ebasis)) {
-			return expand_mul(ex_to<mul>(ebasis), *num_exponent);
+			return expand_mul(ex_to<mul>(ebasis), *num_exponent, 0);
 		}
 	
 		// ^(*(...,x;c1),c2) -> *(^(*(...,x;1),c2),c1^c2)  (c1, c2 numeric(), c1>0)
@@ -518,27 +540,26 @@ ex power::evalm() const
 // from mul.cpp
 extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
 
-ex power::subs(const lst & ls, const lst & lr, unsigned options) const
+ex power::subs(const exmap & m, unsigned options) const
 {	
-	const ex &subsed_basis = basis.subs(ls, lr, options);
-	const ex &subsed_exponent = exponent.subs(ls, lr, options);
+	const ex &subsed_basis = basis.subs(m, options);
+	const ex &subsed_exponent = exponent.subs(m, options);
 
 	if (!are_ex_trivially_equal(basis, subsed_basis)
 	 || !are_ex_trivially_equal(exponent, subsed_exponent)) 
-		return power(subsed_basis, subsed_exponent).subs_one_level(ls, lr, options);
+		return power(subsed_basis, subsed_exponent).subs_one_level(m, options);
 
-	if (!(options & subs_options::subs_algebraic))
-		return subs_one_level(ls, lr, options);
+	if (!(options & subs_options::algebraic))
+		return subs_one_level(m, options);
 
-	lst::const_iterator its, itr;
-	for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+	for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
 		int nummatches = std::numeric_limits<int>::max();
 		lst repls;
-		if (tryfactsubs(*this, *its, nummatches, repls))
-			return (ex_to<basic>((*this) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(ls, lr, options);
+		if (tryfactsubs(*this, it->first, nummatches, repls))
+			return (ex_to<basic>((*this) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches))).subs_one_level(m, options);
 	}
 
-	return subs_one_level(ls, lr, options);
+	return subs_one_level(m, options);
 }
 
 ex power::eval_ncmul(const exvector & v) const
@@ -546,6 +567,16 @@ ex power::eval_ncmul(const exvector & v) const
 	return inherited::eval_ncmul(v);
 }
 
+ex power::conjugate() const
+{
+	ex newbasis = basis.conjugate();
+	ex newexponent = exponent.conjugate();
+	if (are_ex_trivially_equal(basis, newbasis) && are_ex_trivially_equal(exponent, newexponent)) {
+		return *this;
+	}
+	return (new power(newbasis, newexponent))->setflag(status_flags::dynallocated);
+}
+
 // protected
 
 /** Implementation of ex::diff() for a power.
@@ -583,7 +614,7 @@ unsigned power::return_type() const
 {
 	return basis.return_type();
 }
-   
+
 unsigned power::return_type_tinfo() const
 {
 	return basis.return_type_tinfo();
@@ -614,7 +645,7 @@ ex power::expand(unsigned options) const
 			const numeric &num_exponent = ex_to<numeric>(a.overall_coeff);
 			int int_exponent = num_exponent.to_int();
 			if (int_exponent > 0 && is_exactly_a<add>(expanded_basis))
-				distrseq.push_back(expand_add(ex_to<add>(expanded_basis), int_exponent));
+				distrseq.push_back(expand_add(ex_to<add>(expanded_basis), int_exponent, options));
 			else
 				distrseq.push_back(power(expanded_basis, a.overall_coeff));
 		} else
@@ -622,7 +653,7 @@ ex power::expand(unsigned options) const
 		
 		// Make sure that e.g. (x+y)^(1+a) -> x*(x+y)^a + y*(x+y)^a
 		ex r = (new mul(distrseq))->setflag(status_flags::dynallocated);
-		return r.expand();
+		return r.expand(options);
 	}
 	
 	if (!is_exactly_a<numeric>(expanded_exponent) ||
@@ -640,11 +671,11 @@ ex power::expand(unsigned options) const
 	
 	// (x+y)^n, n>0
 	if (int_exponent > 0 && is_exactly_a<add>(expanded_basis))
-		return expand_add(ex_to<add>(expanded_basis), int_exponent);
+		return expand_add(ex_to<add>(expanded_basis), int_exponent, options);
 	
 	// (x*y)^n -> x^n * y^n
 	if (is_exactly_a<mul>(expanded_basis))
-		return expand_mul(ex_to<mul>(expanded_basis), num_exponent);
+		return expand_mul(ex_to<mul>(expanded_basis), num_exponent, options, true);
 	
 	// cannot expand further
 	if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent))
@@ -665,10 +696,10 @@ ex power::expand(unsigned options) const
 
 /** expand a^n where a is an add and n is a positive integer.
  *  @see power::expand */
-ex power::expand_add(const add & a, int n) const
+ex power::expand_add(const add & a, int n, unsigned options) const
 {
 	if (n==2)
-		return expand_add_2(a);
+		return expand_add_2(a, options);
 
 	const size_t m = a.nops();
 	exvector result;
@@ -701,7 +732,7 @@ ex power::expand_add(const add & a, int n) const
 			             !is_exactly_a<mul>(ex_to<power>(b).basis) ||
 			             !is_exactly_a<power>(ex_to<power>(b).basis));
 			if (is_exactly_a<mul>(b))
-				term.push_back(expand_mul(ex_to<mul>(b),numeric(k[l])));
+				term.push_back(expand_mul(ex_to<mul>(b), numeric(k[l]), options, true));
 			else
 				term.push_back(power(b,k[l]));
 		}
@@ -715,7 +746,7 @@ ex power::expand_add(const add & a, int n) const
 		             !is_exactly_a<mul>(ex_to<power>(b).basis) ||
 		             !is_exactly_a<power>(ex_to<power>(b).basis));
 		if (is_exactly_a<mul>(b))
-			term.push_back(expand_mul(ex_to<mul>(b),numeric(n-k_cum[m-2])));
+			term.push_back(expand_mul(ex_to<mul>(b), numeric(n-k_cum[m-2]), options, true));
 		else
 			term.push_back(power(b,n-k_cum[m-2]));
 
@@ -725,7 +756,7 @@ ex power::expand_add(const add & a, int n) const
 
 		term.push_back(f);
 
-		result.push_back((new mul(term))->setflag(status_flags::dynallocated));
+		result.push_back(ex((new mul(term))->setflag(status_flags::dynallocated)).expand(options));
 
 		// increment k[]
 		l = m-2;
@@ -752,7 +783,7 @@ ex power::expand_add(const add & a, int n) const
 
 /** Special case of power::expand_add. Expands a^2 where a is an add.
  *  @see power::expand_add */
-ex power::expand_add_2(const add & a) const
+ex power::expand_add_2(const add & a, unsigned options) const
 {
 	epvector sum;
 	size_t a_nops = a.nops();
@@ -775,7 +806,7 @@ ex power::expand_add_2(const add & a) const
 		
 		if (c.is_equal(_ex1)) {
 			if (is_exactly_a<mul>(r)) {
-				sum.push_back(expair(expand_mul(ex_to<mul>(r),_num2),
+				sum.push_back(expair(expand_mul(ex_to<mul>(r), _num2, options, true),
 				                     _ex1));
 			} else {
 				sum.push_back(expair((new power(r,_ex2))->setflag(status_flags::dynallocated),
@@ -783,7 +814,7 @@ ex power::expand_add_2(const add & a) const
 			}
 		} else {
 			if (is_exactly_a<mul>(r)) {
-				sum.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r),_num2),
+				sum.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r), _num2, options, true),
 				                     ex_to<numeric>(c).power_dyn(_num2)));
 			} else {
 				sum.push_back(a.combine_ex_with_coeff_to_pair((new power(r,_ex2))->setflag(status_flags::dynallocated),
@@ -818,7 +849,7 @@ ex power::expand_add_2(const add & a) const
 
 /** Expand factors of m in m^n where m is a mul and n is and integer.
  *  @see power::expand */
-ex power::expand_mul(const mul & m, const numeric & n) const
+ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool from_expand) const
 {
 	GINAC_ASSERT(n.is_integer());
 
@@ -827,6 +858,8 @@ ex power::expand_mul(const mul & m, const numeric & n) const
 
 	epvector distrseq;
 	distrseq.reserve(m.seq.size());
+	bool need_reexpand = false;
+
 	epvector::const_iterator last = m.seq.end();
 	epvector::const_iterator cit = m.seq.begin();
 	while (cit!=last) {
@@ -835,11 +868,23 @@ ex power::expand_mul(const mul & m, const numeric & n) const
 		} else {
 			// it is safe not to call mul::combine_pair_with_coeff_to_pair()
 			// since n is an integer
-			distrseq.push_back(expair(cit->rest, ex_to<numeric>(cit->coeff).mul(n)));
+			numeric new_coeff = ex_to<numeric>(cit->coeff).mul(n);
+			if (from_expand && is_exactly_a<add>(cit->rest) && new_coeff.is_pos_integer()) {
+				// this happens when e.g. (a+b)^(1/2) gets squared and
+				// the resulting product needs to be reexpanded
+				need_reexpand = true;
+			}
+			distrseq.push_back(expair(cit->rest, new_coeff));
 		}
 		++cit;
 	}
-	return (new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated);
+
+	const mul & result = static_cast<const mul &>((new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated));
+	if (need_reexpand)
+		return ex(result).expand(options);
+	if (from_expand)
+		return result.setflag(status_flags::expanded);
+	return result;
 }
 
 } // namespace GiNaC