<pre>
--------- CPU symmdq
Symmetric minimum degree algorithm by quotient graph

</pre>
425046
public static int cusolverSpXcsrsymmdqHost(
        cusolverSpHandle handle, 
        int n, 
        int nnzA, 
        cusparseMatDescr descrA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer p)
    {
        return checkResult(cusolverSpXcsrsymmdqHostNative(handle, n, nnzA, descrA, csrRowPtrA, csrColIndA, p));
    }
571382
private void bfs(int v, int[] cc, int id) {
        cc[v] = id;
        Queue<Integer> queue = new LinkedList<>();
        queue.offer(v);
        while (!queue.isEmpty()) {
            int t = queue.poll();
            for (Edge edge : graph[t]) {
                int i = edge.v2;
                if (!digraph && i == t) {
                    i = edge.v1;
                }

                if (cc[i] == -1) {
                    queue.offer(i);
                    cc[i] = id;
                }
            }
        }
    }
247389
def calculate_epsilon_effective(rho, epsilon, seed_degree, neighbor_degrees, mean_degree):
    
    # Calculate a weighted neighborhood degree average.
    # neighborhood_degree = rho*seed_degree + (1-rho)*neighbor_degrees.mean()
    neighborhood_degree = neighbor_degrees.mean()

    # Calculate the seed neighborhood normalized effective epsilon.
    epsilon_effective = (epsilon*np.log(1 + seed_degree))/np.log(1 + neighborhood_degree)

    # Calculate the maximum epsilon for at least one push on a neighboring node.
    # Also the minimum epsilon for a push on all the neighboring nodes.
    epsilon_effective_maximum = np.max(1/(seed_degree*neighbor_degrees))
    epsilon_effective_minimum = np.min(1/(seed_degree*neighbor_degrees))

    # print(epsilon_effective, epsilon_effective_maximum, epsilon_effective_minimum)

    # The maximum epsilon is absolute, whereas we regularize for the minimum.
    if epsilon_effective > epsilon_effective_maximum:
        epsilon_effective = epsilon_effective_maximum
    elif epsilon_effective < epsilon_effective_minimum:
        epsilon_effective = (epsilon_effective_minimum + epsilon_effective)/2

    return epsilon_effective
1768
def incremental_neighbor_graph(X, precomputed=False, k=None, epsilon=None,
                               weighting='none'):
  ''''''
  assert ((k is not None) or (epsilon is not None)
          ), "Must provide `k` or `epsilon`"
  assert (_issequence(k) ^ _issequence(epsilon)
          ), "Exactly one of `k` or `epsilon` must be a sequence."
  assert weighting in ('binary','none'), "Invalid weighting param: " + weighting
  is_weighted = weighting == 'none'

  if precomputed:
    D = X
  else:
    D = pairwise_distances(X, metric='euclidean')
  # pre-sort for efficiency
  order = np.argsort(D)[:,1:]

  if k is None:
    k = D.shape[0]

  # generate the sequence of graphs
  # TODO: convert the core of these loops to Cython for speed
  W = np.zeros_like(D)
  I = np.arange(D.shape[0])
  if _issequence(k):
    # varied k, fixed epsilon
    if epsilon is not None:
      D[D > epsilon] = 0
    old_k = 0
    for new_k in k:
      idx = order[:, old_k:new_k]
      dist = D[I, idx.T]
      W[I, idx.T] = dist if is_weighted else 1
      yield Graph.from_adj_matrix(W)
      old_k = new_k
  else:
    # varied epsilon, fixed k
    idx = order[:,:k]
    dist = D[I, idx.T].T
    old_i = np.zeros(D.shape[0], dtype=int)
    for eps in epsilon:
      for i, row in enumerate(dist):
        oi = old_i[i]
        ni = oi + np.searchsorted(row[oi:], eps)
        rr = row[oi:ni]
        W[i, idx[i,oi:ni]] = rr if is_weighted else 1
        old_i[i] = ni
      yield Graph.from_adj_matrix(W)
