#Author := Ayoola Jinadu #email := ajinadu@sfu.ca # If you find a bug in the code please email me. Points := proc(n,Tmax,p,u) local j,Q,f,g,W,i; f := rand(1..p-1); do Q := [seq(f(),j=1..n)]; until nops({op(Q)}) = n; W[1] := [seq(Q[i]&^u mod p, i=1..n ) ]; for j from 2 to Tmax do W[j] := W[j-1]*~Q mod p; od; g := [ seq(W[j],j=1..Tmax) ]; return Q,g; end: BEA := proc(a,b,x,p) local f,n,d,i; f := EEA(a,b,x,p); if f=FAIL then return FAIL; fi; n := numer(f); d := denom(f); i := 1/lcoeff(d,x) mod p; n := i*n mod p; d := i*d mod p; return n/d; end: EvalW := proc(BB,P,A,S,n,p,X,Alpha) local alpha1,Aux,k,r,T1,Q,L,Qk,Pk,Sk,alpha,Num,Dem,Ans,j,degA,deg,i,f,r1,m,H; degA := nops(A); m := nops(P)+1; T1 := [ seq(0,i=1..m) ]; for k from 1 to n do for i from 2 to m do T1[i] := Alpha[i-1]*S[k]+P[i-1] mod p; od; for r to degA do T1[1] := A[r]; Aux[r] := BB(T1,p); od: Q[k] := modp1(Interp(A,[ seq(Aux[i], i=1..degA) ],X[1]),p); H := modp1(Sqrfree(Q[k]), p); L := Match(H[2],Multi,p); Q[k] := [seq( modp1(ConvertOut(r1, X[1]), p), r1=L ) ]; od; [ seq(Q[k], k=1..n ) ]; end: EvalU := proc(Q,P,A,S,n,p,X,Deg,TreeB,t) local Qk,Pk,Sk,alpha,Num,Dem,Ans,j,degA,deg,i,f,r1,m; j := 1; for alpha in Deg do Qk := [seq(coeff(Q[i],X[1],alpha), i=1..n)]; Pk := Interp(S,Qk,s) mod p; Ans[j] := BEA(TreeB,Pk,s,p); if Ans[j]=FAIL then return FAIL; fi; j := j+1; od; [ seq(Ans[j], j=1..t ) ]; end: Support := proc(M,Y) local Ans,Degg,i,j,SupN,SupD,t,qr,WX,W,degg,SuppN,SuppD,Nums,Dems,suppN,suppD,ik,B,G,jh,ijk; global AN,DegN; for ik to FacN do #numelems(M) do jh := 1; for ijk in DegN[ik] do qr := coeff(AN[ik],Y[1],ijk); if qr <> 0 then WX[jh] := qr; jh := jh+1; fi; od; W := [ seq(WX[i], i=1..numelems(WX) ) ]; t := nops(W); Nums := Array(1..t); Dems := Array(1..t); SuppN := Array(1..t); suppN := Array(1..t); SuppD := Array(1..t); suppD := Array(1..t); for i to t do coeffs(expand(numer(W[i])),Y,'SupN'): SuppN[i] := [SupN]; Nums[i] := nops(SuppN[i]); coeffs(expand(denom(W[i])),Y,'SupD'): SuppD[i] := [SupD]; Dems[i] := nops(SuppD[i]); od; G[ik] :=[Nums,Dems,SuppN,SuppD,DegN[ik]]; od; return [seq(G[ik], ik=1..FacN ) ]; end: MKron := proc(BB,degA,degW,Y,p,Cr,Index) local W,SupN,N,Nmax,m,A,X,Ans,i,Aux,IPoints,t,Images,beta,suppN,suppD,Q,k,V,VT,fk,S,TreeB,LCM,KN,KD,Kn,Kd,Alpha,j,J,A1,Nd,Dd,iN,iD,sN,sD,A2,ij,f,h,H,q,qr,WX,degg,u,Qr,degs,x,IMQ,Nums,Dems,SuppN,SuppD,Degg,Nmax1,ki,Ax,kk,degS,Kdd; global AN; m := nops(Y); Qr := rand(1..p-1); A := [seq(Qr(),i=1..degA)]; m := m-1; X := Y[2..m+1]; Aux := Array(1..degA); Nmax := max( seq( max([Cr[kk][1],Cr[kk][2]]), kk in Index)); do S := [seq(Qr(),i=1..degW)] until nops({op(S)}) = degW ; u := Qr(); do Alpha := [seq(Qr(), i=1..m)] until nops({op(Alpha)}) = m; J,IPoints := Points(m,Nmax,p,u): IMQ := [ seq( EvalW(BB,IPoints[i],A,S,degW,p,X,Alpha), i=1..Nmax ) ]; ki := 1; for kk in Index do Nums,Dems,SuppN,SuppD,Degg := op(Cr[kk]); Ans := Y[1]^degree(AN[kk],Y[1]); j := 1; t := nops(DegN[kk]); A1 := Array(1..t); A2 := Array(1..t); Nmax1 := max(Nums,Dems); degS := degSY[kk]; S := [seq(S[i],i=1..degS)]; TreeB := MUL(S,s,p); Q := [ seq( [ seq( IMQ[ki][j][kk], j=1..degS) ],ki=1..nops(IMQ) ) ]; Images := [ seq( EvalU(Q[i],IPoints[i],A,S,degS,p,X,Degg,TreeB,t), i=1..Nmax1 ) ]; for i to t do A1[i] := Array( [seq(numer(Images[j][i]),j=1..Nmax1 )] ); A2[i] := Array( [seq(denom(Images[j][i]),j=1..Nmax1 )] ); Nd,iN := SameDeg( A1[i],s); Dd,iD := SameDeg( A2[i],s); if Nd = false then printf("Bad evaluation point at position %d\n",iN); return FAIL; elif Dd = false then printf("Bad evaluation point at position %d\n",iD); return FAIL; fi; od; LCM := 1; Kd := Array(1..t); Kdd := Array(1..t); Kn := Array(1..t); for i to t do sN := 1; sD := 1; suppN[i] := Array(1..Nums[i]); suppD[i] := Array(1..Dems[i]); for j to Nums[i] do f := SuppN[i][j]; degs := [seq(degree(f,x),x in X)]; suppN[i][j] := EVALMOD( f, X, J, degs, p ); od; V := [seq(coeff(numer(Images[j][i]),s,0),j=1..Nums[i])]; KN := VSolve(V,suppN[i],p); sN := [seq( suppN[i][ij]&^u mod p, ij=1..numelems(suppN[i]) ) ]; sN := [seq( 1/sN[ij] mod p, ij=1..numelems(sN) ) ]; KN := sN*~KN mod p; Kn[i] := add( (KN[j]*SuppN[i][j]),j=1..nops(sN) ); for j to Dems[i] do f := SuppD[i][j]; degs := [seq(degree(f,x),x in X)]; suppD[i][j] := EVALMOD( f, X, J, degs, p ); od; VT := [seq(coeff(denom(Images[j][i]),s,0),j=1..Dems[i])]; KD := VSolve(VT,suppD[i],p); sD := [seq( suppD[i][ij]&^u mod p, ij=1..numelems(suppD[i]) ) ]; sD := [seq( 1/sD[ij] mod p, ij=1..numelems(sD) ) ]; KD := sD*~KD mod p; Kd[i] := add( (KD[j]*SuppD[i][j]),j=1..nops(sD) ); LCM := Lcm(Kd[i],LCM) mod p; od; for i to t do Divide(LCM,Kd[i],'q') mod p; Kdd[i] := Expand(Kn[i]*q*Y[1]^Degg[i]) mod p; od; f := Expand(LCM*Y[1]^degree(AN[kk],Y[1])) mod p; Ax[ki] := add ( Kdd[i],i=1..t)+f; # print(nops( Ax[kk] )); ki++; od; [ seq( Ax[k], k=1..nops(Ax) ) ]; end: SparseKron := proc(BB,degA,degS,pp,X) local p1,LCM,degZ,p,Q,t,K,T,G,KT,B,h,f,U,CNT,j,H,k,q,Degg,Nums,Dems,SuppN,SuppD,Cr,i,P,ik,jk,Index,L,Ix; global LTX; LTX := Bad; KT := 0; p := pp; P := pp; Cr := Support(AN,X); #print(Cr); jk := 0; Index := [ seq(i,i=1..nops(AN) ) ]; G := 1; while true do p := nextprime(p); printf("New prime = %d \n",p); Q := MKron(BB,degA,degS,X,p,Cr,Index); q := P*p; for i in Index do L := chrem( [Q[i],LTX[i]],[p,P] ); H := iratrecon(L,q,maxquo); if H = FAIL then printf("IRA FAIL\n"); jk++; Ix[jk] := i; # Create a new index I. # Next time around the loop, we only do sparse interpolation for the square free factor that iratrecon failed on. T[jk] := L; fi; if H <> FAIL then G := G*H; fi; od; if jk = 0 then if Good <> [] then U := mul(ik, ik in Good); else U := 1; fi; # U is the product of the all the square free factors from Dixonres that succeded. G := U*ANT*G; return G; else LTX := [seq( T[j], j=1..jk ) ]; P := p*P; #FacN := jk; AN := [ seq( AG[i], i=1..jk ) ]; Index := [ seq( Ix[j],j=1..jk) ]; fi; jk := 0; od; return; end: