# How to implement fast series multiplication

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!

• Is there a reason you're not simply writing a wrapper for the built-in SeriesData object? That way, you can simply use the multiplication implemented for SeriesData when implementing the multiplication of your wrapper, with negligible overhead – Lukas Lang Oct 5 '20 at 10:54
• If the series coefficient lists are long, you might get better performance using ListConvolve. – Daniel Lichtblau Oct 5 '20 at 12:46
• @LukasLang I need to make slight modifications to it, so I cannot simply use SeriesData. – QuantumDot Oct 5 '20 at 22:02
• @DanielLichtblau Interesting, though I'm not seeing how to use ListConvolve. Can you provide a short example of how I can use it? – QuantumDot Oct 5 '20 at 22:02

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.

• So it seems ListConvolve[c1, c2, {1, -1}, 0] is exactly what I am looking for! Thank you. I still have to better acquaint myself with ListConvolute and friends. I can raise the second element of the list to truncate the desired number of terms. Maybe something like ListConvolve[c1, c2, {1, -Length[c1]}, 0] Do I understand that correctly? Thanks again! – QuantumDot Oct 7 '20 at 5:35
• I think so. I've not checked parameter settings specific to series-style multiplication but what you write looks like it might be correct. – Daniel Lichtblau Oct 7 '20 at 13:49

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 *)