Python sympy.hypersimp函数代码示例

def hypsum(expr, n, start, prec):
    """
    Sum a rapidly convergent infinite hypergeometric series with
    given general term, e.g. e = hypsum(1/factorial(n), n). The
    quotient between successive terms must be a quotient of integer
    polynomials.
    """
    from sympy import hypersimp, lambdify

    if start:
        expr = expr.subs(n, n + start)
    hs = hypersimp(expr, n)
    if hs is None:
        raise NotImplementedError("a hypergeometric series is required")
    num, den = hs.as_numer_denom()

    func1 = lambdify(n, num)
    func2 = lambdify(n, den)

    h, g, p = check_convergence(num, den, n)

    if h < 0:
        raise ValueError("Sum diverges like (n!)^%i" % (-h))

    # Direct summation if geometric or faster
    if h > 0 or (h == 0 and abs(g) > 1):
        term = expr.subs(n, 0)
        term = (MPZ(term.p) << prec) // term.q
        s = term
        k = 1
        while abs(term) > 5:
            term *= MPZ(func1(k - 1))
            term //= MPZ(func2(k - 1))
            s += term
            k += 1
        return from_man_exp(s, -prec)
    else:
        alt = g < 0
        if abs(g) < 1:
            raise ValueError("Sum diverges like (%i)^n" % abs(1 / g))
        if p < 1 or (p == 1 and not alt):
            raise ValueError("Sum diverges like n^%i" % (-p))
        # We have polynomial convergence: use Richardson extrapolation
        # Need to use at least quad precision because a lot of cancellation
        # might occur in the extrapolation process
        prec2 = 4 * prec
        term = expr.subs(n, 0)
        term = (MPZ(term.p) << prec2) // term.q

        def summand(k, _term=[term]):
            if k:
                k = int(k)
                _term[0] *= MPZ(func1(k - 1))
                _term[0] //= MPZ(func2(k - 1))
            return make_mpf(from_man_exp(_term[0], -prec2))

        with workprec(prec):
            v = nsum(summand, [0, mpmath_inf], method="richardson")

        return v._mpf_

def test_hypersimp():
    n, k = symbols('nk', integer=True)

    assert hypersimp(factorial(k), k) == k + 1
    assert hypersimp(factorial(k**2), k) is None

    assert hypersimp(1/factorial(k), k) == 1/(k + 1)

    assert hypersimp(2**k/factorial(k)**2, k) == 2/(k**2+2*k+1)

    assert hypersimp(binomial(n, k), k) == (n-k)/(k+1)
    assert hypersimp(binomial(n+1, k), k) == (n-k+1)/(k+1)

    term = (4*k+1)*factorial(k)/factorial(2*k+1)
    assert hypersimp(term, k) == (4*k + 5)/(6 + 16*k**2 + 28*k)

    term = 1/((2*k-1)*factorial(2*k+1))
    assert hypersimp(term, k) == (2*k-1)/(6 + 22*k + 24*k**2 + 8*k**3)

    term = binomial(n, k)*(-1)**k/factorial(k)
    assert hypersimp(term, k) == (k - n)/(k**2+2*k+1)

