y:=Indeterminate(Rationals,"y");
f:=y^4-y^3+3*y^2+2*y-1;
coef:=CoefficientsOfUnivariatePolynomial(f);

p:=5;
x:=Indeterminate(GF(p),"x");

RedResClass:=function(int,p)
if int>(p-1)/2 then return int-p;else return int;
fi;
end;

ModToInt:=function(pol)
local coef, coefint;
if pol=0*Z(p)*x^0 then return 0;
else
coef:=CoefficientsOfUnivariatePolynomial(pol);
p:=Characteristic(pol);
coefint:=List(coef,c->RedResClass(Int(c),p));
return Sum(List([1..Degree(pol)+1],i->coefint[i]*y^(i-1)));
fi;
end;

IntToMod:=function(pol,p)
local coef, coefint;
if pol=0*y^0 then return 0;
else
coef:=CoefficientsOfUnivariatePolynomial(pol);
coefint:=List(coef,c->c mod p);
return Sum(List([1..Degree(pol)+1],i->coefint[i]*x^(i-1)));
fi;
end;

ZassenhausBound:=function(pol)
local coef,a0,M;
coef:=CoefficientsOfUnivariatePolynomial(pol);
M:=Maximum(List(coef,AbsInt));
a0:=coef[Length(coef)];
return 1+M/AbsInt(a0);
end;

MignotteBound:=function(pol)
local B,k;
B:=ZassenhausBound(pol);
k:=QuoInt(Degree(pol),2);
return 2*Maximum(List([1..k],s->Binomial(k,s)*B^s));
end;
MB:=MignotteBound(f);
fm:=IntToMod(f,p);
deg:=Degree(f);
Q:=List([0..deg-1],i->CoefficientsOfUnivariatePolynomial(x^(i*5) mod fm));
I:=IdentityMat(deg);
D:=Z(p)^0*(Q-I);
Rank(D);#fornisce 2 quindi ho 2=4-2 fattori irriducibili per fred
ND:=NullspaceMat(D);#fornisce come base $(1,x^3)$

VectorToPol:=function(vec)
return Sum(List([1..deg],i->vec[i]*x^(i-1)));
end;

repeat
vec:=Random(ND);
upol:=VectorToPol(vec);

facs:=List(GF(p),s->Gcd(fm,upol-s));
fac:=Filtered(facs,x->Degree(x)>0);
until Length(fac)>1;

am:=fac[1];
bm:=fac[2];
gcdrep:=GcdRepresentation(am,bm);

a:=ModToInt(am);
b:=ModToInt(bm);
i:=1;
um:=gcdrep[1];
vm:=gcdrep[2];

repeat
dm:=IntToMod((f-a*b)/p^i,p);
ahm:=dm*vm mod am;
bhm:=Quotient(dm-ahm*bm,am);
a:=a+p^i*ModToInt(ahm);
b:=b+p^i*ModToInt(bhm);
am:=IntToMod(a,p);
bm:=IntToMod(b,p);
i:=i+1;
until p^i>2*MB;