os.chdir("/home/was/.sage//spyx/factorial"); from factorial_0 import *; os.chdir("/home/was/papers/bernoulli/talk")

def q_zeta(m, N):
    return ZZ(1)/prod([ZZ(1)-ZZ(1)/(p**m) for p in primes(N)])

def r_zeta(m, prec):
    w = RealField(prec)(m)
    return w.zeta()

def bernoulli1(m):
    r"""
    Returns the Bernoulli number $B_m$; do every detail directly.
    """
    if m == ZZ(1):
        return -ZZ(1)/ZZ(2)
    if m % ZZ(2):
        return ZZ(0)

    prec = ZZ(16)*m
    verbose('prec = %s'%prec)
    tm = verbose('bernoulli_python setup')
    R =  RealField(prec)
    tm = verbose('computing pi...', tm)
    pi = R.pi()
    tm = verbose('computing factorial...', tm)
    m_factorial = factorial(m)
    tm = verbose('computing pi pow...', tm)
    K = ZZ(2)*m_factorial/((ZZ(2)*pi)**m)
    tm = verbose('computing P...', tm)
    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)])
    
    tm = verbose('computing N...', tm)
    dK = d*K
    n = str(dK).find('.') + ZZ(1)
    N = float(ZZ(10))**(n/(m-ZZ(1)))
    tm = verbose('N = %s;\ncomputing product...'%N, tm)
    assert N < m
    z = ZZ(1)
    z_prev = z
    Rl = RealField()
    for p in P:
        if p > N:
            break
        z /= ZZ(1) - ZZ(1)/(R(p)**R(m))
        diff = abs(z - z_prev).log()
        print p, Rl(diff)
        if diff < -prec:
            break
        z_prev =z

    print z
    a = long((dK*z).ceil())
    if m % ZZ(4) == ZZ(0):
        a = -a
    return Rational(a)/Rational(d)
    #return K

def bernoulli2(m):
    r"""
    Returns the Bernoulli number $B_m$.  Use PARI and MPFR when
    possible for parts of the calculation.

    OBSERVATIONS:
       * Pi:
         MPFR computes pi *much* more quickly than PARI:
         sage: time s = pari.new_with_bits_prec (pi, 10^5)
         CPU times: user 2.48 s, sys: 0.01 s, total: 2.49 s
         Wall time: 2.67
         sage: time p = RealField(10^5)(pi)
         CPU times: user 0.20 s, sys: 0.00 s, total: 0.20 s

       * ZETA:
         PARI computes zeta *much* more quickly than MPFR
         sage: time m = zeta(RealField(10000)(1000))
         CPU times: user 5.12 s, sys: 0.00 s, total: 5.13 s
         Wall time: 5.22
         sage: time n = pari.new_with_bits_prec(1000, 10000).zeta()
         CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
         Wall time: 0.01
         sage: time n=pari.new_with_bits_prec(1000, 20000).zeta()
         CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
         Wall time: 0.06
         sage: time m=zeta(RealField(20000)(1000))
         CPU times: user 139.82 s, sys: 0.65 s, total: 140.46 s
         Wall time: 164.22

       * Factorials: My Pyrex compiled wrapping of GMP is best:
        sage: time n = factorial(100000)
        CPU times: user 6.93 s, sys: 0.00 s, total: 6.93 s
        Wall time: 6.95

        MAGMA:
        sage: time n = magma('Factorial(100000)')
        CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
        Wall time: 9.02

        GP interpreter
        sage: time gp('m=prod(n=1,100000,n);')
        Wall time: 16.00

        PARI C library:
        sage: time n=pari('prod(n=1,100000,n)')
        CPU times: user 17.64 s, sys: 0.22 s, total: 17.85 s
        Wall time: 18.81
    """
    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('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 zeta... (prec=%s)'%prec)
    #w = pari.new_with_bits_prec(m, prec).zeta()
    w = pari.new_with_prec(m, prec0).zeta()
    verbose('got zeta...', tm)
    #z = R(w)
    tm = verbose("converting...")
    z = R(w)
    verbose('converted', tm)

    tm = verbose('constructing bernoulli number')
    dKz = d*K*z
    #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 ndigits(m):
    """
    Compute the number of digits of the numerator of the m-th
    Bernoulli number without actually computing the m-th Bernoulli
    number.
    """
    P = prime_range(m+ZZ(2))
    d = prod([p for p in P if m % (p-ZZ(1)) == ZZ(0)])
    x = log(ZZ(2)) + sum(math.log(n) for n in range(ZZ(1), m+ZZ(1)))  \
        - m*log(ZZ(2)) - m*log(pi) + log(d)
    y = x/log(ZZ(10))
    return ZZ(int(y.ceil()))

def ndigits2(m):
    """
    Compute the number of digits of the numerator of the m-th
    Bernoulli number without actually computing the m-th Bernoulli
    number.
    """
    P = prime_range(m+ZZ(2))
    d = prod([p for p in P if m % (p-ZZ(1)) == ZZ(0)])
    x = log(ZZ(2)) + sum(math.log(n) for n in range(ZZ(1), m+ZZ(1)))  \
        - m*log(ZZ(2)) - m*log(pi) + log(d)
    y = x/log(ZZ(10))
    return ZZ(int(y.ceil()))




def genber_numerical(n, chi, prec, primes):
    if not chi.is_primitive():
        raise ValueError, "chi (=%s) must be primitive"%chi
    CC = ComplexField(prec)
    K = genber_K(n, chi, CC)
    a = CC.pi()**n
    g = chi.gauss_sum_numerical(prec)
    L = genber_L(n, chi**(-ZZ(1)), CC, primes)
    d = genber_d(n, chi, prec)
    print "d = ", d
    print "K = ", K
    print "L = ", L
    print "a = ", a
    print "g = ", g
    num = d*K*L/(a*g)
    return num

def genber_d(n, chi, prec):
    f = chi.conductor()
    if f == ZZ(1):
        raise NotImplementedError
    F = f.factor()
    assert len(F) >= ZZ(1)
    if len(F) > ZZ(1):
        return ZZ(1)
    elif f == ZZ(4):
        return ZZ(2)
    P = F[0]
    if P[0] == ZZ(2) and P[1] > ZZ(2):
        return ZZ(1)
    elif P[1] == ZZ(1):
        return n*P[0]
    elif P[1] > ZZ(1):
        R = chi.base_ring()
        phi = R.complex_embedding(prec)
        return phi(ZZ(1) - chi(ZZ(1)+p))


def genber_K(n, chi, CC):
    f = chi.modulus()
    i = CC.gen(0)
    z = ZZ(2)*factorial(n)*(f/(ZZ(2)*i))**n
    if n % ZZ(2) == ZZ(0):
        z *= -ZZ(1)
    return z

def genber_L(n, chi, CC, max_prime):
    R = chi.base_ring()
    m = chi.modulus()
    phi = R.complex_embedding(CC.prec())
    n = phi(n)
    z = ZZ(1)
    zero = ZZ(2)**(-CC.prec())
    for p in prime_range(max_prime):
        h = z * phi(chi(p)) / (CC(p)**n)
        z = z - h
    return ZZ(1)/z

def test1(prec, primes):
    e = DirichletGroup(ZZ(3)).gen(0)
    print "group: ", e.parent()
    print "chi = ", e
    B = genber_numerical(ZZ(3), e, prec, primes)
    print "B3 = ", B
    return e, B