Standard Normal constiate using Box-Muller transform. http://stackoverflow.com/a/36481059
485897
function randnBM() {
  const u = 1 - Math.random() // Subtraction to flip [0, 1) to (0, 1].
  const v = 1 - Math.random()
  return Math.sqrt(-2.0 * Math.log(u)) * Math.cos(2.0 * Math.PI * v)
}
574195
function() {

      // Q^-1 := Q' / |Q|^2
      //       = [w / (w^2 + |v|^2), -v / (w^2 + |v|^2)]

      // Proof:
      // Q * Q^-1 = [w, v] * [w / (w^2 + |v|^2), -v / (w^2 + |v|^2)]
      //          = [1, 0]
      // Q^-1 * Q = [w / (w^2 + |v|^2), -v / (w^2 + |v|^2)] * [w, v]
      //          = [1, 0].

      var w = this['w'];
      var x = this['x'];
      var y = this['y'];
      var z = this['z'];

      var normSq = w * w + x * x + y * y + z * z;

      if (normSq === 0) {
        return Quaternion['ZERO']; // TODO: Is the result zero or one when the norm is zero?
      }

      normSq = 1 / normSq;

      return new Quaternion(w * normSq, -x * normSq, -y * normSq, -z * normSq);
    }
883252
function rotr_64(x, n)
	{
		var retVal = null, tmp = new Int_64(x.highOrder, x.lowOrder);

		if (32 >= n)
		{
			retVal = new Int_64(
					(tmp.highOrder >>> n) | ((tmp.lowOrder << (32 - n)) & 0xFFFFFFFF),
					(tmp.lowOrder >>> n) | ((tmp.highOrder << (32 - n)) & 0xFFFFFFFF)
				);
		}
		else
		{
			retVal = new Int_64(
					(tmp.lowOrder >>> (n - 32)) | ((tmp.highOrder << (64 - n)) & 0xFFFFFFFF),
					(tmp.highOrder >>> (n - 32)) | ((tmp.lowOrder << (64 - n)) & 0xFFFFFFFF)
				);
		}

		return retVal;
	}
589944
function zOrder(x, y, minX, minY, size) {
        // coords are transformed into non-negative 15-bit integer range
        x = 32767 * (x - minX) / size;
        y = 32767 * (y - minY) / size;
        x = (x | (x << 8)) & 0x00FF00FF;
        x = (x | (x << 4)) & 0x0F0F0F0F;
        x = (x | (x << 2)) & 0x33333333;
        x = (x | (x << 1)) & 0x55555555;
        y = (y | (y << 8)) & 0x00FF00FF;
        y = (y | (y << 4)) & 0x0F0F0F0F;
        y = (y | (y << 2)) & 0x33333333;
        y = (y | (y << 1)) & 0x55555555;
        return x | (y << 1);
    }
