/* annlog  
a Risa/Asir program for computing annihilators for
  f^s, f^s(lof f), ..., f^s(log f)^m
by Toshinori Oaku
2010 August 18 --- 
last modified July 3, 2012. 
*/

/*
Usage: F is a polynomial in variables x_1,...,x_n

annlog(F,a,m,[x1,...,x_n])
--> the annihilating ideal of F^a(log F)^m
annlogs(F,a,m,[x1,...,x_n])
--> the annihilating module of (F^a, F^a log F,...,F^a(log F)^m)
annfslog(F,m,[x1,...,x_n])
--> the annihilating ideal of F^s(log F)^m with a parameter s
annlogs(F,m,[x1,...,x_n])
--> the annihilating module of (F^s, F^s log F,...,F^s(log F)^m)
with a paremeter s

The variable s is reserved!
*/

#define Schar "s"            /* variables for BS ideals */
#define Tchar "_tt"           /* variables for the Malgrange ideal */
#define Uchar "_uu"           /* variables for V-homogenization */
#define Vchar "_vv"           /* variables for saturation */

/* the annihilator ideal of f^s = f_1^{s1}...f_p^{sp} */

def annfs(FF,XX) {
  P = length(FF);
  N = length(XX);
  NT = N + 3*P;

  DXX = [];
  for (I=N-1; I >= 0; I--) 
    DXX = cons(strtov("d"+rtostr(XX[I])),DXX);
  
  SS = []; TT = []; DTT = []; UU = []; VV = [];
  for (I=P-1; I >= 0; I--) {
    SS = cons(strtov(Schar+rtostr(I+1)),SS);
    TT = cons(strtov(Tchar+rtostr(I+1)),TT);
    DTT = cons(strtov("d"+Tchar+rtostr(I+1)),DTT);
    UU = cons(strtov(Uchar+rtostr(I+1)),UU);
    VV = cons(strtov(Vchar+rtostr(I+1)),VV);
  }
  W = append(XX,TT); W = append(W,UU); W = append(W,VV);

  for (I = NT-1, DW = []; I >= 0; I-- )
    DW = cons(strtov("d"+rtostr(W[I])),DW);
  WDW = append(W,DW);

  GG = [];
  for (J=0; J < P; J++)
    GG = append(GG,[ TT[J] - UU[J]*FF[J], UU[J]*VV[J]-1 ]);
  for (I=0; I < N; I++) { 
    XI = XX[I]; DXI = DXX[I];
    FI = DXI;
    for (J=0; J < P; J++) 
      FI += diff(FF[J],XI)*UU[J]*DTT[J];
    GG = append(GG,[FI]);
  }
  if (BSI_Verbose) { print("The Malgrange ideal:"); print(GG); }

/* Set the weight vector MW for the variables WDW:
 
  x1..xn t1..tp u1..up v1..vp dx1..dxn dt1..dtp du1..dup dv1..dvp 

  0    0 0    0 1    1 1    1 0      0 0      0 0      0 0      0
  1    1 1    1 0    0 0    0 1      1 1      1 1      1 1      1
  -1
                        ...
                                               -1  
*/

  MW = newmat(2*NT+1,2*NT);
  for (J = N+P; J < N+3*P; J++) MW[0][J]=1;
  for (J=0; J < N+P; J++) MW[1][J]=1;
  for (J=N+3*P; J < 2*NT; J++) MW[1][J]=1;
  for (I=2; I < 2*NT+1; I++) MW[I][I-2]=-1;

  if (BSI_Verbose) { 
    print("Step 1: Computing a Groebner base w.r.t. a term order:");
    print(WDW); print(MW); }

/*  GG1 = nd_weyl_gr(GG,WDW,0,MW); */
  GG1 = dp_weyl_gr_main(GG,WDW,0,0,MW); 

  GG2 = elim(GG1,append(UU,VV));
  GGs = GG2;
  for (J=0; J < P; J++) {
    GGs = map(psi,GGs,TT[J],DTT[J]);
    GGs = map(subst,GGs,TT[J],-SS[J]-1);
  }
  return(GGs);
}

def mulv(A,B) {
  for (C = []; B != [];  ) {
    C = cons(A*car(B),C);
    B = cdr(B);
  }
  return(reverse(C));
}

/* annihilator of (F^s,...,F^s(log F)^K) in D[s]^{K+1} */
def annfslogs(F,K,XX) {
  S = strtov(Schar); 
  S1 = strtov(Schar+rtostr(1)); 
  Annfs = subst(annfs([F],XX),S1,S); 
  for (Annfslog = []; Annfs != []; Annfs = cdr(Annfs)) {
    for (I=0; I<=K; ++I) {
      Plogs = newvect(K+1);
      P = car(Annfs);
      for (J=I; J >= 0; J--) {
        Plogs[J] = choose(I,I-J)*P;
        P = diff(P,s);
      }
      Annfslog = cons(vtol(Plogs),Annfslog);
    }
  }
  return(Annfslog);
}
    
