P<x>:=PolynomialRing(Integers());

PRho:=function(n,f,s)
  x:=s;
  y:=s;
  i:=0;
  repeat
    i+:=1;
    x:=Evaluate(f,x) mod n;
    y:=Evaluate(f,y) mod n;
    y:=Evaluate(f,y) mod n;
    g:=GCD(y-x,n);
  until g gt 1;
  "Factorization of",n,"found after",i,"iterations:";
  return [g,n/g];
end function;

BPRho:=function(n,f,s)
  x:=s;
  y:=s;
  i:=1;
  j:=2;
  repeat
    if i eq j then
      j*:=2;
      y:=x;
    end if;
    i+:=1;
    x:=Evaluate(f,x) mod n;
    g:=GCD(x-y,n);
  until g gt 1;
  "Factorization of",n,"found after",i,"iterations:";
  return [g,n/g];
end function;