1208934
function bnpDivRemTo(m, q, r) {
    var pm = m.abs();
    if (pm.t <= 0) return;
    var pt = this.abs();
    if (pt.t < pm.t) {
        if (q != null) q.fromInt(0);
        if (r != null) this.copyTo(r);
        return;
    }
    if (r == null) r = nbi();
    var y = nbi(), ts = this.s, ms = m.s;
    var nsh = this.DB - nbits(pm[pm.t - 1]);	// normalize modulus
    if (nsh > 0) {
        pm.lShiftTo(nsh, y);
        pt.lShiftTo(nsh, r);
    }
    else {
        pm.copyTo(y);
        pt.copyTo(r);
    }
    var ys = y.t;
    var y0 = y[ys - 1];
    if (y0 === 0) return;
    var yt = y0 * (1 << this.F1) + ((ys > 1) ? y[ys - 2] >> this.F2 : 0);
    var d1 = this.FV / yt, d2 = (1 << this.F1) / yt, e = 1 << this.F2;
    var i = r.t, j = i - ys, t = (q == null) ? nbi() : q;
    y.dlShiftTo(j, t);
    if (r.compareTo(t) >= 0) {
        r[r.t++] = 1;
        r.subTo(t, r);
    }
    BigInteger.ONE.dlShiftTo(ys, t);
    t.subTo(y, y);	// "negative" y so we can replace sub with am later
    while (y.t < ys) y[y.t++] = 0;
    while (--j >= 0) {
        // Estimate quotient digit
        var qd = (r[--i] == y0) ? this.DM : Math.floor(r[i] * d1 + (r[i - 1] + e) * d2);
        if ((r[i] += y.am(0, qd, r, j, 0, ys)) < qd) {	// Try it out
            y.dlShiftTo(j, t);
            r.subTo(t, r);
            while (r[i] < --qd) r.subTo(t, r);
        }
    }
    if (q != null) {
        r.drShiftTo(ys, q);
        if (ts != ms) BigInteger.ZERO.subTo(q, q);
    }
    r.t = ys;
    r.clamp();
    if (nsh > 0) r.rShiftTo(nsh, r);	// Denormalize remainder
    if (ts < 0) BigInteger.ZERO.subTo(r, r);
}
1282703
function inverse(m, dst) {
  dst = dst || new MatType(16);

  const m00 = m[0 * 4 + 0];
  const m01 = m[0 * 4 + 1];
  const m02 = m[0 * 4 + 2];
  const m03 = m[0 * 4 + 3];
  const m10 = m[1 * 4 + 0];
  const m11 = m[1 * 4 + 1];
  const m12 = m[1 * 4 + 2];
  const m13 = m[1 * 4 + 3];
  const m20 = m[2 * 4 + 0];
  const m21 = m[2 * 4 + 1];
  const m22 = m[2 * 4 + 2];
  const m23 = m[2 * 4 + 3];
  const m30 = m[3 * 4 + 0];
  const m31 = m[3 * 4 + 1];
  const m32 = m[3 * 4 + 2];
  const m33 = m[3 * 4 + 3];
  const tmp_0  = m22 * m33;
  const tmp_1  = m32 * m23;
  const tmp_2  = m12 * m33;
  const tmp_3  = m32 * m13;
  const tmp_4  = m12 * m23;
  const tmp_5  = m22 * m13;
  const tmp_6  = m02 * m33;
  const tmp_7  = m32 * m03;
  const tmp_8  = m02 * m23;
  const tmp_9  = m22 * m03;
  const tmp_10 = m02 * m13;
  const tmp_11 = m12 * m03;
  const tmp_12 = m20 * m31;
  const tmp_13 = m30 * m21;
  const tmp_14 = m10 * m31;
  const tmp_15 = m30 * m11;
  const tmp_16 = m10 * m21;
  const tmp_17 = m20 * m11;
  const tmp_18 = m00 * m31;
  const tmp_19 = m30 * m01;
  const tmp_20 = m00 * m21;
  const tmp_21 = m20 * m01;
  const tmp_22 = m00 * m11;
  const tmp_23 = m10 * m01;

  const t0 = (tmp_0 * m11 + tmp_3 * m21 + tmp_4 * m31) -
      (tmp_1 * m11 + tmp_2 * m21 + tmp_5 * m31);
  const t1 = (tmp_1 * m01 + tmp_6 * m21 + tmp_9 * m31) -
      (tmp_0 * m01 + tmp_7 * m21 + tmp_8 * m31);
  const t2 = (tmp_2 * m01 + tmp_7 * m11 + tmp_10 * m31) -
      (tmp_3 * m01 + tmp_6 * m11 + tmp_11 * m31);
  const t3 = (tmp_5 * m01 + tmp_8 * m11 + tmp_11 * m21) -
      (tmp_4 * m01 + tmp_9 * m11 + tmp_10 * m21);

  const d = 1.0 / (m00 * t0 + m10 * t1 + m20 * t2 + m30 * t3);

  dst[ 0] = d * t0;
  dst[ 1] = d * t1;
  dst[ 2] = d * t2;
  dst[ 3] = d * t3;
  dst[ 4] = d * ((tmp_1 * m10 + tmp_2 * m20 + tmp_5 * m30) -
          (tmp_0 * m10 + tmp_3 * m20 + tmp_4 * m30));
  dst[ 5] = d * ((tmp_0 * m00 + tmp_7 * m20 + tmp_8 * m30) -
          (tmp_1 * m00 + tmp_6 * m20 + tmp_9 * m30));
  dst[ 6] = d * ((tmp_3 * m00 + tmp_6 * m10 + tmp_11 * m30) -
          (tmp_2 * m00 + tmp_7 * m10 + tmp_10 * m30));
  dst[ 7] = d * ((tmp_4 * m00 + tmp_9 * m10 + tmp_10 * m20) -
          (tmp_5 * m00 + tmp_8 * m10 + tmp_11 * m20));
  dst[ 8] = d * ((tmp_12 * m13 + tmp_15 * m23 + tmp_16 * m33) -
          (tmp_13 * m13 + tmp_14 * m23 + tmp_17 * m33));
  dst[ 9] = d * ((tmp_13 * m03 + tmp_18 * m23 + tmp_21 * m33) -
          (tmp_12 * m03 + tmp_19 * m23 + tmp_20 * m33));
  dst[10] = d * ((tmp_14 * m03 + tmp_19 * m13 + tmp_22 * m33) -
          (tmp_15 * m03 + tmp_18 * m13 + tmp_23 * m33));
  dst[11] = d * ((tmp_17 * m03 + tmp_20 * m13 + tmp_23 * m23) -
          (tmp_16 * m03 + tmp_21 * m13 + tmp_22 * m23));
  dst[12] = d * ((tmp_14 * m22 + tmp_17 * m32 + tmp_13 * m12) -
          (tmp_16 * m32 + tmp_12 * m12 + tmp_15 * m22));
  dst[13] = d * ((tmp_20 * m32 + tmp_12 * m02 + tmp_19 * m22) -
          (tmp_18 * m22 + tmp_21 * m32 + tmp_13 * m02));
  dst[14] = d * ((tmp_18 * m12 + tmp_23 * m32 + tmp_15 * m02) -
          (tmp_22 * m32 + tmp_14 * m02 + tmp_19 * m12));
  dst[15] = d * ((tmp_22 * m22 + tmp_16 * m02 + tmp_21 * m12) -
          (tmp_20 * m12 + tmp_23 * m22 + tmp_17 * m02));

  return dst;
}
735329
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);
  }
