#Author := Ayoola Jinadu # If you find a bug in the code please email me at ajinadu@sfu.ca with(PolynomialTools): with(LinearAlgebra): with(ArrayTools): with(CodeTools[Profiling]):with(NumberTheory): read det; read Bmbot; read nextprime; MUL := proc(A,x,p) local n,m,Q1,Q2,P; n := nops(A); if n=1 then return x-A[1] mod p; else m := floor(n/2); Q1 := MUL(A[1..m],x,p); Q2 := MUL(A[m+1..n],x,p); P := Expand(Q1*Q2) mod p; fi; end: PolyInterp := proc(A,S,Q,Alpha,X,p) local degs,degA,Pol,j,B,Bx,pol1,H,L,k,i,r1; degs := nops(S); degA := nops(A); Pol := Array(1..degs); B := Array(1..degA); for j to degs do for k to degA do # Bx := (Q*~S[j] mod p)+~Alpha mod p; Bx := Q*~S[j] mod p; Bx := Bx +~Alpha mod p; Bx := [A[k], seq(i,i=Bx) ]; B[k] := BB(Bx,p); od; pol1[j] := modp1(Interp(A,convert(B,list),X[1]),p); H := modp1(Sqrfree(pol1[j]), p); L := Match(H[2],Multi,p); # "Match" checks the degrees of the univariate monic square-free factors. # Their multiplicities are also verified. #If their degrees or multiplicities are not the same, it returns fail. Pol[j] := [seq( modp1(ConvertOut(r1, X[1]), p), r1=L ) ]; od; return Pol; end: RatFun := proc(pol2,W,S,p) local Q,m,i,pol1,f,pol; m := nops(W); f := modp1(Interp(S,W,s),p); pol := modp1(ConvertOut(f, s), p); Q := EEA(pol2,pol,s,p); if Q = FAIL then return FAIL; fi; #This check is totally unnecessary. #Subroutine Image will check the degree later. return Q; end: BMStep := proc(A,Rx,X,p,u,ta) local t,Lam,deg_Lam,Fbar,gx; gx := Same(A); if gx = true then return A[1] ; fi; t := numelems(A); # no of points Lam := BMEA(A,p,s); deg_Lam := degree(Lam,s); if 2*deg_Lam >= t then return FAIL fi; Fbar := InvBot(A,Rx,X,p,u,Lam,ta); return Fbar; end: BPol := proc(A::list,p::prime)local Y,C,i,n,B; global Rt; Y,n := op( eval(Rt,1) ); C := EVALMOD(Y,n,A, degs, p ); return C; end: EVAL := proc (f,p,Y,d,Ww,Alpha) local RG, S,j,m,A,ij,W,x,g,Wx; global degs, Rt; Rt := [f,Y]; Wx := map(convert,Ww,list); degs := [seq(degree(f,x),x in Y)]; RG := rand(1..p-1): do S := [seq( RG(), j=1..d+1)]; until nops({op(S)}) = d+1; m := numelems(Wx); for j from 1 to m do for ij from 1 to d+1 do g := (Wx[j]*~S[ij]) mod p +~Alpha mod p; A[ij] := BPol(g,p); od; W[j] := modp1(Interp(S,[seq(A[ij], ij=1..d+1) ],s),p); W[j] := modp1(ConvertOut(W[j], s), p); od; [seq(W[j], j=1..m) ]; end: Mshift := proc(Fbar,Coeffn,Kw,Rx,Alpha,X,Nd,u,p,ta) local Xx,f2,tj,f,Aux,Coef,ij,Ax,j; if whattype(Fbar)=integer then return Fbar; fi; Xx := X[2..nops(X)]; f2 := Fbar; tj := numelems(Coeffn); f := f2; Aux := 0; Coef := Array(1..tj); for ij from Nd-1 to 0 by -1 do if f2 <> 0 then Ax := EVAL(f2,p,Xx,ij+1,Kw,Alpha); Aux := Aux+~Ax mod p; fi; if ij > 0 then for j to tj do Coef[j] := coeff(Coeffn[j],s,ij)-coeff(Aux[j],s,ij) mod p; od; if IsZero(Coef) = true then f2 := 0; else f2 := BMStep(Coef,Rx,X,p,u,ta); if f2=FAIL then return FAIL; fi; fi; else f2 := coeff(Coeffn[1],s,0)-coeff(Aux[1],s,0) mod p; fi; f := Expand(f+f2) mod p; od; return f; end: Image := proc(CoeffN,FFk,Alpha,Rx,Q,u,ta,p) local Fk,Nd,iN,Ax,fk,Zk; Fk := FFk; (Nd,iN) := SameDeg(CoeffN,s); #Checks the degrees of the auxiliaray rational functions. if Nd = false then printf("Bad evaluation point at position %d\n",iN); return FAIL; fi; Ax := map( proc(B,d) coeff(B,s,d) end proc, CoeffN,Nd); if Fk = FAIL then Fk := BMStep(Ax,Rx,X,p,u,ta); fi; if Fk <> FAIL then fk := Mshift(Fk,CoeffN,Q,Rx,Alpha,X,Nd,u,p,ta); if fk <> FAIL then return fk; else Zk := FAIL; fi; else return FAIL; fi; end: DixonRes := proc(BB,X::list) local t,Dgx,degA,DegArr,degS,Deg,Rr,k,i,Rx,p,a,S,pol2,n,Alpha,A,tryBM,u,ta,ta_s,W,K,Ky,ik,ij,Zk,kl,Fk,Gk,dT,Q,H,kt,Im,j,HH,CoeffN,CoeffD,FF,GG,LCM,F,Qx,ls,WX,U,Ans,jj,Degg,f,Abar,gk,Bd,Gd; global AN,FacN,LT,degSY,DegG,Good,Bad,DegN,ANT,Multi; Dgx,degSY,degA,DegArr := TotDeg(BB,X); degS := max (degSY); FacN := nops(DegArr); #Number of monic square free factors. Degg := [seq( [seq( DegArr[k][i][2],i=1..numelems(DegArr[k])-1)],k=1..FacN ) ]; jj := 1; ANT := 1; #Get any monic square free factor of the form x1^k. for i to FacN do if Degg[i] <> [] then DegG[jj] := Degg[i] ; jj++; else FacN := FacN-1; f := X[1]^(Dgx[i]); ANT := ANT*f; #LTT := ANT; fi; od; LT := Array(1..FacN); DegG := [ seq( DegG[i], i=1..numelems(DegG) ) ]; Rr := BoundPar(BB,X,DegG,degA); #Max Partial degree bounds Rx := [seq(Rr[i]+1,i=1..nops(Rr)-1 )]; p := 4601552919265804289; printf("Number of Monic Square Factors = %d\n",nops(degSY) ); a := rand(1..p-1); do S := [ seq( a(), i=1..degS) ]; until nops({op(S)}) = degS; pol2 := MUL(S,s,p); n := nops(X)-1; do Alpha := [ seq( a(), i=1..n) ]; until nops({op(Alpha)}) = n; do A := [ seq( a(), i=1..degA) ]; until nops({op(A)}) = degA; tryBM := 3; i := 0; u := a(); ta := PrimitiveRoot(p, greaterthan = 30000); ta_s := ta &^u mod p; W := Array(1..n); K := Array(1..n); W[1] := ta_s; K[1] := ta; for i from 2 to n do W[i] := W[i-1] &^Rx[i-1] mod p; K[i] := K[i-1]&^Rx[i-1] mod p; od; Ky := K; ik := 1; ij := 1; Zk := FAIL; kl := 1; Fk := FAIL; Gk := FAIL; #Initialize Fk(Num)= FAIL and Gk(Denom)=FAIL i := 1; while ik <= nops(DegG) do dT := DegG[ik][ij]; if Zk=FAIL then #compute more evaluation points if i=1 then Q[i] := W; else Q[i] := (W*~Ky mod p); Ky := K*~Ky mod p; fi; H[i] := PolyInterp(A,S,Q[i],Alpha,X,p); #compute an image fi; if i=2*tryBM then for kt from kl to i do Im := [ seq( coeff(H[kt][j][ik],X[1],dT),j=1..degS ) ]; HH[kt] := RatFun(pol2,Im,S,p); #Get the rational functions CoeffN[kt] := numer( HH[kt] ); CoeffD[kt] := denom ( HH[kt] ); od; #Routine Image will check the degrees of these rational functions if Fk = FAIL then Fk := Image(CoeffN,Fk,Alpha,Rx,Q,u,ta,p); #Try to get Fk=Numerator #if Fk= Fail then get more images. if Fk = FAIL then kl := i+1; Zk := FAIL; tryBM := tryBM +iquo(tryBM,3); fi; fi; if Fk <> FAIL and Gk=FAIL then Gk := Image(CoeffD,Gk,Alpha,Rx,Q,u,ta,p); #Try to get Gk=Denom if Gk = FAIL then kl := i+1; Zk := FAIL; tryBM := tryBM +iquo(tryBM,3); else FF[ik][ij] := X[1]^dT*Fk; GG[ik][ij] := Gk; #Fk/Gk is found, get the next rational coeff for the current monic square free factor kl := 1; ij++; #ij is counter for the next rational function coeff. Fk := FAIL; Gk := FAIL; Zk := 1; #Zk=1 because more images might not be needed to get the next coeff fi; fi; fi; if Zk= FAIL then i++; fi; if ij > nops(DegG[ik]) then ij := 1; ik++; fi; #ik is counter to get the next square free factor od; printf("deg x1 =%d,emax = %d,#probes for the first prime = %d\n", degA,degS,CNT); LCM := Array(1..nops(DegG) ); Abar := Array(1..FacN); for i from 1 to nops(DegG) do LCM[i] := Lcm( seq(GG[i][ij], ij=1..nops(DegG[i]) )) mod p; if LCM[i] <> 1 then F := 0; Qx := 0; for ij to nops(DegG[i]) do ls := Divide(LCM[i],GG[i][ij],'qx') mod p; if ls = FAIL then return print(ls=FAIL); FAIL; fi; Qx := qx*FF[i][ij]; F := F+Qx; od; else F := add( FF[i][ij], ij=1..nops(DegG[i]) ); fi; WX := LCM[i]*X[1]^(degree(H[1][1][i]))+ F; LT[i] := Expand( WX ) mod p; Abar[i] := X[1]^(degree(H[1][1][i])) + add( FF[i][ij]/GG[i][ij], ij=1..nops(DegG[i]) ); od; # it is possible that iratrecon fails on 1 or more of the monic square free factor, seperate the failed ones from the successful one. # We only perform sparse interpolation on the failed ones. i := 0; j := 0; ij := 1; for Ans in LT do gk := iratrecon(Ans,p,maxquo); if gk = FAIL then i := i+1; Bd[i] := Ans; AN[i] := Abar[ij]; DegN[i] := DegG[ij]; else j := j+1; Gd[j] := gk; fi; ij := ij+1; od; if j >= 1 then Good := [ seq( Gd[ik], ik=1..j ) ]; else Good := []; fi; if i >= 1 then Bad := [ seq( Bd[ik], ik=1..i ) ]; printf("%d of the Monic square free factors needs more prime\n",i); fi; if j=FacN then U := mul(ik , ik in Good ); U := ANT*U; return U; fi; #ANT is any factor of the form x1^k or it could be one as initialized. FacN := i; AN := [ seq( AN[ik], ik=1..i ) ]; DegN := [ seq( DegN[ik],ik=1..i) ]; SparseKron(BB,degA,degS,p,X); #Sparse Interpolation using the support. end: