def mpf_assert_allclose(res, std, atol=0, rtol=1e-17):
    try:
        len(res)
    except TypeError:
        res = list(res)

    n = len(std)
    if len(res) != n:
        raise AssertionError("Lengths of inputs not equal.")

    failures = []
    for k in range(n):
        try:
            assert_(mp.fabs(res[k] - std[k]) <= atol + rtol*mp.fabs(std[k]))
        except AssertionError:
            failures.append(k)

    ndigits = int(abs(np.log10(rtol)))
    msg = [""]
    msg.append("Bad results ({} out of {}) for the following points:"
               .format(len(failures), n))
    for k in failures:
        resrep = mp.nstr(res[k], ndigits, min_fixed=0, max_fixed=0)
        stdrep = mp.nstr(std[k], ndigits, min_fixed=0, max_fixed=0)
        if std[k] == 0:
            rdiff = "inf"
        else:
            rdiff = mp.fabs((res[k] - std[k])/std[k])
            rdiff = mp.nstr(rdiff, 3)
        msg.append("{}: {} != {} (rdiff {})".format(k, resrep, stdrep, rdiff))
    if failures:
        assert_(False, "\n".join(msg))

 def add_gwb(self, gwb, dist, InjectionFile="None", verbose=-1):
     """ Add GW background on simulated TOAs using a GWB object 
     from libstempo and the pulsar distance in kpc."""
     if verbose==-1 : verbose = self.verbose
     
     if verbose!=0 : print "Add GWB ..."
     #### Making libstempo.tempopulsar save in parfile and timfile
     localverbose=0
     if verbose>1 : localverbose=1
     self.savepar("TmpIdeal", verbose=localverbose)
     self.savetim("TmpIdeal",IdealTOAs=True, verbose=localverbose)
     psr = libstempo.tempopulsar(parfile="TmpIdeal.par",timfile="TmpIdeal.tim",dofit=False)
     #### Creating GWB data
     GWBval = gwb.add_gwb(psr,dist)
     #### Adding data
     fOut = None
     if InjectionFile!="None":
         fOut = open(InjectionFile,'w')
         fOut.write("#TOAIdeal GWB TOARelBeforeInjection TOARelAfterInjection DiffTOAAft-Bef\n")
     for itoa in xrange(len(self.timTOAs)) :
         TOABefInj = self.timTOAs[itoa][0]
         self.timTOAs[itoa][0] += mp.mpf(np.float64(GWBval[itoa]),n=self.prec)
         if fOut!=None :
             fOut.write(mp.nstr(self.timTOAsIdeal[itoa][0],n=self.prec)+" "+repr(GWBval[itoa])+" "\
                        +mp.nstr(TOABefInj-self.timTOAsIdeal[0][0],n=self.prec)+" "\
                        +mp.nstr(self.timTOAs[itoa][0]-self.timTOAsIdeal[0][0],n=self.prec)+" "\
                        +mp.nstr(self.timTOAs[itoa][0]-TOABefInj,n=self.prec)+"\n")
     if fOut!=None :
         fOut.close()

 def savetim(self, basename, sysname='all', multifile=False, IdealTOAs=False, verbose=-1):
     "Save TOA in tim file. It is strictly the original tim, except the column 3 and 4 corresponding to the new toa and error."
     if verbose==-1 : verbose = self.verbose
     if verbose!=0 : print "Save TOAs in",basename+".tim ..."
     
     fOut = open(basename+".tim",'w')
     for line in self.timHead :
         for iw in xrange(len(line)):
             if iw!=0:
                 fOut.write(" ")
             fOut.write(line[iw])
         fOut.write("\n")
     if multifile :
         os.system("mkdir -p "+basename)
         for xsys in self.syslist :
             fOut.write("INCLUDE "+os.path.basename(basename)+"/"+xsys+".tim\n")
             self.savetim(basename+"/"+xsys, sysname=xsys, multifile=False, IdealTOAs=IdealTOAs)
     else :
         tTOAs = self.timTOAs
         if IdealTOAs :
             tTOAs = self.timTOAsIdeal
         for tl in tTOAs :
             if sysname=='all' or sysname==tl[2] :
                 tl[3][2] = "%s" % (mp.nstr(tl[0],n=self.prec)) # update the word corresponding to the TOAs
                 tl[3][3] = "%s" % (mp.nstr(tl[1],n=self.prec)) # update the word corresponding to the errors 
                 ## Write tim line as a series of words
                 for iw in xrange(len(tl[3])):
                     if iw!=0 :
                         fOut.write(" ")
                     fOut.write(tl[3][iw])
                 fOut.write("\n")
     fOut.close()

 def output(self, C):
     """Given the node <-> coordinates mapping, output the mapping
     on the screen. The coordinates are complex numbers."""
     for n in sorted(C.keys()):
         cr = mpmath.nstr(C[n]['c'].real, 50)
         ci = mpmath.nstr(C[n]['c'].imag, 50)
         print "%s %s %s" % (n, cr, ci)
     pass

