How to extract specific terms from the polynomial?

Question

Consider the following polynomial:

expr = (1 + a[x] + c[x, aa] + c[x, aa] da[x, aa] + 
 da[x, aa]*da[x, aa] + 
 dc[x, aa, bb]*da[x, cc]*da[x, bb]*eps[aa, bb, cc, dd])^3 //

Expand

Among all the terms, I want to leave only those which contain exactly three a[x] or da[x,yy_] and one c or dc[x,yy_,zz_]; it is okay if they have eps[aa,bb,cc,dd]. At least one of the terms from expr matching these requirements is 3 a[x]^2 c[x, aa] da[x, aa].

Could you please tell me if this is possible in a more or less efficient way in Mathematica (in my realistic example, I deal with thousands of terms)?

Edit

This is my ugly realization. It uses the following answer:

ruleFixedPower[expr_, field_, dfield_, n_] := (Normal@
     Series[expr /. {factor : 
          field[x] | field[x, a_] | dfield[x, a_] | 
           dfield[x, a_, b_] :> factor*$T}, {$T, 0, n}] /. {$T -> 
      1}) - (Normal@
     Series[expr /. {factor : 
          field[x] | field[x, a_] | dfield[x, a_] | 
           dfield[x, a_, b_] :> factor*$T}, {$T, 0, 
       Max[0, n - 1]}] /. {$T -> 1})

Then,

expr1 = ruleFixedPower[expr, a, da, 3];
expr2 = ruleFixedPower[expr1, c, dc, 1]

6 a(x) da(x,bb) da(x,cc) dc(x,aa,bb) eps(aa,bb,cc,dd)+6 a(x) c(x,aa) da(x,aa)^2+3 a(x)^2 c(x,aa) da(x,aa)+6 c(x,aa) da(x,aa)^3

Unfortunately, it is slow as hell for my realistic tasks with many thousands of terms and more than two different structures, the powers of which have to be fixed. I may pre-select to have only the terms with the total degree equal to n, where n is the total degree of the desired output, but still, it is slow.

The derivative expansion would work much faster than this construction with Series, but it obviously introduces wrong expansion coefficients.

To illustrate the slowness, consider the following expansion:

expr = (1 + a[x] + c[x, aa] + a[x]*b[aa] + c[x, aa]* da[x, aa] + 
      dc[x, bb] dx[x, dd] + da[x, aa]*da[x, aa] + 
      dc[x, aa, bb]*da[x, cc]*da[x, bb]*eps[aa, bb, cc, dd])^16 // 
   Expand;
expr // Length

167637

(*Selecting terms not with the common power not higher than n*)
ruleCommonPower[expr_, 
  n_] := (Normal@
    Series[expr /. {factor : 
         a[x] | ca[x, a_] | da[x, a_] | dc[x, a_, b_] :> 
        factor*$T}, {$T, 0, 14}] /. {$T -> 1})
expr0 = ruleCommonPower[expr, 14]; // AbsoluteTiming
expr1 = ruleFixedPower[expr0, a, da, 10]; // AbsoluteTiming
expr2 = ruleFixedPower[expr1, c, dc, 4] // AbsoluteTiming

This code gets stuck.

Please edit your question to include the expression for the expected result. — Bob Hanlon, Aug 10 '23 at 20:57
@BobHanlon : I have provided the desired output together with my (very slow) solution for this problem. — John Taylor, Aug 10 '23 at 21:23
It is not clear how the terms of the shown desired output satisfies the requirement of "contain exactly three a[x] or da[x,yy_] and one dc[x,yy_,zz_]" — Bob Hanlon, Aug 10 '23 at 22:50
For ruleFixedPower, you ought to be able to use SeriesCoefficient and avoid subtracting two series. — Michael E2, Aug 10 '23 at 23:28
This is very confusing. In the description there is mention of a[x], da[x,yy_] , dc[x,yy_,zz_] but the putative term then given, 3 a[x]^2 c[x, aa] da[x, aa], has no dc[...], contrary to the requirement that there be exactly one (I assume this counting is by multiplicity). Before tackling this I'd want to see the actual objective in unambiguous terms. — Daniel Lichtblau, Aug 11 '23 at 00:14
@DanielLichtblau : excuse me, I have modified the description. — John Taylor, Aug 11 '23 at 06:15

