docs: private void hqr2() {
      // FIXME: does this fail on NaN/inf values?
  
      // This is derived from the Algol procedure hqr2,
      // by Martin and Wilkinson, Handbook for Auto. Comp.,
      // Vol.ii-Linear Algebra, and the corresponding
      // Fortran subroutine in EISPACK.
  
      // Initialize
      final int nn = this.n, low = 0, high = nn - 1;
  
      // Store roots isolated by balanc and compute matrix norm
      double norm = 0.;
      for(int i = 0; i < nn; i++) {
        // If eventually allowing low,high to be set, use this:
        // if(i < low || i > high) { d[i] = H[i][i]; e[i] = 0.; }
        for(int j = (i > 0 ? i - 1 : 0); j < nn; j++) {
          norm += abs(H[i][j]);
        }
      }
  
      // Outer loop over eigenvalue index
      {
        double exshift = 0.0;
        // double p = 0, q = 0, r = 0; // , s = 0, z = 0;
        int iter = 0;
        for(int n = nn - 1; n >= low;) {
          // Look for single small sub-diagonal element
          int l = n;
          for(; l > low; --l) {
            double s = abs(H[l - 1][l - 1]) + abs(H[l][l]);
            if(abs(H[l][l - 1]) < 0x1P-52 * (s == 0. ? norm : s)) {
              break;
            }
          }
  
          // Check for convergence
          if(l == n) {
            // One root found
            d[n] = (H[n][n] += exshift);
            e[n] = 0.;
            n--;
            iter = 0;
            continue;
          }
          final double[] Hn = H[n], Hnm1 = H[n - 1];
          if(l == n - 1) {
            // Two roots found
            double w = Hn[n - 1] * Hnm1[n];
            double p = (Hnm1[n - 1] - Hn[n]) * 0.5, q = p * p + w;
            double x = (Hn[n] += exshift), z = FastMath.sqrt(abs(q));
            Hnm1[n - 1] += exshift;
  
            if(q >= 0) {
              // Real pair
              z = (p >= 0) ? p + z : p - z;
              d[n - 1] = x + z;
              d[n] = (z != 0.) ? x - w / z : d[n - 1];
              e[n] = e[n - 1] = 0.;
              x = Hn[n - 1];
              double s = abs(x) + abs(z);
              p = x / s;
              q = z / s;
              double r = FastMath.hypot(p, q);
              p /= r;
              q /= r;
  
              // Row modification
              for(int j = n - 1; j < nn; j++) {
                double tmp = Hnm1[j];
                Hnm1[j] = q * tmp + p * Hn[j];
                Hn[j] = q * Hn[j] - p * tmp;
              }
  
              // Column modification
              for(int i = 0; i <= n; i++) {
                modifyQP(H[i], n, p, q);
              }
  
              // Accumulate transformations
              for(int i = low; i <= high; i++) {
                modifyQP(V[i], n, p, q);
              }
            }
            else {
              // Complex pair
              d[n] = d[n - 1] = x + p;
              e[n] = -(e[n - 1] = z);
            }
            n -= 2;
            iter = 0;
            continue;
          }
          // No convergence yet
  
          // Form shift
          double x = Hn[n], y = 0., w = 0.;
          y = Hnm1[n - 1];
          w = Hn[n - 1] * Hnm1[n];
  
          // Wilkinson's original ad hoc shift
          if(iter == 10) {
            exshift += x;
            for(int i = low; i <= n; i++) {
              H[i][i] -= x;
            }
            double s = abs(Hn[n - 1]) + abs(Hnm1[n - 2]);
            x = y = 0.75 * s;
            w = -0.4375 * s * s;
          }
  
          // MATLAB's new ad hoc shift
          if(iter == 30) {
            double s = (y - x) * 0.5;
            s = s * s + w;
            if(s > 0) {
              s = FastMath.sqrt(s);
              s = (y < x) ? -s : s;
              s = x - w / ((y - x) * 0.5 + s);
              for(int i = low; i <= n; i++) {
                H[i][i] -= s;
              }
              exshift += s;
              x = y = w = 0.964;
            }
          }
  
          ++iter; // (Could check iteration count here.)
  
          double p = 0., q = 0., r = 0.;
          // Look for two consecutive small sub-diagonal elements
          int m = n - 2;
          for(; m >= l; m--) {
            final double[] Hm = H[m], Hmp1 = H[m + 1];
            final double z = Hm[m], xz = x - z;
            p = (xz * (y - z) - w) / Hmp1[m] + Hm[m + 1];
            q = Hmp1[m + 1] - xz - y;
            r = H[m + 2][m + 1];
            double tmp = abs(p) + abs(q) + abs(r);
            p /= tmp;
            q /= tmp;
            r /= tmp;
            if(m == l || abs(Hm[m - 1]) * (abs(q) + abs(r)) < 0x1P-52 * (abs(p) * (abs(H[m - 1][m - 1]) + abs(z) + abs(Hmp1[m + 1])))) {
              break;
            }
          }
  
          for(int i = m + 2; i <= n; i++) {
            final double[] Hi = H[i];
            Hi[i - 2] = 0.;
            if(i > m + 2) {
              Hi[i - 3] = 0.;
            }
          }
  
          // Double QR step involving rows l:n and columns m:n
          for(int k = m, last = n - 1; k <= last; k++) {
            boolean notlast = (k != last);
            final double[] Hk = H[k], Hkp1 = H[k + 1],
                Hkp2 = notlast ? H[k + 2] : null;
            if(k != m) {
              p = Hk[k - 1];
              q = Hkp1[k - 1];
              r = notlast ? Hkp2[k - 1] : 0.;
              x = abs(p) + abs(q) + abs(r);
              if(x == 0.0) {
                continue; // Jama 1.0.3 fix
              }
              p /= x;
              q /= x;
              r /= x;
            }
            double s = FastMath.hypot(p, q, r);
            s = (p < 0) ? -s : s;
            if(s != 0) {
              if(k != m) {
                Hk[k - 1] = -s * x;
              }
              else if(l != m) {
                Hk[k - 1] = -Hk[k - 1];
              }
              p += s;
              x = p / s;
              y = q / s;
              double z = r / s;
              q /= p;
              r /= p;
  
              // Row modification
              for(int j = k; j < nn; j++) {
                double tmp = Hk[j] + q * Hkp1[j];
                if(notlast) {
                  tmp += r * Hkp2[j];
                  Hkp2[j] -= tmp * z;
                }
                Hk[j] -= tmp * x;
                Hkp1[j] -= tmp * y;
              }
  
              // Column modification
              for(int i = 0; i <= Math.min(n, k + 3); i++) {
                modifyQR(H[i], k, notlast, q, r, x, y, z);
              }
  
              // Accumulate transformations
              for(int i = low; i <= high; i++) {
                modifyQR(V[i], k, notlast, q, r, x, y, z);
              }
            } // (s != 0)
          } // k loop
          // check convergence
        } // while (n >= low)
      }
  
      if(norm == 0.0) {
        return;
      }
  
      // Backsubstitute to find vectors of upper triangular form
      for(int n = nn - 1; n >= 0; n--) {
        double p = d[n], q = e[n];
        if(q == 0) {
          hqr2BacksubstituteReal(n, p, norm);
        }
        else if(q < 0) {
          hqr2BacksubstituteComplex(n, p, q, norm);
        }
      }
  
      // Vectors of isolated roots
      // Only matters if the user can modify low, high:
      // for(int i = 0; i < nn; i++) {
      // if(i < low || i > high) { System.arraycopy(H[i], i, V[i], i, nn - i); }
      // }
  
      // Back transformation to get eigenvectors of original matrix
      hqr2BackTransformation(nn, low, high);
    }