Lecture 2: October  13, 2022

IntegerGCD[f_,g_]:=
Module[{r1,r2},
{r1,r2}={f,g};
LBasic={r1,r2};
While[r2=!=0,
Q=Quotient[r1,r2];
{r2,r1}={r1-Q*r2,r2};
LBasic=Append[LBasic,r2];
];
r1
];

IntegerGCD[-123332,-3232]

LBasic

BasicPolyGCD[f_,g_,var_]:=
Module[{r1,r2},
{r1,r2}={f,g};
LBasic={r1,r2};
While[r2=!=0,
Q=PolynomialQuotient[r1,r2,var];
{r2,r1}={Collect[r1-Q*r2,var,Together],r2};
LBasic=Append[LBasic,r2];
];
r1
];

f=(1+x+2 x^2+x^3)(2x+1)//Expand

g=(1+x+2 x^2)(2x+1)//Expand

BasicPolyGCD[f,g,x]

LBasic

f=Sum[RandomInteger[1234]x^i,{i,0,10}](3x+42)//Expand;

g=Sum[RandomInteger[1234]x^i,{i,0,10}](3x+42)//Expand;

f

res=BasicPolyGCD[f,g,x]//Timing

res/Coefficient[res[[2]],x]//Expand

LBasic

myNormalize[f_,var_]:=
Module[{lc},
If[f===0,Return[0]];
lc=Coefficient[f,var,Exponent[f,var]];
Collect[f/lc,var,Together]
]

NormalizedPolyGCD[f_,g_,var_]:=
Module[{r1,r2},
{r1,r2}={myNormalize[f,var],myNormalize[g,var]};
LMonic={r1,r2};
While[r2=!=0,
Q=PolynomialQuotient[r1,r2,var];
{r2,r1}={Collect[r1-Q*r2,var,Together],r2};
r2=myNormalize[r2,var];
LMonic=Append[LMonic,r2];
];
r1
];

NormalizedPolyGCD[f,g,x]

LMonic

RatLength[f_]:=If[IntegerQ[f],IntegerLength[f],IntegerLength[Numerator[f]]+IntegerLength[Denominator[f]]]

PolySize[f_]:=Max[Map[RatLength,Cases[f,_Rational|_Integer,Infinity]]]

Map[PolySize,LBasic]

Map[PolySize,LMonic]

f=Sum[RandomInteger[1234]x^i,{i,0,100}](3x+42)//Expand;

g=Sum[RandomInteger[1234]x^i,{i,0,100}](3x+42)//Expand;

BasicPolyGCD[f,g,x]//Timing

NormalizedPolyGCD[f,g,x]//Timing

Map[PolySize,LBasic]

Map[PolySize,LMonic]

Lecture 3: October  20, 2022

Inverses in F_p

LU[f_]:=
Module[{lu},
lu=Sign[f];
If[lu===0,lu=1];
lu
];

EEA[f_,g_]:=
Module[{r1,s1,t1,r2,s2,t2,rho,Q},
rho=LU[f];
{r1,s1,t1}={f/rho,1/rho,0};
rho=LU[g];
{r2,s2,t2}={g/rho,0,1/rho};

While[r2=!=0,
Q=Quotient[r1,r2];
{r2,r1}={r1-Q*r2,r2};
rho=LU[r2];
r2=r2/rho;
{s2,s1}={(s1-Q*s2)/rho,s2};
{t2,t1}={(t1-Q*t2)/rho,t2};
];
{r1,s1,t1}
];

{a,p}={11,19};

{r,s,t}=EEA[a,p]

s*a+t*p==r

Mod[s*11,p]

Inverses in F[x]/f

LC[f_,var_]:=
Module[{lc},
lc=Coefficient[f,var,Exponent[f,var]];
If[lc===0,lc=1];
lc
];

myNormalize[h_,rho_,var_]:=Collect[h/rho,var,Together]

EEAPoly[f_,g_,var_]:=
Module[{r1,s1,t1,r2,s2,t2,rho,Q},
rho=LC[f,var];
{r1,s1,t1}={myNormalize[f,rho,var],1/rho,0};
rho=LC[g,var];
{r2,s2,t2}={myNormalize[g,rho,var],0,1/rho};


While[r2=!=0,
Q=PolynomialQuotient[r1,r2,var];
{r2,r1}={Collect[r1-Q*r2,var,Together],r2};
rho=LC[r2,x];
r2=myNormalize[r2,rho,var];
{s2,s1}={myNormalize[s1-Q*s2,rho,var],s2};
{t2,t1}={myNormalize[t1-Q*t2,rho,var],t2};
];
{r1,s1,t1}
];