1281631
function _mmult(a, m) {
        m = m.slice();

        var m0 = m[0];
        var m1 = m[1];
        var m2 = m[2];
        var m3 = m[3];
        var m4 = m[4];
        var m5 = m[5];
        var m6 = m[6];
        var m7 = m[7];
        var m8 = m[8];

        m[0] = a[0] * m0 + a[1] * m3;
        m[1] = a[0] * m1 + a[1] * m4;
        m[3] = a[3] * m0 + a[4] * m3;
        m[4] = a[3] * m1 + a[4] * m4;
        m[6] = a[6] * m0 + a[7] * m3 + m6;
        m[7] = a[6] * m1 + a[7] * m4 + m7;

        return transform(m);
    }
116516
function mpn_dbl(r, n, l) {
	var c = 0;
	for (var i = 0; i < l; ++i) {
		r[i] = (c = c >>> 31 | n[i] << 1) & 0x7FFFFFFF;
	}
	return c >>> 31;
}
370157
private static void hqr(DenseMatrix A, double[] d, double[] e) {
        int n = A.nrows();
        int nn, m, l, k, j, its, i, mmin;
        double z, y, x, w, v, u, t, s, r = 0.0, q = 0.0, p = 0.0, anorm = 0.0;

        for (i = 0; i < n; i++) {
            for (j = Math.max(i - 1, 0); j < n; j++) {
                anorm += Math.abs(A.get(i, j));
            }
        }
        nn = n - 1;
        t = 0.0;
        while (nn >= 0) {
            its = 0;
            do {
                for (l = nn; l > 0; l--) {
                    s = Math.abs(A.get(l - 1, l - 1) + Math.abs(A.get(l, l)));
                    if (s == 0.0) {
                        s = anorm;
                    }
                    if (Math.abs(A.get(l, l - 1)) <= Math.EPSILON * s) {
                        A.set(l, l - 1, 0.0);
                        break;
                    }
                }
                x = A.get(nn, nn);
                if (l == nn) {
                    d[nn--] = x + t;
                } else {
                    y = A.get(nn - 1, nn - 1);
                    w = A.get(nn, nn - 1) * A.get(nn - 1, nn);
                    if (l == nn - 1) {
                        p = 0.5 * (y - x);
                        q = p * p + w;
                        z = Math.sqrt(Math.abs(q));
                        x += t;
                        if (q >= 0.0) {
                            z = p + Math.copySign(z, p);
                            d[nn - 1] = d[nn] = x + z;
                            if (z != 0.0) {
                                d[nn] = x - w / z;
                            }
                        } else {
                            d[nn] = x + p;
                            e[nn] = -z;
                            d[nn - 1] = d[nn];
                            e[nn - 1] = -e[nn];
                        }
                        nn -= 2;
                    } else {
                        if (its == 30) {
                            throw new IllegalStateException("Too many iterations in hqr");
                        }
                        if (its == 10 || its == 20) {
                            t += x;
                            for (i = 0; i < nn + 1; i++) {
                                A.sub(i, i, x);
                            }
                            s = Math.abs(A.get(nn, nn - 1)) + Math.abs(A.get(nn - 1, nn - 2));
                            y = x = 0.75 * s;
                            w = -0.4375 * s * s;
                        }
                        ++its;
                        for (m = nn - 2; m >= l; m--) {
                            z = A.get(m, m);
                            r = x - z;
                            s = y - z;
                            p = (r * s - w) / A.get(m + 1, m) + A.get(m, m + 1);
                            q = A.get(m + 1, m + 1) - z - r - s;
                            r = A.get(m + 2, m + 1);
                            s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                            p /= s;
                            q /= s;
                            r /= s;
                            if (m == l) {
                                break;
                            }
                            u = Math.abs(A.get(m, m - 1)) * (Math.abs(q) + Math.abs(r));
                            v = Math.abs(p) * (Math.abs(A.get(m - 1, m - 1)) + Math.abs(z) + Math.abs(A.get(m + 1, m + 1)));
                            if (u <= Math.EPSILON * v) {
                                break;
                            }
                        }
                        for (i = m; i < nn - 1; i++) {
                            A.set(i + 2, i, 0.0);
                            if (i != m) {
                                A.set(i + 2, i - 1, 0.0);
                            }
                        }
                        for (k = m; k < nn; k++) {
                            if (k != m) {
                                p = A.get(k, k - 1);
                                q = A.get(k + 1, k - 1);
                                r = 0.0;
                                if (k + 1 != nn) {
                                    r = A.get(k + 2, k - 1);
                                }
                                if ((x = Math.abs(p) + Math.abs(q) + Math.abs(r)) != 0.0) {
                                    p /= x;
                                    q /= x;
                                    r /= x;
                                }
                            }
                            if ((s = Math.copySign(Math.sqrt(p * p + q * q + r * r), p)) != 0.0) {
                                if (k == m) {
                                    if (l != m) {
                                        A.set(k, k - 1, -A.get(k, k - 1));
                                    }
                                } else {
                                    A.set(k, k - 1, -s * x);
                                }
                                p += s;
                                x = p / s;
                                y = q / s;
                                z = r / s;
                                q /= p;
                                r /= p;
                                for (j = k; j < nn + 1; j++) {
                                    p = A.get(k, j) + q * A.get(k + 1, j);
                                    if (k + 1 != nn) {
                                        p += r * A.get(k + 2, j);
                                        A.sub(k + 2, j, p * z);
                                    }
                                    A.sub(k + 1, j, p * y);
                                    A.sub(k, j, p * x);
                                }
                                mmin = nn < k + 3 ? nn : k + 3;
                                for (i = l; i < mmin + 1; i++) {
                                    p = x * A.get(i, k) + y * A.get(i, k + 1);
                                    if (k + 1 != nn) {
                                        p += z * A.get(i, k + 2);
                                        A.sub(i, k + 2, p * r);
                                    }
                                    A.sub(i, k + 1, p * q);
                                    A.sub(i, k, p);
                                }
                            }
                        }
                    }
                }
            } while (l + 1 < nn);
        }
    }
