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 _d(n, j, prec, sq23pi, sqrt8):
    """
    Compute the sinh term in the outer sum of the HRR formula.
    The constants sqrt(2/3*pi) and sqrt(8) must be precomputed.
    """
    j = from_int(j)
    pi = mpf_pi(prec)
    a = mpf_div(sq23pi, j, prec)
    b = mpf_sub(from_int(n), from_rational(1, 24, prec), prec)
    c = mpf_sqrt(b, prec)
    ch, sh = mpf_cosh_sinh(mpf_mul(a, c), prec)
    D = mpf_div(mpf_sqrt(j, prec), mpf_mul(mpf_mul(sqrt8, b), pi), prec)
    E = mpf_sub(mpf_mul(a, ch), mpf_div(sh, c, prec), prec)
    return mpf_mul(D, E)

def _create_evalf_table():
    global evalf_table
    evalf_table = {
        C.Symbol: evalf_symbol,
        C.Dummy: evalf_symbol,
        C.Float: lambda x, prec, options: (x._mpf_, None, prec, None),
        C.Rational: lambda x, prec, options: (from_rational(x.p, x.q, prec), None, prec, None),
        C.Integer: lambda x, prec, options: (from_int(x.p, prec), None, prec, None),
        C.Zero: lambda x, prec, options: (None, None, prec, None),
        C.One: lambda x, prec, options: (fone, None, prec, None),
        C.Half: lambda x, prec, options: (fhalf, None, prec, None),
        C.Pi: lambda x, prec, options: (mpf_pi(prec), None, prec, None),
        C.Exp1: lambda x, prec, options: (mpf_e(prec), None, prec, None),
        C.ImaginaryUnit: lambda x, prec, options: (None, fone, None, prec),
        C.NegativeOne: lambda x, prec, options: (fnone, None, prec, None),
        C.NaN: lambda x, prec, options: (fnan, None, prec, None),
        C.exp: lambda x, prec, options: evalf_pow(C.Pow(S.Exp1, x.args[0], evaluate=False), prec, options),
        C.cos: evalf_trig,
        C.sin: evalf_trig,
        C.Add: evalf_add,
        C.Mul: evalf_mul,
        C.Pow: evalf_pow,
        C.log: evalf_log,
        C.atan: evalf_atan,
        C.Abs: evalf_abs,
        C.re: evalf_re,
        C.im: evalf_im,
        C.floor: evalf_floor,
        C.ceiling: evalf_ceiling,
        C.Integral: evalf_integral,
        C.Sum: evalf_sum,
        C.Piecewise: evalf_piecewise,
        C.bernoulli: evalf_bernoulli,
    }

def npartitions(n, verbose=False):
    """
    Calculate the partition function P(n), i.e. the number of ways that
    n can be written as a sum of positive integers.

    P(n) is computed using the Hardy-Ramanujan-Rademacher formula,
    described e.g. at http://mathworld.wolfram.com/PartitionFunctionP.html

    The correctness of this implementation has been tested for 10**n
    up to n = 8.
    """
    n = int(n)
    if n < 0: return 0
    if n <= 5: return [1, 1, 2, 3, 5, 7][n]
    # Estimate number of bits in p(n). This formula could be tidied
    pbits = int((math.pi*(2*n/3.)**0.5-math.log(4*n))/math.log(10)+1)*\
        math.log(10,2)
    prec = p = int(pbits*1.1 + 100)
    s = fzero
    M = max(6, int(0.24*n**0.5+4))
    sq23pi = mpf_mul(mpf_sqrt(from_rational(2,3,p), p), mpf_pi(p), p)
    sqrt8 = mpf_sqrt(from_int(8), p)
    for q in xrange(1, M):
        a = A(n,q,p)
        d = D(n,q,p, sq23pi, sqrt8)
        s = mpf_add(s, mpf_mul(a, d), prec)
        if verbose:
            print "step", q, "of", M, to_str(a, 10), to_str(d, 10)
        # On average, the terms decrease rapidly in magnitude. Dynamically
        # reducing the precision greatly improves performance.
        p = bitcount(abs(to_int(d))) + 50
    np = to_int(mpf_add(s, fhalf, prec))
    return int(np)

 def calc_part(expr, nexpr):
     nint = int(to_int(nexpr, round_nearest))
     expr = C.Add(expr, -nint, evaluate=False)
     x, _, x_acc, _ = evalf(expr, 10, options)
     check_target(expr, (x, None, x_acc, None), 3)
     nint += int(no*(mpf_cmp(x or fzero, fzero) == no))
     nint = from_int(nint)
     return nint, fastlog(nint) + 10

 def _convert_tol(ctx, tol):
     if isinstance(tol, int_types):
         return from_int(tol)
     if isinstance(tol, float):
         return from_float(tol)
     if hasattr(tol, "_mpf_"):
         return tol._mpf_
     prec, rounding = ctx._prec_rounding
     if isinstance(tol, basestring):
         return from_str(tol, prec, rounding)
     raise ValueError("expected a real number, got %s" % tol)

 def calc_part(expr, nexpr):
     nint = int(to_int(nexpr, rnd))
     expr = C.Add(expr, -nint, evaluate=False)
     x, _, x_acc, _ = evalf(expr, 10, options)
     try:
         check_target(expr, (x, None, x_acc, None), 3)
     except PrecisionExhausted:
         if not expr.equals(0):
             raise PrecisionExhausted
         x = fzero
     nint += int(no * (mpf_cmp(x or fzero, fzero) == no))
     nint = from_int(nint)
     return nint, fastlog(nint) + 10

 def _as_mpf_val(self, prec):
     return mlib.from_int(self.p)