{a,f}={1+2x+3x^ 2+x^3,x^3+x+1};

IrreduciblePolynomialQ[f]

{r,s,t}=EEAPoly[a,f,x]

s*a+t*f-r//Expand

PolynomialRemainder[s*a,f,x]

Inverses in F_{19^3} (clumsy/inefficient)

{a,f}={1+2x+3x^ 2+x^3,x^3+x+1};

{r,s,t}=EEAPoly[a,f,x]

Inv17=EEA[17,19][[2]]

Mod[Inv17*17,19]

s

sIn19=PolynomialMod[s*17*Inv17,19]

a

PolynomialRemainder[a*sIn19,f,x,Modulus->19]

Lecture 4 : October  27, 2022

Partial Fraction Decomposition (PFD)

den=(x+2)^3(x^2-x-1)^2//Expand

num=x^6+72 x^4+38x^2+123

f=num/den

Apart[f,x]

f1=(1491/(25 (((2+x))^3)))+(334/(25 (((2+x))^2)))+(2933/(125 (2+x)))//Factor

f2=-((2 (-561+227 x))/(25 (((-1-x+(x^2)))^2)))-((24 (-271+117 x))/(125 (-1-x+(x^2))))//Factor

f1+f2-f//Together

CRA (LA Example)

MyRandom[]:=RandomInteger[{-1000000000,1000000000}]

size=300;

A=Table[Table[MyRandom[],{i,1,size}],{i,1,size}];

b=A.Table[100000000*i,{i,1,size}];

LinearSolve[A,b,Method->OneStepRowReduction]//AbsoluteTiming

p=Prime[40000000000]

LinearSolve[A,b,Modulus->p]//AbsoluteTiming

pL={Prime[400000],Prime[400001]}

Apply[Times,pL]

solL=ParallelTable[LinearSolve[A,b,Modulus->pL[[i]]],{i,1,Length[pL]}]//AbsoluteTiming

solL=solL[[2]]//Transpose;

