read "CMBBSHLToeplitz.mpl":
read "KY.mpl":

# This program computes the factors of det(Tn), where Tn is a symmetric Toeplitz matrix represented by a black box. 
# The black box BB is a Maple procedure which calls an external C program to compute the determinant of a symmetric 
# Toeplitz matrix in Zp using Bareiss' O(n^2) algorithm. 

# Copyright: Tian Chen and Michael Monagan, 2024

  
# ----- Initial Settings -----
# Select to use Maple or C code for each of the four subroutines: 
# 1. Probe the black box (Maple option not available, need to use the C code),
# 2. Evaluate the factors f_{i,j-1} at (x2=beta2^k,..., x_{j-1}=beta_{j-1}^k), 
# 3. Bivariate Hensel lifting (BHL),
# 4. Vandermonde solves (Select Maple for the fast Vandermonde solver).  
MapleCode := [C,C,C,C]: 
C := 0: Maple := 1:

BBInput := BB: # The black box which uses the C code to compute Det(Tn).  
Var_Perm := 0: # variable permutation not available for Toeplitz det compt.

N := 10: # Number of variables, also the size of Tn
p := prevprime(2^62-1): # A large prime 
Ntilde := 100003: # Ntilde < p
# ----- End of Initial Settings -----


randzp := rand(p):  
randzp1 := rand(1..Ntilde): 
X := [seq(x[i],i=1..N)];
alpha := Array([seq(randzp1(),i=1..N)],datatype=integer[8]);

# Pre-comting the degrees
tdegbd := N; # Total degree bound for symmetric Toeplitz matrices 
CNT := 0: tBBdet := 0: 
st := time():
for j to N do
    deg[j] := degB(BBInput,N,j,p,tdegbd); 
od:
et := time()-st: printf("time for degree computation = %f\n", et):
degA := [seq(deg[i],i=1..N)]; 
printf("no. of probes for degree compt = %d\n", CNT):

# Interpolate a(x1,alpha) and factor it
a1 := Interp1varZ_chrem( BBInput, X, alpha, degA, N, 1, p ); 
if lcoeff(a1,X[1]) <> 1 then printf("a1 should be monic in %a\n",X[1]); stop; fi: 
a1fac := factor(a1): 
a1f := unpack(a1fac): #printf("a1f = %a\n", a1f);
r := nops(a1f); 
for i to r do ei := a1f[i][2]; 
    if ei <> 1 then printf("a should be square-free\n"); stop; fi: 
od: 
f1 := convert(a1fac,list); # list of factors of a1 
f1 := map(modp, f1, p): # list of factors of a1, in positive mod range 

# Start Algorithm CMBBSHL 
# trace(CMSHL_BB); trace(CMSHL_BBstepj); 
CNT := 0: 
st := time():
fN := CMSHL_BB( f1, N, X, alpha, degA, N, p, Var_Perm, MapleCode ):
et := time()-st: printf("Total time for CMBBSHL = %f\n", et):
printf("Total no.of probes for CMBBSHL = %d\n", CNT):

fN := map( mods, fN, p ): # list of recovered factors in N variables 
nfN := map(nops, fN): # number of terms in each factor
printf("no.of terms in factors = %a;\n", nfN);

filename := cat("factors_Toeplitz_n",N); 
subs(FILENAME=filename, proc() save nfN,fN, FILENAME; end)();

# st := time():
# dd := LinearAlgebra:-Determinant(AA,method=minor):
# et := time()-st: printf("Det compt time in Maple = %f\n", et);
# divide(dd,op(1,fN)); divide(dd,op(2,fN));
# print(nops(fN));

# st := time(): 
# factor(dd):
# et := time()-st: printf("factorizatin time in Maple = %f\n", et);