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 evalf_trig(v, prec, options):
    """
    This function handles sin and cos of complex arguments.

    TODO: should also handle tan of complex arguments.
    """
    if v.func is C.cos:
        func = mpf_cos
    elif v.func is C.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 C.cos:
            return fone, None, prec, None
        elif v.func is C.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 print_boys_test(maxm):
    l1 = []
    l2 = []

    ms = [S(m) for m in xrange(maxm + 1)]
    ts = [
        S(0),
        S(10) ** (-20),
        S(10) ** (-10),
        S(10) ** (-7),
        S(10) ** (-5),
        S(10) ** (-2),
        S(10) ** (-1),
        5 * S(10) ** (-1),
        S(1),
        Rational(3, 2),
        S(2),
        S(10),
        S(10) ** 2,
        S(10) ** 5,
        S(10) ** 10,
        S(10) ** 20,
    ]

    for m in ms:
        for t in ts:
            r = boys_mpmath(m, t).evalf(30, maxn=1000000, strict=True)
            s = to_str(r._mpf_, 20, strip_zeros=False, min_fixed=0, max_fixed=0)
            # print m, t, s
            l1.append("boys_function(%i,%.1e)" % (m, t))
            l2.append(s)

    print "result = np.array([%s])" % ", ".join(l1)
    print "check = np.array([%s])" % ", ".join(l2)

 def _print_Real(self, expr):
     prec = expr._prec
     if prec < 5:
         dps = 0
     else:
         dps = prec_to_dps(expr._prec)
     if self._settings["full_prec"] == True:
         strip = False
     elif self._settings["full_prec"] == False:
         strip = True
     elif self._settings["full_prec"] == "auto":
         strip = self._print_level > 1
     return mlib.to_str(expr._mpf_, dps, strip_zeros=strip)

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"]))
            re, im = x._eval_evalf(prec).as_real_imag()
            if re.has(re_) or im.has(im_):
                raise NotImplementedError
            if re == 0:
                re = None
                reprec = None
            else:
                re = re._to_mpmath(prec, allow_ints=False)._mpf_
                reprec = prec
            if im == 0:
                im = None
                imprec = None
            else:
                im = im._to_mpmath(prec, allow_ints=False)._mpf_
                imprec = prec
            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_boys_inc(maxm, resolution):
    # Computation of all the boys function data
    boys_fn_data = []
    for m in xrange(maxm + 1):
        l = []
        boys_fn_data.append(l)
        m = sympify(m)
        counter = 0
        while True:
            t = Rational(counter, resolution)
            f = boys_mpmath(m, t)
            f_numer = f.evalf(30, maxn=1000, strict=True)
            f_str = to_str(f_numer._mpf_, 16, strip_zeros=False, min_fixed=0, max_fixed=0)
            l.append(f_str)
            ft = boys_function_tail(m, t)
            ft_numer = ft.evalf(30, maxn=1000, strict=True)
            counter += 1
            error = abs(float(ft) - float(f))
            if counter % 100 == 0 or error == 0:
                print "%2i %5i % 20.12e % 20.12e" % (m, counter, f, error)
            if error == 0:
                break

    f = file("boys_inc.cpp", "w")
    print >> f
    print >> f
    print >> f, "#define BOYS_RESOLUTION", resolution
    print >> f, "#define BOYS_MAX_DATA", maxm
    print >> f
    for m in xrange(maxm + 1):
        print >> f, "static double boys_fn_data_%i[%i] = {" % (m, len(boys_fn_data[m]))
        print >> f, wrap_strings(boys_fn_data[m], 22, 10)
        print >> f, "};"
        print >> f

    print >> f, "static double *boys_fn_data[BOYS_MAX_DATA+1] = {"
    print >> f, wrap_strings(["boys_fn_data_%i" % m for m in xrange(maxm + 1)], 16, 4)
    print >> f, "};"
    print >> f

    print >> f, "static long boys_sizes[BOYS_MAX_DATA+1] = {"
    print >> f, wrap_strings([str(len(boys_fn_data[m])) for m in xrange(maxm + 1)], 6, 8)
    print >> f, "};"
    print >> f

    f.close()

 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"] == True:
         strip = False
     elif self._settings["full_prec"] == 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 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:
        #r = finalize_complex(x._eval_evalf(prec)._mpf_, fzero, prec)
        try:
            # Fall back to ordinary evalf if possible
            if 'subs' in options:
                x = x.subs(evalf_subs(prec, options['subs']))
            re, im = x._eval_evalf(prec).as_real_imag()
            if re.has(re_) or im.has(im_):
                raise NotImplementedError
            if re == 0:
                re = None
                reprec = None
            else:
                re = re._to_mpmath(prec, allow_ints=False)._mpf_
                reprec = prec
            if im == 0:
                im = None
                imprec = None
            else:
                im = im._to_mpmath(prec, allow_ints=False)._mpf_
                imprec = prec
            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
    if options.get("chop"):
        r = chop_parts(r, prec)
    if options.get("strict"):
        check_target(x, r, prec)
    return r

def evalf(x, prec, options):
    try:
        rf = evalf_table[x.func]
        r = rf(x, prec, options)
    except KeyError:
        #r = finalize_complex(x._eval_evalf(prec)._mpf_, fzero, prec)
        try:
            # Fall back to ordinary evalf if possible
            if 'subs' in options:
                x = x.subs(options['subs'])
            r = x._eval_evalf(prec)._mpf_, None, prec, None
        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
    if options.get("chop"):
        r = chop_parts(r, prec)
    if options.get("strict"):
        check_target(x, r, prec)
    return r

    def _print_Float(self, expr):
        # Based off of that in StrPrinter
        dps = prec_to_dps(expr._prec)
        str_real = mlib.to_str(expr._mpf_, dps, strip_zeros=True)

        # Must always have a mul symbol (as 2.5 10^{20} just looks odd)
        separator = r" \times "

        if self._settings['mul_symbol'] is not None:
            separator = self._settings['mul_symbol_latex']

        if 'e' in str_real:
            (mant, exp) = str_real.split('e')

            if exp[0] == '+':
                exp = exp[1:]

            return r"%s%s10^{%s}" % (mant, separator, exp)
        elif str_real == "+inf":
            return r"\infty"
        elif str_real == "-inf":
            return r"- \infty"
        else:
            return str_real

 def _print_Real(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)