The sor.cpp changes being checked into release 20 of yuqk:

...
    // Transform the ray into the surface of revolution space.

    MInvTransPoint(P, ray.Origin, Trans);

    MInvTransDirection(D, ray.Direction, Trans);

    // Test if ray misses object's bounds.

    // Note! When global general bounding is on, the initial bounding tests
    // below are less effective, but still save time.

    // Credit: Bill Walker (a-d) code comments for conditional just below and
    // corrections to overall cylinder intersection just after (Oct 2025).

    // (a) Checks if the ray direction is flat or points up, and the ray origin is
    // higher than the highest sor control point (Height2 = ymax) So it will never
    // intersect the rotated sor spline

    // (b) Checks if the ray direction is flat or points down, and the ray origin
    // is lower than the lowest sor control point (Height1 = ymin) So it will never
    // intersect the rotated sor spline

    // (c) Checks if the ray direction is straight or points right, and the ray
    // origin is to the right of the rightmost sor control point (Radius2 = xmax)
    // So it will never intersect the rotated sor spline

    // (d) Checks if the ray direction is straight or points left, and the ray
    // origin is to the left of the leftmost (rotated) sor control point (-Radius2
    // = -xmax) So it will never intersect the rotated sor spline

#ifdef SOR_EXTRA_STATS
    Thread->Stats()[Sor_Bound_Tests]++;
#endif

    if (((D[Y] >= 0.0) && (P[Y] >  Height2)) ||  // (a)
        ((D[Y] <= 0.0) && (P[Y] <  Height1)) ||  // (b)
        ((D[X] >= 0.0) && (P[X] >  Radius2)) ||  // (c)
        ((D[X] <= 0.0) && (P[X] < -Radius2)))    // (d)
    {
        return(false);
    }

    // Get distance of the ray from rotation axis (= y axis).
    // (Normalization of D vector is not required due r0 /= sqrt(a) below)

    r0 = P[X] * D[Z] - P[Z] * D[X];

    // Test if ray misses object's bounds.

                                        // Insure sqrt() never effectively 0.0.
    if ((a = D[X] * D[X] + D[Z] * D[Z]) > gkMinIsectDepthReturned)
    {
        r0 /= std::sqrt(a); // Actual +- distance to y axis.
        if (fabs(r0) > Radius2)
        {
            return(false);
        }
    }

    // We've delayed the normalization of D (working space ray direction vector)
    // until we need it. Avoiding the cost until initial bounding checks above
    // completed.

    len = D.length();

    D /= len;

    // Test base/cap plane.

...

        // Restored old code due delaying normalization of D. (old use of common 'a' calculation is gone)
        x[1] = D[Y] * D[Y] * (3.0 * Entry->A * P[Y] + Entry->B) - D[X] * D[X] - D[Z] * D[Z];

...