Roland F · Answer 1 · 2023-08-11T07:09:19.633

I am not sure, if the following selection meets your demand, but its the way to implement such sieves

   Select[expr, ( Count[#, _a | _da] == 3 && Count[#, _d | _dc] == 1 &)]
 6 a[x] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, cc, dd] + 
     6 a[x] c[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] 
       eps[aa, bb, cc,dd] +
     6 c[x, aa] da[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] 
   eps[aa, bb, cc, dd] + 
     6 c[x, aa]^2 da[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, cc, dd] + 
     6 a[x] da[x, aa]^2 da[x, bb] da[x, cc] dc[x, aa, bb]  bb, cc, eps[aa,dd]

Eventually one has to check for the count of additional features as integer multiples, integer powers or

    Count[#,pattern,depth]

or specifying the variables pattern inside the given functions a by

 Count[#, a[x]|a[y],Infinity]

score 2 · Answer 2 · answered Aug 11 '23 at 07:31

You could try with pattern matching. The pattern will be complicated and needs to be protected against evaluation by "HoldPattern".

Then we will first search for terms with 3 a or da by the pattern: "___ (a[x] | da[x, _]) ___ (a[x] | da[x, _]) ___ (a[x] | da[x, _]) ___":

tmp = Cases[expr, 
  HoldPattern[___   (a[x] | da[x, _])  ___  (a[x] | 
      da[x, _])  ___  (a[x] | da[x, _])  ___]];

Then with search the results for c or dc by the pattern: "___ (c | dc[x, _, _]) ___":

Cases[tmp, HoldPattern[___ (c | dc[x, _, _]) ___]]
{6 a[x] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, cc, dd], 
 6 a[x] c[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, cc, 
   dd], 6 c[x, aa] da[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa,
    bb, cc, dd], 
 6 a[x] c[x, aa] da[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, 
   bb, cc, dd], 
 6 c[x, aa]^2 da[x, aa] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, 
   cc, dd], 
 6 a[x] da[x, aa]^2 da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, cc, 
   dd]}

score 1 · Accepted Answer · answered Aug 11 '23 at 23:05

One can create artificial powers in a new variable and then extract only those of interest from the expanded form. For this example we want to have exactly 3 of one type of variable and exactly 1 of another. We can take use E in the exponents of the first type, and integers for the second. Then we are seeking powers of 1+3*E.

expr = (1 + a[x] + c[x, aa] + c[x, aa] da[x, aa] + 
     da[x, aa]*da[x, aa] + 
     dc[x, aa, bb]*da[x, cc]*da[x, bb]*eps[aa, bb, cc, dd])^3;
powers = {E, E, 1, 1};
reps = Thread[{a[x], da[x, yy_], dc[x, yy_, zz_], c[x, yy_]} -> 
    t^powers*{a[x], da[x, yy], dc[x, yy, zz], c[x, yy]}];
e2 = Expand[expr] /. reps;
terms = Cases[e2, t^(1 + 3E)aa_ :> aa]
(* Out[213]= {3 a[x]^2 c[x, aa] da[x, aa], 6 a[x] c[x, aa] da[x, aa]^2, 
 6 c[x, aa] da[x, aa]^3, 
 6 a[x] da[x, bb] da[x, cc] dc[x, aa, bb] eps[aa, bb, cc, dd]} *)

Thanks! Looks very elegant, and it may be extended to an arbitrary sets of variables. — John Taylor, Aug 12 '23 at 12:28

How to extract specific terms from the polynomial?

3 Answers3