%D.

    Formula:
    %D = SMA(%K, 3)
546094
def percent_d(data, period):
    
    p_k = percent_k(data, period)
    percent_d = sma(p_k, 3)
    return percent_d
826201
private static double scalbn(double x, int n) {
        int hx = hi(x);
        int lx = lo(x);
        int k = (hx & 0x7ff00000) >> 20;         // extract exponent
        if (k == 0) {                          // 0 or subnormal x
            if ((lx | (hx & 0x7fffffff)) == 0) {
                return x; // +-0
            }
            x *= two54;
            hx = hi(x);
            k = ((hx & 0x7ff00000) >> 20) - 54;
            if (n < -50000) {
                return tiny * x;                   // underflow
            }
        }
        if (k == 0x7ff) {
            return x + x;                        // NaN or Inf
        }
        k = k + n;
        if (k > 0x7fe) {
            return huge * copysign(huge, x);      // overflow
        }
        if (k > 0) {
            // normal result
            return fromhilo((hx & 0x800fffff) | (k << 20), lo(x));
        }
        if (k <= -54) {
            if (n > 50000) {
                // in case integer overflow in n+k
                return huge * copysign(huge, x);  // overflow
            }
            else {
                return tiny * copysign(tiny, x);    // underflow
            }
        }
        k += 54;                             // subnormal result
        x = fromhilo((hx & 0x800fffff) | (k << 20), lo(x));
        return x * twom54;
    }
229210
public Vector4d fma(Vector4dc a, Vector4dc b, Vector4d dest) {
        dest.x = x + a.x() * b.x();
        dest.y = y + a.y() * b.y();
        dest.z = z + a.z() * b.z();
        dest.w = w + a.w() * b.w();
        return dest;
    }
1162734
def ccmodmd_xstep(k):
    

    YU0 = mp_D_Y0 - mp_D_U0[k]
    YU1 = mp_D_Y1[k] + mp_S[k] - mp_D_U1[k]
    b = sl.rfftn(YU0, None, mp_cri.axisN) + \
      np.conj(mp_Zf[k]) * sl.rfftn(YU1, None, mp_cri.axisN)
    Xf = sl.solvedbi_sm(mp_Zf[k], 1.0, b, axis=mp_cri.axisM)
    mp_D_X[k] = sl.irfftn(Xf, mp_cri.Nv, mp_cri.axisN)
    mp_DX[k] = sl.irfftn(sl.inner(Xf, mp_Zf[k]), mp_cri.Nv, mp_cri.axisN)
