the destructor, copy constructor, and assignment operator (which were the

[ginac.git] / ginac / pseries.cpp

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index bed134ed10abfab7506c7f8cbecd8ede4662191b..d447aaa64de9d98e135beb67ec7e32a03509cc62 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -4,7 +4,7 @@
  *  methods for series expansion. */
 
 /*
- *  GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2001 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
@@ -50,28 +50,6 @@ pseries::pseries() : basic(TINFO_pseries)
        debugmsg("pseries default constructor", LOGLEVEL_CONSTRUCT);
 }
 
-pseries::~pseries()
-{
-       debugmsg("pseries destructor", LOGLEVEL_DESTRUCT);
-       destroy(false);
-}
-
-pseries::pseries(const pseries &other)
-{
-       debugmsg("pseries copy constructor", LOGLEVEL_CONSTRUCT);
-       copy(other);
-}
-
-const pseries &pseries::operator=(const pseries & other)
-{
-       debugmsg("pseries operator=", LOGLEVEL_ASSIGNMENT);
-       if (this != &other) {
-               destroy(true);
-               copy(other);
-       }
-       return *this;
-}
-
 void pseries::copy(const pseries &other)
 {
        inherited::copy(other);
@@ -755,17 +733,34 @@ ex mul::series(const relational & r, int order, unsigned options) const
  *  @param deg  truncation order of series calculation */
 ex pseries::power_const(const numeric &p, int deg) const
 {
-       int i;
+       // method:
+       // let A(x) be this series and for the time being let it start with a
+       // constant (later we'll generalize):
+       //     A(x) = a_0 + a_1*x + a_2*x^2 + ...
+       // We want to compute
+       //     C(x) = A(x)^p
+       //     C(x) = c_0 + c_1*x + c_2*x^2 + ...
+       // Taking the derivative on both sides and multiplying with A(x) one
+       // immediately arrives at
+       //     C'(x)*A(x) = p*C(x)*A'(x)
+       // Multiplying this out and comparing coefficients we get the recurrence
+       // formula
+       //     c_i = (i*p*a_i*c_0 + ((i-1)*p-1)*a_{i-1}*c_1 + ...
+       //                    ... + (p-(i-1))*a_1*c_{i-1})/(a_0*i)
+       // which can easily be solved given the starting value c_0 = (a_0)^p.
+       // For the more general case where the leading coefficient of A(x) is not
+       // a constant, just consider A2(x) = A(x)*x^m, with some integer m and
+       // repeat the above derivation.  The leading power of C2(x) = A2(x)^2 is
+       // then of course x^(p*m) but the recurrence formula still holds.
        const symbol *s = static_cast<symbol *>(var.bp);
        int ldeg = ldegree(*s);
        
-       // Calculate coefficients of powered series
+       // Compute coefficients of the powered series
        exvector co;
        co.reserve(deg);
-       ex co0;
-       co.push_back(co0 = power(coeff(*s, ldeg), p));
+       co.push_back(power(coeff(*s, ldeg), p));
        bool all_sums_zero = true;
-       for (i=1; i<deg; ++i) {
+       for (int i=1; i<deg; ++i) {
                ex sum = _ex0();
                for (int j=1; j<=i; ++j) {
                        ex c = coeff(*s, j + ldeg);
@@ -777,13 +772,13 @@ ex pseries::power_const(const numeric &p, int deg) const
                }
                if (!sum.is_zero())
                        all_sums_zero = false;
-               co.push_back(co0 * sum / numeric(i));
+               co.push_back(sum / coeff(*s, ldeg) / numeric(i));
        }
        
        // Construct new series (of non-zero coefficients)
        epvector new_seq;
        bool higher_order = false;
-       for (i=0; i<deg; ++i) {
+       for (int i=0; i<deg; ++i) {
                if (!co[i].is_zero())
                        new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
                if (is_order_function(co[i])) {
@@ -814,12 +809,20 @@ ex power::series(const relational & r, int order, unsigned options) const
 {
        ex e;
        if (!is_ex_exactly_of_type(basis, pseries)) {
-               // Basis is not a series, may there be a singulary?
-               if (!exponent.info(info_flags::negint))
+               // Basis is not a series, may there be a singularity?
+               bool must_expand_basis = false;
+               try {
+                       basis.subs(r);
+               } catch (pole_error) {
+                       must_expand_basis = true;
+               }
+               
+               // Is the expression of type something^(-int)?
+               if (!must_expand_basis && !exponent.info(info_flags::negint))
                        return basic::series(r, order, options);
                
-               // Expression is of type something^(-int), check for singularity
-               if (!basis.subs(r).is_zero())
+               // Is the expression of type 0^something?
+               if (!must_expand_basis && !basis.subs(r).is_zero())
                        return basic::series(r, order, options);
                
                // Singularity encountered, expand basis into series
@@ -901,13 +904,6 @@ ex ex::series(const ex & r, int order, unsigned options) const
 
 unsigned pseries::precedence = 38;  // for clarity just below add::precedence
 
-//////////
-// global constants
-//////////
-
-const pseries some_pseries;
-const std::type_info & typeid_pseries = typeid(some_pseries);
-
 #ifndef NO_NAMESPACE_GINAC
 } // namespace GiNaC
 #endif // ndef NO_NAMESPACE_GINAC