/*
divide a differential operator P by Q w.r.t. revlex ordering
for variables XX and DXX and returns the quotient U 
such that P = UQ + R
*/
def divD(P,Q,XX) {
  XDX = xdx(XX);
  dp_ord(0); 
  LT=dp_dtop(dp_hm(dp_ptod(Q,XDX)),XDX);
  U=0; /* qutotient */  
  R=P; /* remainder　*/
  while (( R != 0) && dp_redble(dp_ptod(R,XDX),dp_ptod(Q,XDX))) { 
    S=red(dp_dtop(dp_hm(dp_ptod(R,XDX)),XDX)/LT);
    U=U+S;
    R=R-dp_dtop(dp_weyl_mul(dp_ptod(S,XDX),dp_ptod(Q,XDX)),XDX);
/*  print([U,R]); */
  }
/*  return ([U,R]); */
  return(U); 
} 

def choose(N,K) {
  C = 1;
  for (I=0; I<K; ++I) C = C*(N-I);
  return(C/fac(K));
}

/* annihilator of F^s(log F)^M in D[s] 
   by elimination in D[s]^(M+1)
*/
def annfslog(F,M,XX) {
  N = length(XX);
  DXX = [];
  for (I=N-1; I >= 0; I--) 
    DXX = cons(strtov("d"+rtostr(XX[I])),DXX);
  XSDXS = append([s],XX);
  XSDXS = append(XSDXS,[ds]);
  XSDXS = append(XSDXS,DXX);
  Ann = annfslogs(F,M,XX);
  T0 = time();
  AnnG = nd_weyl_gr(Ann,XSDXS,0,[1,0]);
  Ann0 = last_component_only(AnnG);
  T1 = time();
  print("Elimination completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds"); 
  return(Ann0);
}

/* the annihilator of the vector [F^A,F^A(log F),...,F^A(log F)^K]
*/
def annlogs(F,A,K,XX) {
  Annfslogs = annfslogs2(F,K,XX);
  Annfs = car(Annfslogs); 
  Annfslog = Annfslogs[1];
  B = bfminRoot(F,Annfs,A,XX);
  print(["Substituting s =", B]);  
  M = A-B;
  AnnfBlog = subst(Annfslog,s,B);
  print("Computing module quotient:");
  T0 = time();
  AnnfAlog = modulequo(AnnfBlog,F^M,XX);
  T1 = time();
  print("quotient completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds"); 
  return(AnnfAlog);
}

/* the annihilator of F^A(log F)^K
*/
def annlog(F,A,K,XX) {
  AnnFAlogs = annlogs(F,A,K,XX);
  N = length(XX);
  T0 = time();
  AnnG = nd_weyl_gr(AnnFAlogs,xdx(XX),0,[1,0]);
  Ann0 = last_component_only(AnnG);
  T1 = time();
  print("Elimination completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds");  return(Ann0);
}

/* the minimum root of the b-function
of the form A-1, A-2,... 
returns A if none.
*/
def bfminRoot(F,Ann,A,XX) { 
  N = length(XX);
  Min = A;
  for (I=1; A-I >= -N; ++I) {
    Ans = checkRoot(F,Ann,XX,A-I);
    if (Ans == 1) Min = A-I;
  }
  return(Min);
}

def checkRoot(F,Ann,XX,A) {
  Ann = subst(Ann,s,A);
  G = nd_weyl_gr(cons(F,Ann),xdx(XX),0,0);
  Ans = 1;
  if ((length(G) == 1) && (type(car(G)) == 1) && (car(G) != 0)) 
    Ans = 0;
  return(Ans); 
}

/* annihilator of (F^s,...,F^s(log F)^K) in D[s]^{K+1} 
   together with that of F^s 
*/
def annfslogs2(F,K,XX) {
  Annfs = subst(annfs([F],XX),s1,s);
  Ops = Annfs;
  for (Annfslog = []; Ops != []; Ops = cdr(Ops)) {
    for (I=0; I<=K; ++I) {
      Plogs = newvect(K+1);
      P = car(Ops);
      for (J=I; J >= 0; J--) {
        Plogs[J] = choose(I,I-J)*P;
        P = diff(P,s);
      }
      Annfslog = cons(vtol(Plogs),Annfslog);
    }
  }
  return([Annfs,Annfslog]);
}

