Modification of series expansion

diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index 545fcb842dd03c50ac0aa47c46977277e3cac231..6c80ed6db0ed5e32b94d04ed8ed9657bcd5e83f1 100644 (file)
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -193,7 +193,7 @@ static ex tgamma_series(const ex & arg,
        //   tgamma(x) == tgamma(x+1) / x
        // from which follows
        //   series(tgamma(x),x==-m,order) ==
-       //   series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order+1);
+       //   series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order);
        const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
        if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
                throw do_taylor();  // caught by function::series()
@@ -202,7 +202,7 @@ static ex tgamma_series(const ex & arg,
        ex ser_denom = _ex1;
        for (numeric p; p<=m; ++p)
                ser_denom *= arg+p;
-       return (tgamma(arg+m+_ex1)/ser_denom).series(rel, order+1, options);
+       return (tgamma(arg+m+_ex1)/ser_denom).series(rel, order, options);
 }
 
 

diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index 73fd4f48399f451fa3a009da5eab00c603e73208..806ed3aa5d8b51686aef3d52c4a370ba3ac406f0 100644 (file)
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -469,7 +469,7 @@ static ex tan_series(const ex &x,
        if (!(2*x_pt/Pi).info(info_flags::odd))
                throw do_taylor();  // caught by function::series()
        // if we got here we have to care for a simple pole
-       return (sin(x)/cos(x)).series(rel, order+2, options);
+       return (sin(x)/cos(x)).series(rel, order, options);
 }
 
 REGISTER_FUNCTION(tan, eval_func(tan_eval).
@@ -879,7 +879,7 @@ static ex tanh_series(const ex &x,
        if (!(2*I*x_pt/Pi).info(info_flags::odd))
                throw do_taylor();  // caught by function::series()
        // if we got here we have to care for a simple pole
-       return (sinh(x)/cosh(x)).series(rel, order+2, options);
+       return (sinh(x)/cosh(x)).series(rel, order, options);
 }
 
 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index a0250478bee30de79671149e7c434940074a0c6e..dd4f51ab91ba3becc8007706491b30e4fd6afa01 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -21,7 +21,7 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include <iostream>
+#include <numeric>
 #include <stdexcept>
 
 #include "pseries.h"
@@ -747,19 +747,59 @@ ex mul::series(const relational & r, int order, unsigned options) const
 {
        pseries acc; // Series accumulator
 
+       GINAC_ASSERT(is_a<symbol>(r.lhs()));
+       const ex& sym = r.lhs();
+               
+       // holds ldegrees of the series of individual factors
+       std::vector<int> ldegrees;
+       // flag if series have to be re-calculated
+       bool negldegree = false;
+
        // Multiply with remaining terms
        const epvector::const_iterator itbeg = seq.begin();
        const epvector::const_iterator itend = seq.end();
        for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
                ex op = recombine_pair_to_ex(*it).series(r, order, options);
 
+               int ldeg = op.ldegree(sym);
+               ldegrees.push_back(ldeg);
+               if (ldeg < 0) {
+                       negldegree = true;
+               }
+
                // Series multiplication
                if (it==itbeg)
                        acc = ex_to<pseries>(op);
                else
                        acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
        }
-       return acc.mul_const(ex_to<numeric>(overall_coeff));
+
+       if (negldegree && (seq.size() > 1)) {
+
+               // re-calculation of series
+
+               pseries newacc;
+
+               const int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
+
+               // Multiply with remaining terms
+               const epvector::const_reverse_iterator itbeg = seq.rbegin();
+               const epvector::const_reverse_iterator itend = seq.rend();
+               std::vector<int>::const_reverse_iterator itd = ldegrees.rbegin();
+               for (epvector::const_reverse_iterator it=itbeg; it!=itend; ++it, ++itd) {
+
+                       ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options);
+
+                       // Series multiplication
+                       if (it==itbeg)
+                               newacc = ex_to<pseries>(op);
+                       else
+                               newacc = ex_to<pseries>(newacc.mul_series(ex_to<pseries>(op)));
+               }
+               return newacc.mul_const(ex_to<numeric>(overall_coeff));
+       } else {
+               return acc.mul_const(ex_to<numeric>(overall_coeff));
+       }
 }
 
 
@@ -805,6 +845,9 @@ ex pseries::power_const(const numeric &p, int deg) const
        if (!(p*ldeg).is_integer())
                throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
 
+       // adjust number of coefficients
+       deg = deg - p.to_int()*ldeg;
+       
        // O(x^n)^(-m) is undefined
        if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
                throw pole_error("pseries::power_const(): division by zero",1);
@@ -842,6 +885,7 @@ ex pseries::power_const(const numeric &p, int deg) const
        }
        if (!higher_order && !all_sums_zero)
                new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
+
        return pseries(relational(var,point), new_seq);
 }
 
@@ -895,8 +939,13 @@ ex power::series(const relational & r, int order, unsigned options) const
        }
 
        // No, expand basis into series
+       int intexp = ex_to<numeric>(exponent).to_int();
        ex e = basis.series(r, order, options);
-       return ex_to<pseries>(e).power_const(ex_to<numeric>(exponent), order);
+       int ldeg = ex_to<pseries>(e).ldegree(r.lhs());
+       if (intexp * ldeg < 0) {
+               e = basis.series(r, order + ldeg*(1-intexp), options);
+       }
+       return ex_to<pseries>(e).power_const(intexp, order);
 }