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

#.........这里部分代码省略.........
    special = []
    for arg in args:
        arg = evalf(arg, prec, options)
        if arg[0] is None:
            continue
        arg = C.Float._new(arg[0], 1)
        if arg is S.NaN or arg.is_unbounded:
            special.append(arg)
    if special:
        from sympy.core.mul import Mul

        special = Mul(*special)
        return evalf(special, prec + 4, {})

    # With guard digits, multiplication in the real case does not destroy
    # accuracy. This is also true in the complex case when considering the
    # total accuracy; however accuracy for the real or imaginary parts
    # separately may be lower.
    acc = prec

    # XXX: big overestimate
    working_prec = prec + len(args) + 5

    # Empty product is 1
    start = man, exp, bc = MPZ(1), 0, 1

    # First, we multiply all pure real or pure imaginary numbers.
    # direction tells us that the result should be multiplied by
    # I**direction; all other numbers get put into complex_factors
    # to be multiplied out after the first phase.
    last = len(args)
    direction = 0
    args.append(S.One)
    complex_factors = []

    for i, arg in enumerate(args):
        if i != last and pure_complex(arg):
            args[-1] = (args[-1] * arg).expand()
            continue
        elif i == last and arg is S.One:
            continue
        re, im, re_acc, im_acc = evalf(arg, working_prec, options)
        if re and im:
            complex_factors.append((re, im, re_acc, im_acc))
            continue
        elif re:
            (s, m, e, b), w_acc = re, re_acc
        elif im:
            (s, m, e, b), w_acc = im, im_acc
            direction += 1
        else:
            return None, None, None, None
        direction += 2 * s
        man *= m
        exp += e
        bc += b
        if bc > 3 * working_prec:
            man >>= working_prec
            exp += working_prec
        acc = min(acc, w_acc)
    sign = (direction & 2) >> 1
    if not complex_factors:
        v = normalize(sign, man, exp, bitcount(man), prec, rnd)
        # multiply by i
        if direction & 1:
            return None, v, None, acc
        else:
            return v, None, acc, None
    else:
        # initialize with the first term
        if (man, exp, bc) != start:
            # there was a real part; give it an imaginary part
            re, im = (sign, man, exp, bitcount(man)), (0, MPZ(0), 0, 0)
            i0 = 0
        else:
            # there is no real part to start (other than the starting 1)
            wre, wim, wre_acc, wim_acc = complex_factors[0]
            acc = min(acc, complex_accuracy((wre, wim, wre_acc, wim_acc)))
            re = wre
            im = wim
            i0 = 1

        for wre, wim, wre_acc, wim_acc in complex_factors[i0:]:
            # acc is the overall accuracy of the product; we aren't
            # computing exact accuracies of the product.
            acc = min(acc, complex_accuracy((wre, wim, wre_acc, wim_acc)))

            use_prec = working_prec
            A = mpf_mul(re, wre, use_prec)
            B = mpf_mul(mpf_neg(im), wim, use_prec)
            C = mpf_mul(re, wim, use_prec)
            D = mpf_mul(im, wre, use_prec)
            re = mpf_add(A, B, use_prec)
            im = mpf_add(C, D, use_prec)
        if options.get("verbose"):
            print "MUL: wanted", prec, "accurate bits, got", acc
        # multiply by I
        if direction & 1:
            re, im = mpf_neg(im), re
        return re, im, acc, acc

def evalf_mul(v, prec, options):
    args = v.args
    # With guard digits, multiplication in the real case does not destroy
    # accuracy. This is also true in the complex case when considering the
    # total accuracy; however accuracy for the real or imaginary parts
    # separately may be lower.
    acc = prec
    target_prec = prec
    # XXX: big overestimate
    prec = prec + len(args) + 5
    direction = 0
    # Empty product is 1
    man, exp, bc = MPZ(1), 0, 1
    direction = 0
    complex_factors = []
    # First, we multiply all pure real or pure imaginary numbers.
    # direction tells us that the result should be multiplied by
    # i**direction
    for arg in args:
        re, im, re_acc, im_acc = evalf(arg, prec, options)
        if re and im:
            complex_factors.append((re, im, re_acc, im_acc))
            continue
        elif re:
            (s, m, e, b), w_acc = re, re_acc
        elif im:
            (s, m, e, b), w_acc = im, im_acc
            direction += 1
        else:
            return None, None, None, None
        direction += 2*s
        man *= m
        exp += e
        bc += b
        if bc > 3*prec:
            man >>= prec
            exp += prec
        acc = min(acc, w_acc)
    sign = (direction & 2) >> 1
    v = normalize(sign, man, exp, bitcount(man), prec, round_nearest)
    if complex_factors:
        # make existing real scalar look like an imaginary and
        # multiply by the remaining complex numbers
        re, im = v, (0, MPZ(0), 0, 0)
        for wre, wim, wre_acc, wim_acc in complex_factors:
            # acc is the overall accuracy of the product; we aren't
            # computing exact accuracies of the product.
            acc = min(acc,
                      complex_accuracy((wre, wim, wre_acc, wim_acc)))
            A = mpf_mul(re, wre, prec)
            B = mpf_mul(mpf_neg(im), wim, prec)
            C = mpf_mul(re, wim, prec)
            D = mpf_mul(im, wre, prec)
            re, xre_acc = add_terms([(A, acc), (B, acc)], prec, target_prec)
            im, xim_acc = add_terms([(C, acc), (D, acc)], prec, target_prec)

        if options.get('verbose'):
            print "MUL: wanted", target_prec, "accurate bits, got", acc
        # multiply by i
        if direction & 1:
            return mpf_neg(im), re, acc, acc
        else:
            return re, im, acc, acc
    else:
        # multiply by i
        if direction & 1:
            return None, v, None, acc
        else:
            return v, None, acc, None

 def __mul__(self, other):
     if isinstance(other, Number):
         rhs, prec = other._as_mpf_op(self._prec)
         return Real._new(mlib.mpf_mul(self._mpf_, rhs, prec, rnd), prec)
     return Number.__mul__(self, other)

 def _make_tol(ctx):
     hundred = (0, 25, 2, 5)
     eps = (0, MPZ_ONE, 1-ctx.prec, 1)
     return mpf_mul(hundred, eps)