Return to Revisions

I included this code as a friendly version to the excellent code from @anderstood. This script handles the cases with number of variables from 1 to 4. Extension to more variables can be easily done by extending the Switch span.

Clear[ExtractFactors]
ExtractFactors[f_, vars_] := Module[{rf, ef},
  rf = CoefficientRules[f, vars];
  ef = Map[First, rf];
  Return[ef]
]

pol = 1 + a b d + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d;
pol = 1 - a^2 b^2 c d^2;
pol = 1 + a b d + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d;
pol = 1 - a b c d;
pol = 1 - a^2 b^2 c^2 d^2;
pol = 1 - a b c;
pol = 1 - a - a^2;
pol = 1 - a b;
pol = a b + a c + c b - a b c;
pol = 1 - a b c d;
pol = 1 - a + a^2;
pol = 1 - a b c d e;
pol = 1 - a b + a^2 b^2;
pol = 1 + a b d + b + a^2 b + a b^2 - b c + 2 a b c - a^2 b c - 4 a b c d^2;

vars = Variables[pol];
monom = ExtractFactors[pol, vars];
numvars = Length[vars];
mexp = Table[Max[Transpose[monom][[All ;; k]]], {k, 1, numvars}];
table = Switch[Length[mexp],
   1, amax = mexp[[1]]; Table[If[i + j > amax || i + j < 2, 0, 1]*a^i (1 - a)^j, {i, 0, amax}, {j, 0, amax}],
   2, amax = mexp[[1]]; bmax = mexp[[2]]; Table[If[i + j > amax || k + l > bmax || i + j + k + l < 2, 0, 1]*a^i (1 - a)^j*b^k*(1 - b)^l, {i, 0, amax}, {j, 0, amax}, {k, 0, bmax}, {l, 0, bmax}],
   3, amax = mexp[[1]]; bmax = mexp[[2]]; cmax = mexp[[3]];Table[If[i + j > amax || k + l > bmax || m + n > cmax || i + j + k + l + m + n < 2, 0, 1]*a^i (1 - a)^j*b^k*(1 - b)^l*c^m*(1 - c)^n, {i, 0, amax}, {j, 0, amax}, {k, 0, bmax}, {l, 0, bmax}, {m, 0, cmax}, {n, 0, cmax}],
   4, amax = mexp[[1]]; bmax = mexp[[2]]; cmax = mexp[[3]]; dmax = mexp[[4]]; Table[If[i + j > amax || k + l > bmax || m + n > cmax || o + p > dmax || i + j + k + l + m + n + o + p < 2, 0, 1]*a^i (1 - a)^j*b^k*(1 - b)^l*c^m*(1 - c)^n*d^o*(1 - d)^p, {i, 0, amax}, {j, 0, amax}, {k, 0, bmax}, {l, 0, bmax}, {m, 0, cmax}, {n, 0, cmax}, {o, 0, dmax}, {p, 0, dmax}], 
  _, Print["Num. of variables should be less than 5"]];
If[Length[mexp] < 5,
   monomials = Union[Flatten[table]];
   np = Length[monomials];
   A = Table[Subscript[\[Alpha], i], {i, 1, np}];
   restrsA = Thread[A >= 0];
   dif = A.monomials - pol;
   coefs = CoefficientRules[dif, vars];
   rels = Map[Last, coefs];
   equs = Thread[rels == 0];
   sol = Quiet@FindInstance[Join[equs, restrsA], A];
   If[Length[sol] > 0,
      polfound = A.monomials /. sol[[1]]; 
      Print[polfound]; 
      Print[Simplify[polfound - pol /. sol[[1]]]], Print["No match found"]
   ]
]