How to implement fast series multiplication

Question

I need to generalize the SeriesData object for my own purposes. One of the things I need to do is to reimplement code to multiply product of series.

I've made two attempts at this, and both of them are slower than the built-in SeriesData. Is there an algorithm that has better Timing than mine?

(*Two of my implementations*)
multSerList1[lists__] := 
    Array[Plus @@ Times @@@ (MapThread[Part, {{lists}, #}, 1] & /@ 
      Flatten[Permutations /@ IntegerPartitions[#, {3}], 1]) &, Min[Length /@ {lists}], Length[{lists}]];

multSerList2[listFirst_, listRest__] := 
    Fold[Function[{a1, a2}, Array[Inner[Times, Take[a1, #], Reverse[Take[a2, #]], Plus] &, Min[Length /@ 
      {listFirst, listRest}]]], listFirst, {listRest}];

To test this, I try to multiply the following three series together:

realExampleList = {
  List @@ Normal[Series[Exp[y x], {x, 0, 4}]], 
  List @@ Normal[Series[Log[1 + c x], {x, 0, 5}]],
  List @@ Normal[Series[PolyLog[2, -n x], {x, 0, 3}]]}

Then multSerList1@@realExampleList or multSerList2@@ realExampleList both yield

I can apply AbsoluteTiming to time my code. It takes 0.0002 s for first one and 0.0001 s for the second one on my machine. But Multiplying the SeriesData (need to remove the Normal and List)

realExampleListSerData = {
  Series[Exp[y x], {x, 0, 4}], 
  Series[Log[1 + c x], {x, 0, 5}],
  Series[PolyLog[2, -n x], {x, 0, 3}]};

By simply doing Times@@realExampleListSerData, it gets the answer in 0.00002 s, which is five times faster.

I need help implementing the multiplication of series that performs approximately as well as the SeriesData.

Thanks!

Answer 1

Here is the idea for using convolution. I suspect there is similar in the reference page for ListConvolve. I first illustrate in a slightly roundabout way, by creating explicit polynomials and then extracting coefficient lists.

SeedRandom[1234]
s = 10;
n = 8;
p1 = randomPoly[n, s, x]
p2 = randomPoly[n, s, x]
c1 = CoefficientList[p1, x]
c2 = CoefficientList[p2, x]

(* Out[1843]= -9 + 10 x + 10 x^2 - 9 x^3 - 10 x^4 + 6 x^5 + 7 x^6 - 
 6 x^7 + x^8

Out[1844]= -5 + 9 x + 5 x^2 + 6 x^3 + 9 x^4 - 8 x^5 + 2 x^6 - 
 5 x^7 + x^8

Out[1845]= {-9, 10, 10, -9, -10, 6, 7, -6, 1}

Out[1846]= {-5, 9, 5, 6, 9, -8, 2, -5, 1} *)

Now compare expanding and convolution, with the latter having argument settings to achieve Cauchy product convolution.

Expand[p1*p2]
ListConvolve[c1, c2, {1, -1}, 0]

(* Out[1847]= 45 - 131 x - 5 x^2 + 131 x^3 - 2 x^4 + 57 x^5 - 93 x^6 - 
 33 x^7 - 45 x^8 + 97 x^9 + 19 x^10 - 51 x^11 + 31 x^12 - 49 x^13 + 
 39 x^14 - 11 x^15 + x^16

Out[1848]= {45, -131, -5, 131, -2, 57, -93, -33, -45, 97, 19, -51, \
31, -49, 39, -11, 1} *)

Here is a larger example. For Series I have to double the order so as to get all coefficients in the product. You would not actually do that for series products, I'm doing it to show that the results match.

SeedRandom[1234]
bign = 10^4;
bigp1 = randomPoly[bign, s, x];
bigp2 = randomPoly[bign, s, x];
bigs1 = Series[bigp1, {x, 0, 2*bign}];
bigs2 = Series[bigp2, {x, 0, 2*bign}];
bigc1 = CoefficientList[bigp1, x];
bigc2 = CoefficientList[bigp2, x];
Timing[bigsprod = bigs1*bigs2;]
Timing[bigpprod = Expand[bigp1*bigp2];]
Timing[biglc = ListConvolve[bigc1, bigc2, {1, -1}, 0];]
biglc === bigsprod[[3]] === CoefficientList[bigpprod, x]

(* Out[1910]= {42.8077, Null}

Out[1911]= {0.02545, Null}

Out[1912]= {0.000905, Null}

Out[1913]= True *)

One can improve a bit on the convolution timing by setting the second parameter so that higher order terms are omitted.

Answer 2

If you have a fixed maximum degree of polynomials, you can compute the product of two generic polynomials. This is in some sense similar to a naive implementation of symbolic automatic differentiation.

ClearAll[make$mult] ;
make$mult[] := Block[
    {x, a, b, c1, c2, p1, p2, p3, arg, list},
    p1 = a0 + Dot[c1 = ToExpression[Map[StringTemplate["a``"],Range[order]]],x^Range[order]] ;
    p2 = b0 + Dot[c2 = ToExpression[Map[StringTemplate["b``"],Range[order]]],x^Range[order]] ;
    p3 = Collect[p1*p2,x] /. x^w_ /; w > order -> 0 ;
    arg = {a0,c1,b0,c2} // Flatten ;
    list = CoefficientList[p3,x] ;
    Apply[Function,{arg,list}]
] ;

ClearAll[make$jet] ;
make$jet[expression_] := Apply[jet,CoefficientList[Normal[Series[expression,{x,0,order}]] + w*x^order,x] /. w -> 0]

ClearAll[make$poly] ;
make$poly[jet[list__]] := {list}.x^Range[0,Length[{list}]-1] ;

ClearAll[jet] ;
jet /: Times[jet[a__],jet[b__]] := Apply[jet,Apply[mult,{a,b}]] ;

(* example *)
(* set max order *)
order = 4 ;
(* generate mult (one time computation, can be compiled) *)
ClearAll[mult] ;
mult = make$mult[] ;
(* polynomials *)
p1 = Exp[y x] ;
p2 = Log[1 + c x] ;
p3 = PolyLog[2, -n x] ;
(* convert to jets *)
j1 = make$jet[p1] ;
j2 = make$jet[p2] ;
j3 = make$jet[p3] ;
(* mult and convert to polynomial *)
res = make$poly[j1*j2*j3] ; // RepeatedTiming
(* compare *)
realExampleListSerData = {
  Series[Exp[y x], {x, 0, 4}], 
  Series[Log[1 + c x], {x, 0, 5}],
  Series[PolyLog[2, -n x], {x, 0, 3}]};
ref = Normal[Times@@realExampleListSerData] ; // RepeatedTiming
res - ref // Expand
(* {0.0000534, Null} *)
(* {0.0000131, Null} *)
(* 0 *)