def _to_mpmath(self, prec, allow_ints=True):
     # mpmath functions accept ints as input
     errmsg = "cannot convert to mpmath number"
     if allow_ints and self.is_Integer:
         return self.p
     try:
         re, im, _, _ = evalf(self, prec, {})
         if im:
             return make_mpc((re, im))
         else:
             return make_mpf(re)
     except NotImplementedError:
         v = self._eval_evalf(prec)
         if v is None:
             raise ValueError(errmsg)
         if v.is_Float:
             return make_mpf(v._mpf_)
         # Number + Number*I is also fine
         re, im = v.as_real_imag()
         if allow_ints and re.is_Integer:
             re = from_int(re.p)
         elif re.is_Float:
             re = re._mpf_
         else:
             raise ValueError(errmsg)
         if allow_ints and im.is_Integer:
             im = from_int(im.p)
         elif im.is_Float:
             im = im._mpf_
         else:
             raise ValueError(errmsg)
         return make_mpc((re, im))

 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))

 def _mpmath_(self, prec, rnd):
     return mpmath.make_mpf(mlib.from_rational(self.p, self.q, prec, rnd))

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_