[pixelesque]

JFYI: Your inverse ray direction calculation is not NaN-safe: if rays are completely axis-parallel in one dimension, so the direction value is 0.0 for that axis, you'll be doing the val / 0.0 which results in a NaN...

Also, as you're using full double/f64-precision all the time, you're leaving a fair bit of performance on the table: transcendentals (sin(), cos(), etc) in particular - can be a lot slower than when using f32, and generally double precision can be special-cased to particular areas of the renderer that need it (curve, sphere intersection, and some situations where volume scattering produces very small distances).

[deliveryboyman]

What's the proper way to handle a zero in the direction vector when calculating the reciprocal direction? Should it evaluate to infinity?

[pixelesque]

Inverse is still 0.0 technically, but yes, there is a trick you can use with Inf and SIMD to mask them out, so Inf is sometimes used.

However, I'd just condition it for the moment.

so:

invDirX = dirX != 0.0 ? 1.0 / dirX : 0.0; etc, etc for each dimension.

Obviously doing the != 0.0 comparison is not great, as it suffers from potential issues again (especially if you have denormals), but you can generally get away with it I've found in most cases.