Computing Groebner basis for a complicated systems of polynomials I am trying to solve complicated systems of polynomial equations. The first step is to determine maximal sets of independent variables for the solution manifold (ideal) or the number of isolated solutions using Gröbner bases. In some cases, Mathematica, Maple and SymPy do not seem to be able to determine the Gröbner basis (in reasonable time).
Do you recommend some tools for this task?
Let me give a test case for which Mathematica, Maple and SymPy did not finish in reasonable time (~<10h) on my computer. I used grevlex monomial order. Do you know a tool that can compute the Gröbner basis for this case? There are 21 polynomials and 8 variables u, q1, r1, q2, r2, q3, r3, q4 (using Gram-Schmidt, one can reduce to 8 polynomials):
https://pastebin.com/mpqZUQqC
 A: I suggest you first try to work modulo a machine sized prime, e.g. p=2^31-1 for Maple 2016 or newer, or p=65521 for older versions of Maple.  Maple can compute a grevlex basis mod p in about a minute, and you can count solutions assuming the system behaves generically for that prime.  E.g.
infolevel[Groebner] := 4:  # see what is going on
Groebner:-IsZeroDimensional(sys,characteristic=2^31-1);  # false?

A: As you have 8 polynomials in 8 variables then you can generically expect that the number of complex solutions is finite.  In Maple, for example, you can try to find them using fsolve() with random initial data.  It is probably best to use the default precision initially, then you can set Digits to 100 or so and then refine the solutions.  Then you can use identify() to see if the solutions are rational, or algebraic of low degree, or have some other nice form that Maple can understand (if that is not already obvious).  Once you have a good idea of what the solutions are, and any symmetries, then you can use that to guide a more theoretical analysis of the equations.
