Numerical bug in intersecting the equation of ray and torus when the camera is far from the torus

I am trying to ray trace a torus without triangulating the torus and just by intersecting the ray and torus analytic equation. I did that with the following code:

void circularTorusIntersectFunc(const CircularTorus* circularToruses, RTCRay& ray, size_t item)
{
  const CircularTorus& torus = circularToruses[item];

  Vec3fa O = ray.org /*- sphere.p*/;
  Vec3fa Dir = ray.dir;
  O.w = 1.0f;
  Dir.w = 0.0f;
  O = torus.inv_transform.mult(O);
  Dir = torus.inv_transform.mult(Dir);

  // r1: cross section of torus
  // r2: the ring's radius
  //  _____                     ____
  // / r1  \------->r2<--------/    \
  // \_____/                   \____/

  float r2 = sqr(torus.r1);
  float R2 = sqr(torus.r2);

  double a4 = sqr(dot(Dir, Dir));
  double a3 = 4 * dot(Dir, Dir) * dot(O, Dir);
  double a2 = 4 * sqr(dot(O, Dir)) + 2 * dot(Dir, Dir) * (dot(O, O) - r2 - R2) + 4 * R2 * sqr(Dir.z);
  double a1 = 4 * dot(O, Dir) * (dot(O, O) - r2 - R2) + 8 * R2 * O.z * Dir.z;
  double a0 = sqr(dot(O, O) - r2 - R2) + 4 * R2 * sqr(O.z) - 4 * R2 * r2;

  a3 /= a4; a2 /= a4; a1 /= a4; a0 /= a4;

  double roots[4];
  int n_real_roots;
  n_real_roots = SolveP4(roots, a3, a2, a1, a0);

  if (n_real_roots == 0) return;

  Vec3fa intersect_point;
  for (int i = 0; i < n_real_roots; i++)
  {
    float root = static_cast<float>(roots[i]);
    intersect_point = root * Dir + O;

    if ((ray.tnear <= root) && (root <= ray.tfar)) {

      ray.u = 0.0f;
      ray.v = 0.0f;
      ray.tfar = root;
      ray.geomID = torus.geomID;
      ray.primID = item;
      Vec3fa normal(
        4.0 * intersect_point.x * (sqr(intersect_point.x) + sqr(intersect_point.y) + sqr(intersect_point.z) - r2 - R2),
        4.0 * intersect_point.y * (sqr(intersect_point.x) + sqr(intersect_point.y) + sqr(intersect_point.z) - r2 - R2),
        4.0 * intersect_point.z * (sqr(intersect_point.x) + sqr(intersect_point.y) + sqr(intersect_point.z) - r2 - R2) + 8 * R2*intersect_point.z,
        0.0f
        );

      ray.Ng = normalize(torus.transform.mult(normal));
    }
  }
}

The code to solve the equation for SolveP4 function is taken from Solution of cubic and quatric functions.

The problem is when we are looking at the torus closely, it works pretty nice as follows:

But when I zoom out the camera, so camera is looking at the torus far from it, it suddenly gets so noisy and it is shape is not well identified. I tried to use more than 1 samples per pixels but still I have the same problem. It is as follows:

It seems I am facing a numerical problem but I dont know how to solve it. Anyone can help me with that?

Also, it is good to mention that I am raytracing the torus with Intel's Embree Lib.

Update (Single Color):

I think a lot of the problem is using single precision float rather than double precision.

Define two functions

double dsqr(double x) { return x*x; }

double ddot(const Vec3fa &a,Vec3fa &b) {
  double x1 = a.x, y1 = a.y, z1 = a.z;
  double x2 = b.x, y2 = b.y, z2 = b.z;
  return x1*x2 + y1*y2 + z1*z2;
}

to find the square and the dot product but using double precision. Change the calculations of r2 R2 a4 a3 a2 a1 and a0 to use these

double r2 = dsqr(torus.r1);
double R2 = dsqr(torus.r2);

double a4 = dsqr(ddot(Dir, Dir));
double a3 = 4 * ddot(Dir, Dir) * ddot(O, Dir);
double a2 = 4 * dsqr(ddot(O, Dir)) + 2 * ddot(Dir, Dir) * (ddot(O, O) - r2 - R2)
    + 4 * R2 * dsqr(Dir.z);
double a1 = 4 * ddot(O, Dir) * (ddot(O, O) - r2 - R2) + 8 * R2 * O.z * Dir.z;
double a0 = dsqr(ddot(O, O) - r2 - R2) + 4 * R2 * dsqr(O.z) - 4 * R2 * r2;

all the remaining code is the same. In my test this made a fuzzy looking image look perfectly sharp.