Quelle polyrat.gi Sprache: Shell

#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Alexander Hulpke.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file contains functions for polynomials over the rationals
##

#############################################################################
##
#F APolyProd(<a>,,) . . . . . . . . . . polynomial product a*b mod p
##
##return a*b mod p;
InstallGlobalFunction(APolyProd,function(a,b,p)
local ac,bc,i,j,pc,pv,ci,fam;

 ci:=CIUnivPols(a,b);
 fam:=FamilyObj(a);
 a:=CoefficientsOfLaurentPolynomial(a);
 b:=CoefficientsOfLaurentPolynomial(b);

 pv:=a[2]+b[2];
 #pc:=ProductCoeffs(a.coefficients,b.coefficients);
 ac:=List(a[1],i->i mod p);
 bc:=List(b[1],i->i mod p);
 if Length(ac)>Length(bc) then
 pc:=ac;
 ac:=bc;
 bc:=pc;
 fi;
 # prepare result list
 pc:=[];
 for i in [1..Length(ac)+Length(bc)-1] do
 pc[i]:=0;
 od;
 for i in [1..Length(ac)] do
 #pc:=SumCoeffs(pc,i*bc) mod p;
 for j in [1..Length(bc)] do
 pc[i+j-1]:=(pc[i+j-1]+ac[i]*bc[j]) mod p;
 od;
 od;
 pv:=pv+RemoveOuterCoeffs(pc,fam!.zeroCoefficient);
 return LaurentPolynomialByExtRepNC(fam,pc,pv,ci);
end);

#############################################################################
##
#F BPolyProd(<a>,,<m>,) . . . . . . polynomial product a*b mod m mod p
##
## return EuclideanRemainder(PolynomialRing(Rationals),a*b mod p,m) mod p;
InstallGlobalFunction(BPolyProd,function(a,b,m,p)
local ac,bc,mc,i,j,pc,ci,f,fam;

 ci:=CIUnivPols(a,b);
 fam:=FamilyObj(a);
 a:=CoefficientsOfLaurentPolynomial(a);
 b:=CoefficientsOfLaurentPolynomial(b);
 m:=CoefficientsOfLaurentPolynomial(m);
 # we shift as otherwise the mod will mess up valuations (should occur
 # rarely anyhow)
 ac:=List(a[1],i->i mod p);
 ac:=ShiftedCoeffs(ac,a[2]);
 bc:=List(b[1],i->i mod p);
 bc:=ShiftedCoeffs(bc,b[2]);
 mc:=List(m[1],i->i mod p);
 mc:=ShiftedCoeffs(mc,m[2]);
 ReduceCoeffsMod(ac,mc,p);
 ShrinkRowVector(ac);
 ReduceCoeffsMod(bc,mc,p);
 ShrinkRowVector(bc);
 if Length(ac)>Length(bc) then
 pc:=ac;
 ac:=bc;
 bc:=pc;
 fi;
 pc:=[];
 f:=false;
 for i in ac do
 if f then
 # only do it 2nd time (here to avoin doing once too often)
 bc:=ShiftedCoeffs(bc,1);
 ReduceCoeffsMod(bc,mc,p);
 ShrinkRowVector(bc);
 else
 f:=true;
 fi;
 #pc:=SumCoeffs(pc,i*bc) mod p;
 for j in [Length(pc)+1..Length(bc)] do
 pc[j]:=0;
 od;
 for j in [1..Length(bc)] do
 pc[j]:=(pc[j]+i*bc[j] mod p) mod p;
 od;
 ShrinkRowVector(pc);
 ReduceCoeffsMod(pc,mc,p);
 ShrinkRowVector(pc);
 od;
 p:=RemoveOuterCoeffs(pc,fam!.zeroCoefficient);
 return LaurentPolynomialByExtRepNC(fam,pc,p,ci);
end);

#############################################################################
##
#F ApproxRational: approximativ k"urzen
##
BindGlobal("ApproxRational",function(r,s)
local n,d,u;
 n:=NumeratorRat(r);
 d:=DenominatorRat(r);
 u:=LogInt(d,10)-s+1;
 if u>0 then
 u:=10^u;
 return QuoInt(n,u)/QuoInt(d,u);
 else
 return r;
 fi;
end);

#############################################################################
##
#F ApproximateRoot(<num>,<n>[,<digits>]) . . approximate th n-th root of num
## numerically with a denominator of 'digits' digits.
##
BIND_GLOBAL( "APPROXROOTS", NEW_SORTED_CACHE(false) );

BindGlobal("ApproximateRoot",function(arg)
local r,e,f,store,maker;
 r:=arg[1];
 e:=arg[2];
 if Length(arg)>2 then
 f:=arg[3];
 else
 f:=10;
 fi;

 store:= e<=10 and IsInt(r) and 0<=r and r<=100 and f=10;

 maker := function()
 local x,nf,lf,c,letzt;

 x:=RootInt(NumeratorRat(r),e)/RootInt(DenominatorRat(r),e);
 nf:=r;
 c:=0;
 letzt:=[];
 repeat
 lf:=nf;
 Add(letzt,x);
 x:=ApproxRational(1/e*((e-1)*x+r/(x^(e-1))),f+6);
 nf:=AbsoluteValue(x^e-r);
 if nf=0 then
 c:=6;
 else
 if nf>lf then
 lf:=nf/lf;
 else
 lf:=lf/nf;
 fi;
 if lf<2 then
 c:=c+1;
 else
 c:=0;
 fi;
 fi;
 # until 3 times no improvement
 until c>2 or x in letzt;
 return x;
 end;

 if store then
 return GET_FROM_SORTED_CACHE(APPROXROOTS, [e,r], maker);
 fi;
 return maker();
end);

#############################################################################
##
#F ApproxRootBound(f) Numerical approximation of RootBound (better, but
## may fail)
##
BindGlobal("ApproxRootBound",function(f)
local x,p,tp,diff,app,d,scheit,v,nkon;

 x:=IndeterminateNumberOfLaurentPolynomial(f);
 p:=CoefficientsOfLaurentPolynomial(f);
 if p[2]<0 or not ForAll(p[1],IsRat) then
 # avoid complex conjugation etc.
 Error("only yet implemented for rational polynomials");
 fi;

 # eliminate valuation
 f:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,p[1],x);
 x:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,[0,1],x);

 # probably first test, whether polynomial should be inverted. However,
 # we expect roots larger than one.
 d:=DegreeOfLaurentPolynomial(f);
 f:=Value(f,1/x)*x^d;
 app:=1/2;
 diff:=1/4;
 nkon:=true;
 repeat
 # pol, whose roots are the 1/app of the roots of f
 tp:=Value(f,x*app);
 tp:=CoefficientsOfLaurentPolynomial(tp)[1];
 tp:=tp/tp[1];
 tp:=List(tp,i->ApproxRational(i,10));
 # now check, by using the Lehmer/Schur method, whether tp has a root
 # in the unit circle, i.e. f has a root in the app-circle
 repeat
 scheit:=false;
 p:=tp;
 repeat
 d:=Length(p);
 # compute T[p]=\bar a_n p-a_0 p*, everything rational.
 p:=p[1]*p-p[d]*Reversed(p);
 p:=List(p,i->ApproxRational(i,10));
 d:=Length(p);
 while d>1 and p[d]=0 do
 Unbind(p[d]);
 d:=d-1;
 od;
 v:=p[1];
 if v=0 then
 scheit:=nkon;
 fi;
 nkon:=ForAny(p{[2..Length(p)]},i->i<>0);
 until v<=0;
 if scheit then
 # we fail due to rounding errors
 return fail;
 else
 if v<0 then
 # zero in the unit circle, app smaller
 app:=app-diff;
 else
 # no zero in the unit circle, app larger
 app:=app+diff;
 fi;
 fi;
 until not scheit;
 diff:=diff/2;
 # until good circle found, which does not contain roots.
 until v=0 and (1-app/(app+diff))<1/40;

 # revert last enlargement and add accuracy to be secure
 app:=app-2*diff;
 return 1/app+1/20;
end);

#############################################################################
##
#F RootBound(<f>) . . . . bound for absolute value of (complex) roots of f
##
InstallGlobalFunction(RootBound,function(f)
local a,b,c,d;
 # valuation gives only 0 as zero, this can be neglected
 f:=CoefficientsOfLaurentPolynomial(f)[1];
 # normieren
 f:=f/Last(f);
 f:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,f,1);

 a:=ApproxRootBound(f);
 # did the numerical part fail?
 if a=fail then
 c:=CoefficientsOfLaurentPolynomial(f)[1];
 c:=List(c,AbsInt);
 d:=Length(c);
 a:=Maximum(1,Sum(c{[1..d-1]}));
 b:=1+Maximum(c);
 if b<a then
 a:=b;
 fi;
 b:=Maximum(List([1..d-1],i->RootInt(d*Int(AbsInt(c[d-i])+1/2),i)+1));
 if b<a then
 a:=b;
 fi;
 if ForAll(c,i->i<>0) then
 b:=List([3..d],i->2*AbsInt(c[i-1]/c[i]));
 Add(b,AbsInt(c[1]/c[2]));
 b:=Maximum(b);
 if b<a then
 a:=b;
 fi;
 fi;
 b:=Sum([1..d-1],i->AbsInt(c[i]-c[i+1]))+AbsInt(c[1]);
 if b<a then
 a:=b;
 fi;
 b:=1/20+Maximum(List([1..d-1],
 i->RootInt(Int(AbsInt(c[d-i]/Binomial(d-1,i))+1/2),i)+1))
 /(ApproximateRoot(2,d-1)-1)+10^(-10);
 if b<a then
 a:=b;
 fi;

 fi;
 return a;
