//
//          Symbolic-numerical polynomial solving
// -------------------------------------------------------
LIB "solve.lib";
option(redSB);    // option for computing always reduced GB's


// -------------------------------------------------------
//               A simple example
// -------------------------------------------------------
ring  r = 0,(x,y,z),lp;
ideal I = x2+y+z-1,
          x+y2+z-1,
          x+y+z2-1;


//=====       a lexicographical Groebner basis     =====
ideal j=std(I);
j;

//=====    computing the triangular sets           =====


// Step 1: compute fi_tilde 
ring  r1 = 0,(y,z),lp;
poly f2_t = y2-y-z2+z;
poly f3_t = 2yz2+z4-z2;
poly f4_t = z6-4z4+4z3-z2;


// Step 2: compute triangular sets for the fi_tilde
ideal j1 = f2_t, f3_t, f4_t;
j1 = std(j1);
j1;

// not triangular


// Recursion: Step 2a
ring  r2 = 0,(z),lp;
poly g2_t = 2z2;
poly g3_t = z6-4z4+4z3-z2;
std(g2_t,g3_t);

// is triangular: Triangular set for the fi_tilde's:
//  --> 1)   (z2, y2-y-z2+z)


// Recursion: Step 2b
setring r1;
ideal j2 = std( sat( j1, 2z2 )[1] );
j2;

// is triangular: Triangular set for the fi_tilde's:
//  --> 2)   (z4-4z2+4z-1, 2y+z2-1)

//
// g3_t is multiple of g2_t, hence: end of recursion
//
// Output:  Triangular sets for I:
//
// --> 1)   (z2, y2-y-z2+z, x+y+z2-1)
// --> 2)   (z4-4z2+4z-1, 2y+z2-1, x+y+z2-1)


// Step 3: all fi_t are already in j, hence: nothing to do!


//  Direct computation using SINGULAR command:
setring r;
list tr = triangMH(j);
tr;

// correspond to reduced GB's of the above, e.g.:
ideal it = z4-4z2+4z-1, 2y+z2-1, x+y+z2-1 ;
std(it);

// === apply numerical solving up to 10 digits ===
triang_solve(tr,10);   // use triangular sets
rlist;

basering;   // basering has changed to rC
 
kill rC;


// === apply algebraic solving ===
setring r;
tr;

// for first triangular set, obviously (1,0,0) and (0,1,0) solutions 
// (both with multiplicity 2), hence, consider the second one:
//

factorize(tr[2][1]);

// hence: (0,0,1) is solution with multiplicity 2
//
// other solutions are not rational!
// ==> goto field extension

ring r_ext = (32003,a),(x,y,z),lp;   // finite char since factorization over
                                      // Q(a) not yet implemented in Singular 
minpoly = a2+2a-1;

factorize(z4-4z2+4z-1);

// compute the x and y components:
(1-a2)/2;
(1-(a+2)^2)/2;

// hence: solutions are (a,a,a) and (-a-2,-a-2,-a-2)
//
// note that  a = -1+sqrt(2) 

kill r;


// -------------------------------------------------------
//                 The ISSAC Challenge
// -------------------------------------------------------

ring r = 0,(w,x,y,z),lp; 
ideal I = -2w2+9wx+8x2+9wy+9xy+6y2-7wz-3xz-7yz-6z2-4w+8x+4y+8z+2,
           3w2-5wx+4x2-3wy+2xy+9y2-6wz-2xz+6yz+7z2+9w+7x+5y+7z+5,
           7w2+5wx+2x2+3wy+9xy-4y2-5wz-7xz-5yz-4z2-5w+4x+6y-9z+2,
           8w2+5wx+5x2-4wy+2xy+7y2+2wz-7xz-8yz+7z2+3w-7x-7y-8z+8;

ideal J = std(I);
J;

// check whether the ideal is zero-dimensional:
degree(J);

// we compute a triangular system with Moeller's method (improved
// by Hillebrand)
//
printlevel = 6;   // just to see what happens during computation

list tr = triangMH(J,2);	


// tr is a list of 1 triangular system
//
// Let's look at the most important data: first the univariate polynomial
tr[1][1];


// then the leading monomials of tr[1][2], tr[1][3]:
tr[1][2][1];

tr[1][3][1];

// we solve the triangular system  using recursive solving of univariate
// polynomials (using Laguerre's method) up to 30 digits
//
triang_solve(tr,30);

size(rlist);

//Hence we have 16 different complex solutions 
//
// ------------------------------------------------

kill r;
printlevel=0;


// -------------------------------------------------------
//           Another famous example: "cyclic_5"
// -------------------------------------------------------

ring r = 0,(a,b,c,d,e),lp;
ideal I = a+b+c+d+e, 
          ab+ae+bc+cd+de, 
          abc+abe+ade+bcd+cde, 
          abcd+abce+abde+acde+bcde, 
          abcde-1; 

ideal J=groebner(I);

degree(J);

// we compute a triangular system with Moeller's method (improved
// by Hillebrand)
//
list tr = triangMH(J,2);	

size(tr);

// tr is a list of 12 triangular system of degrees, e.g.:
tr[1];

// we solve the triangular system using recursive solving of univariate
// polynomials (using Laguerre's method) up to 30 digits
//
triang_solve(tr,30);
size(rlist);

// Hence we have 70 different complex solutions, e.g.
rlist[3];

quit;
*/