Como implementar a multiplicação rápida de séries

Question

Preciso generalizar o SeriesDataobjeto para meus próprios objetivos. Uma das coisas que preciso fazer é reimplementar o código para multiplicar o produto da série.

Eu fiz duas tentativas para fazer isso, e ambas são mais lentas do que o integrado SeriesData. Existe algum algoritmo melhor Timingque o meu?

(*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}];

Para testar isso, tento multiplicar as três séries a seguir:

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}]]}

Então multSerList1@@realExampleListou multSerList2@@ realExampleListambos rendem

Posso aplicar AbsoluteTimingpara cronometrar meu código. Leva 0,0002 s para o primeiro e 0,0001 s para o segundo na minha máquina. Mas Multiplicando o SeriesData(é necessário remover o Normale 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}]};

Simplesmente fazendo Times@@realExampleListSerData, ele obtém a resposta em 0,00002 s, que é cinco vezes mais rápido.

Preciso de ajuda para implementar a multiplicação de séries que tem um desempenho tão bom quanto o SeriesData.

Obrigado!

Answer 1

Aqui está a ideia para usar a convolução. Suspeito que haja algo semelhante na página de referência para ListConvolve. Primeiro, ilustro de uma forma ligeiramente indireta, criando polinômios explícitos e, em seguida, extraindo listas de coeficientes.

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

Agora compare a expansão e a convolução, com a última tendo configurações de argumento para obter a convolução do produto de Cauchy.

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

Aqui está um exemplo maior. Pois Serieseu tenho que dobrar o pedido para obter todos os coeficientes no produto. Você não faria isso para produtos em série, estou fazendo isso para mostrar que os resultados correspondem.

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

Pode-se melhorar um pouco no tempo de convolução definindo o segundo parâmetro de forma que os termos de ordem superior sejam omitidos.

Answer 2

Se você tiver um grau máximo fixo de polinômios, poderá calcular o produto de dois polinômios genéricos. Em certo sentido, isso é semelhante a uma implementação ingênua de diferenciação automática simbólica.

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