Map[ChineseRemainder[#,pL]&,solL]//AbsoluteTiming

Lecture 5: November  3, 2022

Pade approximation

LC[f_,var_]:=
Module[{lc},
lc=Coefficient[f,var,Exponent[f,var]];
If[lc===0,lc=1];
lc
];

myNormalize[h_,rho_,var_]:=Collect[h/rho,var,Together]

RatRecon[f_,g_,var_,k_]:=
Module[{r1,s1,t1,r2,s2,t2,rho,Q},
rho=LC[f,var];
{r1,t1}={myNormalize[f,rho,var],0};
rho=LC[g,var];
{r2,t2}={myNormalize[g,rho,var],1/rho};

While[Exponent[r2,x]>=k,
Q=PolynomialQuotient[r1,r2,var];
{r2,r1}={Collect[r1-Q*r2,var,Together],r2};
rho=LC[r2,x];
r2=myNormalize[r2,rho,var];
{t2,t1}={myNormalize[t1-Q*t2,rho,var],t2};
];
rho=LC[t2,x];
{myNormalize[r2,rho,var],myNormalize[t2,rho,var]}
];

rat=x/(x^2-x-1)/(x-1)/(1+x^2)

{h,f}={x^7,Series[rat,{x,0,6}]//Normal}

{r,t}=RatRecon[h,f,x,2]

ratReconstruct=r/t//Factor

ratReconstruct-rat

Lecture 6: November  10, 2022

Rational number reconstruction (not complete version)

LU[f_]:=
Module[{lu},
lu=Sign[f];
If[lu===0,lu=1];
lu
];

RatNumber[f_,g_,k_]:=
Module[{r1,s1,t1,r2,s2,t2,rho,Q},
rho=LU[f];
{r1,t1}={f/rho,0};
rho=LU[g];
{r2,t2}={g/rho,1/rho};
While[r2>=k,
Q=Quotient[r1,r2];
{r2,r1}={r1-Q*r2,r2};
rho=LU[r2];
r2=r2/rho;
{t2,t1}={(t1-Q*t2)/rho,t2};
];
rho=LU[t2];
{r2,t2}/rho
];

{A,b}={{{2,2,3},{1,3,1},{1,1,-9}},{1,1,1}}

sol=LinearSolve[A,b]

m=659

k=20;
{k,m/k}//N

solHom=LinearSolve[A,b,Modulus->m]

solHomRat=Map[RatNumber[m,#,k]&,solHom]

solHomRatRecon=Map[#[[1]]/#[[2]]&,solHomRat]

sol-solHomRatRecon

Lecture 9: December  1, 2022

Exit[]

LC[f_,var_]:=
Module[{lc},
lc=Coefficient[f,var,Exponent[f,var]];
If[lc===0,lc=1];
lc
];

myNormalize[h_,rho_,var_]:=Collect[h,var,Together[#/rho]&]

GetInverse[f_,g_,var_]:=
Module[{r1,s1,r2,s2,rho,Q},
rho=LC[f,var];
{r1,s1}={myNormalize[f,rho,var],1/rho};
rho=LC[g,var];
{r2,s2}={myNormalize[g,rho,var],0};
While[r2=!=0,
Q=PolynomialQuotient[r1,r2,var];
{r2,r1}={Collect[r1-Q*r2,var,Together],r2};
rho=LC[r2,var];
r2=myNormalize[r2,rho,var];
{s2,s1}={myNormalize[s1-Q*s2,rho,var],s2};
];
s1
];

s=GetInverse[x^2+1,x^2+x+1,x]

PolynomialRemainder[s*(x^2+1),x^2+x+1,x]

s=GetInverse[1+4 y+7 (y^2)+6 (y^3)+7 (y^4)+8 (y^5)+10 (y^6),1+4 y+7 (y^2)+6 (y^3)+7 (y^4)+8 (y^5)+10 (y^6)+(y^7),y]

PolynomialRemainder[s*(1+4 y+7 (y^2)+6 (y^3)+7 (y^4)+8 (y^5)+10 (y^6)),1+4 y+7 (y^2)+6 (y^3)+7 (y^4)+8 (y^5)+10 (y^6)+(y^7),y]

ReducePoly[fIn_,rho_,p_,x_,y_]:=
Collect[fIn,x,Together[PolynomialRemainder[#*rho,p,y]]&];

ReducePoly[f_,p_,x_,y_]:=ReducePoly[f,1,p,x,y]

EEAField[f_,g_,p_,x_,y_]:=
Module[{r1,r2,rho,Q},
rho=GetInverse[LC[f,x],p,y];
r1=ReducePoly[f,rho,p,x,y];
rho=GetInverse[LC[g,x],p,y];
r2=ReducePoly[g,rho,p,x,y];
While[r2=!=0,
(*since r2 is monic, we can use the naive version...*)
{r2,r1}={ReducePoly[PolynomialRemainder[r1,r2,x],p,x,y],r2};
rho=GetInverse[LC[r2,x],p,y];
r2=ReducePoly[r2,rho,p,x,y];
];
r1
];

PP[f_,x_,y_]:=
Module[{coeffL,gcd},
coeffL=CoefficientList[f,x];
gcd=Apply[PolynomialGCD,CoefficientList[f,x]];
gcd=gcd/LC[gcd,y];
coeffL=Together[coeffL/gcd];
coeffL.Table[x^i,{i,0,Length[coeffL]-1}]
]

RandomIrrPoly[y_,d_]:=
Module[{f},
f=y^d;
While[!IrreduciblePolynomialQ[f],
f=Sum[RandomInteger[{1,10}]y^i,{i,0,d-1}]+y^d
];
f
]

PolyGCDBigPrime[f_,g_,x_,y_]:=
Module[{b,done,d,yDegBound,p,fTilde,gTilde,w,v,fStar,gStar,i},
b=PolynomialGCD[LC[f,x],LC[g,x]];

d=Max[Exponent[f,y],Exponent[g,y]];
yDegBound=d+1+Exponent[b,y];

done=False;
While[done===False,
p=RandomIrrPoly[y,yDegBound];
{fTilde,gTilde}={ReducePoly[f,p,x,y],ReducePoly[g,p,x,y]};
v=EEAField[fTilde,gTilde,p,x,y];
w=ReducePoly[b*v,p,x,y];
fStar=ReducePoly[PolynomialQuotient[fTilde,v,x],p,x,y];
gStar=ReducePoly[PolynomialQuotient[gTilde,v,x],p,x,y];
done=Exponent[ReducePoly[fStar*w,p,x,y],y]==Exponent[ReducePoly[b*fTilde,p,x,y],y]&&Exponent[ReducePoly[gStar*w,p,x,y],y]==Exponent[ReducePoly[b*gTilde,p,x,y],y];
If[done===False,Print[Redo]]
];
PP[w,x,y]
]

gcd=(1+x+y^3*x^2+(y^7+y^3+12345)x^4)//Expand;
f=(x+x y^3+1)gcd//Expand;
g=(x^2+y^3+x y +1)gcd//Expand;

res=PolyGCDBigPrime[f,g,x,y]//Timing

res[[2]]-gcd//Together

Do[PolyGCDBigPrime[f,g,x,y],{i,0,100}]

Lecture 10: December  15, 2022

GCD

Exit[]

LC[f_,var_]:=
Module[{lc},
lc=Coefficient[f,var,Exponent[f,var]];
If[lc===0,lc=1];
lc
];

NormalPolyGCD[f_,g_,x_]:=
Module[{res},
res=PolynomialGCD[f,g];
res=Collect[res/LC[res,x],x,Together];
res
]

LaunchKernels[5];

DistributeDefinitions[{LC,NormalPolyGCD}];

RandomPoints[no_,b_,var_]:=
Module[{pointL,point},
pointL={};
While[Length[pointL]=!=no,
point=RandomInteger[{-1000,1000}];
If[Together[b/.var->point]=!=0,
pointL=Union[pointL,{point}]
];
];
pointL
]

PP[f_,x_,y_]:=
Module[{coeffL,gcd},
coeffL=CoefficientList[f,x];
gcd=Apply[PolynomialGCD,CoefficientList[f,x]];
gcd=gcd/LC[gcd,y];
coeffL=Together[coeffL/gcd];
coeffL.Table[x^i,{i,0,Length[coeffL]-1}]
]

PolyGCDSmallPrime[f_,g_,x_,y_]:=
Module[{b,done,d,yDegBound,pL,eL,e,gcdL,fTilde,gTilde,w,v,fStar,gStar,i,gL,gcdLL,bL,fL},
Print["Compute GCD of LCs"];
b=PolynomialGCD[LC[f,x],LC[g,x]];

d=Max[Exponent[f,y],Exponent[g,y]];
yDegBound=d+1+Exponent[b,y];

yDegBound=Floor[yDegBound*2/3];

Print["Number of points: ",yDegBound];

done=False;
While[done===False,
pL=RandomPoints[yDegBound,b,y];

Print["Compute GCDs at evaluation points"];
gcdL=ParallelMap[NormalPolyGCD[f/.y->#,g/.y->#,x]&,pL];

eL=Map[Exponent[#,x]&,gcdL];
e=Max[eL];
eL={eL,pL,gcdL}//Transpose;
eL=Select[eL,#[[1]]===e&];
If[Length[eL]<yDegBound,Continue[]];
{eL,pL,gcdL}=Take[eL,yDegBound]//Transpose;

Print["Interpolation to get the possible GCD"];
gcdLL=Map[CoefficientList[#,x]&,gcdL]//Transpose;
bL=Table[b/.y->pL[[i]],{i,Length[pL]}]//Together;
w=ParallelTable[InterpolatingPolynomial[{pL,gcdLL[[i]]*bL}//Transpose,y],{i,Length[gcdLL]}];
w=w.Table[x^i,{i,0,e}];

Print[PolynomialQ[f*b/w//Together]&&PolynomialQ[g*b/w//Together]];



Print["Compute the potential cofactor for the first polynomial (up to a content)"];
fL=ParallelTable[PolynomialQuotient[(f/.y->pL[[i]]),gcdL[[i]],x],{i,Length[pL]}];
fL=Map[CoefficientList[#,x]&,fL]//Transpose;



fStar=ParallelTable[InterpolatingPolynomial[{pL,fL[[i]]}//Transpose,y],{i,Length[fL]}];
fStar=fStar.Table[x^i,{i,0,Length[fL]-1}];

Print["Compute the potential cofactor for the second polynomial (up to a content)"];
gL=ParallelTable[PolynomialQuotient[(g/.y->pL[[i]]),gcdL[[i]],x],{i,Length[pL]}];
gL=Map[CoefficientList[#,x]&,gL]//Transpose;
gStar=ParallelTable[InterpolatingPolynomial[{pL,gL[[i]]}//Transpose,y],{i,Length[gL]}];
gStar=gStar.Table[x^i,{i,0,Length[gL]-1}];


done=Exponent[fStar,y]+Exponent[w,y]==Exponent[b,y]+Exponent[f,y]&&Exponent[gStar,y]+Exponent[w,y]==Exponent[b,y]+Exponent[g,y];

If[done===False,Print[Redo]];
];
PP[w,x,y]
]

gcd=(1+x+y^3*x^2+(y^7+y^3+12345)x^4)//Expand;
f=(x+x y^3+1)gcd//Expand;
g=(x^2+y^3+x y +1)gcd//Expand;

res=PolyGCDSmallPrime[f,g,x,y]//Timing

res=PolyGCDSmallPrime[f,g,x,y]//Timing

MyRandom[]:=RandomInteger[{-100000000,1000000000}]

gcd=Sum[MyRandom[]x^i*y^i,{i,0,100},{j,0,100}]//Expand;
f=Sum[MyRandom[]x^i*y^i,{i,0,150},{j,0,200}]gcd//Expand;
g=Sum[MyRandom[]x^i*y^i,{i,0,150},{j,0,200}]gcd//Expand;

res=PolyGCDSmallPrime[f,g,x,y]//AbsoluteTiming

PolyGCDSmallPrimeFaster[f_,g_,x_,y_]:=
Module[{b,done,d,yDegBound,pL,eL,gcdL,fTilde,gTilde,w,v,fStar,gStar,i},
Print["Compute GCD of LCs"];
b=PolynomialGCD[LC[f,x],LC[g,x]];

d=Max[Exponent[f,y],Exponent[g,y]];

yDegBound=d+1+Exponent[b,y];


yDegBound=Floor[yDegBound 2/3];

Print["Number of points: ",yDegBound];

done=False;
While[done===False,
pL=RandomPoints[yDegBound,b,y];

Print["Compute GCDs at evaluation points"];
gcdL=ParallelMap[NormalPolyGCD[f/.y->#,g/.y->#,x]&,pL];

eL=Map[Exponent[#,x]&,gcdL];
e=Max[eL];
eL={eL,pL,gcdL}//Transpose;
eL=Select[eL,#[[1]]===e&];
If[Length[eL]<yDegBound,Continue[]];
{eL,pL,gcdL}=Take[eL,yDegBound]//Transpose;

Print["Interpolation to get the possible GCD"];
gcdLL=Map[CoefficientList[#,x]&,gcdL]//Transpose;
bL=Table[b/.y->pL[[i]],{i,Length[pL]}]//Together;
w=ParallelTable[InterpolatingPolynomial[{pL,gcdLL[[i]]*bL}//Transpose,y],{i,Length[gcdLL]}];


w=w.Table[x^i,{i,0,e}];

done=PolynomialQ[f*b/w//Together]&&PolynomialQ[g*b/w//Together];

If[done===False,Print[Redo]];
];
PP[w,x,y]
]

res=PolyGCDSmallPrimeFaster[f,g,x,y]//AbsoluteTiming

Groebner basis

f=3x^7 y^4+2x^4 y^10-3x^3y^12

monL=MonomialList[f]

monL==MonomialList[f,"Lexicographic"]

W={{1,0},{0,1}};
W//MatrixForm

monL==MonomialList[f,W]

monL=MonomialList[f,"DegreeLexicographic"]

W={{1,1},{1,0}};

monL==MonomialList[f,W]

monL=MonomialList[f,"DegreeReverseLexicographic"]

W={{1,1},{0,-1}};
W//MatrixForm

monL==MonomialList[f,W]

Lecture 11: January 12 ,2023

Monomial order

f=3x^7 y^4+2x^4 y^10-3x^3y^12

monL=MonomialList[f]

monL==MonomialList[f,"Lexicographic"]

monL=MonomialList[f,"DegreeLexicographic"]

Polynomial reduction

f=3 (x^4) (y^2)+6 (x^7) (y^4)+2 (x^4) (y^10)-3 (x^3) (y^12)

S={y+(x^3) (y^3),(x^4) (y^6),x (y^3)};

res=PolynomialReduce[f,S,{x,y}]

f-(res[[1,1]]S[[1]]+res[[1,2]]S[[2]]+res[[1,3]]S[[3]])//Together

red1=f-S[[1]]*3 (x^4) y//Expand

PolynomialReduce[red1,S,{x,y}]

Lecture 14: January 26, 2023

Applications of Groebner basis

RandomPolynomial[varL_,degree_,size_,no_]:=
Module[{},
Sum[RandomInteger[{-size,size}]Apply[Times,varL^Table[RandomInteger[degree],{Length[varL]}]],{no}]
]

fL=Table[RandomPolynomial[{x,y,z},4,3,3],{3}]
G=GroebnerBasis[fL,{x,y}]

Ideal member ship

Example

fL={(x^3) (z^2)+2 (x^2) (z^3)-(y^3) (z^3),-(y^2)+(z^2)-x y (z^2)-(x^3) (y^3) (z^3),-2 x (y^2) z-(y^2) (z^2)+2 (x^3) y (z^3)+(x^2) (y^3) (z^3),(x^2) y z+(x^2) (y^3) z+2 (z^2)-2 (x^3) (z^2),-2 (x^2) (y^3)+2 (x^2) (y^2) (z^2)-2 (x^3) y (z^3)+2 x (y^3) (z^3)};

f=-2 (x^2) (y^2)+x (y^3) z-(y^5) z+2 (x^2) (z^2)+(x^2) y (z^2)-2 (x^3) y (z^2)-(x^3) (y^2) (z^2)+2 (x^6) (y^2) (z^2)+(x^2) (y^3) (z^2)-(x^3) (y^4) (z^2)+2 (y^5) (z^2)-3 x y (z^3)+4 (y^3) (z^3)-2 x (y^3) (z^3)+(x^4) y (z^4)-2 (x^2) (y^2) (z^4)-(x^5) (y^2) (z^4)-2 (y^3) (z^4)-(x^4) (y^3) (z^4)-4 y (z^5)+(x^2) y (z^5)-2 (x^3) y (z^5)-2 (x^2) (y^2) (z^5)-(x^3) (y^2) (z^5)-(y^3) (z^5)-(y^3) (z^6);

PolynomialReduce[f,fL,{x,y,z}]

G=GroebnerBasis[fL,{x,y,z}]

PolynomialReduce[f,G,{x,y,z}]

Elimination

polys={ -5x^2+y z-x-1, 2x+3 x y+y^2,x-3y+x z-2 z^2};

G=GroebnerBasis[polys,{z,x,y}]
Select[G,FreeQ[#,z]&]

Finding roots

Example 1

fL={-6 (x^5) y-8 x (y^2)-(x^4) (y^3)+x (y^4)+(x^5) (y^5),6 (y^3)+7 (x^5) (y^3)+10 (y^4)+10 (x^5) (y^4)+8 (x^3) (y^5),-7 (x^4)+9 (x^5)+10 y-3 (x^3) (y^2)-8 (x^5) (y^2)}

G=GroebnerBasis[fL,{y,z,x}]

solX=Solve[G[[1]]==0,x]//Union//Flatten

Solve[G[[2]]==0/.solX[[1]],y]

Solve[G[[2]]==0/.solX[[2]],y]

solL={{0,0},{7/9,0}}

fL/.{x->0,y->0}

fL/.{x->7/9,y->0}

Example 2

fL={-1-2 x y,2-x-2 x (y^2)+2 (x^2) (y^3)};

G=GroebnerBasis[fL,{x,y}]

G=GroebnerBasis[fL,{y,x}]

solX=Solve[G[[1]]==0,x]//Flatten//Union

Solve[G[[2]]==0/.solX[[1]],y]

Solve[G[[2]]==0/.solX[[2]],y]

solL={{1/2,-1},{3/2,-1/3}}

Table[fL/.{x->solL[[i,1]],y->solL[[i,2]]},{i,1,2}]