459354
def permute_graph(G, order):
  '''
  '''
  adj = G.matrix('dense')
  adj = adj[np.ix_(order, order)]
  return Graph.from_adj_matrix(adj)
1247897
public static void symmLowerToFull( DMatrixRMaj A )
    {
        if( A.numRows != A.numCols )
            throw new MatrixDimensionException("Must be a square matrix");

        final int cols = A.numCols;

        for (int row = 0; row < A.numRows; row++) {
            for (int col = row+1; col < cols; col++) {
                A.data[row*cols+col] = A.data[col*cols+row];
            }
        }
    }
383279
def graph_simplify(G):
    
    spurs = []
    path_nodes = []
    for k, d in G.degree().iteritems():
        if d == 1:
            spurs.append(k)
        elif d == 2:
            path_nodes.append(k)

    logging.debug("Remove {0} spurs.".format(len(spurs)))
    G.remove_nodes_from(spurs)

    SG = G.subgraph(path_nodes)
    cc = nx.connected_components(SG)
    for c in cc:
        if len(c) == 1:
            continue
        c = set(c)
        neighbors = set()
        for x in c:
            neighbors |= set(G.neighbors(x))
        neighbors -= c
        newtag = list(c)[0] + "*"
        for n in neighbors:
            G.add_edge(newtag, n)
        G.remove_nodes_from(c)
    logging.debug("Contract {0} path nodes into {1} nodes.".\
                    format(len(path_nodes), len(cc)))
729125
private void dfs(int v, int[] cc, int id) {
        cc[v] = id;
        for (Edge edge : graph[v]) {
            int t = edge.v2;
            if (!digraph && t == v) {
                t = edge.v1;
            }

            if (cc[t] == -1) {
                dfs(t, cc, id);
            }
        }
    }
817087
public static void symmLowerToFull( DMatrixSparseCSC A , DMatrixSparseCSC B , @Nullable IGrowArray gw )
    {
        ImplCommonOps_DSCC.symmLowerToFull(A, B, gw);
    }
1367122
def parametrized_bottleneck(x, hparams):
  """"""
  if hparams.bottleneck_kind == "tanh_discrete":
    d, _ = tanh_discrete_bottleneck(
        x, hparams.bottleneck_bits, hparams.bottleneck_noise * 0.5,
        hparams.discretize_warmup_steps, hparams.mode)
    return d, 0.0
  if hparams.bottleneck_kind == "isemhash":
    return isemhash_bottleneck(
        x, hparams.bottleneck_bits, hparams.bottleneck_noise * 0.5,
        hparams.discretize_warmup_steps, hparams.mode,
        hparams.isemhash_noise_dev, hparams.isemhash_mix_prob)
  if hparams.bottleneck_kind == "vq":
    return vq_discrete_bottleneck(x, hparams.bottleneck_bits, hparams.vq_beta,
                                  hparams.vq_decay, hparams.vq_epsilon)
  if hparams.bottleneck_kind == "em":
    return vq_discrete_bottleneck(
        x,
        hparams.bottleneck_bits,
        hparams.vq_beta,
        hparams.vq_decay,
        hparams.vq_epsilon,
        soft_em=True,
        num_samples=hparams.vq_num_samples)
  if hparams.bottleneck_kind == "gumbel_softmax":
    return gumbel_softmax_discrete_bottleneck(
        x,
        hparams.bottleneck_bits,
        hparams.vq_beta,
        hparams.vq_decay,
        hparams.vq_epsilon,
        hparams.temperature_warmup_steps,
        hard=False,
        summary=True)

  raise ValueError(
      "Unsupported hparams.bottleneck_kind %s" % hparams.bottleneck_kind)