1295614
private boolean smo(double epsgr) {
        SupportVector i = svmax;
        int ii = gmaxindex;
        double old_alpha_i = i.alpha[ii];

        if (i.kcache == null) {
            gram(i);
        }

        SupportVector j = svmin;        
        int jj = gminindex;        
        double old_alpha_j = j.alpha[jj];
        
        // Second order working set selection.
        double best = 0.0;
        double gi = ii == 0 ? -i.g[0] : i.g[1];
        for (SupportVector v : sv) {
            double curv = i.k + v.k - 2 * kernel.k(i.x, v.x);
            if (curv <= 0.0) curv = TAU;
            
            double gj = -v.g[0];
            if (v.alpha[0] > 0.0 && gj < gi) {
                double gain = -Math.sqr(gi - gj) / curv;
                if (gain < best) {
                    best = gain;
                    j = v;
                    jj = 0;
                    old_alpha_j = j.alpha[0];
                }
            }
            
            gj = v.g[1];
            if (v.alpha[1] < v.C && gj < gi) {
                double gain = -Math.sqr(gi - gj) / curv;
                if (gain < best) {
                    best = gain;
                    j = v;
                    jj = 1;
                    old_alpha_j = j.alpha[1];
                }                
            }
        }
        
        if (j.kcache == null) {
            gram(j);
        }
        
        // Determine curvature
        double curv = i.k + j.k - 2 * kernel.k(i.x, j.x);
        if (curv <= 0.0) curv = TAU;

        if (ii != jj) {
            double delta = (-i.g[ii] - j.g[jj]) / curv;
            double diff = i.alpha[ii] - j.alpha[jj];
            i.alpha[ii] += delta;
            j.alpha[jj] += delta;

            if (diff > 0.0) {
                // Region III
                if (j.alpha[jj] < 0.0) {
                    j.alpha[jj] = 0.0;
                    i.alpha[ii] = diff;
                }
            } else {
                // Region IV
                if (i.alpha[ii] < 0.0) {
                    i.alpha[ii] = 0.0;
                    j.alpha[jj] = -diff;
                }
            }
            
            if (diff > i.C - j.C) {
                // Region I
                if (i.alpha[ii] > i.C) {
                    i.alpha[ii] = i.C;
                    j.alpha[jj] = i.C - diff;
                }
            } else {
                // Region II
                if (j.alpha[jj] > j.C) {
                    j.alpha[jj] = j.C;
                    i.alpha[ii] = j.C + diff;
                }
            }
        } else {

            double delta = (i.g[ii] - j.g[jj]) / curv;
            double sum = i.alpha[ii] + j.alpha[jj];
            i.alpha[ii] -= delta;
            j.alpha[jj] += delta;

            if (sum > i.C) {
                if (i.alpha[ii] > i.C) {
                    i.alpha[ii] = i.C;
                    j.alpha[jj] = sum - i.C;
                }
            } else {
                if (j.alpha[jj] < 0) {
                    j.alpha[jj] = 0;
                    i.alpha[ii] = sum;
                }
            }

            if (sum > j.C) {
                if (j.alpha[jj] > j.C) {
                    j.alpha[jj] = j.C;
                    i.alpha[ii] = sum - j.C;
                }
            } else {
                if (i.alpha[ii] < 0) {
                    i.alpha[ii] = 0.0;
                    j.alpha[jj] = sum;
                }
            }
        }
        
        double delta_alpha_i = i.alpha[ii] - old_alpha_i;
        double delta_alpha_j = j.alpha[jj] - old_alpha_j;

        int si = 2 * ii - 1;
        int sj = 2 * jj - 1;
        for (int k = 0; k < sv.size(); k++) {
            SupportVector v = sv.get(k);
            v.g[0] -= si * i.kcache.get(k) * delta_alpha_i + sj * j.kcache.get(k) * delta_alpha_j;
            v.g[1] += si * i.kcache.get(k) * delta_alpha_i + sj * j.kcache.get(k) * delta_alpha_j;
        }

        // optimality test
        minmax();
        b = -(gmax + gmin) / 2;
        if (gmax - gmin < epsgr) {
            return false;
        }

        return true;
    }
852580
function(d) {
    return d.wm || (d.wm =
      (later.D.val(d) +
      (later.dw.val(later.M.start(d)) - 1) + (7 - later.dw.val(d))) / 7);
  }
800600
function diagonal(s, d) {
      return `M ${s.x} ${s.y}
              C ${(s.x)} ${(s.y + d.y)/2},
                ${(s.x)} ${(s.y + d.y)/2},
                ${d.x} ${d.y}`;
    }
1186035
public static void multAddOuter( double alpha , DMatrix6x6 A , double beta , DMatrix6 u , DMatrix6 v , DMatrix6x6 C ) {
        C.a11 = alpha*A.a11 + beta*u.a1*v.a1;
        C.a12 = alpha*A.a12 + beta*u.a1*v.a2;
        C.a13 = alpha*A.a13 + beta*u.a1*v.a3;
        C.a14 = alpha*A.a14 + beta*u.a1*v.a4;
        C.a15 = alpha*A.a15 + beta*u.a1*v.a5;
        C.a16 = alpha*A.a16 + beta*u.a1*v.a6;
        C.a21 = alpha*A.a21 + beta*u.a2*v.a1;
        C.a22 = alpha*A.a22 + beta*u.a2*v.a2;
        C.a23 = alpha*A.a23 + beta*u.a2*v.a3;
        C.a24 = alpha*A.a24 + beta*u.a2*v.a4;
        C.a25 = alpha*A.a25 + beta*u.a2*v.a5;
        C.a26 = alpha*A.a26 + beta*u.a2*v.a6;
        C.a31 = alpha*A.a31 + beta*u.a3*v.a1;
        C.a32 = alpha*A.a32 + beta*u.a3*v.a2;
        C.a33 = alpha*A.a33 + beta*u.a3*v.a3;
        C.a34 = alpha*A.a34 + beta*u.a3*v.a4;
        C.a35 = alpha*A.a35 + beta*u.a3*v.a5;
        C.a36 = alpha*A.a36 + beta*u.a3*v.a6;
        C.a41 = alpha*A.a41 + beta*u.a4*v.a1;
        C.a42 = alpha*A.a42 + beta*u.a4*v.a2;
        C.a43 = alpha*A.a43 + beta*u.a4*v.a3;
        C.a44 = alpha*A.a44 + beta*u.a4*v.a4;
        C.a45 = alpha*A.a45 + beta*u.a4*v.a5;
        C.a46 = alpha*A.a46 + beta*u.a4*v.a6;
        C.a51 = alpha*A.a51 + beta*u.a5*v.a1;
        C.a52 = alpha*A.a52 + beta*u.a5*v.a2;
        C.a53 = alpha*A.a53 + beta*u.a5*v.a3;
        C.a54 = alpha*A.a54 + beta*u.a5*v.a4;
        C.a55 = alpha*A.a55 + beta*u.a5*v.a5;
        C.a56 = alpha*A.a56 + beta*u.a5*v.a6;
        C.a61 = alpha*A.a61 + beta*u.a6*v.a1;
        C.a62 = alpha*A.a62 + beta*u.a6*v.a2;
        C.a63 = alpha*A.a63 + beta*u.a6*v.a3;
        C.a64 = alpha*A.a64 + beta*u.a6*v.a4;
        C.a65 = alpha*A.a65 + beta*u.a6*v.a5;
        C.a66 = alpha*A.a66 + beta*u.a6*v.a6;
    }