end);

#############################################################################
##
#F BombieriNorm(<pol>) . . . . . . . . . . . . compute weighted Norm [pol]_2
##
InstallGlobalFunction(BombieriNorm,function(f)
local c,i,n,s;
 c:=CoefficientsOfLaurentPolynomial(f);
 c:=ShiftedCoeffs(c[1],c[2]);
 n:=Length(c)-1;
 s:=0;
 for i in [0..n] do
 s:=s+AbsInt(c[i+1])^2/Binomial(n,i);
 od;
 return ApproximateRoot(s,2);
end);

#############################################################################
##
#F MinimizedBombieriNorm(<pol>) . . . . . . . minimize weighted Norm [pol]_2
## by shifting roots
##
InstallMethod(MinimizedBombieriNorm,true,
 [IsPolynomial and IsRationalFunctionsFamilyElement],0,
function(f)
local bn,bnf,a,b,c,d,bb,bf,bd,step,x,cnt;

 step:=1;
 bb:=infinity;
 bf:=f;
 bd:=0;
 bn:=[];
 # evaluation of norm, including storing it (avoids expensive double evals)
 bnf := function(dis)
 local p,g;
 p:=Filtered(bn,i->i[1]=dis);
 if p=[] then
 g:=Value(f,x+dis);
 p:=[dis,BombieriNorm(g)];
 Add(bn,p);
 if bb>p[2] then
 # note record
 bf:=g;
 bb:=p[2];
 bd:=dis;
 fi;
 return p[2];
 else
 return p[1][2];
 fi;
 end;

 x:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,[0,1],
 IndeterminateNumberOfLaurentPolynomial(f));
 d:=0;
 cnt:=0;
 repeat
 cnt:=cnt+1;
 Info(InfoPoly,2,"Minimizing BombieriNorm, x->x+(",d,")");
 # local parabola approximation
 a:=bnf(d-step);
 b:=bnf(d);
 c:=bnf(d+step);
 if ab and c>b) and (a+c<>2*b) then
 a:=-(c-a)/2/(a+c-2*b)*step;
 # stets aufrunden (wir wollen weg)
 a:=step*Int(AbsInt(a)/step+1)*SignInt(a);
 if a=0 then
 Error("???");
 else
 d:=d+a;
 fi;
 fi;
 until (a>b and c>b) # no better can be reached
 or cnt>6
 or ForAll([d-1,d,d+1],i->Filtered(bn,j->j[1]=i)<>[]); #or loop
 # best value
 return [bf,bd];
end);

#############################################################################
##
#F BeauzamyBound(<pol>) . . . . . Beauzamy's Bound for Factors Coefficients
## cf. JSC 13 (1992), 463-472
##
BindGlobal("BeauzamyBound",function(f)
local n;
 n:=DegreeOfLaurentPolynomial(f);
 return Int(
 # the strange number in the next line is an (upper) rational approximation
 # for 3^{3/4}/2/\sqrt(\pi)
 643038/1000000*ApproximateRoot(3^n,2)/ApproximateRoot(n,2)*BombieriNorm(f))+1;
end);

#############################################################################
##
#F OneFactorBound(<pol>) . . . . . . . . . . . . Bound for one factor of pol
##
InstallGlobalFunction(OneFactorBound,function(f)
local d,n;
 n:=DegreeOfLaurentPolynomial(f);
 if n>=3 then
 # Single factor bound of Beauzamy, Trevisan and Wang (1993)
 return Int(10912/10000*(ApproximateRoot(2^n,2)/ApproximateRoot(n^3,8)
 *(ApproximateRoot(BombieriNorm(f),2))))+1;
 else
 # Mignotte's single factor bound
 d:=QuoInt(n,2);
 return
 Binomial(d,QuoInt(d,2))
 *(1+RootInt(Sum(CoefficientsOfLaurentPolynomial(f)[1],
 i->i^2),2));
 fi;
end);

#############################################################################
##
#F PrimitivePolynomial(<f>) . . . remove denominator and coefficients gcd
##
InstallMethod(PrimitivePolynomial,"univariate polynomial",true,
 [IsUnivariatePolynomial],0,
function(f)
local lcm, c, fc,fac;

 fc:=CoefficientsOfLaurentPolynomial(f)[1];
 # compute lcm of denominator
 lcm := 1;
 for c in fc do
 lcm := LcmInt(lcm,DenominatorRat(c));
 od;

 # remove all denominators
 f := f*lcm;
 fac:=1/lcm;
 fc:=CoefficientsOfLaurentPolynomial(f)[1];

 # remove gcd of coefficients
 if Length(fc)>0 then
 fc:=Gcd(fc);
 else
 fc:=1;
 fi;
 fac:=fac*fc;
 return [f*(1/fc),fac];

end);

BindGlobal("PrimitiveFacExtRepRatPol",function(e)
local d,lcm,i,fac;
 d:=e{[2,4..Length(e)]};
 if not ForAll(d,IsRat) then
 TryNextMethod();
 fi;
 lcm:=1;
 for i in d do
 lcm := LcmInt(lcm,DenominatorRat(i));
 od;
 fac:=1/lcm;
 d:=d*lcm;
 if Length(d)>0 then
 fac:=fac*Gcd(d);
 fi;
 return fac;
end);

InstallMethod(PrimitivePolynomial,"rational polynomial",true,
 [IsPolynomial],0,
function(f)
local e,fac;
 e:=ExtRepPolynomialRatFun(f);
 fac:=PrimitiveFacExtRepRatPol(e);
 return [f/fac,fac];
end);

#############################################################################
##
#F BeauzamyBoundGcd(<f>,<g>) . . . . . Beauzamy's Bound for Gcd Coefficients
##
## cf. JSC 13 (1992),463-472
##
BindGlobal("BeauzamyBoundGcd",function(f,g)
local n, A, B,lf,lg;

 lf:=LeadingCoefficient(f);
 if not IsOne(lf) then
 f:=f/lf;
 fi;

 lg:=LeadingCoefficient(f);
 if not IsOne(lg) then
 g:=g/lg;
 fi;
 n := DegreeOfLaurentPolynomial(f);
 # the strange number in the next line is an (upper) rational
 # approximation for 3^{3/4}/2/\sqrt(\pi)
 A := Int(643038/1000000
 * ApproximateRoot(3^n,2)/ApproximateRoot(n,2)
 * BombieriNorm(f))+1;

 # the strange number in the next line is an (upper) rational
 # approximation for 3^{3/4}/2/\sqrt(\pi)
 n := DegreeOfLaurentPolynomial(g);
 B := Int(643038/1000000
 * ApproximateRoot(3^n,2)/ApproximateRoot(n,2)
 * BombieriNorm(g))+1;
 B:=Minimum(A,B);
 if not (IsOne(lf) or IsOne(lg)) then
 B:=B*GcdInt(lf,lg);
 fi;
 return B;

end);

#############################################################################
##
#F RPGcdModPrime(<f>,<g>,,<a>,<brci>) . . gcd mod 
##
BindGlobal("RPGcdModPrime",function(f,g,p,a,brci)
local fam,gcd, u, v, w, val, r, s;

 fam:=CoefficientsFamily(FamilyObj(f));
 f:=CoefficientsOfLaurentPolynomial(f);
 g:=CoefficientsOfLaurentPolynomial(g);
 # compute in the finite field F_
 val := Minimum(f[2], g[2]);
 s := ShiftedCoeffs(f[1],f[2]-val);
 r := ShiftedCoeffs(g[1],g[2]-val);
 ReduceCoeffsMod(s,p); ShrinkRowVector(s);
 ReduceCoeffsMod(r,p); ShrinkRowVector(r);

 # compute the gcd
 u := r;
 v := s;
 while 0 < Length(v) do
 w := v;
 ReduceCoeffsMod(u,v,p);
 ShrinkRowVector(u);
 v := u;
 u := w;
 od;
 #gcd := u * (a/Last(u));
 gcd:=u;
 MultVector(gcd,a/Last(u));
 ReduceCoeffsMod(gcd,p);

 # and return the polynomial
 return LaurentPolynomialByCoefficients(fam,gcd,val,brci);

end);