215618
def _chien_search_fast(self, sigma):
        ''''''
        # TODO: doesn't work when fcr is different than 1 (X values are incorrectly "shifted"...)
        # TODO: try to mix this approach with the optimized walk on only interesting values, implemented in _chien_search_faster()
        X = []
        j = []
        p = GF2int(self.generator)
        if not hasattr(self, 'const_poly'): self.const_poly = [GF2int(self.generator)**(i+self.fcr) for i in _range(self.gf2_charac, -1, -1)] # constant polynomial that will allow us to update the previous polynomial evaluation to get the next one
        const_poly = self.const_poly # caching for more efficiency since it never changes
        ev_poly, ev = sigma.evaluate_array( p**1 ) # compute the first polynomial evaluation
        # Try for each possible location
        for l in _range(1, self.gf2_charac+1): # range 1:256 is important: if you use range 0:255, if the last byte of the ecc symbols is corrupted, it won't be correctable! You need to use the range 1,256 to include this last byte.
            #l = (i+self.fcr)

            # Check if it's a root for the polynomial
            if ev == 0: # If it's 0, then bingo! It's an error location
                # Compute the error location polynomial X (will be directly used to compute the errors magnitudes inside the Forney algorithm)
                X.append( p**(-l) )
                # Compute the coefficient position (not the error position, it's actually the reverse: we compute the degree of the term where the error is located. To get the error position, just compute n-1-j).
                # This is different than the notes, I think the notes were in error
                # Notes said j values were just l, when it's actually 255-l
                j.append(self.gf2_charac - l)

            # Update the polynomial evaluation for the next iteration
            # we simply multiply each term[k] with alpha^k (where here alpha = p = GF2int(generator)).
            # For more info, see the presentation by Andrew Brown, or this one: http://web.ntpu.edu.tw/~yshan/BCH_decoding.pdf
            # TODO: parallelize this loop
            for i in _range(1, len(ev_poly)+1): # TODO: maybe the fcr != 1 fix should be put here?
                ev_poly[-i] *= const_poly[-i]
            # Compute the new evaluation by just summing
            ev = sum(ev_poly)

        return X, j