def intersection_module(II,JJ,XX) {
  T = strtov(Tchar);

  for (IIT=[]; II != []; II = cdr(II)) 
    IIT = cons(mulv(T,car(II)),IIT);
  IIT=reverse(IIT);
  for (JJT=[]; JJ != []; JJ = cdr(JJ)) 
    JJT = cons(mulv(1-T,car(JJ)),JJT);
  JJT=reverse(JJT);
  IJT=append(IIT,JJT);

  N = length(XX);
  XT = cons(T,XX);
  XTDXT = xdx(XT);

  M=newmat(2*N+3,2*N+2);
  M[0][0]=1;
  for (J=0; J < 2*N+2; ++J) M[1][J]=1;
  for (I=2; I < 2*N+3; ++I) M[I][2*N+3-I]=-1;

  G=nd_weyl_gr(IJT,XTDXT,0,[0,M]);
  GT=eli_module(G,T);  
  return(GT);
}

def eli_module(X,T) {
  L=length(X);
  Y = [];
  for (I=0; I<L; I++) {
    if (deg_module(X[I],T)<1) {
      Y = cons(X[I],Y);
    }   
  }
  Y = reverse(Y);
  return(Y);
}

def deg_module(P,T) {
  Deg = -1;
  M = length(P);
  for (I=0; I < M; I++) 
    if (deg(P[I],T) > Deg) Deg = deg(P[I],T);
  return(Deg);
}

/* quotient module R:Q = {(P1,...,Pr) | (P1Q1,...,PrQr) in R}
of the module R by a vector Q */
def quotient_module(R,Q,XX) {
  M = length(car(R));
  if (type(Q) == 2) {
    Qlist = [];
    for (I=0; I<M; ++I)
      Qlist = cons(Q,Qlist);
    Q = Qlist;
  }
  Qvecs = [];
  for (I=0; I<M; ++I) {
    Qvec = newvect(M);
    Qvec[I]=Q[I];
    Qvec = vtol(Qvec);
    Qvecs = cons(Qvec,Qvecs);
  }
  Qvecs = reverse(Qvecs);

  G = intersection_module(Qvecs,R,XX);

  Quo = [];
  N = length(G);
  for (J=0; J < N; J++) {
    QuoElement = newvect(M);
    for (I=0; I < M; I++)    
      QuoElement[I] = divD(G[J][I],Q[I],XX);
    Quo = cons(vtol(QuoElement),Quo);
  }
  return(reverse(Quo));
}

/* returns [x,y,..,dx,dy,..] for X = [x,y,..] */
def xdx(XX) {
  N = length(XX);
  DXX = [];
  for (I=N-1; I >= 0; I--) 
    DXX = cons(strtov("d"+rtostr(XX[I])),DXX);
  XXDXX = append(XX,DXX);
  return(XXDXX);
}

/* returns lists of a list L of lists free from variables V */
def elim_list(L,V) {
	for (LL=[]; L != []; L = cdr(L)) {
		F = car(L); N = length(F);
		G = elim(F,V);
		if (length(G) == N) LL = cons(F,LL);
	}
	return reverse(LL);
}

/* module quotient II:F for a left D-module II and a polynomial F */

def modulequo(II,F,X) {
  R = length(car(II));
  T = strtov(Tchar);
  TX = cons(T,X);
  N = length(X);
  DTX = [];
  for (I=N; I >= 0; I--)
    DTX = cons(strtov("d"+rtostr(TX[I])),DTX);
  W = newmat(2*(N+1)+1, 2*(N+1));
  W[0][0] = 1;
  for (J=0; J < 2*(N+1); J++) W[1][J] = 1;
  for (I=2; I <= 2*(N+1); I++) W[I][I-2] = -1;
  for (JJ=[]; II != []; II = cdr(II)) 
    JJ = cons(map(mul,car(II),T),JJ);
  for (I=0; I < R; ++I) {
    FI = newvect(R); 
    FI[I] = (1-T)*F;
    FI = vtol(FI);
    JJ = cons(FI,JJ);
  }
  Gr = nd_weyl_gr(JJ,append(TX,DTX),0,[0,W]);
  JJ0 = elim_list(Gr,[T]);
  for (JJ1 = []; JJ0 != []; JJ0 = cdr(JJ0)) 
    JJ1 = cons(map(divD,car(JJ0),F,X),JJ1);
  return(JJ1);
}

def mul(A,B) {
return A*B;
}

/* returns elements of L free from variables V */
def elim(L,V) {
	for (LL=[]; L != []; L = cdr(L)) {
		F = car(L);
		E = 1;
		W = V;
		for (; W != []; W = cdr(W)) 	
			if (deg(F,car(W)) > 0) {
				E = 0; break; 
			}
		if (E == 1) LL = cons(F,LL);
	}
	return LL;
}

/* returns lists of form [0,...,0,*] */
def last_component_only(L) {
  if (L == []) return([]);
  L0 = car(L); 
  N = length(L0);
  for (LL = [] ; L != []; L = cdr(L) ) {
    L0 = car(L);
    OK = 1;
    for (I=0; I < N-1; I++)
      if (L0[I] != 0) OK = 0;
    if (OK == 1) LL = cons(L0[N-1],LL);
  }
  return(reverse(LL));
}

end$