BindGlobal("RPGcdCRT",function(f,p,g,q,ci)
local min, cf, lf, cg, lg, i, P, m, r, fam;

 fam := CoefficientsFamily(FamilyObj(f));
 f:=CoefficientsOfLaurentPolynomial(f);
 g:=CoefficientsOfLaurentPolynomial(g);
 # remove valuation
 min := Minimum(f[2],g[2]);
 if f[2] <> min then
 cf := ShiftedCoeffs(f[1],f[2] - min);
 else
 cf := ShallowCopy(f[1]);
 fi;
 lf := Length(cf);
 if g[2] <> min then
 cg := ShiftedCoeffs(g[1],g[2] - min);
 else
 cg := ShallowCopy(g[1]);
 fi;
 lg := Length(cg);

 # use chinese remainder
 r := [ p,q ];
 P := p * q;
 m := P/2;
 for i in [ 1 .. Minimum(lf,lg) ] do
 cf[i] := ChineseRem(r,[ cf[i],cg[i] ]);
 if m < cf[i] then cf[i] := cf[i] - P; fi;
 od;
 if lf < lg then
 for i in [ lf+1 .. lg ] do
 cf[i] := ChineseRem(r,[ 0,cg[i] ]);
 if m < cf[i] then cf[i] := cf[i] - P; fi;
 od;
 elif lg < lf then
 for i in [ lg+1 .. lf ] do
 cf[i] := ChineseRem(r,[ cf[i],0 ]);
 if m < cf[i] then cf[i] := cf[i] - P; fi;
 od;
 fi;

 # return the polynomial
 return LaurentPolynomialByCoefficients(fam,cf,min,ci);

end);

BindGlobal("RPGcd1",function(t,a,f,g)
local G, P, l, m, i;

 # will hold the product of primes use so far
 t.modulo := t.prime;

 # <G> will hold the approximation of the gcd
 G := t.gcd;

 # use next prime until we reach the Beauzamy bound
 while t.modulo < t.bound do
 repeat t.prime := NextPrimeInt(t.prime); until a mod t.prime <> 0;

 # compute modular gcd
 t.gcd := RPGcdModPrime(f,g,t.prime,a,t.brci);
 Info(InfoPoly,3,"gcd mod ",t.prime," = ",t.gcd);

 # if the degree of <C> is smaller we started with wrong 
 if DegreeOfLaurentPolynomial(t.gcd)
 < DegreeOfLaurentPolynomial(G)
 then
 Info(InfoPoly,3,"found lower degree,restarting");
 return false;
 fi;

 # if the degrees of <C> and <G> are equal use chinese remainder
 if DegreeOfLaurentPolynomial(t.gcd)
 = DegreeOfLaurentPolynomial(G)
 then
 P := G;
 G := RPGcdCRT(P,t.modulo,t.gcd,t.prime,t.brci);
 t.modulo := t.modulo * t.prime;
 Info(InfoPoly,3,"gcd mod ",t.modulo," = ",G);
 if G = P then
 t.correct := IsZero(f mod G) and IsZero(g mod G);
 if t.correct then
 Info(InfoPoly,3,"found correct gcd");
 t.gcd := G;
 return true;
 fi;
 fi;
 fi;
 od;

 # get <G> into the -<t.modulo>/2 to +<t.modulo> range
 G:=CoefficientsOfLaurentPolynomial(G);
 l := [];
 m := t.modulo/2;
 for i in [ 1 .. Length(G[1]) ] do
 if m < G[1][i] then
 l[i] := G[1][i] - t.modulo;
 else
 l[i] := G[1][i];
 fi;
 od;
 G := LaurentPolynomialByExtRepNC(FamilyObj(f),l,G[2],t.brci);
 Info(InfoPoly,3,"gcd mod ",t.modulo," = ",G);

 # check if <G> is correct but return 'true' in any case
 t.correct := IsZero(f mod G) and IsZero(g mod G);
 t.gcd := G;
 return true;

end);

BindGlobal("RPIGcd", function(f,g)
local a,t;

 # special case zero:
 if IsZero(f) then
 return g;
 elif IsZero(g) then
 return f;
 fi;
 # compute the Beauzamy bound for the gcd
 t := rec(prime := 1000);
 t.brci:=CIUnivPols(f,g);

 t.bound := 2 * Int(BeauzamyBoundGcd(f,g)+1);
 Info(InfoPoly,3,"Beauzamy bound = ",t.bound/2);

 # avoid gcd of leading coefficients
 a := GcdInt(LeadingCoefficient(f),LeadingCoefficient(g));
 repeat

 # start with first prime avoiding gcd of leading coefficients
 repeat t.prime := NextPrimeInt(t.prime); until a mod t.prime <> 0;

 # compute modular gcd with leading coefficient <a>
 t.gcd := RPGcdModPrime(f,g,t.prime,a,t.brci);
 Info(InfoPoly,3,"gcd mod ",t.prime," = ",t.gcd);

 # loop until we have success
 repeat
 if 0 = DegreeOfLaurentPolynomial(t.gcd) then
 Info(InfoPoly,3,"<f> and <g> are relative prime");
 return One(f);
 fi;
 until RPGcd1(t,a,f,g);
 until t.correct;

 # return the gcd
 return t.gcd;

end);

#############################################################################
##
#F GcdOp( <R>, <f>, <g> ) . . . . . . . for rational univariate polynomials
##
InstallRingAgnosticGcdMethod("rational univariate polynomials",
 IsCollsElmsElms,IsIdenticalObj,
 [IsRationalsPolynomialRing and IsEuclideanRing,
 IsUnivariatePolynomial,IsUnivariatePolynomial],0,
function(f,g)
local brci,gcd,fam,fc,gc;

 fam:=FamilyObj(f);

 if not (IsIdenticalObj(CoefficientsFamily(fam),CyclotomicsFamily)
 and ForAll(CoefficientsOfLaurentPolynomial(f)[1],IsRat)
 and ForAll(CoefficientsOfLaurentPolynomial(g)[1],IsRat)) then
 TryNextMethod(); # not applicable as not rational
 fi;

 brci:=CIUnivPols(f,g);
 if brci=fail then TryNextMethod();fi;
 # check trivial cases
 if -infinity = DegreeOfLaurentPolynomial(f) then
 return g;
 elif -infinity = DegreeOfLaurentPolynomial(g) then
 return f;
 elif 0 = DegreeOfLaurentPolynomial(f)
 or 0 = DegreeOfLaurentPolynomial(g)
 then
 return One(f);
 fi;

 # convert polynomials into integer polynomials
 f := PrimitivePolynomial(f)[1];
 g := PrimitivePolynomial(g)[1];
 Info(InfoPoly,3,"<f> = ",f);
 Info(InfoPoly,3,"<g> = ",g);

 fc:=CoefficientsOfLaurentPolynomial(f);
 gc:=CoefficientsOfLaurentPolynomial(g);

 # try heuristic method:
 gcd:=HeuGcdIntPolsCoeffs(fc[1],gc[1]);
 if gcd=fail then
 # fall back to the original version:
 gcd:=RPIGcd(f,g);
 return StandardAssociate(gcd);

 fi;
 fc:=Minimum(fc[2],gc[2]);
 fc:=fc+RemoveOuterCoeffs(gcd,fam!.zeroCoefficient);
 if Length(gcd)>0 and not IsOne(Last(gcd)) then
 gcd:=gcd/Last(gcd);
 fi;
 return LaurentPolynomialByExtRepNC(fam,gcd,fc,brci);
end);

InstallMethod(\mod,"reduction of univariate rational polynomial at a prime",
 true,[IsUnivariatePolynomial,IsInt],0,
function(f,p)
local c;
 c:=CoefficientsOfLaurentPolynomial(f);
 if Length(c[1])>0 and
 ForAny(c[1],i->not (IsRat(i) or IsAlgebraicElement(i))) then
 TryNextMethod();
 fi;
 return LaurentPolynomialByCoefficients(
 CoefficientsFamily(FamilyObj(f)),List(c[1],i->i mod p),c[2],
 IndeterminateNumberOfLaurentPolynomial(f));
end);

InstallMethod(\mod,"reduction of general rational polynomial at a prime",
 true,[IsPolynomial,IsInt],0,
function(f,p)
local c,d,i,m;
 c:=ExtRepPolynomialRatFun(f);
 d:=[];
 for i in [2,4..Length(c)] do
 if not (IsRat(c[i]) or IsAlgebraicElement(c[i])) then
 TryNextMethod();
 fi;
 m:=c[i] mod p;
 if m<>0 then
 Add(d,c[i-1]);
 Add(d,m);
 fi;
 od;
 return PolynomialByExtRepNC(FamilyObj(f),d);
end);

#############################################################################
##
#F RPQuotientModPrime(<f>,<g>,) . . . quotient
##
BindGlobal("RPQuotientModPrime",function(f,g,p)
local m, n, i, k, c, q, val, fc,gc,brci,fam;

 # get base ring
 brci:=CIUnivPols(f,g);
 fam:=FamilyObj(f);
 # reduce <f> and <g> mod 
 f := f mod p;
 g := g mod p;

 fc:=CoefficientsOfLaurentPolynomial(f);
 gc:=CoefficientsOfLaurentPolynomial(g);

 # if <f> is zero return it
 if 0 = Length(fc[1]) then
 return f;
 fi;

 # check the value of the valuation of <f> and <g>
 if fc[2] < gc[2] then
 return false;
 fi;
 val := fc[2]-gc[2];

 # Try to divide <f> by <g>,compute mod 
 q := [];
 n := Length(gc[1]);
 m := Length(fc[1]) - n;
 gc:=gc[1];
 f := ShallowCopy(fc[1]);
 for i in [0..m] do
 c := f[m-i+n]/gc[n] mod p;
 for k in [1..n] do
 f[m-i+k] := (f[m-i+k] - c*gc[k]) mod p;
 od;
 q[m-i+1] := c;
 od;

 # Did the division work?

 for i in [ 1 .. m+n ] do
 if f[i] <> fam!.zeroCoefficient then
 return false;
 fi;
 od;
 val:=val+RemoveOuterCoeffs(q,fam!.zeroCoefficient);
 return LaurentPolynomialByExtRepNC(fam,q,val,brci);

end);

#############################################################################
##
#F RPGcdRepresentationModPrime(<f>,<g>,) . gcd
##
BindGlobal("RPGcdRepresentationModPrime",function(f,g,p)

 local val, # the minimal valuation of <f> and <g>
 s, sx, # first line of gcd algorithm
 t, tx, # second line of gcd algorithm
 h, hx, # temp for swapping lines
 q, # quotient
 n,m,r,c, # used in quotient
 brci,
 i,k; # loops

 Info(InfoPoly,3,"f=",f,"g=",g);
 # get base ring
 brci:=CIUnivPols(f,g);
 brci:=[CoefficientsFamily(FamilyObj(f)),brci];

 # remove common x^i term
 f:=CoefficientsOfLaurentPolynomial(f);
 g:=CoefficientsOfLaurentPolynomial(g);

 val:=Minimum(f[2],g[2]);
 f :=ShiftedCoeffs(f[1],f[2]-val);
 g :=ShiftedCoeffs(g[1],g[2]-val);
 ReduceCoeffsMod(f,p); ShrinkRowVector(f);
 ReduceCoeffsMod(g,p); ShrinkRowVector(g);

 # compute the gcd and representation mod 
 s := ShallowCopy(f); sx := [ One(brci[1]) ];
 t := ShallowCopy(g); tx := [];
 while 0 < Length(t) do
 Info(InfoPoly,3,"<s> = ",s,", <sx> = ",sx,"\n",
 "#I <t> = ",t,", <tx> = ",tx);

 # compute the euclidean quotient of <s> by <t>
 q := [];
 n := Length(t);
 m := Length(s) - n;
 r := ShallowCopy(s);
 for i in [ 0 .. m ] do
 c := r[m-i+n] / t[n] mod p;
 for k in [ 1 .. n ] do
 r[m-i+k] := (r[m-i+k] - c * t[k]) mod p;
 od;
 q[m-i+1] := c;
 od;
 Info(InfoPoly,3,"<q> = ",q);

 # update representation
 h := t;
 hx := tx;
 t := s;
 AddCoeffs(t,ProductCoeffs(q,h),-1);
 ReduceCoeffsMod(t,p);
 ShrinkRowVector(t);
 tx := sx;
 AddCoeffs(tx,ProductCoeffs(q,hx),-1);
 ReduceCoeffsMod(tx,p);
 ShrinkRowVector(tx);
 s := h;
 sx := hx;
 od;
 Info(InfoPoly,3,"<s> = ",s,", <sx> = ",sx);

 # compute conversion for standard associate
 q := (1/Last(s)) mod p;

 # convert <s> and <x> back into polynomials
 if 0 = Length(g) then
 #sx := q * sx;
 MultVector(sx,q);
 ReduceCoeffsMod(sx,p);
 return [ LaurentPolynomialByCoefficients(brci[1],sx,0,brci[2]),
 Zero(brci[1]) ];
 else
 #hx := q * sx;
 hx:=ShallowCopy(sx);
 MultVector(hx,q);
 ReduceCoeffsMod(hx,p);
 hx := LaurentPolynomialByCoefficients(brci[1],hx,0,brci[2]);
 AddCoeffs(s,ProductCoeffs(sx,f),-1);
 #s := q * s;
 MultVector(s,q);
 ReduceCoeffsMod(s,p);
 s := LaurentPolynomialByCoefficients(brci[1],s,0,brci[2]);
 g := LaurentPolynomialByCoefficients(brci[1],g,0,brci[2]);
 q := RPQuotientModPrime(s,g,p);
 return [ hx,q ];
 fi;

end);

#############################################################################
##
#F HenselBound(<pol>,[<minpol>,<den>]) . . . Bounds for Factor coefficients
## if the computation takes place over an algebraic extension, then
## minpol and denominator must be given
##
InstallGlobalFunction(HenselBound,function(arg)
local pol,n,nalpha,d,dis,rb,bound,a,i,j,k,l,w,bin,lm,bea,polc,ro,rbpow;

 pol:=arg[1];
 if Length(arg)>1 then
 n:=arg[2];
 d:=arg[3];

 dis:=Discriminant(n);
 nalpha:=RootBound(n); # bound for norm of \alpha.

 polc:=CoefficientsOfLaurentPolynomial(pol)[1];
 if not ForAll(polc,IsRat) then
 # now try to bound the roots of f accordingly. As in all estimates by
 # RootBound only the absolute value of the coefficients is used, we will
 # estimate these first, and replace f by the polynomial
 # x^n-b_{n-1}x^(n-1)-...-b_0 whose roots are certainly larger
 a:=[];
 for i in polc do
 # bound for coefficients of pol
 if IsRat(i) then
 Add(a,AbsInt(i));
 else
 Add(a,Sum(ExtRepOfObj(i),AbsInt)*nalpha);
 fi;
 od;
 a:=-a;
 a[Length(a)]:=-a[Length(a)];
 pol:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,a,1);
 else
 pol:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,polc,1);
 fi;
 n:=DegreeOfLaurentPolynomial(n);
 else
 n:=1;
 fi;

 bound:=[];
 rb:=0;
 #BeauzamyBound
 bea:=BeauzamyBound(pol);
 # compute Landau-Mignotte bound for absolute values of
 # coefficients of any factor
 polc:=CoefficientsOfLaurentPolynomial(pol)[1];
 w:=Sum(polc,i->i^2);
 # we want an upper bound of the root, RootInt will give a lower
 # bound. So we compute the root of w-1 (in case w is a perfect square)
 # and add 1. As we nowhere selected a specific galois representative,
 # this bound (which is rational!) will bound all conjugactes as well.
 lm:=(RootInt(Int(w)-1,2)+1);
 for k in [1..DegreeOfLaurentPolynomial(pol)] do

 l:=2^k*lm;
 if l<bea then
 w:=l;
 else
 w:=bea;
 fi;

 #if lb>10^30 or n>1 then
 if bea>10^200 or n>1 then
 if rb=0 then
 rb:=RootBound(pol);
 if rb>1000 and not IsInt(rb) then
 rb:=Int(rb+1);
 fi;
 rbpow:=[rb];
 a:=rb;
 fi;
 # now try factor deg k
 bin:=1;
 for j in [1..k] do
 bin:=bin*(k-j+1)/j;
 if not IsBound(rbpow[j]) then
 rbpow[j]:=rbpow[j-1]*rb;
 if rbpow[j]>10 and not IsInt(rbpow[j]) then
 rbpow[j]:=Int(rbpow[j]+1);
 fi;
 fi;
 w:=bin*rbpow[j];
 if w>a then
 a:=w;
 fi;
 od;

 # select the better bound
 if a<l then
 w:=a;
 else
 w:=l;
 fi;
 fi;

 if n>1 then
 # algebraic Extension case
 # finally we have to bound (again) the coefficients of \alpha when
 # writing the coefficients of the factor as \sum c_i/d\alpha^i.

 ro:=AbsInt(dis);
 ro:=RootInt(NumeratorRat(ro))/(1+RootInt(DenominatorRat(ro)-1));
 w:=Int(d*w*Factorial(n)/ro*nalpha^(n*(n-1)/2))+1;
 fi;

 bound[k]:=Int(w)+1;
 od;
 return bound;
end);

