How to export Bone matrices for animation?
After another few days of frustration, I finally was able to make it work. Basically, I gave up on exporting euler angles, head and tail, instead, I just exported the bone matrix itself like this:

mat = global_matrix * armature.matrix_world * bone.matrix_local


And in the game's code, I use this exported matrix for computing inverse and for relative to parent matrices (by multiplying bone's matrix by inverse of parent's).