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

$\endgroup$

2 Answers 2

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.

$\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$ 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$ Dec 31, 2016 at 0:43
0
$\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]

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.