- 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 92cb325fe04c00db04e0875c55774c9ed7856f38..17faae2c60244e76a2cdfd8f998d1294facf68fa 100644 (file)
--- 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 1e036627c68c8d39bbc37c25e944490c827e7e2c..eae04860af7cf0e409de264c86d9b41a25524a8f 100644 (file)
--- 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 90b9c1452baa6828788f51764275e161db2fbcd3..1990ebcba938bd4b5fedf81f1ac74b5724d5729f 100644 (file)
--- 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 446d6d100b0ced4cab01ede4fe0bbfdf2e17d392..1b5f2cf3d5bdc391a4d0f08ff20746d159f18828 100644 (file)
--- 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 f13c96572e2e14b4903c4a2c6cafb1588e48440b..549f066e097cc99107c14a3139d691e22248ec88 100755 (executable)
--- 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 499229e09072e4c972a61bf691457cab0bf76b0b..84b05d887dab363b35f1350b4e9ded0e193ad5b4 100644 (file)
--- 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 9421570a2bfc6c3cb35a76ed6b010f72f0c5d56d..22669e55bd01b569b9bc449197097ae11928144a 100644 (file)
--- 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 3851a5783101d03679d1619c1886ef10d6d77bf9..be5e014cb536d0f8feebdc1e756ad8a8c436e124 100644 (file)
--- 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 7081eb5b333ebe22f7fd365e80a0a513d60e0603..de5522ad119ef9e2b30a392a8122bb8e74683f1c 100644 (file)
--- 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 96c256077f92f0964c4f7fa8360f407e61e8405d..3d170ff8df5435d332971e9731e3795c3cbaefc4 100644 (file)
--- 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 6c22ec08dbca76c800fea1654292a3c23c4ae9b0..46df9595bd16d0f7efd259fd08cc9bd4bea93b92 100644 (file)
--- 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 9eab1a594cb85a0e28eebece3c7cab9339eebfe5..5e29e4820171d431c83d49560307d0ab6922d214 100644 (file)
--- 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 d47c89e2bb049c36db8cb0ccd1abd3d6b94fa85c..5767d4cfb924f91e4ec5adc7c43a59cecc409f05 100644 (file)
--- 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 0e7660b969fcc28d0d905725ac99bf55fad4162a..08f99edea92015177dd14beaed3182ffb3ef9b3e 100644 (file)
--- 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 8e641caf2ddec551ccfd6f86308649b5690b9d06..a4bb802ec9475c6b8f87ac4becbeaaf22b13cea5 100644 (file)
--- 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 d5a73ccf27ea1f188ff6e347a37fe07a768ba9f7..2808db99addaf9312ad29a9a6d5b06ce9f484193 100644 (file)
--- 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 8a0cc1220311d5eb7c0070fb76dd552e1ba28813..43ded2110572f258ce180103e80125c23275aab7 100644 (file)
--- 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 f3fc1019dcb5df152cafb8ca1034971ab3a0ac5d..86204484c099fd9e9f37fdc9b8cff3cc95928aa8 100644 (file)
--- 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 9d4b0c86c7a749bc7f38e8106209138a1a97e4c2..11fcbffeebbfc3a5ae2180ef7156c67625e365f7 100644 (file)
--- 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 8f7ac50142933bf154cf27e128c04265e416a54b..00f02b1f09f1dd7b3b706c839dc7ad1f3aa2f60d 100644 (file)
--- 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 7f6337d2907d3a282693a90b44bd4a6940528439..c4b3ec07bfc4733e813b0501c5556297a6ced740 100644 (file)
--- 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 c3b223b5eb0f494f722ffc25f20a10ce72de5431..8f4efed9ad555d9fc402a7c125cd12e685695b24 100644 (file)
--- 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;