# bernoulli numbers of higher order

from PyM import *


def vprod(x,y,d=''):
    if d=='': d = len(x)+len(y)
    if d==0: return []
    z = d*[0]
    for j in range(min(d,len(y))): z[j]=y[j]
    for j in range(min(d,len(x))):
        z[j] += x[j]
        for k in range(min(d-1-j,len(y))):
            z[j+k+1] += x[j]*y[k]
    return vector(z)

show('vprod examples')    
x = [1,2,3]
y = [7,8,9,10]

show(vprod(x,y))
show(vprod(x,y,10))
show(vprod(x,y,5))

show(vpower(x,2))
show(vpower(x,3,12))
show(len(vpower(y,100)))

def vpower(x,m,d=''):
    if d=='': d=len(x)
    if m==0: return []
    if m==1: return pad(x,d)
    z = m*[0]
    while m>0:
        if m%2: z = vprod(z,x,d)
        x = vprod(x,x,d)
        m = m//2
    return vector(z)
    
# Given list or vector [c0,c1,...,c(r-1)], 
# computes list or vector  
# [s0,s1,..,s(r-1)s such that 
# 1+s0*t+s1*t**2+...+s(r-1)t**r is inverse of
# 1+c0*t+c1*t**2+...+c(r-1)t**r mod t**(r+1)
def invert_vector(c, r=''): 
    if isinstance(c,Vector_type): 
        return vector(invert_vector(list(c),r))
    if len(c)==0: return invert_vector([0],r)
    if r==0: return []
    if r=='': r=len(c)
    c = pad(c,r)
    s = [-c[0]]
    #if r==1: return s
    for k in range(1,r):
        s += [-c[k]-convolution(s,c,k-1)]
    return vector(s) 

show('Invert vector examples')
x = vector([3,-1,1])
y = invert_vector(x)

show('y =',y)

show(vprod(x,invert_vector(x,15),10))
show(vprod(x,invert_vector(x,15),18))

def bernoulli_numbers(N):
    u = 1>>Q_
    if N==0: return [u]
    B = invert_vector([u/factorial(j+1) for j in range(1,N+1)])
    return [u]+[factorial(k+1)*B[k] for k in range(N)]
#
BV_ = bernoulli_numbers

def bernoulli_number(N):
    if N==0: return 1>>Q_
    elif N==1: return -1/2>>Q_
    else:
        if odd(N): return 0>>Q_
        else: return bernoulli_numbers(N)[-1]
#
B_=bernoulli_number

show('Bernoulli numbers examples')
show('BV_(20) =',bernoulli_numbers(20))

show('B_(100) =',B_(100))
 
# Taylor coeffs of x/(e^x - 1) at 0
# -1/2, 1/12, 0, -1/720, 0, 1/30240, 0, -1/1209600, 
# 0, 1/47900160, 0, -691/1307674368000, 0, 1/74724249600, 
# 0, -3617/10670622842880000, 0, 43867/5109094217170944000, 
# 0, -174611/802857662698291200000
def todd_numbers(n):
    u = 1>>Q_
    return invert_vector([u/factorial(j+1) for j in range(1,n+1)])

show('Todd Number 21:',todd_numbers(21))


def high_bernoulli_numbers(N,k):
    u = 1>>Q_
    if N==0: return [u]
    B = invert_vector([u/factorial(j+1) for j in range(1,N+1)])
    B = vpower(B,k,N)
    return [u]+[factorial(j+1)*B[j] for j in range(N)]
#    
HBs_=high_bernoulli_numbers

show('High Bernoulli Numbers examples')
show(HB_(5,4))

N=10
for p in [2*j+1 for j in range(1,N+1,2) if is_prime(2*j+1)]:
    show(p**3,HB_(p+2,p+1))

show(HB_(4,3))

show(HB_(3,3)-3**2/2,HB_(5,5)-5**2/2)

def bernoulli_polynomial(n,x='x'):
    u = 1>>Q_
    [Qx,x] = polynomial_ring(Q_,x)
    P = x**n+n*B_(1)*x**(n-1)
    ntop = n
    if even(ntop): ntop+=1
    for k in range(2,ntop,2):
        P += binom(n,k)*B_(k)*x**(n-k)
    return P
#
BP_ = bernoulli_polynomial

k = 6;
B = bernoulli_polynomial(k,'x')
show('B6(x) = ',B)
show('B(6) =',B_(k))

# B(N,k) 
def high_bernoulli_number(N,k):
    if N==0: return 1>>Q_
    return high_bernoulli_numbers(N,k)[-1]
#
HB_=high_bernoulli_number
    
# I checked that B(p+2,p+1) is divisible by p**3 
# for all primes p between 5 and 101 (Carlitz theorem, 1951)
# For p=3, B(5,4)=-9 is divisible by 9.
# On the the otehr hand, B(2+2,2+1) = 19/10.
# I have also checked that B(p,p)-p**2/2 is divisible by p**3 for p>=3
# in a few cases, and other similar relations in the same paper
# There are very many.

def zhe(n,j):
    if igcd(j,2*n)!=1: return 'zhe: coprime condition not satisfied'
    k=0
    for _ in range(j):
        k=2**k
        #print(k)
    #print(k)
    return HB_(n,k)/(2**(2*n))

show('Zhe examples')
N = [n for n in range(1,20) if igcd(n,3)==1 and igcd(n,5)==1]

for n in N: print(abs(zhe(n,3)))
print()
for n in N: print(abs(zhe(n,5)))

from math import pi,log
euler_gamma=0.577215664901532
kappa = pi/euler_gamma/log(2)

def atiyah(n):
    def t(j):
        x = (j*log(j)-j+(log(j)+1)/j)/log(2)
        return (1-x)/2**(j+1)
    s = sum([t(j) for j in range(1,n+1)])
    return s*kappa

print()
show('Atiyah examples')
print(atiyah(50))
print(atiyah(100))
print(atiyah(200))