https://en.wikipedia.org/wiki/Spouge%27s_approximation
773959
public static BigDecimal factorialUsingSpougeCached(BigDecimal x, MathContext mathContext) {
		MathContext mc = new MathContext(mathContext.getPrecision() * 2, mathContext.getRoundingMode());

		int a = mathContext.getPrecision() * 13 / 10;
		List<BigDecimal> constants = getSpougeFactorialConstants(a);

		BigDecimal bigA = BigDecimal.valueOf(a);

		boolean negative = false;
		BigDecimal factor = constants.get(0);
		for (int k = 1; k < a; k++) {
			BigDecimal bigK = BigDecimal.valueOf(k);
			factor = factor.add(constants.get(k).divide(x.add(bigK), mc), mc);
			negative = !negative;
		}

		BigDecimal result = pow(x.add(bigA, mc), x.add(BigDecimal.valueOf(0.5), mc), mc);
		result = result.multiply(exp(x.negate().subtract(bigA, mc), mc), mc);
		result = result.multiply(factor, mc);

		return result.round(mathContext);
	}
39917
def _basis_spline_factory(coef, degree, knots, der, ext):
        """"""
        return functools.partial(interpolate.splev, tck=(knots, coef, degree), der=der, ext=ext)
102025
def ρ(self, e):
        
        d = self._dichalcogenide
        at, λ = d.at, d.λ
        return abs(2 * e - λ) * (2 * at**2)**(-1)
1372707
def approximate_split(x, num_splits, axis=0):
  
  size = shape_list(x)[axis]
  size_splits = [tf.div(size + i, num_splits) for i in range(num_splits)]
  return tf.split(x, size_splits, axis=axis)
561948
function spfa(graph, s) {
  const distance = {};
  const previous = {};
  const queue = {};
  const isInQue = {};
  const cnt = {};
  let head = 0;
  let tail = 1;
  // initialize
  distance[s] = 0;
  queue[0] = s;
  isInQue[s] = true;
  cnt[s] = 1;
  graph.vertices.forEach(v => {
    if (v !== s) {
      distance[v] = Infinity;
      isInQue[v] = false;
      cnt[v] = 0;
    }
  });

  let currNode;
  while (head !== tail) {
    currNode = queue[head++];
    isInQue[currNode] = false;
    const neighbors = graph.neighbors(currNode);
    for (let i = 0; i < neighbors.length; i++) {
      const v = neighbors[i];
      // relaxation
      const newDistance = distance[currNode] + graph.edge(currNode, v);
      if (newDistance < distance[v]) {
        distance[v] = newDistance;
        previous[v] = currNode;
        if (!isInQue[v]) {
          queue[tail++] = v;
          isInQue[v] = true;
          cnt[v]++;
          if (cnt[v] > graph.vertices.size) {
            // indicates negative-weighted cycle
            return {
              distance: {}
            };
          }
        }
      }
    }
  }

  return {
    distance,
    previous
  };
}
1316092
def bound_perturbed_gmres(pseudo, p, epsilon, deltas):
    '''
    '''
    if not numpy.all(numpy.array(deltas) > epsilon):
        raise ArgumentError('all deltas have to be greater than epsilon')

    bound = []
    for delta in deltas:
        # get boundary paths
        paths = pseudo.contour_paths(delta)

        # get vertices on boundary
        vertices = paths.vertices()

        # evaluate polynomial
        supremum = numpy.max(numpy.abs(p(vertices)))

        # compute bound
        bound.append(epsilon/(delta-epsilon)
                     * paths.length()/(2*numpy.pi*delta)
                     * supremum)
    return bound
417260
function approximate_unit_arc(theta1, delta_theta) {
  var alpha = 4/3 * Math.tan(delta_theta/4);

  var x1 = Math.cos(theta1);
  var y1 = Math.sin(theta1);
  var x2 = Math.cos(theta1 + delta_theta);
  var y2 = Math.sin(theta1 + delta_theta);

  return [ x1, y1, x1 - y1*alpha, y1 + x1*alpha, x2 + y2*alpha, y2 - x2*alpha, x2, y2 ];
}
1174634
private ApproximationLine approximateLowerHull(ConvexHull convexHull, double[] log_k, double sum_log_k, double sum_log_k2, double[] log_kDist, double sum_log_kDist, double sum_log_k_kDist) {
    // StringBuilder msg = new StringBuilder(1000);
    int[] lowerHull = convexHull.getLowerHull();
    int l = convexHull.getNumberOfPointsInLowerHull();
    int k_0 = settings.kmax - lowerHull.length + 1;

    // linear search on all line segments on the lower convex hull
    // msg.append("lower hull l = ").append(l).append('\n');
    double low_error = Double.MAX_VALUE;
    double low_m = 0.0;
    double low_t = 0.0;

    for(int i = 1; i < l; i++) {
      double cur_m = (log_kDist[lowerHull[i]] - log_kDist[lowerHull[i - 1]]) / (log_k[lowerHull[i]] - log_k[lowerHull[i - 1]]);
      double cur_t = log_kDist[lowerHull[i]] - cur_m * log_k[lowerHull[i]];
      double cur_error = ssqerr(k_0, settings.kmax, log_k, log_kDist, cur_m, cur_t);
      // msg.append(" Segment = ").append(i).append(" m =
      // ").append(cur_m).append(" t = ").append(cur_t).append(" lowerror =
      // ").append(cur_error).append('\n');
      if(cur_error < low_error) {
        low_error = cur_error;
        low_m = cur_m;
        low_t = cur_t;
      }
    }

    // linear search on all points of the lower convex hull
    boolean is_right = true; // NEEDED FOR PROOF CHECK
    for(int i = 0; i < l; i++) {
      double cur_m = optimize(k_0, settings.kmax, sum_log_k, sum_log_k2, log_k[lowerHull[i]], log_kDist[lowerHull[i]], sum_log_k_kDist, sum_log_kDist);
      double cur_t = log_kDist[lowerHull[i]] - cur_m * log_k[lowerHull[i]];
      // only valid if both neighboring points are underneath y=mx+t
      if((i == 0 || log_kDist[lowerHull[i - 1]] >= log_kDist[lowerHull[i]] - cur_m * (log_k[lowerHull[i]] - log_k[lowerHull[i - 1]])) && (i == l - 1 || log_kDist[lowerHull[i + 1]] >= log_kDist[lowerHull[i]] + cur_m * (log_k[lowerHull[i + 1]] - log_k[lowerHull[i]]))) {
        double cur_error = ssqerr(k_0, settings.kmax, log_k, log_kDist, cur_m, cur_t);
        if(cur_error < low_error) {
          low_error = cur_error;
          low_m = cur_m;
          low_t = cur_t;
        }
      }

      // check proof of bisection search
      if(!(i > 0 && log_kDist[lowerHull[i - 1]] < log_kDist[lowerHull[i]] - cur_m * (log_k[lowerHull[i]] - log_k[lowerHull[i - 1]])) && !is_right) {
        LOG.warning("ERROR lower: The bisection search will not work properly!");
        if(!(i < l - 1 && log_kDist[lowerHull[i + 1]] < log_kDist[lowerHull[i]] + cur_m * (log_k[lowerHull[i + 1]] - log_k[lowerHull[i]]))) {
          is_right = false;
        }
      }
    }
    return new ApproximationLine(k_0, low_m, low_t);
  }
104439
def FZStaeckel(u,v,pot,delta): #pragma: no cover because unused
    
    R,z= bovy_coords.uv_to_Rz(u,v,delta=delta)
    return _evaluatezforces(pot,R,z)
648614
def AcosTriAngles(v, f, normalize):
    

    if normalize:
        nm = lambda x : NormalizedNx3(x)
    else:
        nm = lambda x : x

    return Ch(lambda v :
        Sum3xN(NormalizedNx3(TriEdges(f, 1, 0, nm(v))) * NormalizedNx3(TriEdges(f, 2, 0, nm(v)))) &
        Sum3xN(NormalizedNx3(TriEdges(f, 2, 1, nm(v))) * NormalizedNx3(TriEdges(f, 0, 1, nm(v)))) &
        Sum3xN(NormalizedNx3(TriEdges(f, 0, 2, nm(v))) * NormalizedNx3(TriEdges(f, 1, 2, nm(v)))))
662862
function maybeScaleToFit (G, nested) {
    if (scale !== null) return
    const maxValue = sum(nested.bandValues)
    if (maxValue <= 0) {
      scale = 1
    } else {
      scale = (y1 - y0) / maxValue
      if (whitespace !== 1) scale *= (1 - whitespace)
    }
  }
418781
def _approx_equal(x, y, eps):
    
    if x.space != y.space:
        return False

    if x is y:
        return True

    try:
        return x.dist(y) <= eps
    except NotImplementedError:
        try:
            return x == y
        except NotImplementedError:
            return False
82162
def von_mises_strain(self):
        
        eps = self - 1/3 * np.trace(self) * np.identity(3)

        return np.sqrt(np.sum(eps * eps) * 2/3)
1167095
def g_reuss(self):
        
        return 15. / (8. * self.compliance_tensor.voigt[:3, :3].trace() -
                      4. * np.triu(self.compliance_tensor.voigt[:3, :3]).sum() +
                      3. * self.compliance_tensor.voigt[3:, 3:].trace())
988236
def _laplace_approximation_newton(self, tol=1e-6, stepsz=0.9, verbose=False):
        
        from pylds.util import solve_symm_block_tridiag, scipy_solve_symm_block_tridiag
        scale = self.T * self.D_emission

        def newton_step(x, stepsz):
            assert 0 <= stepsz <= 1
            g = self.gradient_log_joint(x)
            H_diag, H_upper_diag = self.sparse_hessian_log_joint(x)
            Hinv_g = -scipy_solve_symm_block_tridiag(-H_diag / scale,
                                               -H_upper_diag / scale,
                                               g / scale)
            return x - stepsz * Hinv_g

        if verbose:
            print("Fitting Laplace approximation")

        itr = [0]
        def cbk(x):
            print("Iteration: ", itr[0],
                  "\tObjective: ", (self.log_joint(x) / scale).round(4),
                  "\tAvg Grad: ", (self.gradient_log_joint(x).mean() / scale).round(4))
            itr[0] += 1

        # Solve for optimal x with Newton's method
        x = self.gaussian_states
        dx = np.inf
        while dx >= tol:
            xnew = newton_step(x, stepsz)
            dx = np.mean(abs(xnew - x))
            x = xnew

            if verbose:
                cbk(x)

        assert np.all(np.isfinite(x))
        if verbose:
            print("Done")

        return x
1227524
private double hermiteInterpolate(final double Y0, final double Y1, final double Y2, final double Y3,
                                      final double MU, final double TENSION, final double BIAS) {
        double m0;
        double m1;
        final double MU2;
        final double Mu3;
        final double A0;
        final double A1;
        final double A2;
        final double A3;

        MU2 = MU * MU;
        Mu3 = MU2 * MU;
        m0 = (Y1 - Y0) * (1 + BIAS) * (1 - TENSION) / 2;
        m0 += (Y2 - Y1) * (1 - BIAS) * (1 - TENSION) / 2;
        m1 = (Y2 - Y1) * (1 + BIAS) * (1 - TENSION) / 2;
        m1 += (Y3 - Y2) * (1 - BIAS) * (1 - TENSION) / 2;
        A0 = (2 * Mu3) - (3 * MU2) + 1;
        A1 = Mu3 - (2 * MU2) + MU;
        A2 = Mu3 - MU2;
        A3 = (-2 * Mu3) + (3 * MU2);

        return ((A0 * Y1) + (A1 * m0) + (A2 * m1) + (A3 * Y2));
    }
553865
def _get_assumptions(t):
        

        if t.op in ('__le__', '__lt__', 'ULE', 'ULT'):
            return [ t.args[0] >= 0 ]
        elif t.op in ('__ge__', '__gt__', 'UGE', 'UGT'):
            return [ t.args[0] <= 2**len(t.args[0])-1 ]
        elif t.op in ('SLE', 'SLT'):
            return [ _all_operations.SGE(t.args[0], -(1 << (len(t.args[0])-1))) ]
        elif t.op in ('SGE', 'SGT'):
            return [ _all_operations.SLE(t.args[0], (1 << (len(t.args[0])-1)) - 1) ]
        else:
            return [ ]