import math


def product_decompose_rotation(rot):
  """Decomposes a rotation matrix into Davenport chained rotations."""
  dim = rot.shape[0]
  factors = []
  work_rot = rot.copy()
  for axis0 in range(dim):
    for axis1 in sorted(
        range(1 + axis0, dim),
        # Resolve largest-to-be-zeroed-out
        # index first.
        key=lambda n: -abs(work_rot[axis0, n])):  # pylint:disable=cell-var-from-loop
      ea0 = work_rot[axis0, axis0]
      ea1 = work_rot[axis1, axis0]
      angle = math.atan2(ea1, ea0)
      if angle != 0:
        factors.append((axis0, axis1, angle))
      ca = math.cos(angle)
      sa = math.sin(angle)
      row0 = work_rot[axis0].copy()
      row1 = work_rot[axis1].copy()
      work_rot[axis0, :] = +ca * row0 + sa * row1
      work_rot[axis1, :] = -sa * row0 + ca * row1
  return list(reversed(factors))