#############################################################################
##
#F TrialQuotientRPF(<f>,<g>,) . . . . . . f/g if coeffbounds are given by b
##
InstallGlobalFunction(TrialQuotientRPF,function(f,g,b)
local fc,gc,a,m, n, i, k, c, q, val, brci,fam;
 brci:=CIUnivPols(f,g);
 a:=DegreeOfLaurentPolynomial(f)-DegreeOfLaurentPolynomial(g);
 fam:=FamilyObj(f);
 fc:=CoefficientsOfLaurentPolynomial(f);
 gc:=CoefficientsOfLaurentPolynomial(g);
 if a=0 then
 # Special case (that has to return 0)
 a:=b[1];
 else
 a:=b[a];
 fi;
 if 0=Length(fc[1]) then
 return f;
 fi;
 if fc[2]<gc[2] then
 return fail;
 fi;
 val:=fc[2]-gc[2];
 q:=[];
 n:=Length(gc[1]);
 m:=Length(fc[1])-n;
 f:=ShallowCopy(ShiftedCoeffs(fc[1],fc[2]));
 gc:=gc[1];
# if IsField(R) then
 for i in [0..m] do
 c:=f[(m-i+n)]/gc[n];
 if MaxNumeratorCoeffAlgElm(c)>a then
 Info(InfoPoly,3,"early break");
 return fail;
 fi;
 for k in [1..n] do
 f[m-i+k]:=f[m-i+k]-c*gc[k];
 od;
 q[m-i+1]:=c;
 od;
# else
# for i in [0..m] do
# c:=Quotient(R,f[m-i+n],gc[n]);
# if c=fail then
# return fail;
# fi;
# if MaxNumeratorCoeffAlgElm(c)>a then
# Info(InfoPoly,3,"early break\n");
# return fail;
# fi;
# for k in [1..n] do
# f[m-i+k]:=f[m-i+k]-c*g.coefficients[k];
# if MaxNumeratorCoeffAlgElm(f[m-i+k])>b then
# Info(InfoPoly,3,"early break\n");
# return fail;
# fi;
# od;
# q[m-i+1]:=c;
# od;
# fi;
 k:=Zero(f[1]);
 for i in [1..m+n] do
 if f[i]<>k then
 return fail;
 fi;
 od;
 val:=val+RemoveOuterCoeffs(q,fam!.zeroCoefficient);
 return LaurentPolynomialByExtRepNC(fam,q,val,brci);
end);

#############################################################################
##
#F TryCombinations(<f>,...) . . . . . . . . . . . . . . . . try factors
##
InstallGlobalFunction(TryCombinations,function(f,lc,l,p,t,bounds,opt,split,useonefacbound)
local p2, res, j, i,ii,o,d,b,lco,degs, step, cnew, sel, deli,
 degf, good, act, da, prd, cof, q, combi,mind,binoli,alldegs;

 alldegs:=t.degrees;
 # <res> contains the irr/reducible factors and the remaining ones
 res:=rec(irreducibles:=[],
 irrFactors :=[],
 reducibles :=[],
 redFactors :=[],
 remaining :=[ 1 .. Length(l) ]);

 # coefficients should be in -/2 and /2
 p2 :=p/2;
 deli:=List(l,DegreeOfLaurentPolynomial);

 # sel are the still selected indices
 sel:=[ 1 .. Length(l) ];

 # create List of binomial coefficients to speed up the 'Combinations' process
 binoli:=[];
 for i in [0..Length(l)-1] do
 binoli[i+1]:=List([0..i],j->Binomial(i,j));
 od;

 step:=0;
 act :=1;
 repeat

 # factors of larger than half remaining degree we will find as
 # final cofactor. This cannot be used if we factor only using the one
 # factor bound!

 if not useonefacbound then
 degf:=DegreeOfLaurentPolynomial(f);
 degs:=Filtered(alldegs,i -> 2*i<=degf);
 else
 degs:=alldegs;
 fi;

 if IsBound(opt.onlydegs) then
 degs:=Intersection(degs,opt.onlydegs);
 Info(InfoPoly,3,"degs=",degs);
 fi;

 if act in sel then

 # search all combinations of Length step+1 containing the act-th
 # factor,that are allowed
 good:=true;
 da:=List(degs,i -> i-deli[act]);

 # check,whether any combination will be of suitable degree

 cnew:=Set(deli{Filtered(sel,i->i>act)});

 if ForAny(da,i->NrRestrictedPartitions(i,cnew,step)>0) then
 # as we have all combinations including < <act>,we can skip them
 Info(InfoPoly,2,"trying length ",step+1," containing ",act);
 cnew:=Filtered(sel,i -> i > act);
 else
 Info(InfoPoly,2,"length ",step+1," containing ",act," not feasible");
 cnew:=[];
 fi;

 mind:=Sum(deli); # the maximum of the possible degrees. We surely
 # will find something smaller

 lco:=Binomial(Length(cnew),step);
 if 0 = lco then
 # fix mind to make sure,we don't erroneously eliminate the factor
 mind:=0;
 else
 Info(InfoPoly,2,lco," combinations");
 i:=1;
 while good and i<=lco do

 # try combination number i
 # combi:=CombinationNr(cnew,step,i);

 q:=i;
 d:=Length(cnew); # the remaining Length
 o:=0;
 combi:=[];
 for ii in [step-1,step-2..0] do
 j:=1;
 b:=binoli[d][ii+1];
 while q>b do
 q:=q-b;
 # compute b:=Binomial(d-(j+1),ii);
 b:=b*(d-j-ii)/(d-j);
 j:=j+1;
 od;
 o:=j+o;
 d:=d-j;
 Add(combi,cnew[o]);
 od;

 # check whether this yields a minimal degree
 d:=Sum(deli{combi});
 if d<mind then
 mind:=d;
 fi;

 if d in da then
 AddSet(combi,act); # add the 'always' factor

 # make sure that the quotient has a chance,compute the
 # extremal coefficient of the product:
 q:=(ProductMod(List(l{combi},
 i->CoefficientsOfLaurentPolynomial(i)[1][1]),
 p) * lc) mod p;
 if p2 < q then
 q:=q - p;
 fi;

 # As we don't know yet the gcd of all the products
 # coefficients (to make it primitive),we do a slightly
 # weaker test: (test of leading coeffs is first in
 # 'TrialQuotientRPF') this just should reduce the number of
 # 'ProductMod' necessary. the absolute part of the
 # product must divide the absolute part of f up to a
 # divisor of <lc>
 q:=CoefficientsOfLaurentPolynomial(f)[1][1] / q * lc;
 if not IsInt(q) then
 Info(InfoPoly,3,"ignoring combination ",combi);
 q:=fail;
 else
 Info(InfoPoly,2,"testing combination ",combi);

 # compute the product and reduce
 prd:=ProductMod(l{combi},p);
 prd:=CoefficientsOfUnivariatePolynomial(prd);
 cof:=[];
 for j in [ 1 .. Length(prd) ] do
 cof[j]:=(lc*prd[j]) mod p;
 if p2 < cof[j] then
 cof[j]:=cof[j] - p;
 fi;
 od;

 # make the product primitive
 cof:=cof * (1/Gcd(cof));
 prd:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,
 cof,t.ind);
 q:=TrialQuotientRPF(f,prd,bounds);
 fi;

 if q <> fail then
 f:=q;
 Info(InfoPoly,2,"found true factor of degree ",
 DegreeOfLaurentPolynomial(prd));
 if Length(combi)=1 or split then
 q:=0;
 else
 q:=2*lc*OneFactorBound(prd);
 if q <= p then
 Info(InfoPoly,2,"proven irreducible by 'OneFactorBound'");
 fi;
 fi;

 # for some reason,we know,the factor is irred.
 if q <= p then
 Append(res.irreducibles,combi);
 Add(res.irrFactors,prd);

 if IsBound(opt.stopdegs)
 and DegreeOfLaurentPolynomial(prd) in opt.stopdegs then
 Info(InfoPoly,2,"hit stopdegree");
 Add(res.redFactors,f);
 res.stop:=true;
 return res;
 fi;

 else
 Add(res.reducibles,combi);
 Add(res.redFactors,prd);
 fi;
 SubtractSet(res.remaining,combi);
 good:=false;
 SubtractSet(sel,combi);
 fi;

 fi;
 i:=i+1;

 od;
 fi;

 # we can forget about the actual factor,as any longer combination
 # is too big
 if Length(degs)>1 and deli[act]+mind >= Maximum(degs) then
 Info(InfoPoly,2,"factor ",act," can be further neglected");
 sel:=Difference(sel,[act]);
 fi;

 fi;

 # consider next factor
 act:=act + 1;
 if 0 < Length(sel) and act>Maximum(sel) then
 step:=step+1;
 act :=sel[1];
 fi;

 # until nothing is left
 until 0 = Length(sel) or Length(sel)<step;

 # if <split> is true we *must* find a complete factorization.
 if split and 0 < Length(res.remaining) and not IsOne(f) then
#and not(IsBound(opt.onlydegs) or IsBound(opt.stopdegs)) then

 # the remaining f must be an irreducible factor,larger than deg/2
 Append(res.irreducibles,res.remaining);
 res.remaining:=[];
 Add(res.irrFactors,f);
 fi;

 # return the result
 return res;

end);