1239796
private int trMedian5(int ISAd, int v1, int v2, int v3, int v4, int v5) {
        int t;
        if (SA[ISAd + SA[v2]] > SA[ISAd + SA[v3]]) {
            t = v2;
            v2 = v3;
            v3 = t;
        }
        if (SA[ISAd + SA[v4]] > SA[ISAd + SA[v5]]) {
            t = v4;
            v4 = v5;
            v5 = t;
        }
        if (SA[ISAd + SA[v2]] > SA[ISAd + SA[v4]]) {
            t = v2;
            v2 = v4;
            v4 = t;
            t = v3;
            v3 = v5;
            v5 = t;
        }
        if (SA[ISAd + SA[v1]] > SA[ISAd + SA[v3]]) {
            t = v1;
            v1 = v3;
            v3 = t;
        }
        if (SA[ISAd + SA[v1]] > SA[ISAd + SA[v4]]) {
            t = v1;
            v1 = v4;
            v4 = t;
            t = v3;
            v3 = v5;
            v5 = t;
        }
        if (SA[ISAd + SA[v3]] > SA[ISAd + SA[v4]]) {
            return v4;
        }
        return v3;
    }
447213
function householderTransform (x) {
  // Based on "Matrix Computations" by Golub, Van Loan
  const n = x.dimension;
  const v = x.clone();
  v.x = 1;
  let sigma = 0;
  for (let i = 1; i < n; ++i) {
    const vi = v.getComponent(i);
    sigma += vi * vi;
  }
  let beta = 0;
  if (sigma !== 0) {
    const x1 = x.x;
    const mu = _Math.sqrt(x1 * x1 + sigma);
    const v1 = (x1 <= 0 ? x1 - mu : -sigma / (x1 + mu));
    v.x = v1;
    beta = 2 * v1 * v1 / (sigma + v1 * v1);
    v.multiplyScalar(1.0 / v1);
  }
  return {
    v,
    beta
  };
}
138937
def normalize(v):
    ''''''
    if 0 == np.linalg.norm(v):
        return v
    return v / np.linalg.norm(v)
