I get this message when I run your code:
NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with
-((92291 x)/87992)-(121001 y)/175984 == 1.
It is omitted from your quoted result. Did you not get it? Obviously, NSolve
adds linear conditions until the system is not underdetermined. (I say, obviously, because that's what the message says.) With that in mind, Frank's response seems to be exactly what is done to find the solutions (at least the Gröbner basis part).
For the following equation, it adds two hyperplanes:
NSolve[x^2 + y^2 + z^2 - 1 == 0, {x, y, z}, Reals]
By the way, you can turn off that method and get parametrized solutions (at least for many algebraic equations):
SetSystemOptions["NSolveOptions" -> "UseSlicingHyperplanes" -> False]
NSolve[x^2 + y^2 + z^2 - 1 == 0, {x, y, z}, Reals]
NSolve::svars: Equations may not give solutions for all "solve" variables.
Out[] = {
{z -> ConditionalExpression[-1.` Sqrt[1.` - 1.` x^2 - 1.` y^2],
-1.` < x < 1.` && -1.` Sqrt[1.` - 1.` x^2] <= y <= Sqrt[1.` - 1.` x^2]]},
{z -> ConditionalExpression[Sqrt[1.` - 1.` x^2 - 1.` y^2],
-1.` < x < 1.` && -1.` Sqrt[1.` - 1.` x^2] <= y <= Sqrt[1.` - 1.` x^2]]},
{x -> -1.`, y -> 0, z -> 0},
{x -> 1.`, y -> 0, z -> 0}}