//SINGULAR Example4.3.5

//===================== we need ================================
option(redSB);            //to obtain a reduced normal form 
ring R=0,(x,y),lp;            
ideal a1=x;               //preparation of the example 
ideal a2=y2+2y+1,x-1;
ideal a3=y2-2y+1,x-1;
ideal I=intersect(a1,a2,a3);
poly h=y4-2y2+1;       //the lcm of the leading coefficients
ideal I1=quotient(I,h);

//==============================================================

ideal I2=std(I+ideal(h));
//we compute now the decomposition of I2
indepSet(I2);

list fac=factorize(I2[1]);
fac;

ideal J1=std(I2,(y+1)^2);        //the two candidates 
ideal J2=std(I2,(y-1)^2);        //for primary ideals

J1;
J2;

map phi=R,x,x+y;
map psi=R,x,-x+y;               //and the inverse map

ideal K1=std(phi(J1));
ideal K2=std(phi(J2));
factorize(K1[1]);

ideal K11=std(K1,(y+1)^2);        //the new candidates 
                                  //for primary ideals
ideal K12=std(K1,(y+2)^2);        //coming from K1
factorize(K2[1]);

ideal K21=std(K2,(y-1)^2);         //the new candidates
                                   // for primary ideals
ideal K22=std(K2,y2);              //coming from K2
K11=std(psi(K11));                 //the inverse coordinate
                                   //transformation
K12=std(psi(K12));
K21=std(psi(K21));
K22=std(psi(K22));
K11;                               //the result
K12;
K21;
K22;