[GiNaC-list] Fwd: Performance problem with series

[Frieder Kalisch]

Hello,

Am Freitag 21 Januar 2005 13:10 schrieb Sheplyakov Alexei:
> On Wed, Jan 19, 2005 at 06:31:41PM +0100, Chris Dams wrote:
> > Ad (2): I think that the expansion coefficients should be in expandend
> > form (again in the sense of .expand()).
>
> In ginsh:
>
> series((a+x)^(100), x==1, 3);
> ((1+a)^100)+(100*(1+a)^99)*(-1+x)+(4950*(1+a)^98)*(-1+x)^2+Order((-1+x)^3)
>
> I don't think anyone wants these coefficient to be .expand()'-ed...

Perhaps the coefficients should preserve the "expandedness" of the original expression. The 
patch below implements this. In each step of the composition of a series in pseries::mul_series
and pseries::power_const the coefficients are multiplied term by term. In this way deeply 
nested expressions are avoided.

Another way to implement the same feature would be to add an option to expand to not descend 
into children . Any comments about these two suggestions and the patch? Please note that I am 
not aware of any coding standards within ginac.

Regards

                 Frieder Kalisch


*** ginac/pseries.cpp.save      2005-01-18 10:13:46.000000000 +0100
--- ginac/pseries.cpp   2005-02-03 13:41:45.000000000 +0100
***************
*** 693,698 ****
--- 693,738 ----
  }


+ /** Helper struct for multiply_term_by_term.
+  * map_function, that multiplies its argument by a constant. If the
+  * argument or the constant is a sum, they are multiplied term by term. */
+ struct multiply_by : public map_function {
+     ex c;
+
+     multiply_by(const ex &e) : c(e) {} ;
+
+     ex operator()(const ex &e)
+     {
+         if (is_a<add>(c)) {
+           multiply_by mult_func(e);
+           return c.map(mult_func);
+       }
+
+       return e*c;
+     }
+
+ };
+
+
+ /** Helper function, that uses struct multiply_by.
+  * Multiply two expressions. If either is a sum, they are multiplied
+  * term by term. */
+ ex multiply_term_by_term(const ex &e1, const ex &e2)
+ {
+     if (is_a<add>(e1)) {
+         multiply_by mult_func(e2);
+         return e1.map(mult_func);
+     }
+
+     if (is_a<add>(e2)) {
+         multiply_by mult_func(e1);
+       return e2.map(mult_func);
+     }
+
+     return e1*e2;
+ }
+
+
  /** Multiply a pseries object with a numeric constant, producing a pseries
   *  object that represents the product.
   *
***************
*** 756,762 ****
                        ex a_coeff = coeff(var, i);
                        ex b_coeff = other.coeff(var, cdeg-i);
                        if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
!                               co += a_coeff * b_coeff;
                }
                if (!co.is_zero())
                        new_seq.push_back(expair(co, numeric(cdeg)));
--- 796,802 ----
                        ex a_coeff = coeff(var, i);
                        ex b_coeff = other.coeff(var, cdeg-i);
                        if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
!                               co += multiply_term_by_term(a_coeff, b_coeff);
                }
                if (!co.is_zero())
                        new_seq.push_back(expair(co, numeric(cdeg)));
***************
*** 893,906 ****
        for (int i=1; i<deg; ++i) {
                ex sum = _ex0;
                for (int j=1; j<=i; ++j) {
!                       ex c = coeff(var, j + ldeg);
                        if (is_order_function(c)) {
                                co.push_back(Order(_ex1));
                                break;
                        } else
!                               sum += (p * j - (i - j)) * co[i - j] * c;
                }
!               co.push_back(sum / coeff(var, ldeg) / i);
        }

        // Construct new series (of non-zero coefficients)
--- 933,946 ----
        for (int i=1; i<deg; ++i) {
                ex sum = _ex0;
                for (int j=1; j<=i; ++j) {
!                       ex c = multiply_term_by_term(coeff(var, j + ldeg), _ex1/i/coeff(var, ldeg));
                        if (is_order_function(c)) {
                                co.push_back(Order(_ex1));
                                break;
                        } else
!                               sum += multiply_term_by_term((p * j - (i - j)) * co[i - j], c);
                }
!               co.push_back(sum);
        }

        // Construct new series (of non-zero coefficients)



-- 
Frieder Kalisch                         Institut für theoretische Physik
kalisch at tphys.uni-heidelberg.de         Philosophenweg 19
+49-6221-549-433                        D-69120 Heidelberg