1102730
def _D(self, R):
        """"""
        return (0.3 * self._HNn**2 / self._sin_alpha / R
                + self._HNn + R * self._sin_alpha) / (0.3 * self._HNn + R * self._sin_alpha)
423067
def dcmdottoang_vel(R,Rdot):
    
    w = vee_map(Rdot.dot(R.T))
    Omega = vee_map(R.T.dot(Rdot))

    return (w, Omega)
337141
def Kdiag(self, X, target):
        """"""
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU
        
        k1 = (2*lu+ly)/(lu+ly)**2
        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
172531
function ecp_sign(r, s, p, a, G, n, d, z, l) {
	var R = ecp_new(l), S = ecp_new(l), k = mpn_new(l);
	for (;;) {
		for (var i = 0; i < l; ++i) {
			k[i] = Math.random() * 0x80000000 | 0;
		}
		var r0 = ecp_proj(R, ecp_mul(S, k, G, a, p, l), p, l)[0];
		if (mpn_cmp(r0, n, l) >= 0) {
			mpn_sub(r0, r0, n, l);
		}
		if (!mpn_zero_p(r0, l)) {
			var t0 = _ecp_t0, t1 = _ecp_t1;
			fp_mul(s, fp_inv(t0, k, n, l), fp_add(t1, z, fp_mul(t1, r0, d, n, l), n, l), n, l);
			if (!mpn_zero_p(s, l)) {
				mpn_copyi(r, r0, l);
				break;
			}
		}
	}
}
831455
def d(self,value):
        
        assert value.shape[0]==self._P*self._N, 'd dimension mismatch'
        self._d = value
        self.clear_cache('Yhat','Xhat','Areml','Areml_eigh','Areml_chol','Areml_inv','beta_hat','B_hat',
                         'LRLdiag_Xhat_tens','LRLdiag_Yhat','Areml_grad',
                         'beta_grad','Xstar_beta_grad','Zstar','DLZ')
1208940
def det4D(m):
    '''
    
    '''
    # I just solved this in Mathematica, copy-pasted, and replaced the string '] m' with ']*m':
    # Mathematica code: Det@Table[m[i][j], {i, 0, 3}, {j, 0, 3}]
    return (m[0][3]*m[1][2]*m[2][1]*m[3][0] - m[0][2]*m[1][3]*m[2][1]*m[3][0] -
            m[0][3]*m[1][1]*m[2][2]*m[3][0] + m[0][1]*m[1][3]*m[2][2]*m[3][0] +
            m[0][2]*m[1][1]*m[2][3]*m[3][0] - m[0][1]*m[1][2]*m[2][3]*m[3][0] -
            m[0][3]*m[1][2]*m[2][0]*m[3][1] + m[0][2]*m[1][3]*m[2][0]*m[3][1] +
            m[0][3]*m[1][0]*m[2][2]*m[3][1] - m[0][0]*m[1][3]*m[2][2]*m[3][1] -
            m[0][2]*m[1][0]*m[2][3]*m[3][1] + m[0][0]*m[1][2]*m[2][3]*m[3][1] +
            m[0][3]*m[1][1]*m[2][0]*m[3][2] - m[0][1]*m[1][3]*m[2][0]*m[3][2] -
            m[0][3]*m[1][0]*m[2][1]*m[3][2] + m[0][0]*m[1][3]*m[2][1]*m[3][2] +
            m[0][1]*m[1][0]*m[2][3]*m[3][2] - m[0][0]*m[1][1]*m[2][3]*m[3][2] -
            m[0][2]*m[1][1]*m[2][0]*m[3][3] + m[0][1]*m[1][2]*m[2][0]*m[3][3] +
            m[0][2]*m[1][0]*m[2][1]*m[3][3] - m[0][0]*m[1][2]*m[2][1]*m[3][3] -
            m[0][1]*m[1][0]*m[2][2]*m[3][3] + m[0][0]*m[1][1]*m[2][2]*m[3][3])
