lst now provides (read-only) iterator types, the use of which makes iterating

[ginac.git] / ginac / power.cpp

diff --git a/ginac/power.cpp b/ginac/power.cpp
index ed26a4f7d9291b924317cf80f6bb82c2d98d682a..2cae8a88f4a2e4c867a25c4775bbfc677e5dfbf6 100644 (file)
--- a/ginac/power.cpp
+++ b/ginac/power.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
 
 /*
- *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2003 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
@@ -23,6 +23,7 @@
 #include <vector>
 #include <iostream>
 #include <stdexcept>
+#include <limits>
 
 #include "power.h"
 #include "expairseq.h"
@@ -36,6 +37,7 @@
 #include "matrix.h"
 #include "indexed.h"
 #include "symbol.h"
+#include "lst.h"
 #include "print.h"
 #include "archive.h"
 #include "utils.h"
@@ -71,7 +73,7 @@ DEFAULT_DESTROY(power)
 // archiving
 //////////
 
-power::power(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+power::power(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
 {
        n.find_ex("basis", basis, sym_lst);
        n.find_ex("exponent", exponent, sym_lst);
@@ -226,14 +228,13 @@ bool power::info(unsigned inf) const
        return inherited::info(inf);
 }
 
-unsigned power::nops() const
+size_t power::nops() const
 {
        return 2;
 }
 
-ex & power::let_op(int i)
+ex power::op(size_t i) const
 {
-       GINAC_ASSERT(i>=0);
        GINAC_ASSERT(i<2);
 
        return i==0 ? basis : exponent;
@@ -523,16 +524,31 @@ ex power::evalm(void) const
        return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
 }
 
-ex power::subs(const lst & ls, const lst & lr, bool no_pattern) const
-{
-       const ex &subsed_basis = basis.subs(ls, lr, no_pattern);
-       const ex &subsed_exponent = exponent.subs(ls, lr, no_pattern);
+// from mul.cpp
+extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
 
-       if (are_ex_trivially_equal(basis, subsed_basis)
-        && are_ex_trivially_equal(exponent, subsed_exponent))
-               return basic::subs(ls, lr, no_pattern);
-       else
-               return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, no_pattern);
+ex power::subs(const lst & ls, const lst & lr, unsigned options) const
+{      
+       const ex &subsed_basis = basis.subs(ls, lr, options);
+       const ex &subsed_exponent = exponent.subs(ls, lr, options);
+
+       if (!are_ex_trivially_equal(basis, subsed_basis)
+        || !are_ex_trivially_equal(exponent, subsed_exponent)) 
+               return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, options);
+
+       if(!(options & subs_options::subs_algebraic))
+               return basic::subs(ls, lr, options);
+
+       lst::const_iterator its, itr;
+       for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+               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))).basic::subs(ls, lr, options);
+       }
+
+       ex result=basic::subs(ls, lr, options);
+       return result;
 }
 
 ex power::eval_ncmul(const exvector & v) const
@@ -664,7 +680,7 @@ ex power::expand_add(const add & a, int n) const
        if (n==2)
                return expand_add_2(a);
 
-       const int m = a.nops();
+       const size_t m = a.nops();
        exvector result;
        // The number of terms will be the number of combinatorial compositions,
        // i.e. the number of unordered arrangement of m nonnegative integers
@@ -676,7 +692,7 @@ ex power::expand_add(const add & a, int n) const
        intvector upper_limit(m-1);
        int l;
 
-       for (int l=0; l<m-1; ++l) {
+       for (size_t l=0; l<m-1; ++l) {
                k[l] = 0;
                k_cum[l] = 0;
                upper_limit[l] = n;
@@ -732,10 +748,10 @@ ex power::expand_add(const add & a, int n) const
                // recalc k_cum[] and upper_limit[]
                k_cum[l] = (l==0 ? k[0] : k_cum[l-1]+k[l]);
 
-               for (int i=l+1; i<m-1; ++i)
+               for (size_t i=l+1; i<m-1; ++i)
                        k_cum[i] = k_cum[i-1]+k[i];
 
-               for (int i=l+1; i<m-1; ++i)
+               for (size_t i=l+1; i<m-1; ++i)
                        upper_limit[i] = n-k_cum[i-1];
        }
 
@@ -749,7 +765,7 @@ ex power::expand_add(const add & a, int n) const
 ex power::expand_add_2(const add & a) const
 {
        epvector sum;
-       unsigned a_nops = a.nops();
+       size_t a_nops = a.nops();
        sum.reserve((a_nops*(a_nops+1))/2);
        epvector::const_iterator last = a.seq.end();