from scipy import *

M = 100
J = 6
K = 10

#test
R = 3.
N = 10

def g_k( k, M, J ):
  rng = 1.*arange( 1, M, 1 )
  s1 = add.reduce( (log(rng)**k) / rng )
  s1 = s1+(.5*(log(M)**k))/M - (log(M)**(k+1))/(k+1)

  B = special.bernoulli(2*J)
  s2 = 0.

  P = 1.*zeros( k+1 )
  P[0] = 1.
  P_0_k = poly1d(P)
  P_k = P_0_k.deriv()-P_0_k

  j = 1

  while 1:
    tmp = B[2*j]/(2.*j)/(M**(2.*j))*P_k(log(M))
    s2 = s2-tmp

    if( j>=J ):
      break

    tmp = P_k.deriv()/(2.*j)-P_k
    P_k = tmp.deriv()/(2.*j+1)-tmp

    j=j+1

  result = (-1)**k/special.gamma(k+1.)*(s1+s2)

  return result


g_vec = poly1d([g_k( K-i, M, J ) for i in range( K+1 )])

def zeta( z ):
  return 1./(z-1.)+g_vec(z-1.)

def test(R, N):
  print "z              |zeta(z)        |1.+special.zetac(z)"
  print "==================================================="
  for z in (R*(2.*rand(N)-1.)+1.):
    print "%15.10g|%15.10g|%15.10g"%( z, zeta(z), 1.+special.zetac(z) )


if __name__ == "__main__":
  test( R, N )