From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Fri, 1 Jun 2001 22:50:00 +0000 (+0000)
Subject: - in the output, list delimiters are now { } and matrix delimiters are [ ]
X-Git-Tag: release_0-9-0~20
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=c8ba9c6cf819792cbf88d25b324406b67d5cc49a

- in the output, list delimiters are now { } and matrix delimiters are [ ]
- added evalm() method for evaluating sums and products of matrices
- added map() method for applying a function to subexpressions
- power(noncommut, posint) automatically expands the product, removed ncpow()
---

diff --git a/ginac/add.cpp b/ginac/add.cpp
index 92cb325f..17faae2c 100644
--- a/ginac/add.cpp
+++ b/ginac/add.cpp
@@ -25,6 +25,7 @@
 
 #include "add.h"
 #include "mul.h"
+#include "matrix.h"
 #include "archive.h"
 #include "debugmsg.h"
 #include "utils.h"
@@ -336,6 +337,38 @@ ex add::eval(int level) const
 	return this->hold();
 }
 
+ex add::evalm(void) const
+{
+	// Evaluate children first and add up all matrices. Stop if there's one
+	// term that is not a matrix.
+	epvector *s = new epvector;
+	s->reserve(seq.size());
+
+	bool all_matrices = true;
+	bool first_term = true;
+	matrix sum;
+
+	epvector::const_iterator it = seq.begin(), itend = seq.end();
+	while (it != itend) {
+		const ex &m = recombine_pair_to_ex(*it).evalm();
+		s->push_back(split_ex_to_pair(m));
+		if (is_ex_of_type(m, matrix)) {
+			if (first_term) {
+				sum = ex_to_matrix(m);
+				first_term = false;
+			} else
+				sum = sum.add(ex_to_matrix(m));
+		} else
+			all_matrices = false;
+		it++;
+	}
+
+	if (all_matrices)
+		return sum + overall_coeff;
+	else
+		return (new add(s, overall_coeff))->setflag(status_flags::dynallocated);
+}
+
 ex add::simplify_ncmul(const exvector & v) const
 {
 	if (seq.size()==0) {
diff --git a/ginac/add.h b/ginac/add.h
index 1e036627..eae04860 100644
--- a/ginac/add.h
+++ b/ginac/add.h
@@ -53,6 +53,7 @@ public:
 	int ldegree(const ex & s) const;
 	ex coeff(const ex & s, int n=1) const;
 	ex eval(int level=0) const;
+	ex evalm(void) const;
 	ex series(const relational & r, int order, unsigned options = 0) const;
 	ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
 	numeric integer_content(void) const;
diff --git a/ginac/basic.cpp b/ginac/basic.cpp
index 90b9c145..1990ebcb 100644
--- a/ginac/basic.cpp
+++ b/ginac/basic.cpp
@@ -221,16 +221,32 @@ bool basic::has(const ex & other) const
 {
 	GINAC_ASSERT(other.bp!=0);
 	lst repl_lst;
-	if (match(*other.bp, repl_lst)) return true;
-	if (nops()>0) {
-		for (unsigned i=0; i<nops(); i++)
-			if (op(i).has(other))
-				return true;
-	}
+	if (match(*other.bp, repl_lst))
+		return true;
+	for (unsigned i=0; i<nops(); i++)
+		if (op(i).has(other))
+			return true;
 	
 	return false;
 }
 
+/** Construct new expression by applying the specified function to all
+ *  sub-expressions. */
+ex basic::map(map_func f) const
+{
+	unsigned num = nops();
+	if (num == 0)
+		return *this;
+
+	basic *copy = duplicate();
+	copy->setflag(status_flags::dynallocated);
+	copy->clearflag(status_flags::hash_calculated);
+	ex e(*copy);
+	for (unsigned i=0; i<num; i++)
+		e.let_op(i) = f(e.op(i));
+	return e.eval();
+}
+
 /** Return degree of highest power in object s. */
 int basic::degree(const ex & s) const
 {
@@ -347,6 +363,16 @@ ex basic::evalf(int level) const
 	return *this;
 }
 
+/** Evaluate sums and products of matrices. */
+ex basic::evalm(void) const
+{
+	unsigned num = nops();
+	if (num == 0)
+		return *this;
+	else
+		return map(GiNaC::evalm);
+}
+
 /** Perform automatic symbolic evaluations on indexed expression that
  *  contains this object as the base expression. */
 ex basic::eval_indexed(const basic & i) const
diff --git a/ginac/basic.h b/ginac/basic.h
index 446d6d10..1b5f2cf3 100644
--- a/ginac/basic.h
+++ b/ginac/basic.h
@@ -54,6 +54,8 @@ class print_context;
   typedef std::vector<ex> exvector;
 #endif
 
+typedef ex (*map_func)(const ex & e);
+
 /** This class is the ABC (abstract base class) of GiNaC's class hierarchy.
  *  It is responsible for the reference counting. */
 class basic
@@ -111,12 +113,14 @@ public: // only const functions please (may break reference counting)
 	virtual ex operator[](const ex & index) const;
 	virtual ex operator[](int i) const;
 	virtual bool has(const ex & other) const;
+	virtual ex map(map_func f) const;
 	virtual int degree(const ex & s) const;
 	virtual int ldegree(const ex & s) const;
 	virtual ex coeff(const ex & s, int n = 1) const;
 	virtual ex collect(const ex & s, bool distributed = false) const;
 	virtual ex eval(int level = 0) const;
 	virtual ex evalf(int level = 0) const;
+	virtual ex evalm(void) const;
 	virtual ex series(const relational & r, int order, unsigned options = 0) const;
 	virtual bool match(const ex & pattern, lst & repl_lst) const;
 	virtual ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const;
diff --git a/ginac/container.pl b/ginac/container.pl
index f13c9657..549f066e 100755
--- a/ginac/container.pl
+++ b/ginac/container.pl
@@ -35,8 +35,8 @@ if ($type eq 'exprseq') {
 	$reserve=0;
 	$prepend=1;
 	$let_op=1;
-	$open_bracket='[';
-	$close_bracket=']';
+	$open_bracket='{';
+	$close_bracket='}';
 
 } else {
 	die "invalid type $type";
@@ -210,6 +210,7 @@ public:
 	bool info(unsigned inf) const;
 	unsigned nops() const;
 	ex & let_op(int i);
+	ex map(map_func f) const;
 	ex expand(unsigned options=0) const;
 	ex eval(int level=0) const;
 	ex evalf(int level=0) const;
@@ -430,6 +431,17 @@ unsigned ${CONTAINER}::nops() const
 
 ${LET_OP_IMPLEMENTATION}
 
+ex ${CONTAINER}::map(map_func f) const
+{
+	${STLT} s;
+	RESERVE(s,seq.size());
+	for (${STLT}::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
+		s.push_back(f(*it));
+	}
+
+	return this${CONTAINER}(s);
+}
+
 ex ${CONTAINER}::expand(unsigned options) const
 {
 	${STLT} s;
diff --git a/ginac/ex.h b/ginac/ex.h
index 499229e0..84b05d88 100644
--- a/ginac/ex.h
+++ b/ginac/ex.h
@@ -84,6 +84,7 @@ public:
 	unsigned nops() const { return bp->nops(); }
 	ex expand(unsigned options=0) const;
 	bool has(const ex & other) const { return bp->has(other); }
+	ex map(map_func f) const { return bp->map(f); }
 	int degree(const ex & s) const { return bp->degree(s); }
 	int ldegree(const ex & s) const { return bp->ldegree(s); }
 	ex coeff(const ex & s, int n = 1) const { return bp->coeff(s, n); }
@@ -104,6 +105,7 @@ public:
 	ex collect(const ex & s, bool distributed = false) const { return bp->collect(s, distributed); }
 	ex eval(int level = 0) const { return bp->eval(level); }
 	ex evalf(int level = 0) const { return bp->evalf(level); }
+	ex evalm(void) const { return bp->evalm(); }
 	ex diff(const symbol & s, unsigned nth = 1) const;
 	ex series(const ex & r, int order, unsigned options = 0) const;
 	bool match(const ex & pattern) const;
@@ -374,6 +376,9 @@ inline ex eval(const ex & thisex, int level = 0)
 inline ex evalf(const ex & thisex, int level = 0)
 { return thisex.evalf(level); }
 
+inline ex evalm(const ex & thisex)
+{ return thisex.evalm(); }
+
 inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
 { return thisex.diff(s, nth); }
 
diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp
index 9421570a..22669e55 100644
--- a/ginac/expairseq.cpp
+++ b/ginac/expairseq.cpp
@@ -299,6 +299,20 @@ ex &expairseq::let_op(int i)
 	throw(std::logic_error("let_op not defined for expairseq and derived classes (add,mul,...)"));
 }
 
+ex expairseq::map(map_func f) const
+{
+	epvector *v = new epvector;
+	v->reserve(seq.size());
+
+	epvector::const_iterator cit = seq.begin(), last = seq.end();
+	while (cit != last) {
+		v->push_back(split_ex_to_pair(f(recombine_pair_to_ex(*cit))));
+		cit++;
+	}
+
+	return thisexpairseq(v, f(overall_coeff));
+}
+
 ex expairseq::eval(int level) const
 {
 	if ((level==1) && (flags &status_flags::evaluated))
diff --git a/ginac/expairseq.h b/ginac/expairseq.h
index 3851a578..be5e014c 100644
--- a/ginac/expairseq.h
+++ b/ginac/expairseq.h
@@ -93,6 +93,7 @@ public:
 	unsigned nops() const;
 	ex op(int i) const;
 	ex & let_op(int i);
+	ex map(map_func f) const;
 	ex eval(int level=0) const;
 	ex evalf(int level=0) const;
 	ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 7081eb5b..de5522ad 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -509,20 +509,6 @@ ex lsolve(const ex &eqns, const ex &symbols)
 	return sollist;
 }
 
-/** non-commutative power. */
-ex ncpow(const ex & basis, unsigned exponent)
-{
-	if (exponent == 0)
-		return _ex1();
-
-	exvector v;
-	v.reserve(exponent);
-	for (unsigned i=0; i<exponent; ++i)
-		v.push_back(basis);
-
-	return ncmul(v, true);
-}
-
 // Symmetrize/antisymmetrize over a vector of objects
 static ex symm(const ex & e, exvector::const_iterator first, exvector::const_iterator last, bool asymmetric)
 {
diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 96c25607..3d170ff8 100644
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -133,9 +133,6 @@ DECLARE_FUNCTION_2P(Derivative)
 
 ex lsolve(const ex &eqns, const ex &symbols);
 
-/** Power of non-commutative basis. */
-ex ncpow(const ex & basis, unsigned exponent);
-
 /** Symmetrize expression over a set of objects (symbols, indices). */
 ex symmetrize(const ex & e, exvector::const_iterator first, exvector::const_iterator last);
 
diff --git a/ginac/input_lexer.ll b/ginac/input_lexer.ll
index 6c22ec08..46df9595 100644
--- a/ginac/input_lexer.ll
+++ b/ginac/input_lexer.ll
@@ -101,10 +101,6 @@ Digits			ginac_yylval = (long)Digits; return T_DIGITS;
 "<="			return T_LESSEQ;
 ">="			return T_GREATEREQ;
 
-			/* matrix delimiters */
-\[\[			return T_MATRIX_BEGIN;
-\]\]			return T_MATRIX_END;
-
 			/* numbers */
 {D}+			|
 "#"{D}+"R"{AN}+		|
diff --git a/ginac/input_parser.yy b/ginac/input_parser.yy
index 9eab1a59..5e29e482 100644
--- a/ginac/input_parser.yy
+++ b/ginac/input_parser.yy
@@ -52,7 +52,7 @@ static std::string parser_error;
 
 /* Tokens (T_LITERAL means a literal value returned by the parser, but not
    of class numeric or symbol (e.g. a constant or the FAIL object)) */
-%token T_NUMBER T_SYMBOL T_LITERAL T_DIGITS T_EQUAL T_NOTEQ T_LESSEQ T_GREATEREQ T_MATRIX_BEGIN T_MATRIX_END
+%token T_NUMBER T_SYMBOL T_LITERAL T_DIGITS T_EQUAL T_NOTEQ T_LESSEQ T_GREATEREQ
 
 /* Operator precedence and associativity */
 %right '='
@@ -112,8 +112,8 @@ exp	: T_NUMBER		{$$ = $1;}
 	| exp '^' exp		{$$ = pow($1, $3);}
 	| exp '!'		{$$ = factorial($1);}
 	| '(' exp ')'		{$$ = $2;}
-	| '[' list_or_empty ']'	{$$ = $2;}
-	| T_MATRIX_BEGIN matrix T_MATRIX_END	{$$ = lst_to_matrix(ex_to_lst($2));}
+	| '{' list_or_empty '}'	{$$ = $2;}
+	| '[' matrix ']'	{$$ = lst_to_matrix(ex_to_lst($2));}
 	;
 
 exprseq	: exp			{$$ = exprseq($1);}
@@ -128,8 +128,8 @@ list	: exp			{$$ = lst($1);}
 	| list ',' exp		{lst l(static_cast<lst &>(*($1.bp))); $$ = l.append($3);}
 	;
 
-matrix	: T_MATRIX_BEGIN row T_MATRIX_END		{$$ = lst($2);}
-	| matrix ',' T_MATRIX_BEGIN row	T_MATRIX_END	{lst l(static_cast<lst &>(*($1.bp))); $$ = l.append($4);}
+matrix	: '[' row ']'		{$$ = lst($2);}
+	| matrix ',' '[' row ']' {lst l(static_cast<lst &>(*($1.bp))); $$ = l.append($4);}
 	;
 
 row	: exp			{$$ = lst($1);}
diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index d47c89e2..5767d4cf 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -156,23 +156,23 @@ void matrix::print(const print_context & c, unsigned level) const
 
 	} else {
 
-		c.s << "[[ ";
+		c.s << "[";
 		for (unsigned y=0; y<row-1; ++y) {
-			c.s << "[[";
+			c.s << "[";
 			for (unsigned x=0; x<col-1; ++x) {
 				m[y*col+x].print(c);
 				c.s << ",";
 			}
 			m[col*(y+1)-1].print(c);
-			c.s << "]], ";
+			c.s << "],";
 		}
-		c.s << "[[";
+		c.s << "[";
 		for (unsigned x=0; x<col-1; ++x) {
 			m[(row-1)*col+x].print(c);
 			c.s << ",";
 		}
 		m[row*col-1].print(c);
-		c.s << "]] ]]";
+		c.s << "]]";
 
 	}
 }
@@ -599,6 +599,19 @@ matrix matrix::mul(const numeric & other) const
 }
 
 
+/** Product of matrix and scalar expression. */
+matrix matrix::mul_scalar(const ex & other) const
+{
+	exvector prod(row * col);
+
+	for (unsigned r=0; r<row; ++r)
+		for (unsigned c=0; c<col; ++c)
+			prod[r*col+c] = m[r*col+c] * other;
+
+	return matrix(row, col, prod);
+}
+
+
 /** operator() to access elements.
  *
  *  @param ro row of element
diff --git a/ginac/matrix.h b/ginac/matrix.h
index 0e7660b9..08f99ede 100644
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -49,6 +49,7 @@ public:
 	ex expand(unsigned options=0) const;
 	ex eval(int level=0) const;
 	ex evalf(int level=0) const;
+	ex evalm(void) const {return *this;}
 	ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const;
 	ex eval_indexed(const basic & i) const;
 	ex add_indexed(const ex & self, const ex & other) const;
@@ -69,6 +70,7 @@ public:
 	matrix sub(const matrix & other) const;
 	matrix mul(const matrix & other) const;
 	matrix mul(const numeric & other) const;
+	matrix mul_scalar(const ex & other) const;
 	const ex & operator() (unsigned ro, unsigned co) const;
 	matrix & set(unsigned ro, unsigned co, ex value);
 	matrix transpose(void) const;
diff --git a/ginac/mul.cpp b/ginac/mul.cpp
index 8e641caf..a4bb802e 100644
--- a/ginac/mul.cpp
+++ b/ginac/mul.cpp
@@ -26,6 +26,7 @@
 #include "mul.h"
 #include "add.h"
 #include "power.h"
+#include "matrix.h"
 #include "archive.h"
 #include "debugmsg.h"
 #include "utils.h"
@@ -404,6 +405,45 @@ ex mul::evalf(int level) const
 	return mul(s,overall_coeff.evalf(level));
 }
 
+ex mul::evalm(void) const
+{
+	// numeric*matrix
+	if (seq.size() == 1 && is_ex_of_type(seq[0].rest, matrix))
+		return ex_to_matrix(seq[0].rest).mul(ex_to_numeric(overall_coeff));
+
+	// Evaluate children first, look whether there are any matrices at all
+	// (there can be either no matrices or one matrix; if there were more
+	// than one matrix, it would be a non-commutative product)
+	epvector *s = new epvector;
+	s->reserve(seq.size());
+
+	bool have_matrix = false;
+	epvector::iterator the_matrix;
+
+	epvector::const_iterator it = seq.begin(), itend = seq.end();
+	while (it != itend) {
+		const ex &m = recombine_pair_to_ex(*it).evalm();
+		s->push_back(split_ex_to_pair(m));
+		if (is_ex_of_type(m, matrix)) {
+			have_matrix = true;
+			the_matrix = s->end() - 1;
+		}
+		it++;
+	}
+
+	if (have_matrix) {
+
+		// The product contained a matrix. We will multiply all other factors
+		// into that matrix.
+		matrix m = ex_to_matrix(the_matrix->rest);
+		s->erase(the_matrix);
+		ex scalar = (new mul(s, overall_coeff))->setflag(status_flags::dynallocated);
+		return m.mul_scalar(scalar);
+
+	} else
+		return (new mul(s, overall_coeff))->setflag(status_flags::dynallocated);
+}
+
 ex mul::simplify_ncmul(const exvector & v) const
 {
 	if (seq.size()==0) {
diff --git a/ginac/mul.h b/ginac/mul.h
index d5a73ccf..2808db99 100644
--- a/ginac/mul.h
+++ b/ginac/mul.h
@@ -55,6 +55,7 @@ public:
 	ex coeff(const ex & s, int n = 1) const;
 	ex eval(int level=0) const;
 	ex evalf(int level=0) const;
+	ex evalm(void) const;
 	ex series(const relational & s, int order, unsigned options = 0) const;
 	ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
 	numeric integer_content(void) const;
diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp
index 8a0cc122..43ded211 100644
--- a/ginac/ncmul.cpp
+++ b/ginac/ncmul.cpp
@@ -28,6 +28,7 @@
 #include "ex.h"
 #include "add.h"
 #include "mul.h"
+#include "matrix.h"
 #include "print.h"
 #include "archive.h"
 #include "debugmsg.h"
@@ -432,6 +433,36 @@ ex ncmul::eval(int level) const
 										  status_flags::evaluated);
 }
 
+ex ncmul::evalm(void) const
+{
+	// Evaluate children first
+	exvector *s = new exvector;
+	s->reserve(seq.size());
+	exvector::const_iterator it = seq.begin(), itend = seq.end();
+	while (it != itend) {
+		s->push_back(it->evalm());
+		it++;
+	}
+
+	// If there are only matrices, simply multiply them
+	it = s->begin(); itend = s->end();
+	if (is_ex_of_type(*it, matrix)) {
+		matrix prod(ex_to_matrix(*it));
+		it++;
+		while (it != itend) {
+			if (!is_ex_of_type(*it, matrix))
+				goto no_matrix;
+			prod = prod.mul(ex_to_matrix(*it));
+			it++;
+		}
+		delete s;
+		return prod;
+	}
+
+no_matrix:
+	return (new ncmul(s))->setflag(status_flags::dynallocated);
+}
+
 ex ncmul::thisexprseq(const exvector & v) const
 {
 	return (new ncmul(v))->setflag(status_flags::dynallocated);
@@ -518,9 +549,10 @@ exvector ncmul::expandchildren(unsigned options) const
 {
 	exvector s;
 	s.reserve(seq.size());
-
-	for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
-		s.push_back((*it).expand(options));
+	exvector::const_iterator it = seq.begin(), itend = seq.end();
+	while (it != itend) {
+		s.push_back(it->expand(options));
+		it++;
 	}
 	return s;
 }
diff --git a/ginac/ncmul.h b/ginac/ncmul.h
index f3fc1019..86204484 100644
--- a/ginac/ncmul.h
+++ b/ginac/ncmul.h
@@ -61,6 +61,7 @@ public:
 	ex expand(unsigned options=0) const;
 	ex coeff(const ex & s, int n=1) const;
 	ex eval(int level=0) const;
+	ex evalm(void) const;
 	exvector get_free_indices(void) const;
 	ex thisexprseq(const exvector & v) const;
 	ex thisexprseq(exvector * vp) const;
diff --git a/ginac/power.cpp b/ginac/power.cpp
index 9d4b0c86..11fcbffe 100644
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -28,9 +28,9 @@
 #include "expairseq.h"
 #include "add.h"
 #include "mul.h"
+#include "ncmul.h"
 #include "numeric.h"
 #include "inifcns.h"
-#include "relational.h"
 #include "symbol.h"
 #include "print.h"
 #include "archive.h"
@@ -68,7 +68,6 @@ DEFAULT_DESTROY(power)
 power::power(const ex & lh, const ex & rh) : basic(TINFO_power), basis(lh), exponent(rh)
 {
 	debugmsg("power ctor from ex,ex",LOGLEVEL_CONSTRUCT);
-	GINAC_ASSERT(basis.return_type()==return_types::commutative);
 }
 
 /** Ctor from an ex and a bare numeric.  This is somewhat more efficient than
@@ -76,7 +75,6 @@ power::power(const ex & lh, const ex & rh) : basic(TINFO_power), basis(lh), expo
 power::power(const ex & lh, const numeric & rh) : basic(TINFO_power), basis(lh), exponent(rh)
 {
 	debugmsg("power ctor from ex,numeric",LOGLEVEL_CONSTRUCT);
-	GINAC_ASSERT(basis.return_type()==return_types::commutative);
 }
 
 //////////
@@ -243,6 +241,11 @@ ex & power::let_op(int i)
 	return i==0 ? basis : exponent;
 }
 
+ex power::map(map_func f) const
+{
+	return (new power(f(basis), f(exponent)))->setflag(status_flags::dynallocated);
+}
+
 int power::degree(const ex & s) const
 {
 	if (is_exactly_of_type(*exponent.bp,numeric)) {
@@ -320,17 +323,17 @@ ex power::eval(int level) const
 	const ex & ebasis    = level==1 ? basis    : basis.eval(level-1);
 	const ex & eexponent = level==1 ? exponent : exponent.eval(level-1);
 	
-	bool basis_is_numerical = 0;
-	bool exponent_is_numerical = 0;
+	bool basis_is_numerical = false;
+	bool exponent_is_numerical = false;
 	numeric * num_basis;
 	numeric * num_exponent;
 	
 	if (is_exactly_of_type(*ebasis.bp,numeric)) {
-		basis_is_numerical = 1;
+		basis_is_numerical = true;
 		num_basis = static_cast<numeric *>(ebasis.bp);
 	}
 	if (is_exactly_of_type(*eexponent.bp,numeric)) {
-		exponent_is_numerical = 1;
+		exponent_is_numerical = true;
 		num_exponent = static_cast<numeric *>(eexponent.bp);
 	}
 	
@@ -360,89 +363,97 @@ ex power::eval(int level) const
 	if (ebasis.is_equal(_ex1()))
 		return _ex1();
 	
-	if (basis_is_numerical && exponent_is_numerical) {
+	if (exponent_is_numerical) {
+
 		// ^(c1,c2) -> c1^c2 (c1, c2 numeric(),
 		// except if c1,c2 are rational, but c1^c2 is not)
-		bool basis_is_crational = num_basis->is_crational();
-		bool exponent_is_crational = num_exponent->is_crational();
-		numeric res = num_basis->power(*num_exponent);
+		if (basis_is_numerical) {
+			bool basis_is_crational = num_basis->is_crational();
+			bool exponent_is_crational = num_exponent->is_crational();
+			numeric res = num_basis->power(*num_exponent);
 		
-		if ((!basis_is_crational || !exponent_is_crational)
-			|| res.is_crational()) {
-			return res;
-		}
-		GINAC_ASSERT(!num_exponent->is_integer());  // has been handled by now
-		// ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-h)<1, q integer
-		if (basis_is_crational && exponent_is_crational
-			&& num_exponent->is_real()
-			&& !num_exponent->is_integer()) {
-			numeric n = num_exponent->numer();
-			numeric m = num_exponent->denom();
-			numeric r;
-			numeric q = iquo(n, m, r);
-			if (r.is_negative()) {
-				r = r.add(m);
-				q = q.sub(_num1());
+			if ((!basis_is_crational || !exponent_is_crational)
+				|| res.is_crational()) {
+				return res;
 			}
-			if (q.is_zero())  // the exponent was in the allowed range 0<(n/m)<1
-				return this->hold();
-			else {
-				epvector res;
-				res.push_back(expair(ebasis,r.div(m)));
-				return (new mul(res,ex(num_basis->power_dyn(q))))->setflag(status_flags::dynallocated | status_flags::evaluated);
+			GINAC_ASSERT(!num_exponent->is_integer());  // has been handled by now
+
+			// ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-h)<1, q integer
+			if (basis_is_crational && exponent_is_crational
+				&& num_exponent->is_real()
+				&& !num_exponent->is_integer()) {
+				numeric n = num_exponent->numer();
+				numeric m = num_exponent->denom();
+				numeric r;
+				numeric q = iquo(n, m, r);
+				if (r.is_negative()) {
+					r = r.add(m);
+					q = q.sub(_num1());
+				}
+				if (q.is_zero())  // the exponent was in the allowed range 0<(n/m)<1
+					return this->hold();
+				else {
+					epvector res;
+					res.push_back(expair(ebasis,r.div(m)));
+					return (new mul(res,ex(num_basis->power_dyn(q))))->setflag(status_flags::dynallocated | status_flags::evaluated);
+				}
 			}
 		}
-	}
 	
-	// ^(^(x,c1),c2) -> ^(x,c1*c2)
-	// (c1, c2 numeric(), c2 integer or -1 < c1 <= 1,
-	// case c1==1 should not happen, see below!)
-	if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,power)) {
-		const power & sub_power = ex_to_power(ebasis);
-		const ex & sub_basis = sub_power.basis;
-		const ex & sub_exponent = sub_power.exponent;
-		if (is_ex_exactly_of_type(sub_exponent,numeric)) {
-			const numeric & num_sub_exponent = ex_to_numeric(sub_exponent);
-			GINAC_ASSERT(num_sub_exponent!=numeric(1));
-			if (num_exponent->is_integer() || abs(num_sub_exponent)<1)
-				return power(sub_basis,num_sub_exponent.mul(*num_exponent));
+		// ^(^(x,c1),c2) -> ^(x,c1*c2)
+		// (c1, c2 numeric(), c2 integer or -1 < c1 <= 1,
+		// case c1==1 should not happen, see below!)
+		if (is_ex_exactly_of_type(ebasis,power)) {
+			const power & sub_power = ex_to_power(ebasis);
+			const ex & sub_basis = sub_power.basis;
+			const ex & sub_exponent = sub_power.exponent;
+			if (is_ex_exactly_of_type(sub_exponent,numeric)) {
+				const numeric & num_sub_exponent = ex_to_numeric(sub_exponent);
+				GINAC_ASSERT(num_sub_exponent!=numeric(1));
+				if (num_exponent->is_integer() || (abs(num_sub_exponent) - _num1()).is_negative())
+					return power(sub_basis,num_sub_exponent.mul(*num_exponent));
+			}
 		}
-	}
 	
-	// ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
-	if (exponent_is_numerical && num_exponent->is_integer() &&
-		is_ex_exactly_of_type(ebasis,mul)) {
-		return expand_mul(ex_to_mul(ebasis), *num_exponent);
-	}
+		// ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
+		if (num_exponent->is_integer() && is_ex_exactly_of_type(ebasis,mul)) {
+			return expand_mul(ex_to_mul(ebasis), *num_exponent);
+		}
 	
-	// ^(*(...,x;c1),c2) -> ^(*(...,x;1),c2)*c1^c2 (c1, c2 numeric(), c1>0)
-	// ^(*(...,x,c1),c2) -> ^(*(...,x;-1),c2)*(-c1)^c2 (c1, c2 numeric(), c1<0)
-	if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,mul)) {
-		GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above
-		const mul & mulref = ex_to_mul(ebasis);
-		if (!mulref.overall_coeff.is_equal(_ex1())) {
-			const numeric & num_coeff = ex_to_numeric(mulref.overall_coeff);
-			if (num_coeff.is_real()) {
-				if (num_coeff.is_positive()) {
-					mul * mulp = new mul(mulref);
-					mulp->overall_coeff = _ex1();
-					mulp->clearflag(status_flags::evaluated);
-					mulp->clearflag(status_flags::hash_calculated);
-					return (new mul(power(*mulp,exponent),
-					                power(num_coeff,*num_exponent)))->setflag(status_flags::dynallocated);
-				} else {
-					GINAC_ASSERT(num_coeff.compare(_num0())<0);
-					if (num_coeff.compare(_num_1())!=0) {
+		// ^(*(...,x;c1),c2) -> ^(*(...,x;1),c2)*c1^c2 (c1, c2 numeric(), c1>0)
+		// ^(*(...,x,c1),c2) -> ^(*(...,x;-1),c2)*(-c1)^c2 (c1, c2 numeric(), c1<0)
+		if (is_ex_exactly_of_type(ebasis,mul)) {
+			GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above
+			const mul & mulref = ex_to_mul(ebasis);
+			if (!mulref.overall_coeff.is_equal(_ex1())) {
+				const numeric & num_coeff = ex_to_numeric(mulref.overall_coeff);
+				if (num_coeff.is_real()) {
+					if (num_coeff.is_positive()) {
 						mul * mulp = new mul(mulref);
-						mulp->overall_coeff = _ex_1();
+						mulp->overall_coeff = _ex1();
 						mulp->clearflag(status_flags::evaluated);
 						mulp->clearflag(status_flags::hash_calculated);
 						return (new mul(power(*mulp,exponent),
-						                power(abs(num_coeff),*num_exponent)))->setflag(status_flags::dynallocated);
+						                power(num_coeff,*num_exponent)))->setflag(status_flags::dynallocated);
+					} else {
+						GINAC_ASSERT(num_coeff.compare(_num0())<0);
+						if (num_coeff.compare(_num_1())!=0) {
+							mul * mulp = new mul(mulref);
+							mulp->overall_coeff = _ex_1();
+							mulp->clearflag(status_flags::evaluated);
+							mulp->clearflag(status_flags::hash_calculated);
+							return (new mul(power(*mulp,exponent),
+							                power(abs(num_coeff),*num_exponent)))->setflag(status_flags::dynallocated);
+						}
 					}
 				}
 			}
 		}
+
+		// ^(nc,c1) -> ncmul(nc,nc,...) (c1 positive integer)
+		if (num_exponent->is_pos_integer() && ebasis.return_type() != return_types::commutative) {
+			return ncmul(exvector(num_exponent->to_int(), ebasis), true);
+		}
 	}
 	
 	if (are_ex_trivially_equal(ebasis,basis) &&
diff --git a/ginac/power.h b/ginac/power.h
index 8f7ac501..00f02b1f 100644
--- a/ginac/power.h
+++ b/ginac/power.h
@@ -53,6 +53,7 @@ public:
 	bool info(unsigned inf) const;
 	unsigned nops() const;
 	ex & let_op(int i);
+	ex map(map_func f) const;
 	int degree(const ex & s) const;
 	int ldegree(const ex & s) const;
 	ex coeff(const ex & s, int n = 1) const;
diff --git a/ginac/relational.cpp b/ginac/relational.cpp
index 7f6337d2..c4b3ec07 100644
--- a/ginac/relational.cpp
+++ b/ginac/relational.cpp
@@ -171,6 +171,11 @@ ex & relational::let_op(int i)
 
 	return i==0 ? lh : rh;
 }
+
+ex relational::map(map_func f) const
+{
+	return (new relational(f(lh), f(rh), o))->setflag(status_flags::dynallocated);
+}
 	
 ex relational::eval(int level) const
 {
diff --git a/ginac/relational.h b/ginac/relational.h
index c3b223b5..8f4efed9 100644
--- a/ginac/relational.h
+++ b/ginac/relational.h
@@ -56,6 +56,7 @@ public:
 	bool info(unsigned inf) const;
 	unsigned nops() const;
 	ex & let_op(int i);
+	ex map(map_func f) const;
 	ex eval(int level=0) const;
 	ex evalf(int level=0) const;
 	ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;