def hypsum(expr, n, start, prec):
    """
    Sum a rapidly convergent infinite hypergeometric series with
    given general term, e.g. e = hypsum(1/factorial(n), n). The
    quotient between successive terms must be a quotient of integer
    polynomials.
    """
    from sympy import Float, hypersimp, lambdify

    if prec == float('inf'):
        raise NotImplementedError('does not support inf prec')

    if start:
        expr = expr.subs(n, n + start)
    hs = hypersimp(expr, n)
    if hs is None:
        raise NotImplementedError("a hypergeometric series is required")
    num, den = hs.as_numer_denom()

    func1 = lambdify(n, num)
    func2 = lambdify(n, den)

    h, g, p = check_convergence(num, den, n)

    if h < 0:
        raise ValueError("Sum diverges like (n!)^%i" % (-h))

    term = expr.subs(n, 0)
    if not term.is_Rational:
        raise NotImplementedError("Non rational term functionality is not implemented.")

    # Direct summation if geometric or faster
    if h > 0 or (h == 0 and abs(g) > 1):
        term = (MPZ(term.p) << prec) // term.q
        s = term
        k = 1
        while abs(term) > 5:
            term *= MPZ(func1(k - 1))
            term //= MPZ(func2(k - 1))
            s += term
            k += 1
        return from_man_exp(s, -prec)
    else:
        alt = g < 0
        if abs(g) < 1:
            raise ValueError("Sum diverges like (%i)^n" % abs(1/g))
        if p < 1 or (p == 1 and not alt):
            raise ValueError("Sum diverges like n^%i" % (-p))
        # We have polynomial convergence: use Richardson extrapolation
        vold = None
        ndig = prec_to_dps(prec)
        while True:
            # Need to use at least quad precision because a lot of cancellation
            # might occur in the extrapolation process; we check the answer to
            # make sure that the desired precision has been reached, too.
            prec2 = 4*prec
            term0 = (MPZ(term.p) << prec2) // term.q

            def summand(k, _term=[term0]):
                if k:
                    k = int(k)
                    _term[0] *= MPZ(func1(k - 1))
                    _term[0] //= MPZ(func2(k - 1))
                return make_mpf(from_man_exp(_term[0], -prec2))

            with workprec(prec):
                v = nsum(summand, [0, mpmath_inf], method='richardson')
            vf = Float(v, ndig)
            if vold is not None and vold == vf:
                break
            prec += prec  # double precision each time
            vold = vf

        return v._mpf_

def hypsum(expr, n, start, prec):
    """
    Sum a rapidly convergent infinite hypergeometric series with
    given general term, e.g. e = hypsum(1/factorial(n), n). The
    quotient between successive terms must be a quotient of integer
    polynomials.
    """
    from sympy import hypersimp, lambdify

    if start:
        expr = expr.subs(n, n+start)
    hs = hypersimp(expr, n)
    if hs is None:
        raise NotImplementedError("a hypergeometric series is required")
    num, den = hs.as_numer_denom()
    func1 = lambdify(n, num)
    func2 = lambdify(n, den)

    h, g, p = check_convergence(num, den, n)

    if h < 0:
        raise ValueError("Sum diverges like (n!)^%i" % (-h))

    # Direct summation if geometric or faster
    if h > 0 or (h == 0 and abs(g) > 1):
        one = MP_BASE(1) << prec
        term = expr.subs(n, 0)
        term = (MP_BASE(term.p) << prec) // term.q
        s = term
        k = 1
        while abs(term) > 5:
            term *= MP_BASE(func1(k-1))
            term //= MP_BASE(func2(k-1))
            s += term
            k += 1
        return from_man_exp(s, -prec)
    else:
        alt = g < 0
        if abs(g) < 1:
            raise ValueError("Sum diverges like (%i)^n" % abs(1/g))
        if p < 1 or (p == 1 and not alt):
            raise ValueError("Sum diverges like n^%i" % (-p))
        # We have polynomial convergence:
        # Use Shanks extrapolation for alternating series,
        # Richardson extrapolation for nonalternating series
        if alt:
            # XXX: better parameters for Shanks transformation
            # This tends to get bad somewhere > 50 digits
            N = 5 + int(prec*0.36)
            M = 2 + N//3
            NTERMS = M + N + 2
        else:
            N = 3 + int(prec*0.15)
            M = 2*N
            NTERMS = M + N + 2
        # Need to use at least double precision because a lot of cancellation
        # might occur in the extrapolation process
        prec2 = 2*prec
        one = MP_BASE(1) << prec2
        term = expr.subs(n, 0)
        term = (MP_BASE(term.p) << prec2) // term.q
        s = term
        table = [make_mpf(from_man_exp(s, -prec2))]
        for k in xrange(1, NTERMS):
            term *= MP_BASE(func1(k-1))
            term //= MP_BASE(func2(k-1))
            s += term
            table.append(make_mpf(from_man_exp(s, -prec2)))
            k += 1
        orig = mp.prec
        try:
            mp.prec = prec
            if alt:
                v = shanks_extrapolation(table, N, M)
            else:
                v = richardson_extrapolation(table, N, M)
        finally:
            mp.prec = orig
        return v._mpf_