/*  bsi ... version 20090504
Risa/Asir programs for Bernstein-Sato ideals
by Toshinori Oaku

Usage: f_1,...,f_p are polynomials in variables x_1,...,x_n

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

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]
via primary decomposition
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) 

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

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 ann_{D[s]}f^s \cap Q[x,s]

bsilocalSE([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]
with one-by-one elimination of differentiations 
% SE for successive elimination: the efficiency depends on 
the order of x_1,...,x_n.
Example:
bsilocalSE([z,x^4+y^4+2*z*x^2*y^2],[z,y,x],[0,0,0]);

bsiglobalSE([f_1,...,f_p],[x_1,...,x_n]);
--> the global Bernstein-Sato ideal of [f_1,...,f_p] with one-by-one 
elimination of differentiations 

bsi_stratifySE([f_1,...,f_p],[x_1,...,x_n]);
--> stratified data for local Bernstein-Sato ideals 
with one-by-one elimination

bsi_timer(B); --> print timing data for each step if B != 0
bsi_verbose(B); --> print intermediate results if B != 0
*****************************************************************
Variables s1,tt1,uu1,vv1,... are reserved.
You may rename them by editing the following 4 lines:
*****************************************************************/
#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 */

extern  BSI_Verbose;         /* output intermediate results if not 0 */
extern  BSI_Timer;           /* output timing data for each step if not 0 */
/****************************************************************
The functions defined in this file:
annfs, bsilocal, bsilocalSE, bsiglobal, bsiglobalSE, bsi_stratify, 
bsi_stratifySE, bsi_step2, bsi_step2SE, bsi_stratify_data, 
bsi_global, bsi_local, bsi_timer, bsi_verbose, bsi_eliminate,
bsi_elim, bsi_intersect, bsi_mul, bsi_remove0, bsi_toSingular
*****************************************************************
Dependencies...the following functions in other library files are called:
psi in "bfct", primadec in "primdec"
*****************************************************************/
if (!module_definedp("bfct")) load("bfct")$ else{ }$  

def bf0(F,XX,AA) {
  BF = bsilocal([F],XX,AA);
  B = car(BF);
  Bs = subst(B,s1,s);
  return(Bs);
}


