/**   bsi  version20090326
Risa/Asir programs for Bernstein-Sato ideals
by Toshinori Oaku 
***********************************************/

/* Usage:

annfs([f_1,...,f_p],[x_1,...,x_n]);
--> the annihilator ideal of f_1^{s1}...f_p^{sp}
by the algorithm of Oaku-Takayama (JPAA 139 (1999), 201--233)

bsiglobal([f_1,...,f_p],[x_1,...,x_n]);
--> the global Bernstein-Sato ideal of [f_1,...,f_p]

bsilocal([f_1,...,f_p],[x_1,...,x_n],[a_1,...,a_n]);
--> the local Bernstein-Sato ideal of [f_1,...,f_p] at [a_1,...,a_n]

Example:
bsilocal([z,x^4+y^4+2*z*x^2*y^2],[x,y,z],[0,0,0]);
an example of non-principal BS ideal by Briancon-Maynadier
(J. Math. Kyoto Univ. 39 (1999), 215--232) 

bsiglobalHE([f_1,...,f_p],[x_1,...,x_n]);
--> the global Bernstein-Sato ideal of [f_1,...,f_p] by one-by-one 
elimination of differentiations (the efficiency depends on the order
of x_1,...,x_n)   % HE for heuristic elimination

bsilocalHE([f_1,...,f_p],[x_1,...,x_n],[a_1,...,a_n]);
--> the local Bernstein-Sato ideal of [f_1,...,f_p] at [a_1,...,a_n]
by one-by-one elimination of differentiations

bsi_stratify([f_1,...,f_p],[x_1,...,x_n]);
--> stratified data for local Bernstein-Sato ideals in a list of
[ Q_i cap Q[s], the radical of Q_i cap Q[x] ],
where Q_i are the primary components of the output ideal of Step 2 
of the algorithm of Bahloul-Oaku. 

*********************************************************************
Variables s1,tt1,uu1,vv1,... are reserved.
You may replace them by other strings by editing the next 4 lines:
*********************************************************************/
#define Schar "s"            /* variables for BS ideals */
#define Tchar "tt"           /* variables for the Malgrange ideals */
#define Uchar "uu"           /* variables for V-homogenization */
#define Vchar "vv"           /* variables for saturation */

extern  Verbose;             /* output intermediate results if not 0 */
extern  Timer;               /* output timing data for each step if not 0 */

if (!module_definedp("bfct")) load("bfct")$ else{ }$


/* the annihilator ideal for 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 = DW[I];
    FI = DXI;
    for (J=0; J < P; J++) 
      FI += diff(FF[J],XI)*UU[J]*DTT[J];
    GG = append(GG,[FI]);
  }
  if (Verbose==1) print(["GG",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 (Verbose==1) print(WDW);
  if (Verbose==1) print(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);
}


/* local Bernstein-Sato ideal of FF in XX at AA 
usage: bsilocal([f_1,...,f_p],[x_1,...,x_n],[a_1,...,a_n]);
*/

