Quelle polyrat.gi Sprache: unbekannt

#############################################################################
##
## 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);

Quelle polyrat.gi Sprache: unbekannt

[ Dauer der Verarbeitung: 0.54 Sekunden (vorverarbeitet) ]