766162
private void tql2() {
    // This is derived from the Algol procedures tql2, by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.

    System.arraycopy(e, 1, e, 0, n - 1);
    e[n - 1] = 0.;

    double f = 0., tst1 = 0.;
    for(int l = 0; l < n; l++) {
      // Find small subdiagonal element
      tst1 = max(tst1, abs(d[l]) + abs(e[l]));
      int m = l;
      for(; m < n; ++m) {
        if(abs(e[m]) <= 0x1P-52 * tst1) {
          break;
        }
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.
      if(m > l) {
        do {
          // Compute implicit shift
          f += tql2ComputeImplicitShift(l);
          // Implicit QL transformation.
          tql2ImplicitQL(l, m, d[l + 1]);
          // Check for convergence.
          // TODO: Iteration limit?
        }
        while(abs(e[l]) > 0x1P-52 * tst1);
      }
      d[l] += f;
      e[l] = 0.;
    }
  }
116888
public static void symmRankTranA_sub(DMatrixRMaj a , DMatrixRMaj c ,
                                         int startIndexC )
    {
        // TODO update so that it doesn't modify/read parts that it shouldn't
        final double dataA[] = a.data;
        final double dataC[] = c.data;

//        for( int i = 0; i < a.numCols; i++ ) {
//            for( int k = 0; k < a.numRows; k++ ) {
//                double valA = dataA[k*a.numCols+i];
//
//                for( int j = i; j < a.numCols; j++ ) {
//                    dataC[startIndexC+i*c.numCols+j] -= valA * dataA[k*a.numCols+j];
//                }
//            }
//        }

        final int strideC = c.numCols + 1;
        for( int i = 0; i < a.numCols; i++ ) {
            int indexA = i;
            int endR = a.numCols;

            for( int k = 0; k < a.numRows; k++ , indexA += a.numCols , endR += a.numCols) {
                int indexC = startIndexC;
                final double valA = dataA[indexA];
                int indexR = indexA;

                while( indexR < endR ) {
                    dataC[indexC++] -= valA * dataA[indexR++];
                }
            }
            startIndexC += strideC;
        }

    }
1029768
private void _nonsymetricHqr2 ()
  {
    // 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 = m_nDim;
    int n = nn - 1;
    final int low = 0;
    final int high = nn - 1;
    double exshift = 0;
    double p = 0;
    double q = 0;
    double r = 0;
    double s = 0;
    double z = 0;
    double t;
    double w;
    double x;
    double y;

    // Store roots isolated by balanc and compute matrix norm
    double norm = 0.0;
    for (int i = 0; i < nn; i++)
    {
      if (i < low || i > high)
      {
        m_aEVd[i] = m_aHessenBerg[i][i];
        m_aEVe[i] = 0.0;
      }
      for (int j = Math.max (i - 1, 0); j < nn; j++)
      {
        norm += MathHelper.abs (m_aHessenBerg[i][j]);
      }
    }

    // Outer loop over eigenvalue index
    int iter = 0;
    while (n >= low)
    {
      // Look for single small sub-diagonal element
      int l = n;
      while (l > low)
      {
        s = MathHelper.abs (m_aHessenBerg[l - 1][l - 1]) + MathHelper.abs (m_aHessenBerg[l][l]);
        if (s == 0.0)
          s = norm;
        if (MathHelper.abs (m_aHessenBerg[l][l - 1]) < EPSILON * s)
          break;
        l--;
      }

      // Check for convergence
      // One root found
      if (l == n)
      {
        m_aHessenBerg[n][n] = m_aHessenBerg[n][n] + exshift;
        m_aEVd[n] = m_aHessenBerg[n][n];
        m_aEVe[n] = 0.0;
        n--;
        iter = 0;
        // Two roots found
      }
      else
        if (l == n - 1)
        {
          w = m_aHessenBerg[n][n - 1] * m_aHessenBerg[n - 1][n];
          p = (m_aHessenBerg[n - 1][n - 1] - m_aHessenBerg[n][n]) / 2.0;
          q = p * p + w;
          z = Math.sqrt (MathHelper.abs (q));
          m_aHessenBerg[n][n] = m_aHessenBerg[n][n] + exshift;
          m_aHessenBerg[n - 1][n - 1] = m_aHessenBerg[n - 1][n - 1] + exshift;
          x = m_aHessenBerg[n][n];

          // Real pair

          if (q >= 0)
          {
            if (p >= 0)
              z = p + z;
            else
              z = p - z;
            m_aEVd[n - 1] = x + z;
            m_aEVd[n] = m_aEVd[n - 1];
            if (z != 0.0)
              m_aEVd[n] = x - w / z;
            m_aEVe[n - 1] = 0.0;
            m_aEVe[n] = 0.0;
            x = m_aHessenBerg[n][n - 1];
            s = MathHelper.abs (x) + MathHelper.abs (z);
            p = x / s;
            q = z / s;
            r = Math.sqrt (p * p + q * q);
            p /= r;
            q /= r;

            // Row modification
            for (int j = n - 1; j < nn; j++)
            {
              z = m_aHessenBerg[n - 1][j];
              m_aHessenBerg[n - 1][j] = q * z + p * m_aHessenBerg[n][j];
              m_aHessenBerg[n][j] = q * m_aHessenBerg[n][j] - p * z;
            }

            // Column modification
            for (int i = 0; i <= n; i++)
            {
              z = m_aHessenBerg[i][n - 1];
              m_aHessenBerg[i][n - 1] = q * z + p * m_aHessenBerg[i][n];
              m_aHessenBerg[i][n] = q * m_aHessenBerg[i][n] - p * z;
            }

            // Accumulate transformations
            for (int i = low; i <= high; i++)
            {
              z = m_aEigenVector[i][n - 1];
              m_aEigenVector[i][n - 1] = q * z + p * m_aEigenVector[i][n];
              m_aEigenVector[i][n] = q * m_aEigenVector[i][n] - p * z;
            }

            // Complex pair
          }
          else
          {
            m_aEVd[n - 1] = x + p;
            m_aEVd[n] = x + p;
            m_aEVe[n - 1] = z;
            m_aEVe[n] = -z;
          }
          n -= 2;
          iter = 0;

          // No convergence yet
        }
        else
        {
          // Form shift
          x = m_aHessenBerg[n][n];
          y = 0.0;
          w = 0.0;
          if (l < n)
          {
            y = m_aHessenBerg[n - 1][n - 1];
            w = m_aHessenBerg[n][n - 1] * m_aHessenBerg[n - 1][n];
          }

          // Wilkinson's original ad hoc shift
          if (iter == 10)
          {
            exshift += x;
            for (int i = low; i <= n; i++)
            {
              m_aHessenBerg[i][i] -= x;
            }
            s = MathHelper.abs (m_aHessenBerg[n][n - 1]) + MathHelper.abs (m_aHessenBerg[n - 1][n - 2]);
            x = 0.75 * s;
            y = x;
            w = -0.4375 * s * s;
          }

          // MATLAB's new ad hoc shift
          if (iter == 30)
          {
            s = (y - x) / 2.0;
            s = s * s + w;
            if (s > 0)
            {
              s = Math.sqrt (s);
              if (y < x)
                s = -s;
              s = x - w / ((y - x) / 2.0 + s);
              for (int i = low; i <= n; i++)
              {
                m_aHessenBerg[i][i] -= s;
              }
              exshift += s;
              x = 0.964;
              y = x;
              w = x;
            }
          }

          iter++; // (Could check iteration count here.)

          // Look for two consecutive small sub-diagonal elements
          int m = n - 2;
          while (m >= l)
          {
            z = m_aHessenBerg[m][m];
            r = x - z;
            s = y - z;
            p = (r * s - w) / m_aHessenBerg[m + 1][m] + m_aHessenBerg[m][m + 1];
            q = m_aHessenBerg[m + 1][m + 1] - z - r - s;
            r = m_aHessenBerg[m + 2][m + 1];
            s = MathHelper.abs (p) + MathHelper.abs (q) + MathHelper.abs (r);
            p /= s;
            q /= s;
            r /= s;
            if (m == l)
              break;
            if (MathHelper.abs (m_aHessenBerg[m][m - 1]) *
                (MathHelper.abs (q) + MathHelper.abs (r)) < EPSILON *
                                                            (MathHelper.abs (p) *
                                                             (MathHelper.abs (m_aHessenBerg[m - 1][m - 1]) +
                                                              MathHelper.abs (z) +
                                                              MathHelper.abs (m_aHessenBerg[m + 1][m + 1]))))
            {
              break;
            }
            m--;
          }

          for (int i = m + 2; i <= n; i++)
          {
            m_aHessenBerg[i][i - 2] = 0.0;
            if (i > m + 2)
            {
              m_aHessenBerg[i][i - 3] = 0.0;
            }
          }

          // Double QR step involving rows l:n and columns m:n

          for (int k = m; k <= n - 1; k++)
          {
            final boolean notlast = (k != n - 1);
            if (k != m)
            {
              p = m_aHessenBerg[k][k - 1];
              q = m_aHessenBerg[k + 1][k - 1];
              r = (notlast ? m_aHessenBerg[k + 2][k - 1] : 0.0);
              x = MathHelper.abs (p) + MathHelper.abs (q) + MathHelper.abs (r);
              if (x == 0.0)
              {
                continue;
              }
              p /= x;
              q /= x;
              r /= x;
            }

            s = Math.sqrt (p * p + q * q + r * r);
            if (p < 0)
            {
              s = -s;
            }
            if (s != 0)
            {
              if (k != m)
              {
                m_aHessenBerg[k][k - 1] = -s * x;
              }
              else
                if (l != m)
                {
                  m_aHessenBerg[k][k - 1] = -m_aHessenBerg[k][k - 1];
                }
              p += s;
              x = p / s;
              y = q / s;
              z = r / s;
              q /= p;
              r /= p;

              // Row modification

              for (int j = k; j < nn; j++)
              {
                p = m_aHessenBerg[k][j] + q * m_aHessenBerg[k + 1][j];
                if (notlast)
                {
                  p += r * m_aHessenBerg[k + 2][j];
                  m_aHessenBerg[k + 2][j] = m_aHessenBerg[k + 2][j] - p * z;
                }
                m_aHessenBerg[k][j] = m_aHessenBerg[k][j] - p * x;
                m_aHessenBerg[k + 1][j] = m_aHessenBerg[k + 1][j] - p * y;
              }

              // Column modification

              for (int i = 0; i <= Math.min (n, k + 3); i++)
              {
                p = x * m_aHessenBerg[i][k] + y * m_aHessenBerg[i][k + 1];
                if (notlast)
                {
                  p += z * m_aHessenBerg[i][k + 2];
                  m_aHessenBerg[i][k + 2] = m_aHessenBerg[i][k + 2] - p * r;
                }
                m_aHessenBerg[i][k] = m_aHessenBerg[i][k] - p;
                m_aHessenBerg[i][k + 1] = m_aHessenBerg[i][k + 1] - p * q;
              }

              // Accumulate transformations

              for (int i = low; i <= high; i++)
              {
                p = x * m_aEigenVector[i][k] + y * m_aEigenVector[i][k + 1];
                if (notlast)
                {
                  p += z * m_aEigenVector[i][k + 2];
                  m_aEigenVector[i][k + 2] = m_aEigenVector[i][k + 2] - p * r;
                }
                m_aEigenVector[i][k] = m_aEigenVector[i][k] - p;
                m_aEigenVector[i][k + 1] = m_aEigenVector[i][k + 1] - p * q;
              }
            } // (s != 0)
          } // k loop
        } // check convergence
    } // while (n >= low)

    // Backsubstitute to find vectors of upper triangular form

    if (norm == 0.0)
    {
      return;
    }

    for (n = nn - 1; n >= 0; n--)
    {
      p = m_aEVd[n];
      q = m_aEVe[n];

      // Real vector

      if (q == 0)
      {
        int l = n;
        m_aHessenBerg[n][n] = 1.0;
        for (int i = n - 1; i >= 0; i--)
        {
          w = m_aHessenBerg[i][i] - p;
          r = 0.0;
          for (int j = l; j <= n; j++)
          {
            r += m_aHessenBerg[i][j] * m_aHessenBerg[j][n];
          }
          if (m_aEVe[i] < 0.0)
          {
            z = w;
            s = r;
          }
          else
          {
            l = i;
            if (m_aEVe[i] == 0.0)
            {
              if (w != 0.0)
              {
                m_aHessenBerg[i][n] = -r / w;
              }
              else
              {
                m_aHessenBerg[i][n] = -r / (EPSILON * norm);
              }

              // Solve real equations

            }
            else
            {
              x = m_aHessenBerg[i][i + 1];
              y = m_aHessenBerg[i + 1][i];
              q = (m_aEVd[i] - p) * (m_aEVd[i] - p) + m_aEVe[i] * m_aEVe[i];
              t = (x * s - z * r) / q;
              m_aHessenBerg[i][n] = t;
              if (MathHelper.abs (x) > MathHelper.abs (z))
              {
                m_aHessenBerg[i + 1][n] = (-r - w * t) / x;
              }
              else
              {
                m_aHessenBerg[i + 1][n] = (-s - y * t) / z;
              }
            }

            // Overflow control

            t = MathHelper.abs (m_aHessenBerg[i][n]);
            if ((EPSILON * t) * t > 1)
            {
              for (int j = i; j <= n; j++)
              {
                m_aHessenBerg[j][n] = m_aHessenBerg[j][n] / t;
              }
            }
          }
        }

        // Complex vector

      }
      else
        if (q < 0)
        {
          int l = n - 1;

          // Last vector component imaginary so matrix is triangular

          if (MathHelper.abs (m_aHessenBerg[n][n - 1]) > MathHelper.abs (m_aHessenBerg[n - 1][n]))
          {
            m_aHessenBerg[n - 1][n - 1] = q / m_aHessenBerg[n][n - 1];
            m_aHessenBerg[n - 1][n] = -(m_aHessenBerg[n][n] - p) / m_aHessenBerg[n][n - 1];
          }
          else
          {
            _cdiv (0.0, -m_aHessenBerg[n - 1][n], m_aHessenBerg[n - 1][n - 1] - p, q);
            m_aHessenBerg[n - 1][n - 1] = m_dCdivr;
            m_aHessenBerg[n - 1][n] = m_dCdivi;
          }
          m_aHessenBerg[n][n - 1] = 0.0;
          m_aHessenBerg[n][n] = 1.0;
          for (int i = n - 2; i >= 0; i--)
          {
            double ra = 0.0;
            double sa = 0.0;
            double vr;
            double vi;
            for (int j = l; j <= n; j++)
            {
              ra += m_aHessenBerg[i][j] * m_aHessenBerg[j][n - 1];
              sa += m_aHessenBerg[i][j] * m_aHessenBerg[j][n];
            }
            w = m_aHessenBerg[i][i] - p;

            if (m_aEVe[i] < 0.0)
            {
              z = w;
              r = ra;
              s = sa;
            }
            else
            {
              l = i;
              if (m_aEVe[i] == 0)
              {
                _cdiv (-ra, -sa, w, q);
                m_aHessenBerg[i][n - 1] = m_dCdivr;
                m_aHessenBerg[i][n] = m_dCdivi;
              }
              else
              {

                // Solve complex equations

                x = m_aHessenBerg[i][i + 1];
                y = m_aHessenBerg[i + 1][i];
                vr = (m_aEVd[i] - p) * (m_aEVd[i] - p) + m_aEVe[i] * m_aEVe[i] - q * q;
                vi = (m_aEVd[i] - p) * 2.0 * q;
                if (vr == 0.0 && vi == 0.0)
                {
                  vr = EPSILON *
                       norm *
                       (MathHelper.abs (w) +
                        MathHelper.abs (q) +
                        MathHelper.abs (x) +
                        MathHelper.abs (y) +
                        MathHelper.abs (z));
                }
                _cdiv (x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
                m_aHessenBerg[i][n - 1] = m_dCdivr;
                m_aHessenBerg[i][n] = m_dCdivi;
                if (MathHelper.abs (x) > (MathHelper.abs (z) + MathHelper.abs (q)))
                {
                  m_aHessenBerg[i + 1][n - 1] = (-ra - w * m_aHessenBerg[i][n - 1] + q * m_aHessenBerg[i][n]) / x;
                  m_aHessenBerg[i + 1][n] = (-sa - w * m_aHessenBerg[i][n] - q * m_aHessenBerg[i][n - 1]) / x;
                }
                else
                {
                  _cdiv (-r - y * m_aHessenBerg[i][n - 1], -s - y * m_aHessenBerg[i][n], z, q);
                  m_aHessenBerg[i + 1][n - 1] = m_dCdivr;
                  m_aHessenBerg[i + 1][n] = m_dCdivi;
                }
              }

              // Overflow control

              t = Math.max (MathHelper.abs (m_aHessenBerg[i][n - 1]), MathHelper.abs (m_aHessenBerg[i][n]));
              if ((EPSILON * t) * t > 1)
              {
                for (int j = i; j <= n; j++)
                {
                  m_aHessenBerg[j][n - 1] = m_aHessenBerg[j][n - 1] / t;
                  m_aHessenBerg[j][n] = m_aHessenBerg[j][n] / t;
                }
              }
            }
          }
        }
    }

    // Vectors of isolated roots

    for (int i = 0; i < nn; i++)
    {
      if (i < low || i > high)
      {
        for (int j = i; j < nn; j++)
        {
          m_aEigenVector[i][j] = m_aHessenBerg[i][j];
        }
      }
    }

    // Back transformation to get eigenvectors of original matrix

    for (int j = nn - 1; j >= low; j--)
    {
      for (int i = low; i <= high; i++)
      {
        z = 0.0;
        for (int k = low; k <= Math.min (j, high); k++)
        {
          z += m_aEigenVector[i][k] * m_aHessenBerg[k][j];
        }
        m_aEigenVector[i][j] = z;
      }
    }
  }
600915
function brcomp(a, b, x, y) {
      var lambda;

      lambda = a - (a + b) * x;
      if (Math.min(a, b) < 8) {
         return a * Math.log(x) + b * Math.log(y) - lbeta(a, b);
      }
      return 0.5 * (Math.log(a) + Math.log(b) - C.log2pi - Math.log(a + b)) -
             a * phi(1 - lambda / a) - b * phi(1 + lambda / b) -
             bcorr(a, b);
   }
459961
static void computePPXP( DMatrixRMaj P , DMatrixRMaj P_plus, DMatrixRMaj X , int offsetX,  DMatrixRMaj output ) {
		final int N = P.numRows;
		output.reshape(N,2);

		// diag(X) * P <-- N x 2
		for (int i = 0, index=0; i < N; i++, index += 2) {
			double x = -X.data[index+offsetX];
			output.data[index] = x*P.data[index];
			output.data[index+1] = x*P.data[index+1];
		}

		// A = P_plus*( diag(x)*P ) <-- 2 x 2
		double a00=0,a01=0,a10=0,a11=0;
		for (int i = 0, index=0; i < N; i++, index += 2) {
			a00 += P_plus.data[i]*output.data[index];
			a01 += P_plus.data[i]*output.data[index+1];
			a10 += P_plus.data[i+N]*output.data[index];
			a11 += P_plus.data[i+N]*output.data[index+1];
		}

		// P*A <-- N x 2
		for (int i = 0, index=0; i < N; i++, index += 2) {
			output.data[index]   = P.data[index]*a00 + P.data[index+1]*a10;
			output.data[index+1] = P.data[index]*a01+ P.data[index+1]*a11;
		}
	}
1134930
function deterministicGenerateK(curve, hash, d, checkSig, nonce) {
    enforceType("Buffer", hash);
    enforceType(BigInteger, d);

    if (nonce) {
        hash = sha256(Buffer.concat([hash, Buffer.alloc(nonce)]));
    }

    // sanity check
    assert.equal(hash.length, 32, "Hash must be 256 bit");

    var x = d.toBuffer(32);
    var k = Buffer.alloc(32);
    var v = Buffer.alloc(32);

    // Step B
    v.fill(1);

    // Step C
    k.fill(0);

    // Step D
    k = HmacSHA256(Buffer.concat([v, new Buffer([0]), x, hash]), k);

    // Step E
    v = HmacSHA256(v, k);

    // Step F
    k = HmacSHA256(Buffer.concat([v, new Buffer([1]), x, hash]), k);

    // Step G
    v = HmacSHA256(v, k);

    // Step H1/H2a, ignored as tlen === qlen (256 bit)
    // Step H2b
    v = HmacSHA256(v, k);

    var T = BigInteger.fromBuffer(v);

    // Step H3, repeat until T is within the interval [1, n - 1]
    while (T.signum() <= 0 || T.compareTo(curve.n) >= 0 || !checkSig(T)) {
        k = HmacSHA256(Buffer.concat([v, new Buffer([0])]), k);
        v = HmacSHA256(v, k);

        // Step H1/H2a, again, ignored as tlen === qlen (256 bit)
        // Step H2b again
        v = HmacSHA256(v, k);

        T = BigInteger.fromBuffer(v);
    }

    return T;
}