1
$\begingroup$

I have an expression which consists of terms with undefined function calls a[n]:

example = 1 - c^2 + c a[1] a[2] + 1/2 c^2 a[1]^2 a[2]^2 + c a[1] a[3]

Now I want to transform each term with individual as to a different function v[m1,m2,m3], such that v[[m_i]] is the exponent of a[[i]]. In addition, I want that the coefficient of v[[m_i]] is multiplied by the (exponent of a[[i]] + 1):

fct[example]
(* example -> (1 - c^2)*(0+1)*(0+1)*(0+1)*v[0,0,0] + c*(1+1)*(1+1)*(0+1)*v[1,1,0] +
              + 1/2 c^2*(2+1)*(2+1)*(0+1)*v[2,2,0] + c*(1+1)*(0+1)*(1+1)*v[1,0,1]

            = (1 - c^2)*v[0,0,0] + c*4*v[1,1,0] +
              + 1/2 c^2*9*v[2,2,0] + c*4*v[1,0,1]
*)

I found one solution, but it is annoyingly slow, and I was hoping that somebody finds faster methods. The idea is to multiply example with a[1]^2 * a[2]^2 * a[3]^2 (such that every term has the form of coeff * a[1]^n1*a[2]^n2*a[3]^n3), then use a Replace-Rule:

fct[expr_] := 
 Expand[expr*a[1]^2*a[2]^2*a[3]^2] /. 
    {a[1]^n1_*a[2]^n2_*a[3]^n3_ -> (n1-1)*(n2-1)*(n3-1)*v[n1-2, n2-2, n3-2]}

fct[example] (* v[0, 0, 0] - c^2 v[0, 0, 0] + 4 c v[1, 0, 1] + 4 c v[1, 1, 0] + 
                9/2 c^2 v[2, 2, 0] *)

My specific questions:

  1. How can I make the algorithm general for arbitrary a[n] with n>3? (For the way I did it here, I dont know how to access the exponents more general)

  2. How can I make the algorithm faster? (It takes already ~5sec for n=8, and roughly 80 terms of a[1]^n1*a[2]^n2*a[3]^n3 in example.)

Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

You may use CoefficientRules and may find the Polynomial Algebra guide useful.

ClearAll[fct];

fct[expr_, vars_] :=
 Total@*KeyValueMap[#2 Times@@(#1 + 1) v@@#1 &]@*Association@@CoefficientRules[expr, vars]

Then following return immediately.

fct[example, Array[a, 3]]

(1 - c^2) v[0, 0, 0] + 4 c v[1, 0, 1] + 4 c v[1, 1, 0] + 9/2 c^2 v[2, 2, 0]

and

fct[example* a[1]^2*a[2]^2*a[3]^2, Array[a, 3]]

27 (1 - c^2) v[2, 2, 2] + 48 c v[3, 2, 3] + 48 c v[3, 3, 2] + 75/2 c^2 v[4, 4, 2]

Hope this helps.

Improve this answer
$\endgroup$
2
  • $\begingroup$ That looks really nice, CoefficientRules seems to be perfect for my question. Unfortunatly, my M10 does not recognize KeyValueMap and Association, so am installing M11 now. Thank you $\endgroup$ Commented Dec 31, 2016 at 0:01
  • $\begingroup$ Wow that is incredible. Your method is more than 1000 times faster in some relevant tests. And I was also learning some nice new insturction :-) Thank you, that helped very much! $\endgroup$ Commented Dec 31, 2016 at 0:43
1
$\begingroup$ 
v /: v[x__] v[y__] := Apply[v, {x} + {y}]

transform[expr_] :=
 Module[{max, temp, free},
  max = Max@Cases[expr, a[i_] :> i, Infinity];
  temp = expr /. a[i_]^p_. :> (p + 1)*(v @@ UnitVector[max, i]);
  free = Cases[temp, x_ /; FreeQ[x, v]];
  With[{t = Total[free]},
   temp - t + (v @@ ConstantArray[0, max])*t]]

transform[example]

(1 - c^2) v[0, 0, 0] + 4 c v[1, 0, 1] + 4 c v[1, 1, 0] + 9/2 c^2 v[1, 1, 0]

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.