#############################################################################
##
#F RPSquareHensel(<f>,<t>,<opt>)
##
BindGlobal("RPSquareHensel",function(f,t,opt)

 local p, # prime
 q, # current modulus
 q1, # last modulus
 l, # factorization mod <q>
 lc, # leading coefficient of <f>
 bounds, # Bounds for Factor Coefficients
 ofb, # OneFactorBound
 k, # Lift boundary
 prd, # product of <l>
 rep, # lifted representation of gcd(<lp>)
 fcn, # index of true factor in <l>
 dis, # distance of <f> and <l>
 cor, # correction
 max, # maximum absolute coefficient of <f>
 res, # result
 gcd, # used in gcd representation
 i, j; # loop

 # get <l> and 
 l:=t.factors;
 p:=t.prime;

 # get the leading coefficient of <f>
 lc:=LeadingCoefficient(f);

 # and maximal coefficient
 max:=Maximum(List(CoefficientsOfUnivariatePolynomial(f),AbsInt));

 # compute the factor coefficient bounds
 ofb:=2*AbsInt(lc)*OneFactorBound(f);
 Info(InfoPoly,2,"One factor bound = ",ofb);
 bounds:=2*AbsInt(lc)*HenselBound(f);
 k:=bounds[Maximum(Filtered(t.degrees,
 i-> 2*i<=DegreeOfLaurentPolynomial(f)))];
 Info(InfoPoly,2,"Hensel bound = ",k);

 # compute a representation of the 1 mod 
 Info(InfoPoly,2,"computing gcd representation: ",Runtime());
 prd:=(1/lc * f) mod p;
 gcd:=RPQuotientModPrime(prd,l[1],p);
 rep:=[ One(f) ];
 for i in [ 2 .. Length(l) ] do
 dis:=RPQuotientModPrime(prd,l[i],p);
 cor:=RPGcdRepresentationModPrime(gcd,dis,p);
 gcd:=(cor[1] * gcd + cor[2] * dis) mod p;
 rep:=List(rep,z -> z * cor[1] mod p);
 Add(rep,cor[2]);
 od;
 Info(InfoPoly,2,"representation computed: ",Runtime());

 # <res> will hold our result
 res:=rec(irrFactors:=[], redFactors:=[], remaining:=[],
 bounds:=bounds);

 # start Hensel until <q> is greater than k
 q :=p^2;
 q1 :=p;
 while q1 < k do
 Info(InfoPoly,2,"computing mod ",q);

 for i in [ 1 .. Length(l) ] do
 dis:=List(CoefficientsOfUnivariatePolynomial(f),i->i/lc mod q);
 ReduceCoeffsMod(dis,CoefficientsOfUnivariatePolynomial(l[i]),q);
 #dis:=EuclideanRemainder(PolynomialRing(Rationals),dis,l[i]) mod q;
 #dis:=CoefficientsOfUnivariatePolynomial(dis);
 dis:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,dis,t.ind);
 l[i]:=l[i] + BPolyProd(rep[i],dis,l[i],q);
 od;

 # NOCH: Assert und leading coeff
 #if not ForAll(CoefficientsOfLaurentPolynomial(Product(l)-f)[1],
 # i->i mod q=0) then
 # Error("not product modulo q");
 #fi;

 # if this is not the last step update <rep> and check for factors
 if q < k then

 # correct the inverses
 for i in [ 1 .. Length(l) ] do
 if Length(l)=1 then
 dis:=l[1]^0;
 else
 dis:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,
 [1],t.ind);
 for j in Difference([1..Length(l)],[i]) do
 dis:=BPolyProd(dis,l[j],l[i],q);
 od;
 fi;
 rep[i]:=BPolyProd(rep[i],
 (2-APolyProd(rep[i],dis,q)),
 l[i],
 q);
 od;

 # try to find true factors
 if max <= q then
 Info(InfoPoly,2,"searching for factors: ",Runtime());
 fcn:=TryCombinations(f,lc,l,q,t,bounds,opt,false,false);
 elif ofb < q then
 Info(InfoPoly,2,"searching for factors: ",Runtime());
 fcn:=TryCombinations(f,lc,l,q,t,bounds,opt,false,true);

 #Info(InfoPoly,2,"finishing search: ",Runtime());
 else
 fcn:=rec(irreducibles:=[], reducibles:=[]);
 fi;

 # if we have found a true factor update everything
 if 0 < Length(fcn.irreducibles)+Length(fcn.reducibles) then
 # append irreducible factors to <res>.irrFactors
 Append(res.irrFactors,fcn.irrFactors);

 # append reducible factors to <res>.redFactors
 Append(res.redFactors,fcn.redFactors);

 # compute new <f>
 prd:=Product(fcn.redFactors) * Product(fcn.irrFactors);

 f :=f/prd;
 if IsBound(fcn.stop) then
 res.stop:=true;
 return res;
 fi;

 lc :=LeadingCoefficient(f);
 ofb:=2*AbsInt(lc)*OneFactorBound(f);
 Info(InfoPoly,2,"new one factor bound = ",ofb);

 # degree arguments or OFB arguments prove f irreducible
 if (ForAll(t.degrees,i->i=0 or 2*i>=DegreeOfLaurentPolynomial(f))
 or ofb<q) and DegreeOfLaurentPolynomial(f)>0 then
 Add(fcn.irrFactors,f);
 Add(res.irrFactors,f);
 f:=f^0;
 fi;

 # if <f> is trivial return
 if DegreeOfLaurentPolynomial(f) < 1 then
 Info(InfoPoly,2,"found non-trivial factorization");
 return res;
 fi;

 # compute the factor coefficient bounds
 k:=HenselBound(f);
 bounds:=List([ 1 .. Length(k) ],
 i -> Minimum(bounds[i],k[i]));
 k:=2 * AbsInt(lc)
 * bounds[Maximum(Filtered(t.degrees,
 i-> 2*i<=DegreeOfLaurentPolynomial(f)))];
 Info(InfoPoly,2,"new Hensel bound = ",k);

 # remove true factors from <l> and corresponding <rep>
 prd:=(1/LeadingCoefficient(prd)) * prd mod q;
 l :=l{fcn.remaining};
 rep:=List(rep{fcn.remaining},x -> prd * x);

 # reduce <rep>[i] mod <l>[i]
 for i in [ 1 .. Length(l) ] do
 #rep[i]:=rep[i] mod l[i] mod q;
 rep[i]:=CoefficientsOfLaurentPolynomial(rep[i]);
 rep[i]:=ShiftedCoeffs(rep[i][1],rep[i][2]);
 j:=CoefficientsOfLaurentPolynomial(l[i]);
 j:=ReduceCoeffsMod(rep[i],ShiftedCoeffs(j[1],j[2]),q);
 # shrink the list rep[i], according to the 'j' value
 rep[i]:=rep[i]{[1..j]};
 rep[i]:=LaurentPolynomialByCoefficients(
 CyclotomicsFamily,rep[i],0,t.ind);
 od;

 # if there was a factor,we ought to have found it
 elif ofb < q then
 Add(res.irrFactors,f);
 Info(InfoPoly,2,"f irreducible,since one factor would ",
 "have been found now");
 return res;
 fi;
 fi;

 # square modulus
 q1:=q;
 q :=q^2;

 # avoid a modulus too big
 if q > k then
 q:=p^(LogInt(k,p)+1);
 fi;
 od;

 # return the remaining polynomials
 res.remPolynomial:=f;
 res.remaining :=l;
 res.primePower :=q1;
 res.lc :=lc;
 return res;

end);

