I found one solution where I can sort-of get it mathematically by getting the cross product of the vectors in the face, but I think I'm losing accuracy since its all in floats, and rounding the array doesn't really work for what I'm doing.

I'm doing because I need a way to compare the angle of each of the faces to tell if faces are 'parallel' to each other. (not sure if 'parallel' is the right word in this case)

This is the math script that's sort-of working for me

p1 = np.array(co_list[0])
p2 = np.array(co_list[1])
p3 = np.array(co_list[2])

v1 = p3 - p1
v1 = p3 - p1

normal = np.cross(v1, v2)

This is my test geometry

And this is an example of the 'parallel' faces I was talking about (if there's a better term for this please let me know.

I feel like the best way to do this is by getting the direction of the normal for each face and comparing them. Does anyone know of a way to extract the normal value of the face using the blender API?