// This is a Singular file that calculates the solution presented in 
// the section 3 of the article:
// "On the Ernst electro-vacuum equations and ergosurfaces"

//We activate the timer to obtain the time for an individual command.

timer=1;

//the command "option(prot);" is useful for big Groebner basis 
//computations since it shows how the computation proceeds 
//(and let you know that the computation is still going on). 

option(prot);

//The variable tt stores the internal time at this moment. Then the
//command timer-tt; returns the total time used until this moment.

int tt = timer;

int k=5; //defines an integer k and set it to 5

//This creates a polynomial ring over the rational numbers with
//variables a(0),...,a(k),b(0),...,b(k),c(1),...,c(k),d(1),...,d(4),rho0,w,z. 
//w denotes the complex conjugate of the variable z. We can treat z and w as
//independent variables as they are linear independent over 
//the real and hence over the rational numbers.

ring R = 0,(a(1..k),b(1..k),c(1..k),d(1..k),rho0,w,z),dp;

//This defines the polynomials E, F, Phi, Alpha, where F, Alphas
//is conjugate to E, Phi, respectively.

//Remark: the factorial coefficients are written explicitly.

//It follows from Eq. 2.2 that a(1)=b(1)=0 or a(2)=b(2)=0.
//In order to examine second possibility comment out below
//and uncomment the block marked with *.

poly E;

E=a(2)*w+1/2*a(3)*z^2+1/2*a(4)*w^2+a(5)*z*w;

poly F;

F=b(2)*z+1/2*b(3)*w^2+1/2*b(4)*z^2+b(5)*w*z;

//*
//poly E;

//E=a(1)*z+1/2*a(3)*z^2+1/2*a(4)*w^2+a(5)*z*w;

//poly F;

//F=b(1)*w+1/2*b(3)*w^2+1/2*b(4)*z^2+b(5)*w*z;
//*

poly Phi;

Phi=c(1)*z+c(2)*w+1/2*c(3)*z^2+1/2*c(4)*w^2+c(5)*z*w;

poly Alpha;

Alpha=d(1)*w+d(2)*z+1/2*d(3)*w^2+1/2*d(4)*z^2+d(5)*w*z;

//The Ernst-Maxwell equations (difference between left hand side 
//and right hand side). The letter "c" refers to conjugation
//(for example EM1c is conjugated to EM1).

poly EM1 = (E+F+2*Alpha*Phi)*(2*(2*rho0+z+w)*diff(diff(E,w),z)+diff(E,z)+diff(E,w))-4*(2*rho0+z+w)*(diff(E,z)*diff(E,w)+Alpha*(diff(E,z)*diff(Phi,w)+diff(E,w)*diff(Phi,z)));

poly EM1c = (E+F+2*Alpha*Phi)*(2*(2*rho0+z+w)*diff(diff(F,z),w)+diff(F,w)+diff(F,z))-4*(2*rho0+z+w)*(diff(F,w)*diff(F,z)+Phi*(diff(F,w)*diff(Alpha,z)+diff(F,z)*diff(Alpha,w)));

poly EM2 = (E+F+2*Alpha*Phi)*(2*(2*rho0+z+w)*diff(diff(Phi,w),z)+diff(Phi,z)+diff(Phi,w))-2*(2*rho0+z+w)*(diff(E,z)*diff(Phi,w)+diff(E,w)*diff(Phi,z)+4*Alpha*diff(Phi,z)*diff(Phi,w));

poly EM2c = (E+F+2*Alpha*Phi)*(2*(2*rho0+z+w)*diff(diff(Alpha,z),w)+diff(Alpha,w)+diff(Alpha,z))-2*(2*rho0+z+w)*(diff(F,w)*diff(Alpha,z)+diff(F,z)*diff(Alpha,w)+4*Phi*diff(Alpha,w)*diff(Alpha,z));

//This defines an integer vector ("1,1" will refer to "z,w")

intvec v=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1; 

//The Ernst-Maxwell equations are truncated at second order in z, w

EM1=jet(EM1,1,v);
EM1c=jet(EM1c,1,v);
EM2=jet(EM2,1,v);
EM2c=jet(EM2c,1,v);
 
//The matrices store in the first row the monomials in z,w with non-zero
//coefficients and in the second row the coefficient of the corresponding 
//monomial.

matrix C1 = coef(EM1,zw);

matrix C1c = coef(EM1c,zw);

matrix C2 = coef(EM2,zw);

matrix C2c = coef(EM2c,zw);

//We define ideals

ideal Co = C1[2,1..2]; 

ideal C=Co;

//and collect the coefficients

Co = C1c[2,1..2];

C=C,Co;

Co = C2[2,1..3];

C=C,Co;

Co = C2c[2,1..3];

C=C,Co;

//We load the appropriate library 

 LIB "primdec.lib";

//We map into the new ring

ring S = 0,(a(1..k),b(1..k),c(1..k),d(1..k),rho0),dp;   //the new ring
ideal C = imap(R,C);               //we map C to the new ring
ideal gC = groebner(C);            //we compute a Groebner basis of C
"dimension:",dim(gC);              //dim(gC) returns the dimension of the
                                   //complex solution set for C=0

//and compute the minimal prime ideals

 list mpr = minAssGTZ(C);

//This shows the total time and the total memory used by Singular.

"";
"used time:",timer-tt, "sec";
"used memory:", memory(2) / (1024*1024), "MB";
"";

//mpr shows the list of minimal associated prime ideals. 

mpr;

//Check the result. The command size(reduce(C,std(mpr[i])));
//returns 0 if the solutions given the prime ideal mpr[i] are
//also solutions of the original system given by the ideal C.

int l,i;
for (i=1; i<=size(mpr); i++) {
    l= l+size(reduce(C,std(mpr[i])));
}
l;     
//This loop checks the result. The integer l is 0 if and only
//if the solutions given by the prime ideals mpr[i] are also
//solutions of the original system given by the ideal C=0.