ContourPlot3d quality very poor

Here's a roundabout way: Plot in spherical coordinates and transform back. There are some slight imperfections, esp. where the boundaries join.

sph = TransformedField[ "Cartesian" -> "Spherical", 
   singu, {x, y, z} -> {r, θ, ϕ}] // Simplify

(* coordinate and vector field (VertexNormals) transformations *)
cXF = CoordinateTransformData["Spherical" -> "Cartesian", "Mapping"];
vXF[{x_, y_, z_}, {a_, b_, c_}] = 
  TransformedField["Spherical" -> "Cartesian", {a, b, c}, {r, θ, ϕ} -> {x, y, z}];

(* spherical plot *)
cpSPH = ContourPlot3D[
   sph == 0, {r, 0, 1}, {θ, 0, Pi}, {ϕ, 0, 2 Pi}, 
   MeshFunctions -> (* cartesian mesh *)
    Thread[cXF /. HoldPattern[Slot[1][[n_]]] :> Slot[n]]];

(* transform back to cartesian *)
cpCAR = Show[
   cpSPH /. 
    GraphicsComplex[p_, g_, rest___] :> 
     With[{x = Transpose@cXF@Transpose@p},
      GraphicsComplex[
       x,
       g,
       VertexNormals -> Transpose[vXF[Transpose@x, Transpose[VertexNormals /. {rest}]]],
       rest
       ]],
   PlotRange -> 1]

SliceContourPlot3D is not as elegant as ContourPlot3D, but it does display the complexity of the surface, and in less than a minute. It is this complexity that gives ContourPlot3D problems unless PlotPoints is large, in which case the computation crashes, at least on my PC.

s = Collect[singu // Rationalize, Sqrt[_], Simplify];
Show[SliceContourPlot3D[s, {"ZStackedPlanes", {#}}, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, 
    PerformanceGoal -> "Quality", MaxRecursion -> 4, PlotPoints -> 50,
    Contours -> {0}, ContourShading -> Opacity[0], 
    ContourStyle -> Directive[Hue[(# + 1)/2], Thick]] & /@ 
    Range[-.75, .75, .15]]
SwatchLegend[Hue[(# + 1)/2] & /@ Range[-.75, .75, .15], Range[-.75, .75, .15], 
    LegendLayout -> "Row"]

Note that s was computed from signu using Collect with Simplify to improves speed. Multiple SliceContourPlot3D plots were generated and combined to give contours on each slice a different color.