520442
function findx0BigA(a, p, q) {
      var w, s, D, B, u, z, zbar, logpg;

      // snTerm -- see formula 34
      logpg = Math.log(p * gamma(a + 1));
      function snTerm(x) {
         return function(n, v) {
            return n === 0 ? 1 : v * x / (a + n);
         };
      }
      // See formula 34, where `n` is off by 1
      function f(x, n) {
         return Math.exp((logpg + x - Math.log(series(snTerm(x), n + 1))) / a);
      }
      B = q * gamma(a);
      D = Math.max(2, a * (a - 1));
      s = funS(p, q);
      w = a + s * Math.sqrt(a) + (s * s - 1) / 3 +
         (s * s * s - 7 * s) / (36 * Math.sqrt(a)) -
         (3 * s * s * s * s + 7 * s * s - 16) / (810 * a) +
         (9 * s * s * s * s * s + 256 * s * s * s - 433 * s) /
            (38880 * a * Math.sqrt(a));
      if (a >= 500 && Math.abs(1 - w / a) < 1e-6) { return w; }
      if (p > 0.5) {
         if (w < 3 * a) { return w; }
         if (B < Math.pow(10, -D)) { return findx0TinyB(a, -Math.log(B)); }
         u = -Math.log(B) + (a - 1) * Math.log(w) -
               Math.log(1 + (1 - a) / (1 + w));
         // formula 33
         return -Math.log(B) + (a - 1) * Math.log(u) -
               Math.log(1 + (1 - a) / (1 + u));
      }
      // handle p <= 0.5
      z = w > 0.15 * (a + 1) ? w : f(f(f(f(w, 1), 3), 3), 4);
      if (z <= 0.01 * (a + 1) || z > 0.7 * (a + 1)) { return z; }
      /* eslint-disable no-extra-parens */
      return (function() {
         // find N for formula 36
         var N, term;

         N = 0;
         term = 1;
         do {
            N += 1;
            term = snTerm(z)(N, term);
         } while (term > 1e-4);
         zbar = f(z, N + 1);
         return zbar * (1 - (a * Math.log(zbar) - zbar - logpg +
                  Math.log(series(snTerm(z), N + 1))) / (a - zbar));
      }());
      /* eslint-enable no-extra-parens */
   }
27663
def mktk03(terms, seed, G2, G3):
    
# random.seed(n)
    p = 0
    n = seed
    gh = []
    g10, sfact, afact = -18e3, 3.8, 2.4
    g20 = G2 * g10
    g30 = G3 * g10
    alpha = g10/afact
    s1 = s_l(1, alpha)
    s10 = sfact * s1
    gnew = random.normal(g10, s10)
    if p == 1:
        print(1, 0, gnew, 0)
    gh.append(gnew)
    gh.append(random.normal(0, s1))
    gnew = gh[-1]
    gh.append(random.normal(0, s1))
    hnew = gh[-1]
    if p == 1:
        print(1, 1, gnew, hnew)
    for l in range(2, terms + 1):
        for m in range(l + 1):
            OFF = 0.0
            if l == 2 and m == 0:
                OFF = g20
            if l == 3 and m == 0:
                OFF = g30
            s = s_l(l, alpha)
            j = (l - m) % 2
            if j == 1:
                s = s * sfact
            gh.append(random.normal(OFF, s))
            gnew = gh[-1]
            if m == 0:
                hnew = 0
            else:
                gh.append(random.normal(0, s))
                hnew = gh[-1]
            if p == 1:
                print(l, m, gnew, hnew)
    return gh
888467
private int ssMedian3(int Td, int PA, int v1, int v2, int v3) {
        if (T[start + Td + SA[PA + SA[v1]]] > T[start + Td + SA[PA + SA[v2]]]) {
            int t = v1;
            v1 = v2;
            v2 = t;
        }
        if (T[start + Td + SA[PA + SA[v2]]] > T[start + Td + SA[PA + SA[v3]]]) {
            if (T[start + Td + SA[PA + SA[v1]]] > T[start + Td + SA[PA + SA[v3]]]) {
                return v1;
            } else {
                return v3;
            }
        }
        return v2;
    }