Check the VirtuaLight exporter by Jano. There’s a function in it called calcNormal performing the stuff you need. It goes something like this:

```
def calcNormal(face, verts, n):
# If it's a quad, calculate normal from diagonal instead
# Assume verts are sorted in circular fasion (i.e. next to each other)
m = (3 * n) | 1 # when n=0 (tri), m=1; when n=1 (quad), m=3
vectorU = [verts[face[0]][0] - verts[face[2]][0],
verts[face[0]][1] - verts[face[2]][1],
verts[face[0]][2] - verts[face[2]][2]]
vectorV = [verts[face[n]][0] - verts[face[m]][0],
verts[face[n]][1] - verts[face[m]][1],
verts[face[n]][2] - verts[face[m]][2]]
normal = [(vectorU[1] * vectorV[2]) - (vectorU[2] * vectorV[1]),
(vectorU[2] * vectorV[0]) - (vectorU[0] * vectorV[2]),
(vectorU[0] * vectorV[1]) - (vectorU[1] * vectorV[0])]
nLength = sqrt((normal[0] ** 2) + (normal[1] ** 2) + (normal[2] ** 2))
# divide by zero check
if not nLength: return "stray"
else: return (normal[0]/nLength, normal[1]/nLength, normal[2]/nLength)
```