docs: def _vmf_normalize(kappa, dim):
      """Compute normalization constant using built-in numpy/scipy Bessel
      approximations.
  
      Works well on small kappa and mu.
      """
      num = np.power(kappa, dim / 2. - 1.)
  
      if dim / 2. - 1. < 1e-15:
          denom = np.power(2. * np.pi, dim / 2.) * i0(kappa)
      else:
          denom = np.power(2. * np.pi, dim / 2.) * iv(dim / 2. - 1., kappa)
  
      if np.isinf(num):
          raise ValueError("VMF scaling numerator was inf.")
  
      if np.isinf(denom):
          raise ValueError("VMF scaling denominator was inf.")
  
      if np.abs(denom) < 1e-15:
          raise ValueError("VMF scaling denominator was 0.")
  
      return num / denom