a = (x^4) + (x^2) + 1

punkte = {-1, 0, 1}

werte = a /. x -> punkte

divisors[z_] := Module[{d = Divisors[z]}, Join[d, -d]]

divisorliste = 
      Table[divisors[werte[[k]]], {k, Length[werte]}]

tabelle = 
      Partition[Flatten[Outer[List, Apply[Sequence, divisorliste]]], 
        3]

poltabelle = 
      Table[{Expand[
            InterpolatingPolynomial[
              Transpose[{{-1, 0, 1}, tabelle[[k]]}], x]]}, {k, 
          Length[tabelle]}]

reste = 
      Table[PolynomialRemainder[a, poltabelle[[k, 1]], x], {k, 
          Length[poltabelle]}]

pos = Position[reste, 0]

Table[poltabelle[[pos[[k]], 1]], {k, Length[pos]}]