def bsilocal(FF,XX,AA) {
  P = length(FF);
  for (I=P-1,SS=[]; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  QQ = bsi_stratify(FF,XX);  
  BS = bsi_local(QQ,XX,SS,AA);
  return(BS);
}


/* local Bernstein-Sato ideal of FF in XX at AA with one by one elimination */

def bsilocalHE(FF,XX,AA) {
  P = length(FF);
  for (I=P-1,SS=[]; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  QQ = bsi_stratifyHE(FF,XX);  
  BS = bsi_local(QQ,XX,SS,AA);
  return(BS);
}


/* Global Bernstein-Sato ideal of FF in variables XX */

def bsiglobal(FF,XX) {
  P = length(FF);
  N = length(XX);
  SS = [];
  for (I=P-1; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  if (Timer) T0 = time(); 

  GGs = annfs(FF,XX);
  if (Timer) {T1 = time(); print(["Step 1:", T1[0]-T0[0]+T1[1]-T0[1]]);}
  if (Verbose==1) print(["GGs",GGs]);

  JJ = bsi_step2(GGs,FF,XX,SS);
  if (Timer) {T2 = time(); print(["Step 2:", T2[0]-T1[0]+T2[1]-T1[1]]);}
  BS = eliminate(JJ,XX,SS);  
  if (Timer) {T3 = time(); print(["Step 3:", T3[0]-T2[0]+T3[1]-T2[1]]);}
  return(BS);
}

/* Global Bernstein-Sato ideal of FF in variables XX 
with one by one elimination */

def bsiglobalHE(FF,XX) {
  P = length(FF);
  N = length(XX);
  SS = [];
  for (I=P-1; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  if (Timer) T0 = time(); 

  GGs = annfs(FF,XX);
  if (Timer) {T1 = time(); print(["Step 1:", T1[0]-T0[0]+T1[1]-T0[1]]);}
  if (Verbose==1) print(["GGs",GGs]);

  JJ = bsi_step2HE(GGs,FF,XX,SS);
  if (Timer) {T2 = time(); print(["Step 2:", T2[0]-T1[0]+T2[1]-T1[1]]);}
  BS = eliminate(JJ,XX,SS);  
  if (Timer) {T3 = time(); print(["Step 3:", T3[0]-T2[0]+T3[1]-T2[1]]);}
  return(BS);
}


/* Bernstein-Sato ideal with a stratification */

def bsi_stratify(FF,XX) {
  P = length(FF);
  N = length(XX);
  SS = [];
  for (I=P-1; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  if (Timer) { T0 = time(); }

  GGs = annfs(FF,XX);
  if (Timer) {T1 = time(); print(["Step 1:", T1[0]-T0[0]+T1[1]-T0[1]]);}
  if (Verbose==1) print(["GGs",GGs]);

  JJ = bsi_step2(GGs,FF,XX,SS);
  if (Timer) {T2 = time(); print(["Step 2:", T2[0]-T1[0]+T2[1]-T1[1]]);}

  XS = append(XX,SS);
  PP = primadec(JJ,XS);
  if (Verbose==1) print(["primary decomposition of JJ",PP]); 
  QQ = stratify(PP,XX,SS);
  K = length(QQ);
  for (I=0, Gen=[]; I < K; I++) {
    if (Verbose==1) print(PP[I][0]);
    Gen = append(Gen,[length(QQ[I][0])]);
  }
  if (Timer) {T3 = time(); print(["Step 3:", T3[0]-T2[0]+T3[1]-T2[1]]);}
  if (Verbose==1) print(Gen);

  for (I=0, RR=[]; I < K; I++) {
    RR = append(RR,[[PP[I],QQ[I][1],QQ[I][0]]]);
  }
  if (Verbose==11) print(RR);
  return(QQ);
}


/* Bernstein-Sato ideal with a stratification 
with one by one elimination */

def bsi_stratifyHE(FF,XX) {
  P = length(FF);
  N = length(XX);
  SS = [];
  for (I=P-1; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  if (Timer) { T0 = time(); }

  GGs = annfs(FF,XX);
  if (Timer) {T1 = time(); print(["Step 1:", T1[0]-T0[0]+T1[1]-T0[1]]);}
  if (Verbose==1) print(["GGs",GGs]);

  JJ = bsi_step2HE(GGs,FF,XX,SS);
  if (Timer) {T2 = time(); print(["Step 2:", T2[0]-T1[0]+T2[1]-T1[1]]);}

  XS = append(XX,SS);
  PP = primadec(JJ,XS);
  if (Verbose==1) print(["primary decomposition of JJ",PP]); 
  QQ = stratify(PP,XX,SS);
  K = length(QQ);
  for (I=0, Gen=[]; I < K; I++) {
    if (Verbose==1) print(PP[I][0]);
    Gen = append(Gen,[length(QQ[I][0])]);
  }
  if (Timer) {T3 = time(); print(["Step 3:", T3[0]-T2[0]+T3[1]-T2[1]]);}
  if (Verbose==1) print(Gen);

  for (I=0, RR=[]; I < K; I++) {
    RR = append(RR,[[PP[I],QQ[I][1],QQ[I][0]]]);
  }
  if (Verbose==1) print(RR);
  return(QQ);
}


/* Step 2 of the algoritm: outputs an ideal of Q[x,s] */

def bsi_step2(GGs,FF,XX,SS) {
  P = length(SS);
  N = length(XX);

  DXX = [];
  for (I=N-1; I >= 0; I--) 
    DXX = cons(strtov("d"+rtostr(XX[I])),DXX);
  XS = append(XX,SS);
  for (I = N+P-1, DXS = []; I >= 0; I--)
    DXS = cons(strtov("d"+rtostr(XS[I])),DXS);
  XSDXS = append(XS,DXS);

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

  if (Verbose==1) print(XSDXS);
  if (Verbose==1) print(MD);

  F = 1;
  for (J=0; J < P; J++) F = F*FF[J];
  if (Verbose==1) print(["F",F]);

  G2 = cons(F,GGs);
  G2 = dp_weyl_gr_main(G2,XSDXS,0,0,MD);
  G2 = elim(G2,DXX);
  JJ = G2;
  return(JJ);
}


/* Step 2 of the algoritm: outputs an ideal of Q[x,s]  
with one by one elimination */

def bsi_step2HE(GGs,FF,XX,SS) {
  P = length(SS);
  N = length(XX);

  DXX = [];
  for (I=N-1; I >= 0; I--) 
    DXX = cons(strtov("d"+rtostr(XX[I])),DXX);
  XS = append(XX,SS);
  for (I = N+P-1, DXS = []; I >= 0; I--)
    DXS = cons(strtov("d"+rtostr(XS[I])),DXS);
  XSDXS = append(XS,DXS);

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

  if (Verbose==1) print(XSDXS);
  if (Verbose==1) print(MD);

  F = 1;
  for (J=0; J < P; J++) F = F*FF[J];
  if (Verbose==1) print(["F",F]);

  G2 = cons(F,GGs);
  for (II=0; II<N; II++) {
    MD = newmat(2*(N+P)+1,2*(N+P));
    MD[0][N+P+II] = 1;
    for (J=0; J < 2*(N+P); J++) MD[1][J] = 1;
    for (I=2; I < 2*(N+P)+1; I++) MD[I][I-2]= -1;
    if (Verbose==1)  print(MD);
    G2 = dp_weyl_gr_main(G2,XSDXS,0,0,MD);
    G2 = elim(G2,[DXX[II]]);
  }
  JJ = G2;
  return(JJ);
}


/* compute the intersections with X and with S for each primary component 
of Q as a list of [ Q \cap Q[s], \sqrt Q \cap Q[x] ] */

def stratify(Q,X,S) {
  N = length(Q);
  QQ = [];
  for (I=0; I < N; I++)
/*    QQ = append(QQ,[[eliminate(Q[I][0],X,S),
        primedec(eliminate(Q[I][1],S,X),X)]]); 
*/
QQ = append(QQ,[[eliminate(Q[I][0],X,S),eliminate(Q[I][1],S,X)]]); 
  return(QQ); 
}


/* computes the global BS from the output of bsi_stratify() */

def bsi_global(QQ,S) {
  K = length(QQ);
  BS = [1];
  for (I=0; I < K; I++)
    BS = intersect(BS,QQ[I][0],S);
  return(BS);
}

/* compute the local BS at A from the output of bsi_stratify() */

def bsi_local(QQ,X,S,A) {
  N = length(X);
  L = length(QQ);
  BS = [1];
  Flags = [];
  for (I=0; I < L; I++) {
    NI = length(QQ[I][1]);
    Flag = 1;
    for (J=0; J<NI; J++) {
      FIJ = QQ[I][1][J];
      for (K=0; K < N; K++) FIJ = subst(FIJ,X[K],A[K]);
      if (FIJ != 0) Flag = 0;
    }
    Flags = cons(Flag,Flags);
    if (Flag) BS = intersect(BS,QQ[I][0],S);
  }
  return(BS);
}


/************ miscellaneous functions ************************/

/* returns elements of L free from variables V 
by elimination via Groevner base 
L: a generating set of an ideal
V: list of variables to be eliminated
X: list of the other variables
*/

def eliminate(L,V,X) {
	NV = length(V); NX = length(X);
	Weight = newmat(NV+NX+1,NV+NX);
	for (J=0; J < NV; J++) Weight[0][J] = 1;
        for (J=0; J < NV+NX; J++) Weight[1][J] = 1;
	for (I=2; I < NV+NX+1; I++) Weight[I][I-2]= -1;
	L = gr(L,append(V,X),Weight);
	return elim(L,V);
}


/* 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;
}

/* the intersection of the ideals I and J in Q[X] */
def intersect(I,J,X) {
    I = remove0(I); J = remove0(J);
    T = strtov(Tchar+"1");
    if ( I != [] && J != [] ) {
	VX = cons(T,X);
	It = map(mul,I,T);
	Jt = map(mul,J,1-T);
	GIJ = gr(append(It,Jt),VX,[[0,1],[0,length(X)]]);
	IJ = elim(GIJ,[T]);
	return IJ;
    }
    else return [];
}

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


/* remove 0 elements */
def remove0(L) {
    ZERO = 0; 
    for (LL=[]; L != []; L = cdr(L)) {
      F = car(L);
      if (F != ZERO) LL = cons(F,LL);
    }
    return LL;
}

/* convert an ideal to Singular data */ 

def toSingular(JJ,Name) {
  L = length(JJ);
  Str = "ideal " + Name + " = ";
  for (I=0; I<L-1; ++I) 
    Str = Str + rtostr(JJ[I]) + ",\n";
  Str = Str + rtostr(JJ[L-1]) + ";";
  return Str;
}

end$