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 == 1:
        return -1/2
    if m % 2:
        return 0

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

    tm = verbose('computing prec')
    t = log(d) + (m + 0.5) * log(m) - m*(1+log(2*pi)) + 1.712086
    u = t/(log(2)*4)
    prec0 = int(u.ceil()) + 3
    prec = 4*(int(u.ceil()) + 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(2)*m_factorial)/((pari(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 % 4 == 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 == 1:
        return -1/2
    if m % 2:
        return 0

    t0 = verbose('computing bernoulli number %s'%m)
    
    tm = verbose('computing d...')
    P = prime_range(m+2)
    #d = prod([p_minus_1 + 1 for p_minus_1 in divisors(m) if is_prime(p_minus_1+1)])
    d = 6*prod([2*d + 1 for d in divisors(m//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-1) == 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 + 0.5) * log(m) - m*(1+log(2*pi)) + 1.712086
    u = t/(log(2)*4)
    prec0 = int(u.ceil()) + 3
    
    prec = 4*(int(u.ceil()) + 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(2*m_factorial))/((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 % 4 == 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-1)) / (n-1))
    lim = 1 + long(D.ceil())
    prec = prec + 1
    R = RealField(int((prec)))
    print lim, prec
    n = R(n)
    z = 1 - 1/R(2)^n
    for p in prime_range(3, lim+1):
        #print p
        h = z/(R(p)^R(n))
        z = z - h
    return z