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"]
       ]
    ]