Frustum culling: Sphere-Cone test with no squareroot

I crossposted this article at GameDev.

I've found a method for checking if a sphere is partially or totally included inside a cone without using a square root, as the square root has a big overhead on some FPUs, especially on mobile platforms (e.g., on the VFP11 and NEON units of the armv6 and armv7 processors).

Original algorithm

I started with Charles Bloom's algorithm as presented in the second part of this document:

V = sphere.center - cone.apex_location
a = dotProduct(V, cone.direction_normal)
b = a * cone.tan
c = sqrt( dotProduct(V,V) - a² )
d = c - b
e = d * cone.cos

now, if ( e >= sphere.radius ):
    cull the sphere (return 0)
else, if ( e <=-sphere.radius ):
    totally include the sphere (return 1)
else:
    the sphere is partially included (return -1)

What's going on in this cone-sphere test? Basically,
I'm trying to find 'e' which is the shortest distance
from the center of the sphere to the surface of the cone.
You can draw some pictures and see what's going on.
'a' is how far along the ray of the cone the sphere's
center is. 'b' is the radius of the cone at 'a'. 'c' is
the distance from the center of the sphere to the axis
of the cone, and 'd' is the distance from the center of
the sphere to the surface of the cone, along a line
perpendicular to the axis of the cone (which is not the
closest distance).

Simplified form

The algorithm can be simplified and turned into the following:

V = sphere.center - cone.apex_location
a = dotProduct(V, cone.direction_normal)
x = cone.cos * sqrt(dotProduct(V,V) - a²) - a*cone.sin
if abs(x)>sphere.radius:
    if x<0, totally include the sphere (return 1)
    else,   cull the sphere (return 0)
else the sphere is partially included (return -1)

sqrt() must die!

Using squared inequalities, this in turn leads to the following sqrt-free formulation:

V = sphere.center - cone.apex_location
a = dotProduct(V, cone.direction_normal)
p = a*cone_sin
q = cone_cos² * dotProduct(V, V) - a²
r = q - sphere_radius²
if (p<sphere_radius) || (q>0):
    if (r < 2 * sphere_radius * p):
        the sphere is partially included (return -1)
    else if q<0:
        the sphere is totally included (return 1)
    else:
        cull the sphere (return 0)
else:
    if ( -r < 2 * sphere_radius * p):
        the sphere is partially included (return -1)
    else if q<0:
        the sphere is totally included (return 1)
    else:
        cull the sphere (return 0)

Source code

Here is the source of the two methods: Charles Bloom's and mine.