Accurate ray casting

I have been experimenting with ray_cast definition that Blender’s Python API offers. Specifically, I am using definition find in below for casting rays towards an object from a given point.

def intersect_a_ray(p0,p1,ob):
    """
    Definition to create a cylinder in between two points with a certain radius.

    Parameters
    ----------
    p0             : list
                     Starting point of a ray.
    p1             : list
                     Ending point of a ray.
    ob             : blender object
                     Object to be tested for intersection with a ray.

    Returns
    ----------
    result         : list
                     Result of ray intersection test containing state, coalision ray, surface normal.
    loc            : bool/float
                     Location of the intersection. If ray isn't intersecting returns False.
    dist           : float
                     Distance from the starting point of the ray to the intersection point.
    """
    mwi             = ob.matrix_world.inverted()
    ray_begin       = mwi @ Vector((p0[0],p0[1],p0[2]))
    ray_end         = mwi @ Vector((p1[0],p1[1],p1[2]))
    ray_direction   = (ray_end-ray_begin).normalized()
    result          = ob.ray_cast(origin=ray_begin,direction=ray_direction)
    loc             = False
    dist            = 0
    if result[0] == True:
        mw   = ob.matrix_world
        loc  = mw @ result[1]
        dist = ((loc[0]-p0[0])**2+(loc[1]-p0[1])**2+(loc[2]-p0[2])**2)**0.5
    return result,loc,dist

In my observation using either Blender 2.80 or 2.81 or 2.82, the accuracy of ray intersection point with a given object is questionable for the application that I am targeting.

In my test-bed, I have a mesh that represents a flat surface standing at (0,0,0) in the reference coordinate system. By definition, I know that the intersection is at some point at X,Y,0. If I shoot a ray using the definition given at above, the returned location value with loc in intersect_a_ray is sometimes off by orders of 10^-6.

I am speculate that this is a common behavior. If not, may I kindly ask if you see any bugs in intersect_a_ray function?

If it is a common behavior, is there a way to get higher precision intersection points, probably being off by the orders of 10^-7 is good for what I am aiming for.

I think that is intrinsic to the limits of float values… i would not suspect a bug, as the ray triangle algorithms are simple and tested.
what you could do:
As i said triangle intersection is easy, finding the right triangle is the hard part - acceleration structures etc. I never used the raycast (in blender), but if you can get the triangle info, you can write you own little function that does the last step - the intersection in double precision. So doubles instead of floats.

just an idea…

edit: but the meshes are float anyway, so not sure if that gets you anywhere… eg.: smallest number bigger than one is 1.0000001192