docs: def compute_triangle_circumcenters(X, ei_dot_ei, ei_dot_ej):
      """Computes the circumcenters of all given triangles.
      """
      # The input argument are the dot products
      #
      #   <e1, e2>
      #   <e2, e0>
      #   <e0, e1>
      #
      # of the edges
      #
      #   e0: x1->x2,
      #   e1: x2->x0,
      #   e2: x0->x1.
      #
      # Note that edge e_i is opposite of node i and the edges add up to 0.
  
      # The trilinear coordinates of the circumcenter are
      #
      #   cos(alpha0) : cos(alpha1) : cos(alpha2)
      #
      # where alpha_k is the angle at point k, opposite of edge k. The Cartesian
      # coordinates are (see
      # <https://en.wikipedia.org/wiki/Trilinear_coordinates#Between_Cartesian_and_trilinear_coordinates>)
      #
      #     C = sum_i ||e_i|| * cos(alpha_i)/beta * P_i
      #
      # with
      #
      #     beta = sum ||e_i||*cos(alpha_i)
      #
      # Incidentally, the cosines are
      #
      #    cos(alpha0) = <e1, e2> / ||e1|| / ||e2||,
      #
      # so in total
      #
      #    C = <e_0, e0> <e1, e2> / sum_i (<e_i, e_i> <e{i+1}, e{i+2}>) P0
      #      + ... P1
      #      + ... P2.
      #
      # An even nicer formula is given on
      # <https://en.wikipedia.org/wiki/Circumscribed_circle#Barycentric_coordinates>: The
      # barycentric coordinates of the circumcenter are
      #
      #   a^2 (b^2 + c^2 - a^2) : b^2 (c^2 + a^2 - b^2) : c^2 (a^2 + b^2 - c^2).
      #
      # This is only using the squared edge lengths, too!
      #
      alpha = ei_dot_ei * ei_dot_ej
      alpha_sum = alpha[0] + alpha[1] + alpha[2]
      beta = alpha / alpha_sum[None]
      a = X * beta[..., None]
      cc = a[0] + a[1] + a[2]
  
      # alpha = numpy.array([
      #     ei_dot_ei[0] * (ei_dot_ei[1] + ei_dot_ei[2] - ei_dot_ei[0]),
      #     ei_dot_ei[1] * (ei_dot_ei[2] + ei_dot_ei[0] - ei_dot_ei[1]),
      #     ei_dot_ei[2] * (ei_dot_ei[0] + ei_dot_ei[1] - ei_dot_ei[2]),
      # ])
      # alpha /= numpy.sum(alpha, axis=0)
      # cc = (X[0].T * alpha[0] + X[1].T * alpha[1] + X[2].T * alpha[2]).T
      return cc