#####################################################################
#  Computation of Bernoulli numbers
#  Requires SAGE >= 1.2.2
#  (c) William Stein
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#####################################################################

def bern_recur(j):
    """
    Compute Bernoulli numbers $B_n$ for n = 0, 1, ...,j
    using the recurrence.
    """
    B = [1, -1/2]
    for n in xrange(2,j+1):
        if n % 2:
            B.append(0)
        else:
            B.append(-sum( binomial(n+1,k)*B[k] for k in xrange(n) ) / (n+1))


    return B

def bern_denom(n):
    """
    Return the denominator of the Bernoulli number $B_n$.

    EXAMPLES:
        sage: bern_denom(100)
        33330
        sage: bernoulli(100)
        -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330
        sage: bern_denom(1000)
        342999030
        sage: bern_denom(2006)
        6
        sage: bern_denom(10000)
        2338224387510
        sage: bern_denom(100000)
        9355235774427510
        sage: bern_denom(100001)
        2
    """
    if n < 0:
        raise ValueError, "n (=%s) must be >= 0"%n
    if n == 0:
        return 1
    # Write d = p-1.  Then p = d+1 for a divisor d of n.
    # So to compute the p-1 | n we compute the divisors d and
    # check primality of d+1.
    return prod([d+1 for d in divisors(n) if is_prime(d+1)])

def bern_numer_digits(n):
    """
    Return the number of digits of the numbers of $B_n$, computed
    using the formula in Remark~\ref{rem:bernnum} of the book.

    EXAMPLES:
        sage: bern_numer_digits(1000)
        1779
        sage: len(str(abs(bernoulli(1000).numerator())))
        1779
        sage: bern_numer_digits(2006)
        4156
        sage: len(str(abs(bernoulli(2006).numerator())))
        4156
    """
    # Note that easy estimates show that B_n has about n*log(n)
    # decimal digits.
    bits = log(2*n*log(n), 2).floor()+10
    R = RealField(bits)   
    PI = R.pi()

    logd = log(R(bern_denom(n)))

    z = R(n + 1)
    
    # Compute log(Gamma(z)) quickly using Stirling's formula
    log_nfactorial = 1/log(2*PI) + (z - 1/2)*log(z) - z
    m = 1; a = 9999; zero = R(2)^(-bits)
    while abs(a) > zero:
        a = bernoulli(2*m)/(2*m*(2*m-1)*z^(2*m-1))
        log_nfactorial += a
        m += 1
                

    logB = (1 - n)*log(R(2))  -  n*log(PI) + log_nfactorial
    
    return ceiling((logd + logB)/ log(R(10)))
    


def bern_analytic(n, verbose=False):
    r"""
    Compute the Bernoulli number $B_n$ using the analytic algorithm.
    The implementation below uses the bounds from the book, which we
    proved are correct.

    NOTE: The implementation of bernoulli in the PARI C library, which
    is what the \code{bernoulli} command in \sage uses, invokes better
    bounds whose derivation is mysterious and undocumented, but seem
    to work in practice.
    """

    # Special cases
    if n < 0:
        raise ValueError
    if n == 0:
        return 1
    elif n == 1:
        retrn -1/2
    elif n % 2:
        return 0
    
    # Number of bits of precision needed 
    bits = ceiling(1.5*log(10,2)*(bern_numer_digits(n)))
    if verbose: print "bits = ", bits
    R = RealField(bits)
    PI = R.pi()

    # Factorial factor
    
    K = 2 * factorial(n) / (2 * PI)^n
    if verbose: print "K = ", K

    d = bern_denom(n)  # waste -- this is computed twice; also
                       # computed in finding bern_numer_digits
    if verbose: print "d = ", d

    S = RealField(log(2*n*log(n), 2).floor()+10)
    M = ceiling((S(K)*d)^(1/(n-1)))
    
    if verbose: print "M = ", M

    t = cputime()
    one = R(1)
    pp = [R(p^n) for p in prime_range(M+1)]
    v = [one - one/w for w in pp]   
    print "time to compute Euler factors = ", cputime(t)
    t = cputime()
    z = 1/prod( v )
    if verbose:
        print "z = ", z
        print "time for product = ", cputime(t)

    dKz = d * K * z
    print "d*K*z = ", dKz
    a = ceiling(dKz);
    if n % 4 == 0:  a = -a

    return a / d