Start be reading about the Python API for BVH trees.
You would want to build one using FromBMesh.
Then loop through all the faces f of your mesh. For each, find a point p in it (the centroid, say, though that is not exactly right for some weirdly shaped faces). The use the bvhtree’s raycast function using p as origin and trying several directions as I suggested above. I haven’t tried this, but it looks like it will return a tuple of hits on faces with the index being the index of the face hit (I hope - you should try to see). Using that you can either just count hits, or if you want to do the parity calculation, you need to know how the ray direction compares with the normal direction of the hit face (maybe the normal result in the raycast result vector is the normal of the hit face, which would make this easy).
P.S., the reason I suggest using directions that are a little off from (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) is that the way models are typically built, it is often quite likely that rays in those directions will hit the seams between faces, which may result in those hits being reported for both faces or neither (I’m not sure about the latter - depends on whether or not the internal code is using “watertight” ray casting).