Return the value of the second kind Bessel function of order ``n`` at
    ``x``.

    ``n`` should be a Python integer.
426514
def yn(n, x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_yn,
        (n, BigFloat._implicit_convert(x)),
        context,
    )
584334
def get_a(values, n):
        
        m = len(values)-n
        a = numpy.empty((m, n), dtype=float)
        for i in range(m):
            i0 = i-1 if i > 0 else None
            i1 = i+n-1
            a[i] = values[i1:i0:-1]
        return numpy.array(a)
1090324
def eint(x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_eint,
        (BigFloat._implicit_convert(x),),
        context,
    )
1256598
def fib_iter(n):
    

    # precondition
    assert n >= 0, 'n must be positive integer'

    fib_1 = 0
    fib_2 = 1
    sum = 0
    if n <= 1:
        return n
    for _ in range(n-1):
        sum = fib_1 + fib_2
        fib_1 = fib_2
        fib_2 = sum
    return sum
1325310
def sbessely(x, N):
    

    out = np.zeros(N, dtype=np.float64)

    out[0] = -np.cos(x) / x
    out[1] = -np.cos(x) / (x ** 2) - np.sin(x) / x

    for n in xrange(2, N):
        out[n] = ((2.0 * n - 1.0) / x) * out[n - 1] - out[n - 2]

    return out
604458
def set_trapezoidal_integration(self, n):
        
        if not (n >= 2):
            raise ValueError('must have n >= 2; got %r' % (n,))

        self.in_vals[IN_VAL_INTEG_METH] = n + 1
        return self
1161757
def rec_sqrt(x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_rec_sqrt,
        (BigFloat._implicit_convert(x),),
        context,
    )
949301
def get_n_r(self):
        """"""
        return math.floor(self.n / self.eta**self.i + _epsilon), math.floor(self.r * self.eta**self.i + _epsilon)
599086
def function(x, ax, ay):
    '''
    
    '''
    with np.errstate(invalid='ignore'):
        return ay * (x - ax)**0.5
154787
def sample_function(self, f, N):
        
        x = self.interpolation_points(N+1)
        return f(x)
954385
def get_X0(X):
    
    if pandas_available and isinstance(X, pd.DataFrame):
        assert len(X) == 1
        x = np.array(X.iloc[0])
    else:
        x, = X
    return x
43841
def get_random_integer(N, randfunc=None):
    
    if randfunc is None:
        randfunc = Random.new().read

    S = randfunc(N>>3)
    odd_bits = N % 8
    if odd_bits != 0:
        char = ord(randfunc(1)) >> (8-odd_bits)
        S = bchr(char) + S
    value = bytes_to_long(S)
    return value
166040
def _t_of_e(self, a0=None, t_start=None, f0=None, ef=None, t_obs=5.0):
        
        if ef is None:
            ef = np.ones_like(self.e0)*0.0000001

        beta = 64.0/5.0*self.m1*self.m2*(self.m1+self.m2)

        e_vals = np.asarray([np.linspace(ef[i], self.e0[i], self.num_points)
                            for i in range(len(self.e0))])
        integrand = self._find_integrand(e_vals)
        integral = np.asarray([np.trapz(integrand[:, i:], x=e_vals[:, i:])
                              for i in range(e_vals.shape[1])]).T

        if a0 is None and f0 is None:

            a0 = (19./12.*t_start*beta*1/integral[:, 0])**(1./4.) * self._f_e(e_vals[:, -1])

        elif a0 is None:
            a0 = ((self.m1 + self.m2)/self.f0**2)**(1./3.)

        c0 = self._c0_func(a0, self.e0)

        a_vals = c0[:, np.newaxis]*self._f_e(e_vals)

        delta_t = 12./19*c0[:, np.newaxis]**4/beta[:, np.newaxis]*integral

        return e_vals, a_vals, delta_t
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;
    }
650582
def zeta(x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_zeta,
        (BigFloat._implicit_convert(x),),
        context,
    )
902797
def ceilpow2(n):
    """"""
    signif,exponent = frexp(n)
    if (signif < 0):
        return 1;
    if (signif == 0.5):
        exponent -= 1;
    return (1) << exponent;
216075
def _ntz(x):
        
        if x == 0:
            return 0
        y = (~x) & (x - 1)    # There is actually a bug in BAP until 0.8

        def bits(y):
            n = 0
            while y != 0:
                n += 1
                y >>= 1
            return n
        return bits(y)