It seems the method developed by J.J. Van Wijk,
and described in his paper, "Ray tracing objects defined by sweeping a sphere"
(Computers & Graphics, Volume 9, Issue 3, 1985, Pages 283-290) has become the standard approach.
The main idea is to avoid the exact 4th-order computation by performing a pre-test that is only 2nd-order, essentially against the shape produced by sweeping a square rather than a circle around a circle.
Then only for the rays that pierce this bounding shape do you solve the 4th-order equations.
I believe with effort you can find downloadable implementations of ray-torus intersection following
this method. Code written by Han-Wen Nienhuys is available here (but I didn't check to see
if he follows Van Wijk).