#############################################################################
##
#F RPFactorsModPrime(<f>) find suitable prime and factor
##
## <f> must be squarefree. We test 3 "small" and 2 "big" primes.
##
BindGlobal("RPFactorsModPrime",function(f)
local i, # loops
 fc, # f's coeffs
 lc, # leading coefficient of <f>
 p, # current prime
 PR, # polynomial ring over F_
 fp, # <f> in <R>
 lp, # factors of <fp>
 min, # minimal number of factors so far
 P, # best prime so far
 LP, # factorization of <f> mod 
 deg, # possible degrees of factors
 t, # return record
 cof, # new coefficients
 tmp;

 fc:=CoefficientsOfLaurentPolynomial(f)[1];
 # set minimal number of factors to the degree of <f>
 min:=DegreeOfLaurentPolynomial(f)+1;
 lc :=LeadingCoefficient(f);

 # find a suitable prime
 t:=rec(ind:=IndeterminateNumberOfLaurentPolynomial(f));
 p:=1;
 for i in [ 1 .. 5 ] do

 # reset to big prime after first 3 test
 if i = 4 then p:=Maximum(p,1000); fi;

 # find a prime not dividing <lc> and <f>_ squarefree
 repeat
 repeat
 p :=NextPrimeInt(p);
 until lc mod p <> 0 and fc[1] mod p <> 0;
 PR:=PolynomialRing(GF(p));
 tmp:=1/lc mod p;
 fp :=UnivariatePolynomialByCoefficients(FamilyObj(Z(p)),
 List(fc,x->((tmp*x) mod p)* Z(p)^0),1);
 until 0 = DegreeOfLaurentPolynomial(Gcd(PR,fp,Derivative(fp)));

 # factorize <f> modulo 
 Info(InfoPoly,2,"starting factorization mod p: ",Runtime());
 lp:=Factors(PR,fp);

 Info(InfoPoly,2,"finishing factorization mod p: ",Runtime());

 # if <fp> is irreducible so is <f>
 if 1 = Length(lp) then
 Info(InfoPoly,2,"<f> mod ",p," is irreducible");
 t.isIrreducible:=true;
 return t;
 else
 Info(InfoPoly,2,"found ",Length(lp)," factors mod ",p);
 Info(InfoPoly,2,"of degree ",List(lp,DegreeOfLaurentPolynomial));
 fi;

 # choose a maximal prime with minimal number of factors
 if Length(lp) <= min then
 min:=Length(lp);
 P :=p;
 LP :=lp;
 fi;

 # compute the possible degrees
 tmp:=Set(Combinations(List(lp,DegreeOfLaurentPolynomial)),Sum);
 if 1 = i then
 deg:=tmp;
 else
 deg:=Intersection(deg,tmp);
 fi;

 # if there is only one possible degree != 0 then <f> is irreducible
 if 2 = Length(deg) then
 Info(InfoPoly,2,"<f> must be irreducible,only one degree left");
 t.isIrreducible:=true;
 return t;
 fi;

 od;

 # convert factors <LP> back to the integers
 lp:=ShallowCopy(LP);
 for i in [ 1 .. Length(LP) ] do
 cof:=CoefficientsOfUnivariatePolynomial(LP[i]);
 #cof:=IntVecFFE(cof);
 cof:=List(cof,IntFFE);
 LP[i]:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,cof,t.ind);
 od;

 # return the chosen prime
 Info(InfoPoly,2,"choosing prime ",P," with ",Length(LP)," factors");
 Info(InfoPoly,2,"possible degrees: ",deg);
 t.isIrreducible:=false;
 t.prime :=P;
 t.factors :=LP;
 t.degrees :=deg;
 return t;

end);

#############################################################################
##
#M FactorsSquarefree(<R>,<f>,<opt>) . factors of <f>
##
## <f> must be square free and must have a constant term.
##
InstallMethod( FactorsSquarefree, "univariate rational poly", true,
 [IsRationalsPolynomialRing,IsUnivariatePolynomial,IsRecord],0,
function(R,f,opt)
local t, h, fac, g, tmp;

 # find a suitable prime,if <f> is irreducible return
 t:=RPFactorsModPrime(f);
 if t.isIrreducible then return [f]; fi;
 Info(InfoPoly,2,"using prime ",t.prime," for factorization");

 # for easy combining,we want large degree factors first
 Sort(t.factors,
 function(a,b)
 return DegreeOfLaurentPolynomial(a) > DegreeOfLaurentPolynomial(b);
 end);

 # start Hensel
 h:=RPSquareHensel(f,t,opt);

 # combine remaining factors
 fac:=[];

 # first the factors found by hensel
 if 0 < Length(h.remaining) then
 Info(InfoPoly,2,"found ",Length(h.remaining)," remaining terms");
 tmp:=TryCombinations(
 h.remPolynomial,
 h.lc,
 h.remaining,
 h.primePower,
 t,
 h.bounds,
 opt,
 true,false);
 Append(fac,tmp.irrFactors);
 Append(fac,tmp.redFactors);
 else
 tmp:=rec();
 fi;

 # append the irreducible ones
 if 0 < Length(h.irrFactors) then
 Info(InfoPoly,2,"found ",Length(h.irrFactors)," irreducibles");
 Append(fac,h.irrFactors);
 fi;

 # and try to factorize the (possible) reducible ones
 if 0 < Length(h.redFactors) then
 Info(InfoPoly,2,"found ",Length(h.redFactors)," reducibles");

 if not (IsBound(tmp.stop) or IsBound(h.stop)) then
 # the stopping criterion has not yet been reached
 for g in h.redFactors do
 Append(fac,FactorsSquarefree(R,g,opt));
 od;
 else
 Append(fac,h.redFactors);
 fi;
 fi;

 # and return
 return fac;

end);

BindGlobal("RPIFactors",function(R,f,opt)
local fc,ind, v, g, q, s, r, x,shift;

 fc:=CoefficientsOfLaurentPolynomial(f);
 ind:=IndeterminateNumberOfLaurentPolynomial(f);
 # if <f> is trivial return
 Info(InfoPoly,2,"starting integer factorization: ",Runtime());
 if 0 = Length(fc[1]) then
 Info(InfoPoly,2,"<f> is trivial");
 return [ f ];
 fi;

 # remove a valuation
 v:=fc[2];
 f:=UnivariatePolynomialByCoefficients(CyclotomicsFamily,fc[1],ind);
 x:=LaurentPolynomialByCoefficients(CyclotomicsFamily,[1],1,ind);

 # if <f> is almost constant return
 if Length(fc[1])=1 then
 s:=List([1..fc[2]],i->x);
 s[1]:=s[1]*LeadingCoefficient(f);
 return s;
 fi;

 # if <f> is almost linear return
 if 1 = DegreeOfLaurentPolynomial(f) then
 Info(InfoPoly,2,"<f> is almost linear");
 s:=List([1..v],f -> x);
 Add(s,f);
 return s;
 fi;

 # shift the zeros of f if appropriate
 if DegreeOfLaurentPolynomial(f) > 20 then
 g:=MinimizedBombieriNorm(f);
 f:=g[1];
 shift:=-g[2];
 else
 shift:=0;
 fi;

 # make <f> integral,primitive and square free
 g:=Gcd(R,f,Derivative(f));
 q:=PrimitivePolynomial(f/g)[1];
 q:=q * SignInt(LeadingCoefficient(q));
 Info(InfoPoly,3,"factorizing polynomial of degree ",
 DegreeOfLaurentPolynomial(q));

 # and factorize <q>
 if DegreeOfLaurentPolynomial(q) < 2 then
 Info(InfoPoly,2,"<f> is a linear power");
 s:=[ q ];
 else
 # treat zeroes (only one possible)
 s:=CoefficientsOfLaurentPolynomial(q)[2];
 if s>0 then
 s:=[x];
 q:=q/x;
 else
 s:=[];
 fi;
 s:=Concatenation(s,FactorsSquarefree(R,q,opt));
 fi;

 # find factors of <g>
 for r in s do
 if 0 < DegreeOfLaurentPolynomial(g)
 and DegreeOfLaurentPolynomial(g) >= DegreeOfLaurentPolynomial(r) then
 q:=Quotient(R,g,r);
 while 0 < DegreeOfLaurentPolynomial(g) and q <> fail do
 Add(s,r);
 g:=q;
 if DegreeOfLaurentPolynomial(g)>=DegreeOfLaurentPolynomial(r) then
 q:=Quotient(R,g,r);
 else
 q:=fail;
 fi;
 od;
 fi;
 od;

 # reshift
 if shift<>0 then
 Info(InfoPoly,2,"shifting zeros back");
 Apply(s,i->Value(i,x+shift));
 fi;

 # sort the factors
 Append(s,List([1..v],f->x));
 Sort(s);

 if not (IsBound(opt.stopdegs) or IsBound(opt.onlydegs)) and Sum(s,DegreeOfLaurentPolynomial)<>DegreeOfLaurentPolynomial(f)+v then
 Error("degree discrepancy!");
 fi;

 # return the (primitive) factors
 return s;

end);

#############################################################################
##
#M Factors(<R>,<f> ) . . factors of rational polynomial
##
InstallMethod(Factors,"univariate rational polynomial",IsCollsElms,
 [IsRationalsPolynomialRing,IsUnivariatePolynomial],0,
function(R,f)
local r,cr,opt,irf,i;

 opt:=ValueOption("factoroptions");
 PushOptions(rec(factoroptions:=rec())); # options do not hold for
 # subsequent factorizations
 if opt=fail then
 opt:=rec();
 fi;

 cr:=CoefficientsRing(R);
 irf:=IrrFacsPol(f);
 i:=PositionProperty(irf,i->i[1]=cr);
 if i<>fail then
 # if we know the factors,return
 PopOptions();
 return ShallowCopy(irf[i][2]);
 fi;

 # handle trivial case
 if DegreeOfLaurentPolynomial(f) < 2 then
 StoreFactorsPol(cr,f,[f]);
 PopOptions();
 return [f];
 fi;

 # compute the integer factors
 r:=RPIFactors(R,PrimitivePolynomial(f)[1],opt);

 # convert into standard associates and sort
 r:=List(r,x -> StandardAssociate(R,x));
 Sort(r);

 if Length(r)>0 then
 # correct leading term
 r[1]:=r[1]*Quotient(R,f,Product(r));
 fi;

 # and return
 if not IsBound(opt.onlydegs) and not IsBound(opt.stopdegs) then
 StoreFactorsPol(cr,f,r);
 for i in r do
 StoreFactorsPol(cr,i,[i]);
 od;
 fi;
 PopOptions();
 return r;

end);

