Python mpmath.mpf方法代码示例

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def compute_d(K, N):
    """d_{k, n} from DLMF 8.12.12"""
    M = N + 2*K
    d0 = [-mp.mpf(1)/3]
    alpha = compute_alpha(M + 2)
    for n in range(1, M):
        d0.append((n + 2)*alpha[n+2])
    d = [d0]
    g = compute_g(K)
    for k in range(1, K):
        dk = []
        for n in range(M - 2*k):
            dk.append((-1)**k*g[k]*d[0][n] + (n + 2)*d[k-1][n+2])
        d.append(dk)
    for k in range(K):
        d[k] = d[k][:N]
    return d 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def gammainc(a, x, dps=50, maxterms=10**8):
    """Compute gammainc exactly like mpmath does but allow for more
    summands in hypercomb. See

    mpmath/functions/expintegrals.py#L134
    
    in the mpmath github repository.

    """
    with mp.workdps(dps):
        z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
        G = [z]
        negb = mp.fneg(b, exact=True)

        def h(z):
            T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
            return (T1,)

        res = mp.hypercomb(h, [z], maxterms=maxterms)
        return mpf2float(res) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def zpkfreqz(z, p, k, worN=None):
    """
    Frequency response of a filter in zpk format, using mpmath.

    This is the same calculation as scipy.signal.freqz, but the input is in
    zpk format, the calculation is performed using mpath, and the results are
    returned in lists instead of numpy arrays.
    """
    if worN is None or isinstance(worN, int):
        N = worN or 512
        ws = [mpmath.pi * mpmath.mpf(j) / N for j in range(N)]
    else:
        ws = worN

    h = []
    for wk in ws:
        zm1 = mpmath.exp(1j * wk)
        numer = _prod([zm1 - t for t in z])
        denom = _prod([zm1 - t for t in p])
        hk = k * numer / denom
        h.append(hk)
    return ws, h 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def test_rf(self):
        if LooseVersion(mpmath.__version__) >= LooseVersion("1.0.0"):
            # no workarounds needed
            mppoch = mpmath.rf
        else:
            def mppoch(a, m):
                # deal with cases where the result in double precision
                # hits exactly a non-positive integer, but the
                # corresponding extended-precision mpf floats don't
                if float(a + m) == int(a + m) and float(a + m) <= 0:
                    a = mpmath.mpf(a)
                    m = int(a + m) - a
                return mpmath.rf(a, m)

        assert_mpmath_equal(sc.poch,
                            mppoch,
                            [Arg(), Arg()],
                            dps=400) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def test_boxcox(self):

        def mp_boxcox(x, lmbda):
            x = mpmath.mp.mpf(x)
            lmbda = mpmath.mp.mpf(lmbda)
            if lmbda == 0:
                return mpmath.mp.log(x)
            else:
                return mpmath.mp.powm1(x, lmbda) / lmbda

        assert_mpmath_equal(sc.boxcox,
                            exception_to_nan(mp_boxcox),
                            [Arg(a=0, inclusive_a=False), Arg()],
                            n=200,
                            dps=60,
                            rtol=1e-13) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def test_boxcox1p(self):

        def mp_boxcox1p(x, lmbda):
            x = mpmath.mp.mpf(x)
            lmbda = mpmath.mp.mpf(lmbda)
            one = mpmath.mp.mpf(1)
            if lmbda == 0:
                return mpmath.mp.log(one + x)
            else:
                return mpmath.mp.powm1(one + x, lmbda) / lmbda

        assert_mpmath_equal(sc.boxcox1p,
                            exception_to_nan(mp_boxcox1p),
                            [Arg(a=-1, inclusive_a=False), Arg()],
                            n=200,
                            dps=60,
                            rtol=1e-13) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def mpf2float(x):
    """
    Convert an mpf to the nearest floating point number. Just using
    float directly doesn't work because of results like this:

    with mp.workdps(50):
        float(mpf("0.99999999999999999")) = 0.9999999999999999

    """
    return float(mpmath.nstr(x, 17, min_fixed=0, max_fixed=0)) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def mp_assert_allclose(res, std, atol=0, rtol=1e-17):
    """
    Compare lists of mpmath.mpf's or mpmath.mpc's directly so that it
    can be done to higher precision than double.

    """
    try:
        len(res)
    except TypeError:
        res = list(res)

    n = len(std)
    if len(res) != n:
        raise AssertionError("Lengths of inputs not equal.")

    failures = []
    for k in range(n):
        try:
            assert_(mpmath.fabs(res[k] - std[k]) <= atol + rtol*mpmath.fabs(std[k]))
        except AssertionError:
            failures.append(k)

    ndigits = int(abs(np.log10(rtol)))
    msg = [""]
    msg.append("Bad results ({} out of {}) for the following points:"
               .format(len(failures), n))
    for k in failures:
        resrep = mpmath.nstr(res[k], ndigits, min_fixed=0, max_fixed=0)
        stdrep = mpmath.nstr(std[k], ndigits, min_fixed=0, max_fixed=0)
        if std[k] == 0:
            rdiff = "inf"
        else:
            rdiff = mpmath.fabs((res[k] - std[k])/std[k])
            rdiff = mpmath.nstr(rdiff, 3)
        msg.append("{}: {} != {} (rdiff {})".format(k, resrep, stdrep, rdiff))
    if failures:
        assert_(False, "\n".join(msg)) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def compute_a(n):
    """a_k from DLMF 5.11.6"""
    a = [mp.sqrt(2)/2]
    for k in range(1, n):
        ak = a[-1]/k
        for j in range(1, len(a)):
            ak -= a[j]*a[-j]/(j + 1)
        ak /= a[0]*(1 + mp.mpf(1)/(k + 1))
        a.append(ak)
    return a 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def lagrange_inversion(a):
    """Given a series

    f(x) = a[1]*x + a[2]*x**2 + ... + a[n-1]*x**(n - 1),

    use the Lagrange inversion formula to compute a series

    g(x) = b[1]*x + b[2]*x**2 + ... + b[n-1]*x**(n - 1)

    so that f(g(x)) = g(f(x)) = x mod x**n. We must have a[0] = 0, so
    necessarily b[0] = 0 too.

    The algorithm is naive and could be improved, but speed isn't an
    issue here and it's easy to read.

    """
    n = len(a)
    f = sum(a[i]*x**i for i in range(len(a)))
    h = (x/f).series(x, 0, n).removeO()
    hpower = [h**0]
    for k in range(n):
        hpower.append((hpower[-1]*h).expand())
    b = [mp.mpf(0)]
    for k in range(1, n):
        b.append(hpower[k].coeff(x, k - 1)/k)
    b = map(lambda x: mp.mpf(x), b)
    return b 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def _binomial_cdf(k, n, p):
    k, n, p = mpmath.mpf(k), mpmath.mpf(n), mpmath.mpf(p)
    if k <= 0:
        return mpmath.mpf(0)
    elif k >= n:
        return mpmath.mpf(1)

    onemp = mpmath.fsub(1, p, exact=True)
    return mpmath.betainc(n - k, k + 1, x2=onemp, regularized=True) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def _f_cdf(dfn, dfd, x):
    if x < 0:
        return mpmath.mpf(0)
    dfn, dfd, x = mpmath.mpf(dfn), mpmath.mpf(dfd), mpmath.mpf(x)
    ub = dfn*x/(dfn*x + dfd)
    res = mpmath.betainc(dfn/2, dfd/2, x2=ub, regularized=True)
    return res 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def _noncentral_chi_cdf(x, df, nc, dps=None):
    if dps is None:
        dps = mpmath.mp.dps
    x, df, nc = mpmath.mpf(x), mpmath.mpf(df), mpmath.mpf(nc)
    with mpmath.workdps(dps):
        res = mpmath.quad(lambda t: _noncentral_chi_pdf(t, df, nc), [0, x])
        return res 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def test_tklmbda_zero_shape(self):
        # When lmbda = 0 the CDF has a simple closed form
        one = mpmath.mpf(1)
        assert_mpmath_equal(
            lambda x: sp.tklmbda(x, 0),
            lambda x: one/(mpmath.exp(-x) + one),
            [Arg()], rtol=1e-7) 

# 需要导入模块: import mpmath [as 别名]
# 或者: from mpmath import mpf [as 别名]
def test_g():
    # Test data for the g_k. See DLMF 5.11.4.
    with mp.workdps(30):
        g = [mp.mpf(1), mp.mpf(1)/12, mp.mpf(1)/288,
             -mp.mpf(139)/51840, -mp.mpf(571)/2488320,
             mp.mpf(163879)/209018880, mp.mpf(5246819)/75246796800]
        mp_assert_allclose(compute_g(7), g)