Dumb mistake: The mesh has a subdiv modifier which increases the number of verts.
The function closest_point_on_mesh should return the face number in its last argument, but apparently this can be larger than the total number of polygons. How can this be?
The following script computes the closest point on the mesh itself. It works as expected on primitive meshes (the max face number is equal to or slightly smaller than the total number of faces), but for a Genesis 8 Female mesh it is much bigger (download).
The test script only uses a single mesh, but the problem arose when I had two meshes.
import bpy
ob = bpy.context.object
data = [ob.closest_point_on_mesh(v.co) for v in ob.data.vertices]
fnums = [idx for res,loc,nor,idx in data if res]
print("Max fnum = %d" % max(fnums))
print("Total fnums = %d" % len(ob.data.polygons))
Output for G8F:
Max fnum = 65468
Total fnums = 16368
One may note that 4*16368 = 65472. The mesh is pure quad.