# Copyright: Tian Chen and Michael Monagan, 2024

VarPermGen := proc( N::nonnegint, opt::nonnegint, MainVar::posint ) 
    # Input: N is the no.of variables, 
    #        opt=0: no permutation, opt=1: reverse order, opt=2: random combination,
    #        opt>2: swap the chosen MainVar with the first variable.
    #        MainVar: for opt>2 only (1<=MainVar<=N).    
    # Output: VPL = [indicator(0 or 1), a list of new variable ordering].
    local VP,VPL,i;
    if opt = 0 or N = 0 or N = 1 then return [0,[seq(i,i=1..N)]]; fi; 
    if opt = 1 then VPL := [1,[seq(N-i+1,i=1..N)]];
    elif opt = 2 then VPL := [1,combinat:-randperm(N)];
    else   
        if MainVar = 1 then return [0,[seq(i,i=1..N)]]; 
        elif MainVar > N then 
            printf("Chosen Main Variable out of range\n"); return FAIL; 
        fi; 
        VP := Array(1..N,i->i);
        VP[1] := MainVar; VP[MainVar] := 1;
	VPL := [1,convert(VP,list)];
    fi;
    VPL;
end:

MakeBBMPoly := proc( a, VarPerm::{list} ) local X,N,Xnew,i; 
    # This makes a modular black box from a given polynomial a.
    # Input: a in Z[x1,...,xN],
    #        VarPerm is a list of variable permutation.
    # Output: A modular black box Bm(alpha,p) for polynomial a s.t.
    #         Bm(alpha,p) = a(alpha) mod p. 
    X := convert(indets(a),list);
    N := nops(X);
    Xnew := Array(1..N);
    for i to N do Xnew[i] := X[VarPerm[i]]; od; 
    proc( alpha::{Array,list}, p::prime ) global CNT,tBBeval,tBBdet;
        local i; CNT++;
	Eval(a,[seq(Xnew[i]=alpha[i],i=1..N)]) mod p;
    end;
end:

GE:= proc(A::Matrix,p::prime)
local n,m,inv,mu,i,j,k,det;
   n,m := op(1,A);
   if n<>m then error "matrix must be square" fi;
   det := 1;
   for k to n do
       i := k;
       while i<=n and A[i,k]=0 do i := i+1; od;
       if i>n then return 0 fi;
       if i>k then # interchange row i and k
          for j from k to m do A[i,j],A[k,j] := A[k,j],A[i,j] od;
          det := -det;
       fi;
       det := det*A[k,k] mod p;
       inv := 1/A[k,k] mod p;
       for i from k+1 to n do
           if A[i,k]=0 then next fi;
           mu := A[i,k]*inv mod p;
           for j from k+1 to n do
               A[i,j] := A[i,j]-mu*A[k,j] mod p;
           od;
           A[i,k] := 0;
       od;
   od;
   det;
end:

MakeBBdet_Maple := proc( A::Matrix, VarPerm::{list}, opt::nonnegint ) 
   # opt = 2 uses GE to compute det(A)
   local X,N,Xnew,i;
   X := convert(indets(A),list); 
   N := nops(X); 
   Xnew := Array(1..N);
   for i to N do Xnew[i] := X[VarPerm[i]]; od; 
   proc( alpha::Array, p::prime ) global CNT,tBBeval,tBBdet; 
       local i,st,et,A_eval,d; CNT++;
       st := time(); #printf("Xnew = %a, alpha = %a\n",Xnew,alpha);
       A_eval := Eval( A,[seq(Xnew[i]=alpha[i],i=1..N)]) mod p;
       et := time()-st; tBBeval += et;
       if opt = 1 then #_params['opt'] = NULL then 
           st := time(); d := Det(A_eval) mod p;
           et := time()-st; tBBdet += et; 
       elif opt = 2 then  
           st := time(); d := GE(A_eval,p);
           et := time()-st; tBBdet += et; 
       fi; 
       return d;
   end;
end:

MakeBBdet_C := proc( A::Matrix, VarPerm::{list} ) 
    local n,X,N,Xnew,i,AL,AA;
    n := LinearAlgebra:-RowDimension(A);
    X := convert(indets(A),list); 
    N := nops(X); 
    Xnew := Array(1..N);
    for i to N do Xnew[i] := X[VarPerm[i]]; od; 
    Xnew := convert(Xnew,list); 
    AL := convert(A,list); 
    AA := Array(1..n,1..n,datatype=integer[8],order=C_order);
    proc( alpha::Array, p::prime ) global CNT,tBBeval,tBBdet; 
        local st,et,d; CNT++;
        st := time(); #printf("AL=%a,AA=%a,Xnew = %a, alpha=%a,p=%a\n",AL,AA,Xnew,alpha,p); 
        EVALMOD1( AL, Xnew, alpha, AA, p ); 
        et := time()-st; tBBeval += et;
        st := time(); 
        d := Det64s(AA, n, p); 
        et := time()-st; tBBdet += et; 
        return d;
    end;
end: