compute (p_{i|j} + p_{j|i})/(2n)
356698
function (D, perplexity, tol) {
		var Nf = Math.sqrt(D.length); // this better be an integer
		var N = Math.floor(Nf);
		assert(N === Nf, "D should have square number of elements.");
		var Htarget = Math.log(perplexity); // target entropy of distribution
		var P = zeros(N * N); // temporary probability matrix

		var prow = zeros(N); // a temporary storage compartment
		for (var i = 0; i < N; i++) {
			var betamin = -Infinity;
			var betamax = Infinity;
			var beta = 1; // initial value of precision
			var done = false;
			var maxtries = 50;

			// perform binary search to find a suitable precision beta
			// so that the entropy of the distribution is appropriate
			var num = 0;
			while (!done) {
				//debugger;

				// compute entropy and kernel row with beta precision
				var psum = 0.0;
				for (var j = 0; j < N; j++) {
					var pj = Math.exp(-D[i * N + j] * beta);
					if (i === j) {
						pj = 0;
					} // we dont care about diagonals
					prow[j] = pj;
					psum += pj;
				}
				// normalize p and compute entropy
				var Hhere = 0.0;
				for (var j = 0; j < N; j++) {
					var pj = prow[j] / psum;
					prow[j] = pj;
					if (pj > 1e-7) Hhere -= pj * Math.log(pj);
				}

				// adjust beta based on result
				if (Hhere > Htarget) {
					// entropy was too high (distribution too diffuse)
					// so we need to increase the precision for more peaky distribution
					betamin = beta; // move up the bounds
					if (betamax === Infinity) {
						beta = beta * 2;
					}
					else {
						beta = (beta + betamax) / 2;
					}

				} else {
					// converse case. make distrubtion less peaky
					betamax = beta;
					if (betamin === -Infinity) {
						beta = beta / 2;
					}
					else {
						beta = (beta + betamin) / 2;
					}
				}

				// stopping conditions: too many tries or got a good precision
				num++;
				if (Math.abs(Hhere - Htarget) < tol) {
					done = true;
				}
				if (num >= maxtries) {
					done = true;
				}
			}

			// console.log('data point ' + i + ' gets precision ' + beta + ' after ' + num + ' binary search steps.');
			// copy over the final prow to P at row i
			for (var j = 0; j < N; j++) {
				P[i * N + j] = prow[j];
			}

		} // end loop over examples i

		// symmetrize P and normalize it to sum to 1 over all ij
		var Pout = zeros(N * N);
		var N2 = N * 2;
		for (var i = 0; i < N; i++) {
			for (var j = 0; j < N; j++) {
				Pout[i * N + j] = Math.max((P[i * N + j] + P[j * N + i]) / N2, 1e-100);
			}
		}

		return Pout;
	}
1276182
function ( t, p0, p1, p2, p3 ) {

		return - 3 * p0 * ( 1 - t ) * ( 1 - t )  +
			3 * p1 * ( 1 - t ) * ( 1 - t ) - 6 * t * p1 * ( 1 - t ) +
			6 * t *  p2 * ( 1 - t ) - 3 * t * t * p2 +
			3 * t * t * p3;

	}
275521
def sbesselh1(x, N):
    ""
    
    jn = sbesselj(x, N)
    yn = sbessely(x, N)

    return jn + 1j * yn
719738
public static int cs_cumsum(int[] p, int[] c, int n) {
        int i, nz = 0;
        float nz2 = 0;
        if (p == null || c == null)
            return (-1); /* check inputs */
        for (i = 0; i < n; i++) {
            p[i] = nz;
            nz += c[i];
            nz2 += c[i]; /* also in float to avoid int overflow */
            c[i] = p[i]; /* also copy p[0.f.n-1] back into c[0.f.n-1]*/
        }
        p[n] = nz;
        return (int) nz2; /* return sum (c [0.f.n-1]) */
    }
537472
private static void compute(int[] in, int[] out, int width, int height, int radius)
    {
        final int widthMinus1 = width - 1;
        final int tableSize = 2 * radius + 1;
        final int[] divide = new int[Constant.UNSIGNED_BYTE * tableSize];

        for (int i = 0; i < Constant.UNSIGNED_BYTE * tableSize; i++)
        {
            divide[i] = i / tableSize;
        }

        int inIndex = 0;
        for (int y = 0; y < height; y++)
        {
            compute(in, out, y, inIndex, divide, width, widthMinus1, height, radius);
            inIndex += width;
        }
    }
1264657
def pow2(x: int, p: int) -> int:
    """"""
    while p > 0:
        x = x * x % q
        p -= 1
    return x
