LIB "primdec.lib";

// Example of Q[x, y, z]
ring r = 0, (x, y, z), lp;


ideal J = x^2, y^2;
// ideal J is not radical

poly f = x;
// poly f in radical of J
// but g is not in radical of J
poly g = 2x + 3y -4;

// This is how you do radical membership
ideal J1 = x^2, y^2, 1-z*f;
J1 = groebner(J1);
printf("GB of J1 = 1 iff f in radical(J):");
J1;
// GB(J1) should be {1} <=> f in radical(J)


ideal J2 = x^2, y^2, 1-z*g;
J2 = groebner(J2);
printf("GB of J2 should not be 1 as g is not in radical of J");
J2;
// J2 should not be {1}

// And this is how you could compute the radical of J
// Should be <x, y>
ideal J3 = radical(J);
printf("Radical of J is:");
J3;


// Now, apply this to Galois fields
// We know that J_0 is radical, and so is J + J_0, but J may not be radical

printf("Let us see what happens over GF(8)");
ring fq = (2, X), (x, y, z), lp;
minpoly = X^3 + X + 1;

ideal J0 = x8 - x, y8 - y, z8 - z;

ideal J = X*x^2 + (X+1)*y^3 -1, X^5*z^4 + z^3 + X^3*z*y^2 + y*z;

// J is a Groebner basis and so is J0


// J is not radical
ideal radJ = radical(J);
printf("Groebner basis of Ideal J is:");
groebner(J);
printf("Radical of J is:");
radJ;

// J0 is radical
ideal radJ0 = radical(J0);
printf("Radical of J0 is J0 itself:");
radJ0;

// J + J0 is radical
ideal JJ0 = J + J0;
ideal radJJ0 = radical(JJ0);
printf("J + J0 is radical");

printf("The Groebner basis of Radical of J+J0 is:");
groebner(radJJ0);

printf("Which is the same thing as GB of J + J0");
groebner(JJ0);