//***************************** EISENSTEIN *****************************

// Gives 0 if it IS reducible, p if irreducible with p, 1 elsewhise.
// Given polynomial may NOT be constant.

eisenstein := function (f);
  coef := Coefficients(f); c0 := coef[1]; cn := coef[#coef]; fc0 := [];
  if c0 eq 0 then return 0; end if;
  for i in Factorization(c0) do Append(~fc0,i[1]); end for;
  if exists(t){p:p in fc0 | not(IsDivisibleBy(c0,p^2)) and not(IsDivisibleBy(cn,p))
               and forall{m:m in coef[2..#coef-1] | IsDivisibleBy(m,p)}} then
    return t; 
    else return 1; 
  end if;
end function;

// exteis(f,B) tries to find an eisenstein-certificate for f(x+a), where
// |a| =< B.

exteis := function(f,B);
  eis := eisenstein(f); if eis ne 1 then return eis,0; end if;
  a:=1;
  while a le B do
    eisplus := eisenstein (Evaluate(f,x+a));
    eismin  := eisenstein (Evaluate(f,x-a));
    if eisplus  ne 1 then return eisplus,a;
    elif eismin ne 1 then return eismin,-a;
    end if;
    a +:= 1;
  end while;
  return 1,0;
end function;