def bernoulli_pari_ish(m):
    r"""
    Returns the Bernoulli number $B_m$.  Use PARI and MPFR when
    possible for parts of the calculation.
    """
    if m == ZZ(1):
        return -ZZ(1)/ZZ(2)
    if m % ZZ(2):
        return ZZ(0)

    t0 = verbose('computing bernoulli number %s'%m)
    
    tm = verbose('computing d...')
    P = prime_range(m+ZZ(2))
    # IDIOTIC -- should compute by factoring m!
    d = prod([p for p in P if m % (p-ZZ(1)) == ZZ(0)])
    verbose('got d=%s'%d, tm)

    tm = verbose('computing prec')
    t = log(d) + (m + RR('0.5')) * log(m) - m*(ZZ(1)+log(ZZ(2)*pi)) + RR('1.712086')
    u = t/(log(ZZ(2))*ZZ(4))
    prec0 = int(u.ceil()) + ZZ(3)
    prec = ZZ(4)*(int(u.ceil()) + ZZ(12))
    verbose('prec = %s; prec0 = %s'%(prec,prec0), tm)

    #prec = 16*m
    #verbose('prec = %s'%prec)

    tm = verbose('bernoulli_python setup')
    R = RealField(prec)
    tm = verbose('computing pi...')
    PI = R.pi()
    verbose('got pi...', tm)
    tm = verbose('converting to PARI')
    PI = pari(PI)
    verbose('done', tm)
    
    tm = verbose('computing factorial...')
    m_factorial = pari(factorial_gmp(m))
    verbose('got factorial', tm)
    
    tm = verbose('computing K ...')
    K = (pari(ZZ(2))*m_factorial)/((pari(ZZ(2))*PI)**pari(m))
    verbose('got K', tm)


    tm = verbose('computing zeta...')
    z = pari.new_with_prec(m, prec0).zeta()
    verbose('got zeta...', tm)

    tm = verbose('constructing bernoulli number')
    dKz = pari(d)*K*z
    a = long((dKz).round())
    if m % ZZ(4) == ZZ(0):
        a = -a
    B = Rational(a)/Rational(d)
    tm = verbose('done!', tm)
    verbose('total time', t0)
    return B





def bernoulli_invzeta(m):
    r"""
    Returns the Bernoulli number $B_m$.  Use PARI and MPFR when
    possible for parts of the calculation.
    """
    if m == ZZ(1):
        return -ZZ(1)/ZZ(2)
    if m % ZZ(2):
        return ZZ(0)

    t0 = verbose('computing bernoulli number %s'%m)
    
    tm = verbose('computing d...')
    P = prime_range(m+ZZ(2))
    #d = prod([p_minus_1 + 1 for p_minus_1 in divisors(m) if is_prime(p_minus_1+1)])
    d = ZZ(6)*prod([ZZ(2)*d + ZZ(1) for d in divisors(m//ZZ(2))[1:] if is_prime(2*d+1)])
    tm = verbose('got d=%s'%d, tm)
    d0 = prod([p for p in P if m % (p-ZZ(1)) == ZZ(0)])
    assert d==d0, 'd=%s, d0=%s'%(d,d0)
    verbose('got d=%s old way'%d, tm)

    tm = verbose('computing prec')

    t = log(d) + (m + RR('0.5')) * log(m) - m*(ZZ(1)+log(ZZ(2)*pi)) + RR('1.712086')
    u = t/(log(ZZ(2))*ZZ(4))
    prec0 = int(u.ceil()) + ZZ(3)
    
    prec = ZZ(4)*(int(u.ceil()) + ZZ(12))
    verbose('prec = %s; prec0 = %s'%(prec,prec0), tm)

    #prec = 16*m
    #verbose('prec = %s'%prec)

    tm = verbose('bernoulli_python setup')
    R = RealField(prec)
    tm = verbose('computing pi...')
    PI = R.pi()
    verbose('got pi...', tm)
    
    tm = verbose('computing factorial...')
    m_factorial = factorial_gmp(m)
    #m_factorial = R.factorial(m)
    verbose('got factorial', tm)
    
    tm = verbose('computing K ...')
    K = (R(ZZ(2)*m_factorial))/((ZZ(2)*PI)**m)
    verbose('got K', tm)


    tm = verbose('computing inv zeta...')
    invz = inverse_zeta(m, t, prec)
    verbose('got zeta...', tm)

    tm = verbose('constructing bernoulli number')
    dKz = d*K/invz
    #print 'frac = ', str(dKz.frac())[:30]
    a = long((dKz).round())
    if m % ZZ(4) == ZZ(0):
        a = -a
    B = Rational(a)/Rational(d)
    tm = verbose('done!', tm)
    verbose('total time', t0)
    return B

def inverse_zeta(n, t, prec):
    D = exp((t - log(n-ZZ(1))) / (n-ZZ(1)))
    lim = ZZ(1) + long(D.ceil())
    prec = prec + ZZ(1)
    R = RealField(int((prec)))
    print lim, prec
    n = R(n)
    z = ZZ(1) - ZZ(1)/R(ZZ(2))**n
    for p in prime_range(ZZ(3), lim+ZZ(1)):
        #print p
        h = z/R(p)**n
        z = z - h
    return z