docs: def _qr_factor_packed(a, enorm, finfo):
      """Compute the packed pivoting Q-R factorization of a matrix.
  
  Parameters:
  a     - An n-by-m matrix, m >= n. This will be *overwritten*
          by this function as described below!
  enorm - A Euclidian-norm-computing function.
  finfo - A Numpy finfo object.
  
  Returns:
  pmut   - An n-element permutation vector
  rdiag  - An n-element vector of the diagonal of R
  acnorm - An n-element vector of the norms of the rows
           of the input matrix 'a'.
  
  Computes the transposed Q-R factorization of the matrix 'a', with
  pivoting, in a packed form, in-place. The packed information can be
  used to construct matrices Q and R such that
  
    A P = R Q or, in Python,
    np.dot(r, q) = a[pmut]
  
  where q is m-by-m and q q^T = ident and r is n-by-m and is lower
  triangular. The function _qr_factor_full can compute these
  matrices. The packed form of output is all that is used by the main LM
  fitting algorithm.
  
  "Pivoting" refers to permuting the rows of 'a' to have their norms in
  nonincreasing order. The return value 'pmut' maps the unpermuted rows
  of 'a' to permuted rows. That is, the norms of the rows of a[pmut] are
  in nonincreasing order.
  
  The parameter 'a' is overwritten by this function. Its new value
  should still be interpreted as an n-by-m array. It comes in two
  parts. Its strict lower triangular part contains the struct lower
  triangular part of R. (The diagonal of R is returned in 'rdiag' and
  the strict upper trapezoidal part of R is zero.) The upper trapezoidal
  part of 'a' contains Q as factorized into a series of Householder
  transformation vectors. Q can be reconstructed as the matrix product
  of n Householder matrices, where the i'th Householder matrix is
  defined by
  
  H_i = I - 2 (v^T v) / (v v^T)
  
  where 'v' is the pmut[i]'th row of 'a' with its strict lower
  triangular part set to zero. See _qr_factor_full for more information.
  
  'rdiag' contains the diagonal part of the R matrix, taking into
  account the permutation of 'a'. The strict lower triangular part of R
  is stored in 'a' *with permutation*, so that the i'th row of R has
  rdiag[i] as its diagonal and a[pmut[i],:i] as its upper part. See
  _qr_factor_full for more information.
  
  'acnorm' contains the norms of the rows of the original input
  matrix 'a' without permutation.
  
  The form of this transformation and the method of pivoting first
  appeared in Linpack."""
  
      machep = finfo.eps
      n, m = a.shape
  
      if m < n:
          raise ValueError('"a" must be at least as tall as it is wide')
  
      acnorm = np.empty(n, finfo.dtype)
      for j in range(n):
          acnorm[j] = enorm(a[j], finfo)
  
      rdiag = acnorm.copy()
      wa = acnorm.copy()
      pmut = np.arange(n)
  
      for i in range(n):
          # Find the row of a with the i'th largest norm, and note it in
          # the pivot vector.
  
          kmax = rdiag[i:].argmax() + i
  
          if kmax != i:
              temp = pmut[i]
              pmut[i] = pmut[kmax]
              pmut[kmax] = temp
  
              rdiag[kmax] = rdiag[i]
              wa[kmax] = wa[i]
  
              temp = a[i].copy()
              a[i] = a[kmax]
              a[kmax] = temp
  
          # Compute the Householder transformation to reduce the i'th
          # row of A to a multiple of the i'th unit vector.
  
          ainorm = enorm(a[i,i:], finfo)
  
          if ainorm == 0:
              rdiag[i] = 0
              continue
  
          if a[i,i] < 0:
              # Doing this apparently improves FP precision somehow.
              ainorm = -ainorm
  
          a[i,i:] /= ainorm
          a[i,i] += 1
  
          # Apply the transformation to the remaining rows and update
          # the norms.
  
          for j in range(i + 1, n):
              a[j,i:] -= a[i,i:] * np.dot(a[i,i:], a[j,i:]) / a[i,i]
  
              if rdiag[j] != 0:
                  rdiag[j] *= np.sqrt(max(1 - (a[j,i] / rdiag[j])**2, 0))
  
                  if 0.05 * (rdiag[j] / wa[j])**2 <= machep:
                      # What does this do???
                      wa[j] = rdiag[j] = enorm(a[j,i+1:], finfo)
  
          rdiag[i] = -ainorm
  
      return pmut, rdiag, acnorm