//**************************** BUNYAKOVSKY ****************************
r<y> := PolynomialRing(RealField()); print "r is", r;

// Bunya(f) tries to find enough prime or unit evaluations for a 
// Bunyakovsky-certificate.

Bunya := function (f);
  d := Degree(f);
  f0 := Evaluate(f,0);
  f1 := Evaluate(f,1);
  f2 := Evaluate(f,2);
  if   IsUnit(f0) or IsPrime(f0) then ies:=[[0,f0]]; posi:=3; negi:=-3;  
  elif IsUnit(f1) or IsPrime(f1) then ies:=[[1,f1]]; posi:=4; negi:=-2;
  elif IsUnit(f2) or IsPrime(f2) then ies:=[[2,f2]]; posi:=5; negi:=-1;
  else ies := []; posi := 3; negi := -1;
  end if; 
  while #ies le d do
    fposi := Evaluate(f,posi); fnegi := Evaluate(f,negi);
    if IsUnit(fposi) or IsPrime(fposi) then
      Append (~ies,[posi,fposi]);
      posi +:= 3;
    else posi +:= 1;
    end if;
    if IsUnit(fnegi) or IsPrime(fnegi) then 
      Append (~ies,[negi,fnegi]);
      negi -:= 3;
    else negi -:= 1;
    end if;
    end while;
  if #ies eq d+2 then ies := ies[1..#ies-1]; end if;
  return ies;
end function;  


// The m needed for Brillharts certificate can be computed in (at least) three
// different ways. The first one is found very quickly, but it is a bit rough.

computem1 := function(f);
  a0 := LeadingCoefficient(f);
  a1n := Coefficients(f)[1..Degree(f)];
  for i := 1 to #a1n do
    a1n[i] := Ceiling(Abs(a1n[i]/a0));
  end for;
  return 1 + Max(a1n);
end function;

// The second way to find m is more sophisticated and it is a good approximation
// the minimum value that m can take. It is computed fast.

computem2 := function(f);
x := Parent(f).1;
  coefs := Coefficients(f);
  f2 := Abs(coefs[#coefs])*x^(#coefs-1);
  for i := 1 to #coefs-1 do
    f2 -:= Abs(coefs[i]*x^(i-1));
  end for;
  m := 1;
  while true do
    q,r := Quotrem(f2,x-m);
    if r gt 0 and Minimum(Coefficients(q)) ge 0 then return m; end if;
    m +:= 1;
  end while;
end function;

// The third way is to really obtain m by looking at all complex roots of f.
// We now get the lowest possible value for m, but it can take some time to find
// this m.

computem3 := function(f);
  rf := Roots(r!f);
  return Ceiling(Maximum([Modulus(rf[i,1]) : i in [1..#rf]]));
end function;

// Brillhart(f) computes a Brillhart-certificate with m calculated in "computem3".

Brillhart := function(f);
  m := computem3(f);
  x0 := m+1;
  while true do
    fposx := Evaluate(f,x0); fnegx := Evaluate(f,-x0);
    if IsPrime(fposx)  then return  x0,fposx;  end if;
    if IsPrime(fnegx)  then return -x0,fnegx;  end if;
    x0 +:= 1;
  end while;
end function;

// Brillhart2(f) computes a Brillhart-certificate with m as in "computem2", so
// this can be more fast, but sometimes the x0-value found is higher.

Brillhart2 := function(f);
  m := computem2(f);
  x0 := m+1;
  while true do
    fposx := Evaluate(f,x0); fnegx := Evaluate(f,-x0);
    if IsPrime(fposx)  then return  x0,fposx;  end if;
    if IsPrime(fnegx)  then return -x0,fnegx;  end if;
    x0 +:= 1;
  end while;
end function;