BooleanPolynomialSequence.eliminate_linear_variables leaves some linear equations in corner cases

[Bouillaguet]

sage: R.<x,y,z> = BooleanPolynomialRing()
sage: S = Sequence([x*y*z+x*y+z*y+x*z, x+y+z+1, x+y])
sage: S.eliminate_linear_variables(return_reductors=True)
[z + 1], [x + y + z + 1])

Obviously, the first part is linear, and it shouldn't. The result should have been:

( [], [z + 1, x + y] )

[Bouillaguet]

The patch repairs the function, along with the following changes:

• by default, the threshold for the size of the tail changes from 3 to +infinity (3 was completely arbitrary anyway)
• the linear reductors form a reduced groebner basis (i.e. they are in RREF).

[malb]

a bit nitpicky but I think it should read skip is None instead of skip==None. But I think other than that the patch looks good.

[Bouillaguet]

replaced == None by is None

[jdemeyer]

This should be rebased to #13847.

[Bouillaguet]

Done