Use R ::= Q[x[1..4]], DegRevLex;
S ::= Q[x[1..4]], Lex;

UnSet Indentation;

Define ByLT(S,T)   
  Return LT(S)>LT(T);   
EndDefine;

Define Print_GB(X)
  Y:=Untagged(X);   
  PrintLn("Reduced GB w.r.t. lex:");
  Foreach P In Y Do
    PrintLn(P);
    PrintLn;
  EndForeach;   
EndDefine;

//S:=[x[1]^2+x[2]+x[3]-1,x[1]+x[2]^2+x[3]-1,x[1]+x[2]+x[3]^2-1];
S:=[3x[1]^2+2x[2]x[3]-2x[1]x[4],2x[1]x[3]-2x[2]x[4],2x[1]x[2]-2x[3]-2x[3]x[4],x[1]^2+x[2]^2+x[3]^2-1];
I:=Ideal(S);

F:=[];
PrintLn;
For N:=1 To Len(x) Do
  L:=WithoutNth(x,N);
  J:=Gens(Elim(L,I));                    // Generators of elimination ideals
  PrintLn("Gen. of elim. ideal for x[N]: ",J[1]);
  Append(F,J[1]);
  F[N]:=F[N]/GCD(F[N],Der(F[N],x[N]));   // make the generators square-free
  PrintLn("Square-free version: ",F[N]);
  PrintLn;
EndFor;

L:=Concat(S,F);
I:=Ideal(L);                             // Now I is radical

TB:=QuotientBasis(I);
PrintLn("Complement of <LT(I)>: ", TB);
PrintLn("Number of solutions: ",Len(TB));

PrintLn("Degree of the gen. of last elim. ideal: ",Deg(F[Len(x)]));

C:=NewMat(Len(x),1);
For N:=1 To Len(x)-1 Do
  C[N,1]:=-Rand(-10,10)/Rand(1,1);
EndFor;
C[Len(x),1]:=1;
CT:=ConcatLists(List(Mat([x])*C));

PrintLn("Substitute for ",x[Len(x)],": ",CT[1]);

T:=Subst(L,x[Len(x)],CT[1]);             // Transformation of the ideal
J:=Ideal(T);
E:=Elim(WithoutNth(x,Len(x)),J);         // Eliminate all variables except of the last one
G:=Gens(E);
PrintLn("Gen. of elim. ideal for ",x[Len(x)],": ",G[1]);
PrintLn("Degree of the gen. of elim. ideal for ",x[Len(x)],": ",Deg(G[1]));

GD:=GBasis(J);                           // Speeding up by adding a Groebner basis w.r.t. DegRevLex
T:=Concat(T,GD);

Use S;                                   // Change to lex term order
T:=BringIn(T);
J:=Ideal(T);

If Deg(G[1])=Len(TB) Then                // Is the transformed ideal in normal position?
  
  RB:=ReducedGBasis(J);                  // Reduced Groebner basis w.r.t. lex order
  SortBy(RB,Function("ByLT"));
  RBT:=Tagged(RB,"GB");
  Print(RBT);
  
  RR:=RealRoots(RB[Len(RB)],10^(-4));    // Find the real roots of the last polynomial in the GB
  
  GJ:=[];
  For N:=1 To Len(RB)-1 Do
    Append(GJ,x[N]-RB[N]);               // Extract polynomials g_i containing only the last variable
  EndFor;

  For N:=1 To Len(RR) Do                 // for each root of the last polynomial in Groebner basis
    A:=(RR[N].Inf+RR[N].Sup)/2;
    B:=NewMat(1,Len(x)); B[1,Len(x)]:=A;
    For M:=1 To Len(x)-1 Do
      B[1,M]:=Subst(GJ[M],x[Len(x)],A);  // evaluate all the remaining polynomials in GB
    EndFor;
    BT:=B*C;                             // transform the last one: root-c1g1-...-c_{n-1}g_{n-1}
    
    PrintLn;
    PrintLn("Solution is:");
    Print("[");
    For M:=1 To Len(x)-1 Do              // Print all components of the solution corresponding to the root A
      STR:=DecimalStr(B[1,M]);
      Print(STR,",");
    EndFor;
    PrintLn(DecimalStr(BT[1,1]),"]");    
  EndFor;
Else
  PrintLn("Ideal isn't in normal position.");
EndIf;