// Computation of the Lojasiewicz exponent of an ideal.
// 
// Author: Santiago Encinas
// http://www.math.arq.uva.es/sencinas
// 

// You need to load the library desing for Singular
// it is avalaible at
// http://www.risc.uni-linz.ac.at/projects/basic/adjoints/blowup/
LIB "desing.lib";
desinginit();

// The followig are configuation variables for the package desing. You may
// remove or add the comments as needed.
// For more details see the help included with desing package
// 
// Set DCDEBUG to 1 if you want to generate a file "desing.log" with all
// computations involved.
// Change the directory name to an existing directory or set DCDEBUG to
// zero if you do not want the output file
Desing::DCFG[Desing::DCDEBUG] = 1;
Desing::DCFG[Desing::DCDEBUGPATH] = "~/tmp";
// Set DCHTML to 1 if you want to generate html files with the resolution
// tree, each file corresponds to an affine chart
// Change the directory name to an existing directory or set DCFG to
// zero if you do not want the output files
Desing::DCFG[Desing::DCHTML] = 1;
Desing::DCFG[Desing::DCHTMLPATH] = "~/tmp";
// DCTR=1 Compute transforms using NormalForm
// DCTR=0 Compute transforms using extended Grobner basis
Desing::DCFG[Desing::DCTR] = 1;
// DCDOWNGRADE Use local variables locales (0) or global variables (1)
Desing::DCFG[Desing::DCDOWNGRADE] = 1;
// DCNCTEST Perform normal crossing test (1)
// or fit Villamayor algorithm (0)
Desing::DCFG[Desing::DCNCTEST] = 0;
// DCPRINCIPALIZE=0 desingularization
// DCPRINCIPALIZE=1 principalization.
Desing::DCFG[Desing::DCPRINCIPALIZE] = 1;
// Desing::DCFG[Desing::DCHLTMODE]=0;
// Desing::DCFG[Desing::DCGLORD] = 1;
// DCPROGRESS Performs how the chart are being computed
// 0 like a stack, 1 like a queue.
// If 1 computation are finished for each blowing (more or less).
Desing::DCFG[Desing::DCPROGRESS]=1;
// Seconds to generate the tree.
Desing::DCFG[Desing::DCTTIME]=240;

proc ideall(list L)
// Convert a list to an ideal
{
    ideal LI;
    if (size(L))
    {
	LI = ideal(L[1..size(L)]);
    }
    return(LI);
}

proc Inclusion(ideal J, poly f)
// Returns 1 if J is included in the ideal (f), 0 otherwise
{
    int c=1;
    ideal F=std(f);
    for (int i=1; i<=size(J); i=i+1)
    {
	if (NF(J[i],F)<>0)
	{
	    c=0;
	}
    }
    return(c);
}

proc ExpDivide(ideal J, poly f)
// Returns the maximal integer such that J=f^a I
{
    int e=0;
    if (Inclusion(J,f))
    {
	ideal Q=quotient(J,f);
	e=1+ExpDivide(Q,f);
    }
    return(e);
}

proc ExpLojasiewicz(J,I,P)
// Computes the Lojasewicz exponent of the ideal I with respec to J at the
// point P
// P is an ideal defining the point.
// The Lojasewicz exponent is the maximum p/q such that
// J^p\subset\overline{I^q} at P
{
    def RR=basering;
    Desing::DCFG[Desing::DCPRINCIPALIZE] = 1;
    desing(I);
    def R0=Desing::DCR0;
    setring(R0);
    ideal I=imap(RR,I);
    ideal J=imap(RR,J);
    ideal P=imap(RR,P);
    int p=0;
    int q=1;
    for (int i=1; i<=size(Desing::DRC); i=i+1)
    {
	def R=Desing::DCR[Desing::DRC[i]];
	setring(R);
	map f=R0,ideall(DCH[4]);
	ideal P=f(P);
	ideal I=f(I);
	ideal J=f(J);
	ideal Dep=std(ideall(DCH[1]));
	qring RQ=std(Dep);
	setring(RQ);
	def DCH=imap(R,DCH);
	def listaE=DCH[11][1][4];
	ideal I=imap(R,I);
	ideal J=imap(R,J);
	ideal P=imap(R,P);
	poly h; int a; int b;
	for (int j=1; j<=size(listaE); j=j+1)
	{
	    h=DCH[2][listaE[j]];
	    ideal mm=std(ideal(P+h));
	    if (mm[1]<>1)
		{
		a=ExpDivide(I,h);
		b=ExpDivide(J,h);
		// Remove this comments if you want some info on how the maximum is 
		// achieved on charts
		// printf("Chart %s, DCH8: %s, ID: %s, h: %s",i,DCH[8],Desing::DRC[i],h);
		// printf("a= %s, b= %s",a,b);
		if (p*b<a*q)
		    {
		    p=a;
		    q=b;
		    // printf("p= %s, q= %s",p,q); // More print Info
		}
	    }
	}
    }
    setring(RR);
    return(number(p)/number(q)));
}

// The following are some examples. Remove the comments in order to do the
// computation

/* 
 * ring r=0,(x(1..3)),dp;
 * poly f=x(2)^6+x(3)^4+x(1)*(x(1)-3*x(3))^2;
 * poly g1=diff(f,x(1));
 * poly g2=diff(f,x(2));
 * poly g3=diff(f,x(3));
 * ideal I=g1,g2,g3;
 * ideal J=x(1),x(2),x(3);
 * ideal P=x(1),x(2),x(3);
 * print("I");
 * I;
 * print("J");
 * J;
 * print("Point P");
 * P;
 * ExpLojasiewicz(J,I,P);
 * // Result is 5
 */

/* 
 * ring r=0,(x(1..3)),dp;
 * poly g1=x(1)^3+x(2)^2;
 * poly g2=x(1)*x(2)^2+7*x(1)^4+x(3)^3;
 * poly g3=x(3)^2;
 * ideal I=g1,g2,g3;
 * ideal J=x(1),x(2),x(3);
 * ideal P=x(1),x(2),x(3);
 * print("I");
 * I;
 * print("J");
 * J;
 * print("Point P");
 * P;
 * ExpLojasiewicz(J,I,P);
 * // Result is 4
 */

/* 
 * ring r=0,(x(1..3)),dp;
 * poly g1=x(1)^3+x(2)^2;
 * poly g2=x(1)*x(2)^2+x(1)^4+x(3)^3;
 * poly g3=x(1)+x(2)+x(3);
 * ideal I=g1,g2,g3;
 * ideal J=x(1),x(2),x(3);
 * ideal P=x(1),x(2),x(3);
 * print("I");
 * I;
 * print("J");
 * J;
 * print("Point P");
 * P;
 * ExpLojasiewicz(J,I,P);
 * // Result is 3
 */