/* 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 = bsi_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 bsilocalSE(FF,XX,AA) {
  P = length(FF);
  for (I=P-1,SS=[]; I >= 0; I--) 
    SS = cons(strtov("s"+rtostr(I+1)),SS);
  QQ = bsi_stratifySE(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 (BSI_Timer) T0 = time(); 

  GGs = annfs(FF,XX);
  if (BSI_Timer) {
    T1 = time(); 
    print("Step 1 completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds"); 
  }
  if (BSI_Verbose) {print("ann f^s:"); print(GGs);}

  JJ = bsi_step2(GGs,FF,XX,SS);
  if (BSI_Timer) {
    T2 = time(); 
    print("Step 2 completed in "+rtostr(T2[0]-T1[0]+T2[1]-T1[1])+"seconds"); 
  }
  BS = bsi_eliminate(JJ,XX,SS);  
  if (BSI_Timer) {
    T3 = time(); 
    print("Step 3 completed in "+rtostr(T3[0]-T2[0]+T3[1]-T2[1])+"seconds"); 
  }
  return(BS);
}

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

def bsiglobalSE(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 (BSI_Timer) T0 = time(); 

  GGs = annfs(FF,XX);
  if (BSI_Timer) {
    T1 = time(); 
    print("Step 1 completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds"); }
  if (BSI_Verbose) {print("ann f^s:"); print(GGs);}

  JJ = bsi_step2SE(GGs,FF,XX,SS);
  if (BSI_Timer) {
    T2 = time(); 
    print("Step 2 completed in "+rtostr(T2[0]-T1[0]+T2[1]-T1[1])+"seconds"); 
  }
  BS = bsi_eliminate(JJ,XX,SS);  
  if (BSI_Timer) {
    T3 = time(); 
    print("Step 3 completed in "+rtostr(T3[0]-T2[0]+T3[1]-T2[1])+"seconds"); 
  }
  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 (BSI_Timer) { T0 = time(); }

  GGs = annfs(FF,XX);
  if (BSI_Timer) {
    T1 = time(); 
    print("Step 1 completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds"); 
  }
  if (BSI_Verbose) {print("ann f^s:"); print(GGs);}

  JJ = bsi_step2(GGs,FF,XX,SS);
  if (BSI_Timer) {
    T2 = time(); 
    print("Step 2 completed in "+rtostr(T2[0]-T1[0]+T2[1]-T1[1])+"seconds"); 
  }

  XS = append(XX,SS);
  PP = primadec(JJ,XS);
  if (BSI_Verbose) {
    print("Primary decomposition of ann f^s cap Q[x,s]:"); 
    print(PP); 
    print("Stratified data: [rad(Q_i), rad(Q_i cap Q[x]), Q_i cap Q[s])]");
  } 
  QQ = bsi_stratify_data(PP,XX,SS);
  K = length(QQ);
  for (I=0, Gen=[]; I < K; I++) {
    if (BSI_Verbose) {print(rtostr(I+1)+"th component:"); 
      print([ PP[I][0], QQ[I][1], QQ[I][0] ]); 
    }
    Gen = append(Gen,[length(QQ[I][0])]);
  }
  if (BSI_Timer) {
    T3 = time(); 
    print("Step 3 completed in "+rtostr(T3[0]-T2[0]+T3[1]-T2[1])+"seconds"); 
  }
  if (BSI_Verbose) {
    print("The number of generators of each (Q_i cap Q[s]):"); 
    print(Gen); }
  return(QQ);
}


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

def bsi_stratifySE(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 (BSI_Timer) { T0 = time(); }

  GGs = annfs(FF,XX);
  if (BSI_Timer) {
    T1 = time(); 
    print("Step 1 completed in "+rtostr(T1[0]-T0[0]+T1[1]-T0[1])+"seconds"); 
  }
  if (BSI_Verbose) {print("ann f^s:"); print(GGs);}

  JJ = bsi_step2SE(GGs,FF,XX,SS);
  if (BSI_Timer) {
    T2 = time(); 
    print("Step 2 completed in "+rtostr(T2[0]-T1[0]+T2[1]-T1[1])+"seconds"); 
  }
  XS = append(XX,SS);
  PP = primadec(JJ,XS);
  if (BSI_Verbose) {
    print("Primary decomposition of ann f^s cap Q[x,s]:"); 
    print(PP); 
    print("Stratified data: [rad(Q_i), rad(Q_i cap Q[x]), Q_i cap Q[s])]");
  } 
  QQ = bsi_stratify_data(PP,XX,SS);
  K = length(QQ);
  for (I=0, Gen=[]; I < K; I++) {
    if (BSI_Verbose) {print(rtostr(I+1)+"th component:"); 
      print([ PP[I][0], QQ[I][1], QQ[I][0] ]); 
    }
    Gen = append(Gen,[length(QQ[I][0])]);
  }
  if (BSI_Timer) {
    T3 = time(); 
    print("Step 3 completed in "+rtostr(T3[0]-T2[0]+T3[1]-T2[1])+"seconds"); 
  }
  if (BSI_Verbose) {
    print("The number of generators of each (Q_i cap Q[s]):"); 
    print(Gen); }
  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);
  XS = append(XX,SS);
  DXX = [];
  for (I=N-1; I >= 0; I--) 
    DXX = cons(strtov("d"+rtostr(XX[I])),DXX);
  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 (BSI_Verbose) {
    print("Eliminating "+rtostr(DXX)+" with a term order:");
    print(XSDXS); print(MD);
  }
  F = 1;
  for (J=0; J < P; J++) F = F*FF[J];

  G2 = cons(F,GGs);
  G2 = dp_weyl_gr_main(G2,XSDXS,0,0,MD);
  G2 = bsi_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_step2SE(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);
  F = 1;
  for (J=0; J < P; J++) F = F*FF[J];
  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 (BSI_Verbose) {
      print("Eliminating "+rtostr(DXX[II])+" with a term order:");
      print(XSDXS); print(MD);
    }
    G2 = dp_weyl_gr_main(G2,XSDXS,0,0,MD);
    G2 = bsi_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 bsi_stratify_data(Q,X,S) {
  N = length(Q);
  QQ = [];
  for (I=0; I < N; I++)
/*    QQ = append(QQ,[[bsi_eliminate(Q[I][0],X,S),
        primedec(bsi_eliminate(Q[I][1],S,X),X)]]); 
*/
QQ = append(QQ,[[bsi_eliminate(Q[I][0],X,S),bsi_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 = bsi_intersect(BS,QQ[I][0],S);
  return(BS);
}

/* compute the local BS ideal at A from the stratification data 
QQ: the output of bsi_stratify, X: the x-variables
S: the s-variables, A: a point */

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 = bsi_intersect(BS,QQ[I][0],S);
  }
  return(BS);
}

def bsi_timer(X) { BSI_Timer = X; }

def bsi_verbose(X) { BSI_Verbose = X; }

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

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

def bsi_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 bsi_elim(L,V);
}

/* returns elements of L free from variables V */
def bsi_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 bsi_intersect(I,J,X) {
    I = bsi_remove0(I); J = bsi_remove0(J);
    T = strtov(Tchar+"1");
    if ( I != [] && J != [] ) {
	VX = cons(T,X);
	It = map(bsi_mul,I,T);
	Jt = map(bsi_mul,J,1-T);
	GIJ = gr(append(It,Jt),VX,[[0,1],[0,length(X)]]);
	IJ = bsi_elim(GIJ,[T]);
	return IJ;
    }
    else return [];
}

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


/* remove 0 elements */
def bsi_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 from Risa/Asir to Singular 
JJ: an ideal as a list of generators in Asir
Name: the name of the ideal in Singular */

def bsi_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$