4
$\begingroup$

I have one case (Nn=1) for 7 Polynomials as following:

p1= a[0]*b[0] + c[0]*d[0] + e[0]*f[0] + g[0]*h[0];
p2= a[1]*c[1] + b[1]*d[1] + e[1]*g[1] + f[1]*h[1]; 
p3= a[2]*d[2] + b[2]*c[2] + e[2]*h[2] + f[2]*g[2]; 
p4= a[3]*e[3] + b[3]*f[3] + c[3]*g[3] + d[3]*h[3]; 
p5= a[4]*f[4] + b[4]*e[4] + c[4]*h[4] + d[4]*g[4]; 
p6= a[5]*g[5] + b[5]*h[5] + c[5]*e[5] + d[5]*f[5]; 
p7= a[6]*h[6] + b[6]*g[6] + c[6]*f[6] + d[6]*e[6];

and I want to pick up 4 terms from the 7*4=28 terms such as a[0]*b[0], c[0]*d[0] and so on, which gives me the output in the form of a[x1]*b[x2]*c[x3]*d[x4]*e[x5]*f[x6]*g[x7]*h[x8] (we can later write it as FF[x1,x2,x3,x4,x5,x6,x7,x8]).

The way I do is as following:

outputstemp={};
outputs={};
FFnCaseList={a[c1_]*b[c2_]*c[c3_]*d[c4_]*e[c5_]*f[c6_]*g[c7_]*h[c8_]->FF[c1,c2,c3,c4,c5,c6,c7,c8], a[_]->0, b[_]->0, c[_]->0,  d[_]->0, e[_]->0, f[_]->0, g[_]->0, h[_]->0,  x_[_]^n_->0};

AllRows = (p1+p2+p3+p4+p5+p6+p7)^4; 
AppendTo[outputstemp, ExpandAll[AllRows]];

Timing[AppendTo[outputs, outputstemp/.FFnCaseList];]

the time is {10.2344, Null}

So is there any quick way to do this? If I have many different cases (for example Nn=6000) of 7 Polynomials, then I can make a loop that run similar things Nn times (6000*10s). That will takes roughly 16.7 hours!!!

Any comments or suggestions are appreciated! Thank you very much!

$\endgroup$

3 Answers 3

3
$\begingroup$

I would use Coefficient after mopdifying so that terms of the desired form have a unique power in a new variable.

p1 = a[0]*b[0] + c[0]*d[0] + e[0]*f[0] + g[0]*h[0];
p2 = a[1]*c[1] + b[1]*d[1] + e[1]*g[1] + f[1]*h[1];
p3 = a[2]*d[2] + b[2]*c[2] + e[2]*h[2] + f[2]*g[2];
p4 = a[3]*e[3] + b[3]*f[3] + c[3]*g[3] + d[3]*h[3];
p5 = a[4]*f[4] + b[4]*e[4] + c[4]*h[4] + d[4]*g[4];
p6 = a[5]*g[5] + b[5]*h[5] + c[5]*e[5] + d[5]*f[5];
p7 = a[6]*h[6] + b[6]*g[6] + c[6]*f[6] + d[6]*e[6];

Multiply each variable by t^2^j where j is the index of that variable in the list {a,b,c,d,e,f,g,h}.

polys = {p1, p2, p3, p4, p5, p6, p7} /. 
  Thread[{v_a, v_b, v_c, v_d, v_e, v_f, v_g, v_h} -> 
    v*t^(2^Range[0, 7])]

(* Out[13]= {t^3 a[0] b[0] + t^12 c[0] d[0] + t^48 e[0] f[0] + 
  t^192 g[0] h[0], 
 t^5 a[1] c[1] + t^10 b[1] d[1] + t^80 e[1] g[1] + t^160 f[1] h[1], 
 t^6 b[2] c[2] + t^9 a[2] d[2] + t^96 f[2] g[2] + t^144 e[2] h[2], 
 t^17 a[3] e[3] + t^34 b[3] f[3] + t^68 c[3] g[3] + t^136 d[3] h[3], 
 t^18 b[4] e[4] + t^33 a[4] f[4] + t^72 d[4] g[4] + t^132 c[4] h[4], 
 t^20 c[5] e[5] + t^40 d[5] f[5] + t^65 a[5] g[5] + t^130 b[5] h[5], 
     t^24 d[6] e[6] + t^36 c[6] f[6] + t^66 b[6] g[6] + t^129 a[6] h[6]} *)

