#!/usr/bin/env python

# Case of non-standard installation
# python setup.py install --home=~/local
import sys;
import os;
homedir = os.environ["HOME"]
sys.path.append(homedir +  "/local/lib/python")
#print sys.path;

from mpmath import *
from optparse import OptionParser

parser = OptionParser()
parser.add_option("-l", "--l", dest="L",type="int",help="L param", metavar="L",default=4)
parser.add_option("-x", "--x", dest="X",type="float",help="X param", metavar="X",default=0.5)
parser.add_option("-m", "--m", dest="M",type="int",help="M param", metavar="M",default=20)
parser.add_option("-n", "--num", dest="Num",type="int",help="Num param", metavar="Num",default=30)
parser.add_option("-p", "--precision", type="int",dest="Precision",help="Precision",
                  metavar="P",default=60)
parser.add_option("-d", "--d", action="store_true", dest="debug",help="Print debug info", default=False)
parser.add_option("-q", "--quiet",
                  action="store_false", dest="verbose", default=True,
                  help="don't print status messages to stdout")

(options, args) = parser.parse_args()

mp.dps = options.Precision; mp.pretty = True

def g(x, Num):
   mp.dps = options.Precision;
   H=0
   twoSq2X = 2.0*sqrt(2.0)*x;
   N=Num;
   if (x < 1): N=floor(10/x + Num)
   for n in range(1,N):
      twoSq2XN=twoSq2X*n;
      b1 = besselk(1,twoSq2XN)
      b0 = besselk(0,twoSq2XN)
      h = 4*(twoSq2XN*b1 - b0)
      H += h
      if(options.debug):
         print '  n=%d, twoSq2XN=%f, K1=%f, K0=%f, h=%f, H=%f' % (n,twoSq2XN,b1,b0,h,H)

   return exp(-H)

def Xsi(M):
   mp.dps = options.Precision;
   NV=2*M
   sign = -1;
   if ((M % 2) == 0): sign = 1
   vec = []
   for k in range(1,NV+1):
      Vk=mpf('0.0')
      jmin=(k+1)/2
      jmax=k
      if jmax > M: jmax=M
      sign = -sign
      for j in range(jmin,jmax+1):
         Vk += (mpf(j)**mpf((M+1)))/factorial(mpf(M))*binomial(mpf(M),mpf(j))*binomial(mpf(2*j),mpf(j))*binomial(mpf(j),mpf(k-j))
      Vk *= sign
      vec.append(Vk)
   return vec

def invLTGfunction(L,M,Num):
   mp.dps = options.Precision;
   NV=2*M
   X = [i/(float(L)) for i in range(1, L+1)]
   F = []
   for l in range(1,L+1):
      F.append(0.0)
      for k in range(1,NV+1):
         F[l-1]+= xsi[k-1]/mpf(k) * g(sqrt(k*log(2))*X[l-1],Num)
   return F

def invLTGfunctionPointEst(X,M,Num):
   mp.dps = options.Precision;
   NV=2*M
   F = 0.0
   for k in range(1,NV+1):
      if(options.debug): print "k=",k-1
      s = sqrt(k*log(2))*X
      if(options.debug): print "s=",s
      fval = g(sqrt(k*log(2))*X,Num)
      if(options.debug): print "fval=",fval
      xsiNorm = xsi[k-1]/mpf(k)
      if(options.debug): print "xsiNorm=",xsiNorm
      prod = xsiNorm * fval
      if(options.debug): print "prod=",prod
      F += prod
      if(options.debug): print "F=",F
   return F


xsi = Xsi(options.M)

if(options.verbose):
   print "X =",options.X
   print "Precision =",options.Precision
   print "M =",options.M
   print "Num =",options.Num
   print "SumXsi=", sum(xsi)
   print "xsi=", xsi

print invLTGfunctionPointEst(options.X,options.M,options.Num)