1144996
def xors_constraint(v=1):
    
    assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__()
    model, x, y, z = _init()
    r = True
    model.addConsXor([x,y,z], r)
    model.addCons(x==v)
    _optimize("Standard XOR (as boolean)", model)
178538
def normalize(x):
    ''''''
    x = x.astype(float)
    x -= x.mean()
    return x / float(x.std())
1059769
private void similarTransform( int k) {
        double t[] = QT.data;

        // find the largest value in this column
        // this is used to normalize the column and mitigate overflow/underflow
        double max = QrHelperFunctions_ZDRM.computeRowMax(QT,k,k+1,N);

        if( max > 0 ) {
            double gamma = QrHelperFunctions_ZDRM.computeTauGammaAndDivide(k*N+k+1, k*N+N, t, max, tau);
            gammas[k] = gamma;

            // divide u by u_0
            double real_u_0 = t[(k*N+k+1)*2]   + tau.real;
            double imag_u_0 = t[(k*N+k+1)*2+1] + tau.imaginary;
            QrHelperFunctions_ZDRM.divideElements(k+2, N, t, k*N, real_u_0,imag_u_0 );

            // A column is zeroed first.  However a row is being used to store because it reduces
            // cache misses.  Need to compute the conjugate to have the correct householder operation
            for (int i = k+2; i < N; i++) {
                t[(k*N+i)*2+1] = -t[(k*N+i)*2+1];
            }

            t[(k*N+k+1)*2]   = 1.0;
            t[(k*N+k+1)*2+1] = 0;

            // ---------- Specialized householder that takes advantage of the symmetry
//            QrHelperFunctions_ZDRM.rank1UpdateMultR(QT,QT.data,k*N,gamma,k+1,k+1,N,w);
//            QrHelperFunctions_ZDRM.rank1UpdateMultL(QT,QT.data,k*N,gamma,k+1,k+1,N);

            householderSymmetric(k,gamma);

            // since the first element in the householder vector is known to be 1
            // store the full upper hessenberg
            t[(k*N+k+1)*2]   = -tau.real*max;
            t[(k*N+k+1)*2+1] = -tau.imaginary*max;

        } else {
            gammas[k] = 0;
        }
    }
545209
private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
        int n = A.nrows();
        for (int mp = n - 2; mp > 0; mp--) {
            for (int k = mp + 1; k < n; k++) {
                V.set(k, mp, A.get(k, mp - 1));
            }
            int i = perm[mp];
            if (i != mp) {
                for (int j = mp; j < n; j++) {
                    V.set(mp, j, V.get(i, j));
                    V.set(i, j, 0.0);
                }
                V.set(i, mp, 1.0);
            }
        }
    }