def test_nstr():
    m = matrix([[0.75, 0.190940654, -0.0299195971],
                [0.190940654, 0.65625, 0.205663228],
                [-0.0299195971, 0.205663228, 0.64453125e-20]])
    assert nstr(m, 4, min_fixed=-inf) == \
    '''[    0.75  0.1909                    -0.02992]
[  0.1909  0.6563                      0.2057]
[-0.02992  0.2057  0.000000000000000000006445]'''
    assert nstr(m, 4) == \
    '''[    0.75  0.1909   -0.02992]

def pariDilog(z):
    assert isinstance(z, mpmath.mpc)

    pariStr = 'myDilog = dilog((%s) + (%s) * I)' % (
        mpmath.nstr(z.real, mpmath.mp.dps),
        mpmath.nstr(z.imag, mpmath.mp.dps))

    pari_eval(pariStr)
    
    return mpmath.mpc(
        pari_eval('real(myDilog)').replace(' E','E'),
        pari_eval('imag(myDilog)').replace(' E','E'))

def main():
    print(__doc__)
    print()
    stirling_coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0) for x in stirling_series(16)]
    taylor_coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0) for x in taylor_series_at_1(16)]
    print("Stirling series coefficients")
    print("----------------------------")
    print("\n".join(stirling_coeffs))
    print()
    print("Taylor series coefficients")
    print("--------------------------")
    print("\n".join(taylor_coeffs))
    print()

def printRealNumberAsFixed(r):

    assert isinstance(r, mpmath.mpf)
    
    if lessThanMaxErr(r):
        return '0'
    else:
        maxErrDigits = globalsettings.getSetting("maximalErrorDigits")
        s = mpmath.nstr(r,
                        maxErrDigits,
                        min_fixed = -mpmath.inf, max_fixed = mpmath.inf)
        a, b = s.split('.')

        # mpmath chops 1.2000000 to 1.2, so we are adding '0' to b

        alen = len(a)
        if a[0] in '+-':
            alen -= 1
        
        b += (maxErrDigits - alen + len(b)) * '0'

        # b should never be more than maximalErrorDigits

        b = b[:maxErrDigits - alen]

        return a + '.' + b

def generateRealArgument( range = [ 0, 10 ], allowNegative = True ):
    factor = 1

    if allowNegative and getRandomInteger( 2 ) == 1:
        factor = -1

    return nstr( fmul( exp( fmul( getRandomNumber( ), random.uniform( *range ) ) ), factor ) )

def NumberStr(n):
  # Replace spaces
  s = n.replace(' ', '')

  # remove "exactly" for the carbon mass
  s = s.replace('(exactly)', '')

  # if only a number, put it three times
  m = bnum.match(s)
  if m: 
    s = "{:<25} {:<25} {:<25}".format(m.group(1), m.group(1), m.group(1))  

  # if parentheses uncertainty...
  m = buncertain.match(s)
  if m:
    # tricky. duplicate the first part as a string
    s2 = m.group(1)

    # but replace with all zero
    s2 = re.sub(r'\d', '0', s2)

    # now replace last characters
    l = len(m.group(2))
    s2 = s2[:len(s2)-l] + m.group(2)

    # convert to a float
    serr = mp.mpf(s2)
    scenter = mp.mpf(m.group(1))

    s = "{:<25} {:<25} {:<25}".format(mp.nstr(scenter, 18), mp.nstr(scenter-serr, 18), mp.nstr(scenter+serr, 18))

  # Replace bracketed ranges with parentheses
  m = brange.match(s)
  if m:
    slow = mp.mpf(m.group(1))
    shigh = mp.mpf(m.group(2))
    smid = (shigh + slow)/mp.mpf("2.0")
    s = "{:<25} {:<25} {:<25}".format(mp.nstr(smid, 18), mp.nstr(slow, 18), mp.nstr(shigh, 18))

  # just a dash?
  if s == "-":
    s = "{:<25} {:<25} {:<25}".format(0, 0, 0)


  return s

    def cstr(obj):
        try:
            return '{' + ', '.join(cstr(val) for val in obj) + '}'
        except TypeError:
            cls = type(obj)
            if ((cls == complex) or (cls == mpc)):
                return 'dynd::dynd_complex<double>({}, {})'.format(cstr(obj.real), cstr(obj.imag))

            return nstr(obj, pdps)

def bradify():
    """Iteratively calculates the ratios of Lucas sequence numbers.

    Spoiler: It approaches the golden ratio.
    """
    argument_parser = argparse.ArgumentParser(
        description='Watch the ratios of numbers in a Lucas series approach '
                    'the golden ratio (φ).')
    argument_parser.add_argument('-i', '--iterations',
                                 help='Maximum number of iterations.',
                                 type=int, default=1000)
    argument_parser.add_argument('-p', '--precision',
                                 help='Digits of precision to calculate',
                                 type=int, default=100)
    sequences = argument_parser.add_mutually_exclusive_group(required=True)
    sequences.add_argument('-f', '--fibonacci', dest='special_case',
                           action='store_const', const='f',
                           help='The Fibonacci numbers')
    sequences.add_argument('-l', '--lucas', dest='special_case',
                           action='store_const', const='l',
                           help='The Lucas numbers')
    sequences.add_argument('-b', '--brady', dest='special_case',
                           action='store_const', const='b',
                           help='The Brady numbers')
    sequences.add_argument('-a',
                           type=str, dest='sequence_definition',
                           help='A generic sequence')
    args = argument_parser.parse_args()

    if args.special_case:
        if args.special_case == 'f':
            sequence = lucas.fibo()
        elif args.special_case == 'l':
            sequence = lucas.lucas_num()
        elif args.special_case == 'b':
            sequence = lucas.brady()
    else:
        definition = args.sequence_definition.split(',')
        definition[0] = int(definition[0])
        definition[1] = int(definition[1])
        sequence = lucas.Lucas(definition, 5)

    mpmath.mp.dps = args.precision * 2
    previous = mpmath.mpf(sequence[0])
    oldphi=1000
    for i in range(1, args.iterations):
        current = mpmath.mpf(sequence[i])
        try:
            phi = current/previous
        except ZeroDivisionError:
            phi = float('NaN')
        print('φ[%d] = %s' % (i, mpmath.nstr(phi, args.precision)))
        if abs(oldphi - phi) < 10**(-args.precision) and i > 2:
            break
        previous = current
        oldphi = phi

def first_digits_are_pandigital(val):
    val = mpmath.nstr(val,n=10)
    val = val.split('e')
    val = val[0]
    val = val.split('.')
    val = val[0]+val[1]
    val = val[:9]
    if len(set(val)) == 9 and '0' not in val:
        return 1
    return 0

def mpf2float(x):
    """
    Convert an mpf to the nearest floating point number. Just using
    float directly doesn't work because of results like this:

    with mp.workdps(50):
        float(mpf("0.99999999999999999")) = 0.9999999999999999

    """
    return float(mpmath.nstr(x, 17, min_fixed=0, max_fixed=0))

def func(eq, x):
    """ Calculate the answer to f(x)

    Example:
    >>> [func([1, 0, -3, -4], i) for i in np.arange(0.0, 5.0, 1)]
    ['-4.0', '-6.0', '-2.0', '14.0', '48.0']
    """
    ans, x = mp.mpf(0), mp.mpf(str(x))
    index_length = len(eq)
    for i in range(index_length):
        order = index_length - i - 1
        ans = mp.fadd(ans, mp.fmul(eq[i], mp.power(x, order)))
    return mp.nstr(ans)

def doubledouble2string(doubdoub,n=19):
    """
    Convert a double-double tuple into a string with n decimal digits of precision.  Default is n=19,
    which is the maximum precision that this longdouble-based double-double format can provide.
    min_fixed is the position of the leading digit such that any number with a leading
    digit less than or equal to min_fixed will be written in exponential format.
    e.g., 0.00123 has its leading digit in the -3 place, so if min_fixed>=-3
    it will be printed as 1.23e-3, and if min_fixed<=-4 it will be printed as 0.00123.
    """
    mpmath.mp.prec=105
    hi=mpmath.mpmathify(doubdoub[0])
    lo=mpmath.mpmathify(doubdoub[1])
    tot=mpmath.fadd(hi,lo)
    return mpmath.nstr(tot,n)

def list_zeros(N=None,
               t=None,
               limit=None,
               fmt=None,
               download=None):
    if N is None:
        N = request.args.get("N", None, int)
    if t is None:
        t = request.args.get("t", 0, float)
    if limit is None:
        limit = request.args.get("limit", 100, int)
    if fmt is None:
        fmt = request.args.get("format", "plain")
    if download is None:
        download = request.args.get("download", "no")

    if limit < 0:
        limit = 100
    if N is not None:  # None is < 0!! WHAT THE WHAT!
        if N < 0:
            N = 0
    if t < 0:
        t = 0

    if limit > 100000:
        # limit = 100000
        #
        bread = [("L-functions", url_for("l_functions.l_function_top_page")),("Zeros of $\zeta(s)$", url_for(".zetazeros"))]
        return render_template('single.html', title="Too many zeros", bread=bread, kid = "dq.zeros.zeta.toomany")

    if N is not None:
        zeros = zeros_starting_at_N(N, limit)
    else:
        zeros = zeros_starting_at_t(t, limit)

    if fmt == 'plain':
        response = flask.Response(("%d %s\n" % (n, nstr(z,31+floor(log(z,10))+1,strip_zeros=False,min_fixed=-inf,max_fixed=+inf)) for (n, z) in zeros))
        response.headers['content-type'] = 'text/plain'
        if download == "yes":
            response.headers['content-disposition'] = 'attachment; filename=zetazeros'
    else:
        response = str(list(zeros))

    return response

def main():
    print(__doc__)
    K = 25
    N = 25
    with mp.workdps(50):
        d = compute_d(K, N)
    fn = os.path.join(os.path.dirname(__file__), "..", "cephes", "igam.h")
    with open(fn + ".new", "w") as f:
        f.write(header.format(K, N))
        for k, row in enumerate(d):
            row = map(lambda x: mp.nstr(x, 17, min_fixed=0, max_fixed=0), row)
            f.write("{")
            f.write(", ".join(row))
            if k < K - 1:
                f.write("},\n")
            else:
                f.write("}};\n")
        f.write(footer)
    os.rename(fn + ".new", fn)

def generate_cases():
    fmt = "{}{:04d} {} {!r} 0.0 -> {} {!r}"
    for fn in sorted(cases_to_generate.keys()):
        print "-- Additional real values (Jython)"
        count, xlist = cases_to_generate[fn]
        for x in xlist:
            # Compute the function (in the reference library)
            func = getattr(mpmath, fn)
            y = func(x)
            # For the benefit of cmath tests, get the sign of imaginary zero right
            zero = 0.0
            if math.copysign(1.0, x) > 0.0:
                if fn == "cos":
                    zero = -0.0
            else:
                if fn == "cosh":
                    zero = -0.0
            # Output one test case at sufficient precision
            print fmt.format(fn, count, fn, x, mpmath.nstr(y, 20), zero)
            count += 1

    def number_form(n):
        if rational(n):
            p, q = rational(n)
            if q == 1:
                if abs(n) < 2**31-1:
                    return str(n), "i"
                if abs(n) < 2**63-1:
                    return str(n), "j"
            else:
                forms, types = zip(*map(number_form, (p, q)))
                if all([t in ("i", "j") for t in types]):
                    return "%".join(map(encode_item, forms, types)), "S"
            return number_form(p / Decimal(q))

        if isinstance(n, Decimal):
            if abs(n) < 10 ** Decimal(-20):
                return "0", "f"
            q = abs(n).log10().quantize(1, rounding=decimal.ROUND_FLOOR)
            q = n.quantize(10 ** (q - 20)).as_tuple()
            l = len(q.digits) + q.exponent
            q = "-" * q.sign + \
                ''.join(map(str, q.digits[:max(0, l)])) + \
                "." + "0" * max(0, -l) + \
                ''.join(map(str, q.digits[max(0, l):]))
            return q.rstrip("0").rstrip("."), "f"

        if isinstance(n, (mp.mpf, mp.mpc)) and not n.imag:
            return number_form(Decimal(mp.nstr(n.real, mp.mp.dps)))

        if complex_pair and isinstance(n, mp.mpc):
            return list_form([n.real, n.imag])

        if n is None:
            return "0n", "f"

        raise Exception()