Pull out terms that have t to the power 1+2+4+...+128, or 255. In order to speed the conversion to the FF[...] notation I replace sums and products by lists (this speeds otherwise slow pattern match). In a more general setting one might have to canonicalize the ordering using Sort.

Timing[
 allvarsmonoms = 
  Expand[Coefficient[(Apply[Plus, polys])^4, t^Total[2^Range[0, 7]]]];
 res0 = Apply[List, allvarsmonoms];
 res1 = res0 /. {Plus -> List, Times -> List};
 res = Apply[Plus, 
   res1 /. {num_, a[v1_], b[v2_], c[v3_], d[v4_], e[v5_], f[v6_], 
      g[v7_], h[v8_]} -> num*FF[v1, v2, v3, v4, v5, v6, v7, v8]];]
Length[res]

(* Out[71]= {0.0625, Null}

Out[72]= 105 *)
$\endgroup$
3
  • $\begingroup$ Thank you very much! If variable in the list is {a,b,c,d,e,f}, should I use v*t^(2^Range[0, 5])? $\endgroup$
    – Xuemei
    Apr 21, 2020 at 17:25
  • 1
    $\begingroup$ Yes. For cleaner code, I should have defined the variable symbols in a list as vars and then used Length[vars] for that exponent. $\endgroup$ Apr 21, 2020 at 18:04
  • $\begingroup$ Thank you! that speed up a lot. Very interesting! $\endgroup$
    – Xuemei
    Apr 21, 2020 at 18:17
5
$\begingroup$

I am sure there is a less awkward method than the following proposal:

n = 4;
expr = p1 + p2 + p3 + p4 + p5 + p6 + p7;

Nest[Map[If[Signature[Head /@ (List @@ #)] == 0, 0, #] &, Expand[expr * #]] &, 1, n] /. 
(a[c1_] b[c2_] c[c3_] d[c4_] e[c5_] f[c6_] g[c7_] h[c8_]) :>
FF[c1, c2, c3, c4, c5, c6, c7, c8]

which yields 24 a[0] b[0] c[0] d[0] e[0] f[0] g[0] h[0] + ... 24 a[6] b[6] c[6] d[6] e[6] f[6] g[6] h[6] in about a hundredth of the time taken for the OP's original code.

$\endgroup$
5
  • $\begingroup$ wow, that's amazing. Well, i didn't saw hundredth of time from my computer but factor of 20 is quite good! Thank you! $\endgroup$
    – Xuemei
    Apr 20, 2020 at 18:44
  • 1
    $\begingroup$ changing the replacement rule to {Times -> (# FF[##2] &), Alternatives[a, b, c, d, e, f, g, h] -> Identity} seems to double the speed. (+1) $\endgroup$
    – kglr
    Apr 20, 2020 at 19:36
  • $\begingroup$ @kglr, yes, that's cool! Thank you very much! $\endgroup$
    – Xuemei
    Apr 20, 2020 at 19:53
  • $\begingroup$ @Xuemei, your "additional" question should surely be a separate question. $\endgroup$ Apr 21, 2020 at 0:50
  • 1
    $\begingroup$ @J.M., thank you very much and I have posted it as a separate question here. Could you have a look? would be extremely helpful and thank you a lot! $\endgroup$
    – Xuemei
    Apr 21, 2020 at 7:02
3
$\begingroup$
expr  = p1 + p2 + p3 + p4 + p5 + p6 + p7;

res = Total[Multinomial[1, 1, 1, 1] FF @@@ 
       (ReplaceAll[ _[a_] :> a] @ Select[DuplicateFreeQ @* Map[Head]]@
         Apply[Sort@*Join, Subsets[expr /. Plus|Times -> List, {4}],  {-4}])]; // 
           RepeatedTiming // First
 .16

this gives the same result as J.M.'s approach:

res2 = Nest[Map[If[Signature[Head /@ (List @@ #)] == 0, 0, #] &, 
        Expand[expr*#]] &, 1, n] /. 
         (a[c1_] b[c2_] c[c3_] d[c4_] e[c5_] f[c6_] g[c7_] h[c8_]) :> 
           FF[c1, c2, c3, c4, c5, c6, c7, c8]; //  RepeatedTiming // First
0.54
res == res2
True
$\endgroup$
1
  • $\begingroup$ thank you very much and I have posted it as a separate question here. Could you have a look? would be extremely helpful and thank you a lot! $\endgroup$
    – Xuemei
    Apr 21, 2020 at 7:05

Your Answer

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

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