power::series(): handle someg (trivial) singularities of the exponent...
[ginac.git] / ginac / pseries.cpp
index bc04b1868a73663ff0dd3983df25a521133c9798..c290fe0a1f7b09b0c73fa2cebdf4e9cb1f609c38 100644 (file)
@@ -4,7 +4,7 @@
  *  methods for series expansion. */
 
 /*
- *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2010 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
  *
  *  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 <numeric>
-#include <stdexcept>
-
 #include "pseries.h"
 #include "add.h"
 #include "inifcns.h" // for Order function
 #include "archive.h"
 #include "utils.h"
 
+#include <limits>
+#include <numeric>
+#include <stdexcept>
+
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(pseries, basic,
@@ -51,7 +52,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(pseries, basic,
  *  Default constructor
  */
 
-pseries::pseries() : inherited(TINFO_pseries) { }
+pseries::pseries() { }
 
 
 /*
@@ -67,7 +68,7 @@ pseries::pseries() : inherited(TINFO_pseries) { }
  *  @param rel_  expansion variable and point (must hold a relational)
  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
  *  @return newly constructed pseries */
-pseries::pseries(const ex &rel_, const epvector &ops_) : basic(TINFO_pseries), seq(ops_)
+pseries::pseries(const ex &rel_, const epvector &ops_) : seq(ops_)
 {
        GINAC_ASSERT(is_a<relational>(rel_));
        GINAC_ASSERT(is_a<symbol>(rel_.lhs()));
@@ -80,16 +81,22 @@ pseries::pseries(const ex &rel_, const epvector &ops_) : basic(TINFO_pseries), s
  *  Archiving
  */
 
-pseries::pseries(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+void pseries::read_archive(const archive_node &n, lst &sym_lst) 
 {
-       for (unsigned int i=0; true; ++i) {
+       inherited::read_archive(n, sym_lst);
+       archive_node::archive_node_cit first = n.find_first("coeff");
+       archive_node::archive_node_cit last = n.find_last("power");
+       ++last;
+       seq.reserve((last-first)/2);
+
+       for (archive_node::archive_node_cit loc = first; loc < last;) {
                ex rest;
                ex coeff;
-               if (n.find_ex("coeff", rest, sym_lst, i) && n.find_ex("power", coeff, sym_lst, i))
-                       seq.push_back(expair(rest, coeff));
-               else
-                       break;
+               n.find_ex_by_loc(loc++, rest, sym_lst);
+               n.find_ex_by_loc(loc++, coeff, sym_lst);
+               seq.push_back(expair(rest, coeff));
        }
+
        n.find_ex("var", var, sym_lst);
        n.find_ex("point", point, sym_lst);
 }
@@ -107,7 +114,6 @@ void pseries::archive(archive_node &n) const
        n.add_ex("point", point);
 }
 
-DEFAULT_UNARCHIVE(pseries)
 
 //////////
 // functions overriding virtual functions from base classes
@@ -287,7 +293,7 @@ int pseries::degree(const ex &s) const
                epvector::const_iterator it = seq.begin(), itend = seq.end();
                if (it == itend)
                        return 0;
-               int max_pow = INT_MIN;
+               int max_pow = std::numeric_limits<int>::min();
                while (it != itend) {
                        int pow = it->rest.degree(s);
                        if (pow > max_pow)
@@ -315,7 +321,7 @@ int pseries::ldegree(const ex &s) const
                epvector::const_iterator it = seq.begin(), itend = seq.end();
                if (it == itend)
                        return 0;
-               int min_pow = INT_MAX;
+               int min_pow = std::numeric_limits<int>::max();
                while (it != itend) {
                        int pow = it->rest.ldegree(s);
                        if (pow < min_pow)
@@ -412,21 +418,53 @@ ex pseries::evalf(int level) const
 
 ex pseries::conjugate() const
 {
+       if(!var.info(info_flags::real))
+               return conjugate_function(*this).hold();
+
        epvector * newseq = conjugateepvector(seq);
-       ex newvar = var.conjugate();
        ex newpoint = point.conjugate();
 
-       if (!newseq     && are_ex_trivially_equal(newvar, var) && are_ex_trivially_equal(point, newpoint)) {
+       if (!newseq     && are_ex_trivially_equal(point, newpoint)) {
                return *this;
        }
 
-       ex result = (new pseries(newvar==newpoint, newseq ? *newseq : seq))->setflag(status_flags::dynallocated);
+       ex result = (new pseries(var==newpoint, newseq ? *newseq : seq))->setflag(status_flags::dynallocated);
        if (newseq) {
                delete newseq;
        }
        return result;
 }
 
+ex pseries::real_part() const
+{
+       if(!var.info(info_flags::real))
+               return real_part_function(*this).hold();
+       ex newpoint = point.real_part();
+       if(newpoint != point)
+               return real_part_function(*this).hold();
+
+       epvector v;
+       v.reserve(seq.size());
+       for(epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
+               v.push_back(expair((i->rest).real_part(), i->coeff));
+       return (new pseries(var==point, v))->setflag(status_flags::dynallocated);
+}
+
+ex pseries::imag_part() const
+{
+       if(!var.info(info_flags::real))
+               return imag_part_function(*this).hold();
+       ex newpoint = point.real_part();
+       if(newpoint != point)
+               return imag_part_function(*this).hold();
+
+       epvector v;
+       v.reserve(seq.size());
+       for(epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
+               v.push_back(expair((i->rest).imag_part(), i->coeff));
+       return (new pseries(var==point, v))->setflag(status_flags::dynallocated);
+}
+
 ex pseries::eval_integ() const
 {
        epvector *newseq = NULL;
@@ -452,6 +490,34 @@ ex pseries::eval_integ() const
        return *this;
 }
 
+ex pseries::evalm() const
+{
+       // evalm each coefficient
+       epvector newseq;
+       bool something_changed = false;
+       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+               if (something_changed) {
+                       ex newcoeff = i->rest.evalm();
+                       if (!newcoeff.is_zero())
+                               newseq.push_back(expair(newcoeff, i->coeff));
+               }
+               else {
+                       ex newcoeff = i->rest.evalm();
+                       if (!are_ex_trivially_equal(newcoeff, i->rest)) {
+                               something_changed = true;
+                               newseq.reserve(seq.size());
+                               std::copy(seq.begin(), i, std::back_inserter<epvector>(newseq));
+                               if (!newcoeff.is_zero())
+                                       newseq.push_back(expair(newcoeff, i->coeff));
+                       }
+               }
+       }
+       if (something_changed)
+               return (new pseries(var==point, newseq))->setflag(status_flags::dynallocated);
+       else
+               return *this;
+}
+
 ex pseries::subs(const exmap & m, unsigned options) const
 {
        // If expansion variable is being substituted, convert the series to a
@@ -653,7 +719,7 @@ ex pseries::add_series(const pseries &other) const
        epvector::const_iterator b = other.seq.begin();
        epvector::const_iterator a_end = seq.end();
        epvector::const_iterator b_end = other.seq.end();
-       int pow_a = INT_MAX, pow_b = INT_MAX;
+       int pow_a = std::numeric_limits<int>::max(), pow_b = std::numeric_limits<int>::max();
        for (;;) {
                // If a is empty, fill up with elements from b and stop
                if (a == a_end) {
@@ -786,8 +852,8 @@ ex pseries::mul_series(const pseries &other) const
        int cdeg_min = a_min + b_min;
        int cdeg_max = a_max + b_max;
        
-       int higher_order_a = INT_MAX;
-       int higher_order_b = INT_MAX;
+       int higher_order_a = std::numeric_limits<int>::max();
+       int higher_order_b = std::numeric_limits<int>::max();
        if (is_order_function(coeff(var, a_max)))
                higher_order_a = a_max + b_min;
        if (is_order_function(other.coeff(var, b_max)))
@@ -808,7 +874,7 @@ ex pseries::mul_series(const pseries &other) const
                if (!co.is_zero())
                        new_seq.push_back(expair(co, numeric(cdeg)));
        }
-       if (higher_order_c < INT_MAX)
+       if (higher_order_c < std::numeric_limits<int>::max())
                new_seq.push_back(expair(Order(_ex1), numeric(higher_order_c)));
        return pseries(relational(var, point), new_seq);
 }
@@ -826,10 +892,13 @@ ex mul::series(const relational & r, int order, unsigned options) const
                
        // holds ldegrees of the series of individual factors
        std::vector<int> ldegrees;
+       std::vector<bool> ldegree_redo;
 
        // find minimal degrees
        const epvector::const_iterator itbeg = seq.begin();
        const epvector::const_iterator itend = seq.end();
+       // first round: obtain a bound up to which minimal degrees have to be
+       // considered
        for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
 
                ex expon = it->coeff;
@@ -843,19 +912,61 @@ ex mul::series(const relational & r, int order, unsigned options) const
                }
 
                int real_ldegree = 0;
+               bool flag_redo = false;
                try {
                        real_ldegree = buf.expand().ldegree(sym-r.rhs());
                } catch (std::runtime_error) {}
 
                if (real_ldegree == 0) {
+                       if ( factor < 0 ) {
+                               // This case must terminate, otherwise we would have division by
+                               // zero.
+                               int orderloop = 0;
+                               do {
+                                       orderloop++;
+                                       real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
+                               } while (real_ldegree == orderloop);
+                       } else {
+                               // Here it is possible that buf does not have a ldegree, therefore
+                               // check only if ldegree is negative, otherwise reconsider the case
+                               // in the second round.
+                               real_ldegree = buf.series(r, 0, options).ldegree(sym);
+                               if (real_ldegree == 0)
+                                       flag_redo = true;
+                       }
+               }
+
+               ldegrees.push_back(factor * real_ldegree);
+               ldegree_redo.push_back(flag_redo);
+       }
+
+       int degbound = order-std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
+       // Second round: determine the remaining positive ldegrees by the series
+       // method.
+       // here we can ignore ldegrees larger than degbound
+       size_t j = 0;
+       for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
+               if ( ldegree_redo[j] ) {
+                       ex expon = it->coeff;
+                       int factor = 1;
+                       ex buf;
+                       if (expon.info(info_flags::integer)) {
+                               buf = it->rest;
+                               factor = ex_to<numeric>(expon).to_int();
+                       } else {
+                               buf = recombine_pair_to_ex(*it);
+                       }
+                       int real_ldegree = 0;
                        int orderloop = 0;
                        do {
                                orderloop++;
                                real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
-                       } while (real_ldegree == orderloop);
+                       } while ((real_ldegree == orderloop)
+                                       && ( factor*real_ldegree < degbound));
+                       ldegrees[j] = factor * real_ldegree;
+                       degbound -= factor * real_ldegree;
                }
-
-               ldegrees.push_back(factor * real_ldegree);
+               j++;
        }
 
        int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
@@ -1005,6 +1116,29 @@ ex power::series(const relational & r, int order, unsigned options) const
                must_expand_basis = true;
        }
 
+       bool exponent_is_regular = true;
+       try {
+               exponent.subs(r, subs_options::no_pattern);
+       } catch (pole_error) {
+               exponent_is_regular = false;
+       }
+
+       if (!exponent_is_regular) {
+               ex l = exponent*log(basis);
+               // this == exp(l);
+               ex le = l.series(r, order, options);
+               // Note: expanding exp(l) won't help, since that will attempt
+               // Taylor expansion, and fail (because exponent is "singular")
+               // Still l itself might be expanded in Taylor series.
+               // Examples:
+               // sin(x)/x*log(cos(x))
+               // 1/x*log(1 + x)
+               return exp(le).series(r, order, options);
+               // Note: if l happens to have a Laurent expansion (with
+               // negative powers of (var - point)), expanding exp(le)
+               // will barf (which is The Right Thing).
+       }
+
        // Is the expression of type something^(-int)?
        if (!must_expand_basis && !exponent.info(info_flags::negint)
         && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
@@ -1035,7 +1169,10 @@ ex power::series(const relational & r, int order, unsigned options) const
        }
        const ex& sym = r.lhs();
        // find existing minimal degree
-       int real_ldegree = basis.expand().ldegree(sym-r.rhs());
+       ex eb = basis.expand();
+       int real_ldegree = 0;
+       if (eb.info(info_flags::rational_function))
+               real_ldegree = eb.ldegree(sym-r.rhs());
        if (real_ldegree == 0) {
                int orderloop = 0;
                do {
@@ -1100,12 +1237,12 @@ ex integral::series(const relational & r, int order, unsigned options) const
        fexpansion.reserve(fseries.nops());
        for (size_t i=0; i<fseries.nops(); ++i) {
                ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
-               fexpansion.push_back(expair(
-                       currcoeff == Order(_ex1)
-                               ? currcoeff
-                               : integral(x, a.subs(r), b.subs(r), currcoeff),
-                       ex_to<pseries>(fseries).exponop(i)
-               ));
+               currcoeff = (currcoeff == Order(_ex1))
+                       ? currcoeff
+                       : integral(x, a.subs(r), b.subs(r), currcoeff);
+               if (currcoeff != 0)
+                       fexpansion.push_back(
+                               expair(currcoeff, ex_to<pseries>(fseries).exponop(i)));
        }
 
        // Expanding lower boundary
@@ -1118,7 +1255,7 @@ ex integral::series(const relational & r, int order, unsigned options) const
                        break;
                ex currexpon = ex_to<pseries>(fseries).exponop(i);
                int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
-               currcoeff=currcoeff.series(r, orderforf);
+               currcoeff = currcoeff.series(r, orderforf);
                ex term = ex_to<pseries>(aseries).power_const(ex_to<numeric>(currexpon+1),order);
                term = ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
                term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
@@ -1134,7 +1271,7 @@ ex integral::series(const relational & r, int order, unsigned options) const
                        break;
                ex currexpon = ex_to<pseries>(fseries).exponop(i);
                int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
-               currcoeff=currcoeff.series(r, orderforf);
+               currcoeff = currcoeff.series(r, orderforf);
                ex term = ex_to<pseries>(bseries).power_const(ex_to<numeric>(currexpon+1),order);
                term = ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
                term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
@@ -1166,12 +1303,10 @@ ex ex::series(const ex & r, int order, unsigned options) const
        else
                throw (std::logic_error("ex::series(): expansion point has unknown type"));
        
-       try {
-               e = bp->series(rel_, order, options);
-       } catch (std::exception &x) {
-               throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")"));
-       }
+       e = bp->series(rel_, order, options);
        return e;
 }
 
+GINAC_BIND_UNARCHIVER(pseries);
+
 } // namespace GiNaC