#############################################################################
##
#F SymAdic( <x>, ) . . . . . . . . . . symmetric -adic expansion of <x>
#F ( and <x> integers)
##
BindGlobal( "SymAdic", function(x,b)
 local l, b2, r;
 b2:=QuoInt(b,2);
 l:=[];
 while x<>0 do
 r:=x mod b;
 if r>b2 then
 r:=r-b;
 fi;
 Add(l,r);
 x:=(x-r)/b;
 od;
 return l;
end );

#############################################################################
##
#F HeuGcdIntPolsExtRep(<fam>,<a>,)
##
## takes two polynomials ext reps, primitivizes them with integer
## coefficients and tries to
## compute a Gcd expanding Gcds of specializations.
## Source: Geddes,Czapor,Labahn: Algorithm 7.4

#V MAXTRYGCDHEU: defines maximal number of attempts to find Gcd heuristically
MAXTRYGCDHEU:=6;

InstallGlobalFunction(HeuGcdIntPolsExtRep,function(fam,d,e)
local x,xd,er,i,j,k,m,p,sel,xi,gamma,G,g,loop,zero,add;

 # find the indeterminates and their degrees:
 x:=[];
 xd:=[[],[]];
 er:=[d,e];
 for k in [1,2] do
 for i in [1,3..Length(er[k])-1] do
 m:=er[k][i];
 for j in [1,3..Length(m)-1] do
 p:=Position(x,m[j]);
 if p=fail then
 Add(x,m[j]);
 p:=Length(x);
 xd[1][p]:=0;
 xd[2][p]:=0;
 fi;
 xd[k][p]:=Maximum(xd[k][p],m[j+1]);
 od;
 od;
 od;

 # discard the indeterminates which occur in only one of the polynomials
 sel:=[];
 for i in [1..Length(x)] do
 if xd[1][i]>0 and xd[2][i]>0 then
 Add(sel,i);
 fi;
 od;
 x:=x{sel};
 xd:=List(xd,i->i{sel});

 # are the variables disjoint or do we have no variables at all?
 if Length(x)=0 then
 # return the gcd of all coefficients involved
 G:=Gcd(Concatenation(d{[2,4..Length(d)]},e{[2,4..Length(e)]}));
 return G;
 fi;

 # pick the first indeterminate
 x:=x[1];
 xd:=[xd[1][1],xd[2][1]];

 xi:=2*Minimum(Maximum(List(d{[2,4..Length(d)]},AbsInt)),
 Maximum(List(e{[2,4..Length(e)]},AbsInt)))+2;

 for loop in [1..MAXTRYGCDHEU] do
 if LogInt(AbsInt(Int(xi)),10)*Maximum(xd)>5000 then
 return fail;
 fi;
 # specialize both polynomials at x=xi
 # and compute their heuristic Gcd
 gamma:=HeuGcdIntPolsExtRep(fam,
 SpecializedExtRepPol(fam,d,x,xi),
 SpecializedExtRepPol(fam,e,x,xi) );

 if gamma<>fail then
 # generate G from xi-adic expansion
 G:=[];
 i:=0;
 if IsInt(gamma) then
 zero:=Zero(gamma);
 else
 zero:=[];
 fi;
 while gamma<>zero do
 if IsInt(gamma) then
 # gamma is an integer value
 g:=gamma mod xi;
 if g>xi/2 then
 g:=g-xi; # symmetric rep.
 fi;
 gamma:=(gamma-g)/xi;
 if i=0 then
 add:=[[],g];
 else
 add:=[[x,i],g];
 fi;
 else
 # gamma is an ext rep
 g:=[];
 for j in [2,4..Length(gamma)] do
 k:=gamma[j] mod xi;
 if k>xi/2 then
 k:=k - xi; #symmetric rep
 fi;
 if k<>0*k then
 Add(g,gamma[j-1]);
 Add(g,k);
 fi;
 od;
 #gamma:=(gamma-g)/xi in ext rep;
 add:=ShallowCopy(g);
 add{[2,4..Length(add)]}:=-add{[2,4..Length(add)]}; #-g
 gamma:=ZippedSum(gamma,add,0,fam!.zippedSum); # gamma-g
 gamma{[2,4..Length(gamma)]}:=gamma{[2,4..Length(gamma)]}/xi; # /xi
 #add:=g*xp^i; in extrep
 add:=ZippedProduct(g,[[x,i],1],0,fam!.zippedProduct);
 fi;
 # G:=G+add in extrep
 G:=ZippedSum(G,add,0,fam!.zippedSum);
 i:=i+1;
 od;

 if Length(G)>0 and Length(G[1])>0 and
 QuotientPolynomialsExtRep(fam,d,G)<>fail and
 QuotientPolynomialsExtRep(fam,e,G)<>fail then
 return G;
 fi;
 fi;
 xi:=QuoInt(xi*73794,27011); #square of golden ratio
 od;
 return fail;
end);

# and the same in univariate
InstallGlobalFunction(HeuGcdIntPolsCoeffs,function(f,g)
local xi, t, h, i, lf, lg, lh,fr,gr;

 if IsEmpty(f) or IsEmpty(g) then
 return fail;
 fi;
 # first test value for heuristic gcd:
 xi:=2+2*Minimum(Maximum(List(f,AbsInt)),Maximum(List(g,AbsInt)));
 i:=0;
 lf:=Last(f);
 lg:=Last(g);

 # and now the tests:
 while i< MAXTRYGCDHEU do

 # xi-adic expansion of Gcd(f(xi),g(xi)) (regarded as coefficient list)
 h:=Gcd( ValuePol(f,xi),
 ValuePol(g,xi) );
 h:=SymAdic(h,xi);

 # make it primitive:
 t:=Gcd(h);
 if t<>1 then
 h:=h/t;
 fi;
 lh:=Last(h);

 # check if it divides f and g: if yes, ready! if no, try larger xi

 ## this should be done more efficiently !
 if RemInt(lg,lh)=0 and RemInt(lf,lh)=0 then
 fr:=ShallowCopy(f);
 ReduceCoeffs(fr,h);
 gr:=ShallowCopy(g);
 ReduceCoeffs(gr,h);
 t:=Set(Concatenation(fr,gr));
 else
 t:=false;
 fi;
 if t=[0] then
 Info(InfoPoly,4,"GcdHeuPrimitiveList: tried ",i+1," values");
 return h;
 else
 i:=i+1;
 xi:=QuoInt(xi*73794,27011);
 fi;
 od;
 Info(InfoPoly,2,"GcdHeuPrimitiveList: failed after trying ",
 MAXTRYGCDHEU," values");
 return fail;
end);

InstallMethod(HeuristicCancelPolynomialsExtRep,"rationals",true,
 [IsRationalFunctionsFamily,IsList,IsList],0,
function(fam,a,b)
local nf,df,g;
 if not (HasCoefficientsFamily(fam)
 and IsIdenticalObj(CoefficientsFamily(fam),CyclotomicsFamily)
 and ForAll(a{[2,4..Length(a)]},IsRat)
 and ForAll(b{[2,4..Length(b)]},IsRat)) then
 # the coefficients are not all rationals
 TryNextMethod();
 fi;

 # make numerator and denominator primitive with integer coefficients
 nf:=PrimitiveFacExtRepRatPol(a);
 if not IsOne(nf) then
 a:=ShallowCopy(a);
 a{[2,4..Length(a)]}:=a{[2,4..Length(a)]}/nf;
 fi;

 df:=PrimitiveFacExtRepRatPol(b);
 if not IsOne(df) then
 b:=ShallowCopy(b);
 b{[2,4..Length(b)]}:=b{[2,4..Length(b)]}/df;
 fi;

 # the remaining common factor
 nf:=nf/df;
 g:=HeuGcdIntPolsExtRep(fam,a,b);
 if IsList(g) and (Length(g)>2 or (Length(g)=2 and Length(g[1])>0)) then
 # make g primitive
 df:=PrimitiveFacExtRepRatPol(g);
 if not IsOne(df) then
 g:=ShallowCopy(g);
 g{[2,4..Length(g)]}:=g{[2,4..Length(g)]}/df;
 fi;

 Info(InfoPoly,3,"Heuristic integer gcd returns ",g);

 a:=QuotientPolynomialsExtRep(fam,a,g);
 b:=QuotientPolynomialsExtRep(fam,b,g);
 if a<>fail and b<>fail then
 a{[2,4..Length(a)]}:=a{[2,4..Length(a)]}*nf;
 g:=[a,b];
 return g;
 fi;
 fi;
 return fail;
end);

InstallGlobalFunction(PolynomialModP,function(pol,p)
local f, cof;
 f:=GF(p);
 cof:=CoefficientsOfUnivariatePolynomial(pol);
 return UnivariatePolynomial(f,cof*One(f),
 IndeterminateNumberOfUnivariateRationalFunction(pol));
end);

Messung V0.5 in Prozent

¤ Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.0.116Bemerkung: (vorverarbeitet am 2026-10-12) ¤

Wurzel

Suchen

PVS Prover

Isabelle Prover

NIST Cobol Testsuite

Cephes Mathematical Library

Vienna Development Method

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.