docs: private void hqr2BacksubstituteComplex(int n, double p, double q, double norm) {
      final double[] Hn = H[n], Hnm1 = H[n - 1];
      // Last vector component imaginary so matrix is triangular
      final double tmp = Hn[n - 1];
      if(abs(tmp) > abs(Hnm1[n])) {
        Hnm1[n - 1] = q / tmp;
        Hnm1[n] = -(Hn[n] - p) / tmp;
      }
      else {
        cdiv(0., -Hnm1[n], Hnm1[n - 1] - p, q, Hnm1, n - 1);
      }
      Hn[n - 1] = 0.;
      Hn[n] = 1.;
      double r = 0., s = 0., z = 0.;
      for(int i = n - 2, l = n - 1; i >= 0; i--) {
        final double[] Hi = H[i], Hip1 = H[i + 1];
        double w = Hi[i] - p, ra = 0., sa = 0.;
        for(int j = l; j <= n; j++) {
          final double Hij = Hi[j];
          final double[] Hj = H[j];
          ra += Hij * Hj[n - 1];
          sa += Hij * Hj[n];
        }
  
        if(e[i] < 0.) {
          z = w;
          r = ra;
          s = sa;
          continue;
        }
        l = i;
        if(!(e[i] > 0.0)) {
          cdiv(-ra, -sa, w, q, Hi, n - 1);
        }
        else {
          // Solve complex equations
          final double x = Hi[i + 1], y = Hip1[i], dip = d[i] - p;
          double vr = dip * dip + e[i] * e[i] - q * q, vi = dip * 2. * q;
          if(vr == 0. && vi == 0.) {
            vr = 0x1P-52 * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z));
          }
          cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi, Hi, n - 1);
          if(abs(x) > (abs(z) + abs(q))) {
            Hip1[n - 1] = (-ra - w * Hi[n - 1] + q * Hi[n]) / x;
            Hip1[n] = (-sa - w * Hi[n] - q * Hi[n - 1]) / x;
          }
          else {
            cdiv(-r - y * Hi[n - 1], -s - y * Hi[n], z, q, Hip1, n - 1);
          }
        }
  
        // Overflow control
        double t = max(abs(Hi[n - 1]), abs(Hi[n]));
        if((0x1P-52 * t) * t > 1) {
          for(int j = i; j <= n; j++) {
            final double[] Hj = H[j];
            Hj[n - 1] /= t;
            Hj[n] /= t;
          }
        }
      }
    }