622140
function (i, p) {
  var k = Math.floor(Math.log(i) / Math.log(2))
    , r = i - Math.pow(2, k);
  if (r === 0) {
    return Math.pow(2, p - k);
  }
  var nr = (i - 2*r).toString(2).split('').reverse().join('');
  return (parseInt(nr, 2) << p - nr.length) + Math.pow(2, p - k - 1);
}
623990
def IJ(mu, N):
    """"""
    if mu == 0: return 1, 1, 1

    if mu not in range(0, N**2):
        raise ValueError('mu has an invalid value mu='+str(mu)+'.')

    if 1 <= mu <= N-1:
        return mu+1, mu+1, 1
    else:
        m = N-1
        M = N*(N-1)/2
        for jj in range(1, N):
            for ii in range(jj+1, N+1):
                m += 1
                if m == mu or m+M == mu:
                    if mu > N*(N+1)/2 - 1:
                        return ii, jj, -1
                    else:
                        return ii, jj, 1
294129
function mpn_shr1(r, n, l) {
	var c = 0;
	while (l > 0) {
		var w = n[--l] | c;
		c = w << 31;
		r[l] = w >>> 1;
	}
	return c;
}
981499
protected void computeNumericalJacobian( DenseMatrix64F param ,
                                             DenseMatrix64F pt ,
                                             DenseMatrix64F deriv )
    {
        double invDelta = 1.0/DELTA;

        func.compute(param,pt, temp0);

        // compute the jacobian by perturbing the parameters slightly
        // then seeing how it effects the results.
        for( int i = 0; i < param.numRows; i++ ) {
            param.data[i] += DELTA;
            func.compute(param,pt, temp1);
            // compute the difference between the two parameters and divide by the delta
            add(invDelta,temp1,-invDelta,temp0,temp1);
            // copy the results into the jacobian matrix
            System.arraycopy(temp1.data,0,deriv.data,i*pt.numRows,pt.numRows);

            param.data[i] -= DELTA;
        }
    }
113298
def logn2(n, p):
    """"""
    with precision(p):
        extra = 10
        while True:
            with precision(p+extra):
                # use extra precision for intermediate step
                log2upper = log2(n, RoundTowardPositive)
                log2lower = log2(n, RoundTowardNegative)

            lower = div(1, log2upper, RoundTowardNegative)
            upper = div(1, log2lower, RoundTowardPositive)

            # if lower and upper are adjacent (or equal) we're done
            if next_up(lower) == upper:
                return (Fraction(*lower.as_integer_ratio()),
                        Fraction(*upper.as_integer_ratio()))

            # otherwise, increase the precision and try again
            extra += 10
376070
def p_jr(p):
    
    p[4] = Expr.makenode(Container('-', p.lineno(3)), p[4], Expr.makenode(Container(MEMORY.org + 2, p.lineno(1))))
    p[0] = Asm(p.lineno(1), 'JR %s,N' % p[2], p[4])
440898
function ecp_mul(R, n1, N2, a, p, l) {
	var add = mpn_one_p(N2[2], l) ? ecp_add_aff : ecp_add;
	var active = false, swaps = 0, S = _ecp_S, T;
	for (var i = l; i > 0;) {
		var w = n1[--i];
		for (var j = 31; j > 0; --j) {
			if (active) {
				ecp_dbl(S, R, a, p, l);
				T = S, S = R, R = T, ++swaps;
			}
			if ((w <<= 1) < 0) {
				if (active) {
					add(S, R, N2, a, p, l);
					T = S, S = R, R = T, ++swaps;
				}
				else {
					ecp_copy(R, N2, l);
					active = true;
				}
			}
		}
	}
	if (swaps & 1) {
		return ecp_copy(S, R, l);
	}
	return R;
}
25987
private static void mult(int[] src, int srcLen, int v0, int v1, int[] dst) {
        long v = v0 & LONG_MASK;
        long carry = 0;
        for (int j = 0; j < srcLen; j++) {
            long product = v * (src[j] & LONG_MASK) + carry;
            dst[j] = (int) product;
            carry = product >>> 32;
        }
        dst[srcLen] = (int) carry;
        v = v1 & LONG_MASK;
        carry = 0;
        for (int j = 0; j < srcLen; j++) {
            long product = (dst[j + 1] & LONG_MASK) + v * (src[j] & LONG_MASK) + carry;
            dst[j + 1] = (int) product;
            carry = product >>> 32;
        }
        dst[srcLen + 1] = (int) carry;
    }
809123
def computeA(self, j, timeSeries):
        
        seasonLength = self.get_parameter("seasonLength")
        A_j = 0
        for i in range(seasonLength):
            A_j += timeSeries[(seasonLength * (j)) + i][1]
        return A_j / seasonLength
813860
function fp_dbl(r, n, p, l) {
	if (mpn_dbl(r, n, l) || mpn_cmp(r, p, l) >= 0) {
		mpn_sub(r, r, p, l);
	}
	return r;
}
1246992
function(p) {
            var m = this.m;
            return {
                x: m[0] * p.x + m[2] * p.y + m[4],
                y: m[1] * p.x + m[3] * p.y + m[5]
            };
        }