With all due respect:
f = (x^2 + y^2 - 1)^2 + (y^2 + z^2 - 1)^2 + (x^2 + z^2 - 1)^2;
FullSimplify[Reduce[f == 0, {x, y, z}, Reals]]
Reduce[x^2 + y^2 == 1 && z^2 + y^2 == 1 && x^2 + z^2 == 1, {x, y, z}]
Note as expected the last 2 results are equivalent. f=0, is a set of 8 points (vertices of a cube). This is separate issue for surfaces of f=n that can be explored by ContourPlot3D.
This can be seen in many ways:
ir = ImplicitRegion[
x^2 + y^2 == 1 && z^2 + y^2 == 1 && x^2 + z^2 == 1, {x, y, z}];
dr = DiscretizeRegion[ir]
or more instructively,
pts = Tuples[{-(1/Sqrt[2]), 1/Sqrt[2]}, 3];
Show[ParametricPlot3D[{{u, v, u^2 + v^2}, {u^2 + v^2, u, v}, {u,
u^2 + v^2, v}, {u, v, -u^2 - v^2}, {-u^2 - v^2, u,
v}, {u, -u^2 - v^2, v}}, {u, -2, 2}, {v, -2, 2}, Mesh -> False,
PlotStyle ->
Join[Table[Opacity[0.2], {6}], Table[Opacity[0.4], {6}]]],
Graphics3D[{Red, PointSize[0.02], Point[pts]}]]
