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 [1]_.


    The correctness of this implementation has been tested through 10**10.

    Examples
    ========

    >>> from sympy.ntheory import npartitions
    >>> npartitions(25)
    1958

    References
    ==========

    .. [1] http://mathworld.wolfram.com/PartitionFunctionP.html

    """
    n = int(n)
    if n < 0:
        return 0
    if n <= 5:
        return [1, 1, 2, 3, 5, 7][n]
    if '_factor' not in globals():
        _pre()
    # 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))
    if M > 10**5:
        raise ValueError("Input too big") # Corresponds to n > 1.7e11
    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 range(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
    return int(to_int(mpf_add(s, fhalf, prec)))

    def _print_Float(self, expr):
        """Print a Float in C-like scientific notation."""
        prec = expr._prec

        if prec < 5:
            dps = 0
        else:
            dps = prec_to_dps(expr._prec)
        if self._settings["full_prec"] is True:
            strip = False
        elif self._settings["full_prec"] is False:
            strip = True
        elif self._settings["full_prec"] == "auto":
            strip = self._print_level > 1

        rv = to_str(expr._mpf_, dps, strip_zeros=strip, max_fixed=-2, min_fixed=2)

        if rv.startswith('-.0'):
            rv = '-0.' + rv[3:]
        elif rv.startswith('.0'):
            rv = '0.' + rv[2:]

        if self.dtype == np.float32:
            rv = rv + 'F'

        return rv

def evalf_trig(v, prec, options):
    """
    This function handles sin and cos of complex arguments.

    TODO: should also handle tan of complex arguments.
    """
    from sympy import cos, sin
    if v.func is cos:
        func = mpf_cos
    elif v.func is sin:
        func = mpf_sin
    else:
        raise NotImplementedError
    arg = v.args[0]
    # 20 extra bits is possibly overkill. It does make the need
    # to restart very unlikely
    xprec = prec + 20
    re, im, re_acc, im_acc = evalf(arg, xprec, options)
    if im:
        if 'subs' in options:
            v = v.subs(options['subs'])
        return evalf(v._eval_evalf(prec), prec, options)
    if not re:
        if v.func is cos:
            return fone, None, prec, None
        elif v.func is sin:
            return None, None, None, None
        else:
            raise NotImplementedError
    # For trigonometric functions, we are interested in the
    # fixed-point (absolute) accuracy of the argument.
    xsize = fastlog(re)
    # Magnitude <= 1.0. OK to compute directly, because there is no
    # danger of hitting the first root of cos (with sin, magnitude
    # <= 2.0 would actually be ok)
    if xsize < 1:
        return func(re, prec, rnd), None, prec, None
    # Very large
    if xsize >= 10:
        xprec = prec + xsize
        re, im, re_acc, im_acc = evalf(arg, xprec, options)
    # Need to repeat in case the argument is very close to a
    # multiple of pi (or pi/2), hitting close to a root
    while 1:
        y = func(re, prec, rnd)
        ysize = fastlog(y)
        gap = -ysize
        accuracy = (xprec - xsize) - gap
        if accuracy < prec:
            if options.get('verbose'):
                print("SIN/COS", accuracy, "wanted", prec, "gap", gap)
                print(to_str(y, 10))
            if xprec > options.get('maxprec', DEFAULT_MAXPREC):
                return y, None, accuracy, None
            xprec += gap
            re, im, re_acc, im_acc = evalf(arg, xprec, options)
            continue
        else:
            return y, None, prec, None

def npartitions(n):
    """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.

    The correctness of this implementation has been tested for 10**n
    up to n = 8.

    Examples
    ========

    >>> npartitions(25)
    1958

    References
    ==========

    * http://mathworld.wolfram.com/PartitionFunctionP.html

    """
    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 range(1, M):
        a = _a(n, q, p)
        d = _d(n, q, p, sq23pi, sqrt8)
        s = mpf_add(s, mpf_mul(a, d), prec)
        debug("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
    return int(to_int(mpf_add(s, fhalf, prec)))

def evalf(x, prec, options):
    from sympy import re as re_, im as im_
    try:
        rf = evalf_table[x.func]
        r = rf(x, prec, options)
    except KeyError:
        try:
            # Fall back to ordinary evalf if possible
            if 'subs' in options:
                x = x.subs(evalf_subs(prec, options['subs']))
            xe = x._eval_evalf(prec)
            re, im = xe.as_real_imag()
            if re.has(re_) or im.has(im_):
                raise NotImplementedError
            if re == 0:
                re = None
                reprec = None
            elif re.is_number:
                re = re._to_mpmath(prec, allow_ints=False)._mpf_
                reprec = prec
            else:
                raise NotImplementedError
            if im == 0:
                im = None
                imprec = None
            elif im.is_number:
                im = im._to_mpmath(prec, allow_ints=False)._mpf_
                imprec = prec
            else:
                raise NotImplementedError
            r = re, im, reprec, imprec
        except AttributeError:
            raise NotImplementedError
    if options.get("verbose"):
        print("### input", x)
        print("### output", to_str(r[0] or fzero, 50))
        print("### raw", r) # r[0], r[2]
        print()
    chop = options.get('chop', False)
    if chop:
        if chop is True:
            chop_prec = prec
        else:
            # convert (approximately) from given tolerance;
            # the formula here will will make 1e-i rounds to 0 for
            # i in the range +/-27 while 2e-i will not be chopped
            chop_prec = int(round(-3.321*math.log10(chop) + 2.5))
            if chop_prec == 3:
                chop_prec -= 1
        r = chop_parts(r, chop_prec)
    if options.get("strict"):
        check_target(x, r, prec)
    return r

 def _print_Float(self, expr):
     prec = expr._prec
     if prec < 5:
         dps = 0
     else:
         dps = prec_to_dps(expr._prec)
     if self._settings["full_prec"] is True:
         strip = False
     elif self._settings["full_prec"] is False:
         strip = True
     elif self._settings["full_prec"] == "auto":
         strip = self._print_level > 1
     rv = mlib.to_str(expr._mpf_, dps, strip_zeros=strip)
     if rv.startswith('-.0'):
         rv = '-0.' + rv[3:]
     elif rv.startswith('.0'):
         rv = '0.' + rv[2:]
     return rv

 def _print_Float(self, expr):
     """
     override method in StrPrinter
     always printing floating point numbers in scientific notation
     """
     prec = expr._prec
     if prec < 5:
         dps = 0
     else:
         dps = prec_to_dps(expr._prec)
     if self._settings["full_prec"] is True:
         strip = False
     elif self._settings["full_prec"] is False:
         strip = True
     elif self._settings["full_prec"] == "auto":
         strip = self._print_level > 1
     rv = to_str(expr._mpf_, dps, strip_zeros=strip, max_fixed=-2, min_fixed=2)
     if rv.startswith('-.0'):
         rv = '-0.' + rv[3:]
     elif rv.startswith('.0'):
         rv = '0.' + rv[2:]
     return rv

 def _print_Float(self, expr):
     dps = prec_to_dps(expr._prec)
     r = mlib.to_str(expr._mpf_, repr_dps(expr._prec))
     return "%s('%s', prec=%i)" % (expr.__class__.__name__, r, dps)

 def _print_Float(self, expr):
     r = mlib.to_str(expr._mpf_, repr_dps(expr._prec))
     return "%s('%s', precision=%i)" % (expr.__class__.__name__, r, expr._prec)