""" Log-sine integrals package for evaluating log-sine integrals in polylogarithmic terms accompanying the paper "Special values of generalized log-sine integrals" -- Jonathan M. Borwein, University of Newcastle -- Armin Straub, Tulane University -- Version 1.0 (2011/07/13) """ #***************************************************************************** # Copyright (C) 2011 Armin Straub # # Distributed under the terms of the GNU General Public License (GPL) # # This code is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # The full text of the GPL is available at: # # http://www.gnu.org/licenses/ #***************************************************************************** from sage.symbolic.function import BuiltinFunction class ls_Function(BuiltinFunction): r""" The log-sine integral. EXAMPLES:: sage: ls(5,1,pi/3) ls(5, 1, 1/3*pi) AUTHORS: - Armin Straub (2011) REFERENCES: - J. M. Borwein, A. Straub. "Special values of generalized log-sine integrals." Proceedings of ISSAC 2011 (International Symposium on Symbolic and Algebraic Computation), 2011. """ def __init__(self): r""" Constructs an object to represent the log-sine integral. EXAMPLES:: sage: ls(5,1,x) ls(5, 1, x) """ BuiltinFunction.__init__(self, "ls", nargs=3) # BuiltinFunction.__init__(self, "ls", nargs=3, latex_name="\\textrm{Ls}") def _print_latex_(self, n, k, x, **kwds): r""" Return latex expression EXAMPLES:: sage: from sage.misc.latex import latex sage: latex(ls(5,1,pi)) \textrm{Ls}_{5}^{(1)}(\pi) sage: latex(ls(5,0,pi)) \textrm{Ls}_{5}(\pi) """ from sage.misc.latex import latex if k==0: return "\\textrm{Ls}_{%s}(%s)" % (latex(n), latex(x)) return "\\textrm{Ls}_{%s}^{(%s)}(%s)" % (latex(n), latex(k), latex(x)) def _evalf_(self, n, k, s, parent=None): r""" EXAMPLES:: sage: n(ls(5,0,pi)) -24.226558345665453 sage: ls(3,1,pi/2).n() 0.12404565613170646 """ num = numerical_integral(lambda t: t^k * log(abs(2*sin(t/2)))^(n-k-1), [0,s]) return -num[0] # print s.prec() # Doesn't work: eps_abs=2^(-s.prec()), eps_rel=2^(-s.prec()) # E.g.: numerical_integral(x^2,0,1,eps_abs=10^(-100), eps_rel=10^(-100)) class cln_Function(BuiltinFunction): r""" The Clausen function for Nielsen polylogarithms. EXAMPLES:: sage: cl_nielsen(1,2,pi) cl_nielsen(1, 2, pi) """ def __init__(self): r""" Constructs an object to represent the Clausen function. EXAMPLES:: sage: cl_nielsen(1,2,x) cl_nielsen(1, 2, x) """ BuiltinFunction.__init__(self, "cl_nielsen", nargs=3) # BuiltinFunction.__init__(self, "cl_nielsen", nargs=3, latex_name="\\textrm{Cln}") def _print_latex_(self, n, p, x, **kwds): r""" Return latex expression EXAMPLES:: sage: from sage.misc.latex import latex sage: latex(cl_nielsen(3,2,pi)) \textrm{Cl}_{4,1}(\pi) sage: latex(cl_nielsen(5,1,pi)) \textrm{Cl}_{6}(\pi) """ from sage.misc.latex import latex return "\\textrm{Cl}_{%s%s}(%s)" % (latex(n+1), (p-1)*",1", latex(x)) # return "\\textrm{Cln}_{%s,%s}(%s)" % (latex(n), latex(p), latex(x)) class gln_Function(BuiltinFunction): r""" The Glaisher function for Nielsen polylogarithms. EXAMPLES:: sage: gl_nielsen(1,2,pi) gl_nielsen(1, 2, pi) """ def __init__(self): r""" Constructs an object to represent the Glaisher function. EXAMPLES:: sage: gl_nielsen(1,2,x) gl_nielsen(1, 2, x) """ BuiltinFunction.__init__(self, "gl_nielsen", nargs=3) def _print_latex_(self, n, p, x, **kwds): r""" Return latex expression EXAMPLES:: sage: from sage.misc.latex import latex sage: latex(gl_nielsen(3,2,pi)) \textrm{Gl}_{4,1}(\pi) sage: latex(gl_nielsen(5,1,pi)) \textrm{Gl}_{6}(\pi) """ from sage.misc.latex import latex return "\\textrm{Gl}_{%s%s}(%s)" % (latex(n+1), (p-1)*",1", latex(x)) # return "\\textrm{Gln}_{%s,%s}(%s)" % (latex(n), latex(p), latex(x)) ls = ls_Function() cl_nielsen = cln_Function() gl_nielsen = gln_Function() # The value 0 as an expression expr0 = Expression(SR) @cached_function def _nielsen_polylog_at1_a(k): if k==0: return 1 return 1/k * sum([(-1)^(m+1)*zeta(m)*_nielsen_polylog_at1_a(k-m) for m in [2..k]]) @cached_function def _nielsen_polylog_at1_b(k): if k==0: return 1 return -1/k * sum([(-1)^(m+1)*zeta(m)*_nielsen_polylog_at1_b(k-m) for m in [2..k]]) def nielsen_polylog_at1(n, p): r""" Returns the Nielsen polylogarithm S(n,p,1) in terms of zeta values. The computation is based on: K. S. Koelbig, Closed expressions for $\int^1_0 t^{-1} \log^{n-1}t \log^p(1-t) dt$, Mathematics of Computation, 39(160):647--654, 1982. INPUT:: - n -- parameter n - p -- parameter p EXAMPLES:: sage: nielsen_polylog_at1(4,3) -1/72*pi^4*zeta(3) - 1/3*pi^2*zeta(5) + 5*zeta(7) """ result = (-1)^(n+p-1) * sum([ _nielsen_polylog_at1_b(n-m-1) * sum([binomial(m+r+1,r)*_nielsen_polylog_at1_b(p-r)*_nielsen_polylog_at1_a(m+r+1) for r in [0..p]]) for m in [0..n-1]]) return result.expand() def lstoli_step(n, k, t): r""" Helper function which rewrites the log-sine integral ls(n,k,t) into a sum of "simpler" log-sine integrals and polylogarithmic terms. EXAMPLES:: sage: lstoli_step(5,3,pi) pi^3*cl_nielsen(1, 1) + 3*pi^2*cl_nielsen(2, 1) - 6*pi*cl_nielsen(3, 1) + 6*zeta(5) - 6*cl_nielsen(4, 1) """ result = 0 if is_odd(k): result += (-1)^(n+(k-1)/2) * factorial(n-k-1) * factorial(k) * nielsen_polylog_at1(n-k-1,k+1) clgl = cl_nielsen if is_even(n+k) else gl_nielsen # Observe that on the next line Cl is only created with odd depth, and Gl only with even depth. # This means that no general polylogarithmic reductions are necessary. result += (-1)^(n+k) * factorial(n-k-1) * factorial(k) \ * sum([(-t)^r/factorial(r) * (-1)^floor((k+r)/2) * clgl(k-r+1,n-k-1,t) for r in [0..k]]) result -= sum([binomial(n-k-1,r)*binomial(r,m)*(1/2)^r*(-1)^(r/2)*(-pi)^(r-m)*ls(n-r+m,k+m,t) \ for r in range(2,n-k,2) for m in [0..r]]) return result def lsgeneraltoli(n, k, t): r""" Helper function which rewrites the log-sine integral ls(n,k,t) into a sum of polylogarithmic terms. No reductions are performed and it is not checked that the argument t is in [0,2pi]. EXAMPLES:: sage: lstoli_step(5,3,pi) pi^3*cl_nielsen(1, 1) + 3*pi^2*cl_nielsen(2, 1) - 6*pi*cl_nielsen(3, 1) + 6*zeta(5) - 6*cl_nielsen(4, 1) """ result = ls(n,k,t) # Reduce the log-sines level by level. for r in range(0,n-k-1,2): def ls_reduce(n0, k0, t0): if (n-k)-(n0-k0) == r: return lstoli_step(int(n0),int(k0),t0) return ls(n0,k0,t0) result = result.substitute_function(ls,ls_reduce).expand() # Possibly remaining log-sines on the last level are of the form ls(m,m-1). These evaluate trivially. def ls_reduce(n0, k0, t0): if n0==k0+1: return -1/n0 * t0^n0 return ls(n0,k0,t0) result = result.substitute_function(ls,ls_reduce).expand() return result def ls_reduceargument(n, k, t): r""" Helper function which rewrites the log-sine integral ls(n,k,t) into a sum of log-sine integrals with arguments either in [0,2pi] or a multiple of 2pi. EXAMPLES:: sage: ls_reduceargument(5,2,3*pi) 4*pi^2*ls(3, 0, pi) + 4*pi*ls(4, 1, pi) + ls(5, 2, pi) + ls(5, 2, 2*pi) """ m = floor(t/(2*pi)) pm = 1 t0 = t - 2*m*pi if t0 > pi: t0 = 2*pi-t0 pm = -1 m += 1 result = ls(n,k,2*m*pi) + sum([pm^(k-j+1)*(2*m*pi)^j*binomial(k,j)*ls(n-j,k-j,t0) for j in [0..k]]) return result def ls_mult2pito2pi(n, k, m): r""" Helper function which rewrites the log-sine integral ls(n,k,m*2pi) into a sum of log-sine integrals at 2pi. EXAMPLES:: sage: ls_mult2pito2pi(5,2,2) 20*pi^2*ls(3, 0, 2*pi) - 12*pi*ls(4, 1, 2*pi) + 2*ls(5, 2, 2*pi) """ return sum([(-1)^(k-j)*(2*pi)^j*binomial(k,j)*ls(n-j,k-j,2*pi)*sum([r^j for r in [1..m]]) for j in [0..k]]) def ls_expand(expr): r""" Helper function which rewrites all log-sine integrals occuring in expr into sums of log-sine integrals with arguments in [0,2pi]. EXAMPLES:: sage: ls_expand(ls(5,3,-4*pi)) 72*pi^3*ls(2, 0, 2*pi) - 60*pi^2*ls(3, 1, 2*pi) + 18*pi*ls(4, 2, 2*pi) - 2*ls(5, 3, 2*pi) sage: ls_expand(ls(5,3,pi)) ls(5, 3, pi) """ result = expr # Treat negative arguments. def ls_reduce(n0, k0, t0): return (-1)^(k0+1)*ls(n0,k0,-t0) if t0<0 else ls(n0,k0,t0) result = result.substitute_function(ls,ls_reduce) # Reduce arguments to range [0,Pi] and multiples of 2Pi. def ls_reduce(n0, k0, t0): return ls_reduceargument(n0,k0,t0) if t0>pi and not t0/2/pi in ZZ else ls(n0,k0,t0) result = result.substitute_function(ls,ls_reduce) # Reduce values at multiples of 2Pi. def ls_reduce(n0, k0, t0): return ls_mult2pito2pi(n0,k0,t0/2/pi) if t0>2*pi and t0/2/pi in ZZ else ls(n0,k0,t0) result = result.substitute_function(ls,ls_reduce) # Do trivial evaluations. def ls_reduce(n0, k0, t0): return -1/n0 * t0^n0 if n0==k0+1 else ls(n0,k0,t0) result = result.substitute_function(ls,ls_reduce) return result.expand() def lstoli(expr): r""" This is the main function of this package. It converts all log-sine integrals in expr into sums of polylogarithmic terms. Reductions are applied. EXAMPLES:: sage: lstoli(-ls(4,2,pi)) 3/2*pi*zeta(3) sage: lstoli(ls(5,1,pi/3)) 1/12*pi^2*zeta(3) + 3/2*pi*cl_nielsen(3, 1, 1/3*pi) - 19/4*zeta(5) sage: lstoli(ls(6,3,pi/3)-2*ls(6,1,pi/3)) 313/204120*pi^6 sage: assume(0=0, x<=2*pi) # sage: lstoli(ls(5,1,x)) # The trouble is that x==2*pi is very slow... in the above example # the difference is >80 vs <1 seconds. if not t0 in RR: return cl_nielsen(n0,p0,t0) # Reduce Nielsen polylogs at 1. if t0==2*pi: return nielsen_polylog_at1(n0,p0) if is_odd(int(n0+p0)) else expr0 # Reduce Nielsen polylogs at -1. if t0==pi and is_even(int(n0+p0)): return expr0 if t0==pi and p0==1: return -(1-2^(-n0))*zeta(n0+1) if t0 == pi/3 and p0==1 and is_even(int(n0)): return 1/2*(1-2^(-n0))*(1-3^(-n0))*zeta(n0+1) if t0 == 2*pi/3 and p0==1 and is_even(int(n0)): return 2^n0/(1-(-2)^n0) * 1/2*(1-2^(-n0))*(1-3^(-n0))*zeta(n0+1) if t0 == pi/2 and p0==1 and is_even(int(n0)): return 2*(1-2^n0)/2^(2*n0+2)*zeta(n0+1) if t0 == pi/3 and (int(n0),int(p0)) in mcvreductions: return mcvreductions[(int(n0),int(p0))] return cl_nielsen(n0,p0,t0) result = result.substitute_function(cl_nielsen,cl_reduce) def gl_reduce(n0, p0, t0): if not t0 in RR: return gl_nielsen(n0,p0,t0) # Reduce Nielsen polylogs at 1. if t0==2*pi: return nielsen_polylog_at1(n0,p0) if is_even(int(n0+p0)) else expr0 # Reduce Nielsen polylogs at -1. if t0==pi and is_odd(int(n0+p0)): return expr0 if t0==pi and p0==1: return -(1-2^(-n0))*zeta(n0+1) #TODO: t0 in RR misses symbolic variables # even if something like assume(x, 'real') is used if p0==1 and t0 in RR: return -(-1)^floor((n0+1)/2)/factorial(n0+1) \ * 2^n0 * pi^(n0+1) * bernoulli_polynomial(t0/(2*pi),n0+1) if t0==pi/3 and n0==1: return (-1)^(p0+floor((p0+2)/2)) * (pi/3)^(p0+1)/(2*factorial(p0+1)) if t0 == pi/3 and (int(n0),int(p0)) in mgvreductions: return mgvreductions[(int(n0),int(p0))] return gl_nielsen(n0,p0,t0) result = result.substitute_function(gl_nielsen,gl_reduce) return result.expand() mcvreductions = {(1, 3): -1/18*pi^2*cl_nielsen(1, 1, 1/3*pi) - 1/9*pi*zeta(3) + cl_nielsen(3, 1, 1/3*pi), (3, 3): -11/324*pi^3*zeta(3) - 1/18*pi^2*cl_nielsen(3, 1, 1/3*pi) - 29/108*pi*zeta(5) + 3*cl_nielsen(5, 1, 1/3*pi), (1, 5): 1/1944*pi^4*cl_nielsen(1, 1, 1/3*pi) + 1/486*pi^3*zeta(3) - 1/18*pi^2*cl_nielsen(3, 1, 1/3*pi) - 25/162*pi*zeta(5) + cl_nielsen(5, 1, 1/3*pi), (2, 3): 1/108*pi^2*zeta(3) - 1/3*pi*cl_nielsen(3, 1, 1/3*pi) + 29/36*zeta(5), (4, 3): -17/4860*pi^4*zeta(3) - 25/972*pi^2*zeta(5) - 2/3*pi*cl_nielsen(5, 1, 1/3*pi) + 3229/1296*zeta(7), (1, 7): -1/524880*pi^6*cl_nielsen(1, 1, 1/3*pi) - 1/87480*pi^5*zeta(3) + 1/1944*pi^4*cl_nielsen(3, 1, 1/3*pi) + 25/8748*pi^3*zeta(5) - 1/18*pi^2*cl_nielsen(5, 1, 1/3*pi) - 637/3888*pi*zeta(7) + cl_nielsen(7, 1, 1/3*pi), (2, 5): -1/11664*pi^4*zeta(3) + 1/162*pi^3*cl_nielsen(3, 1, 1/3*pi) + 25/648*pi^2*zeta(5) - 2/3*pi*cl_nielsen(5, 1, 1/3*pi) + 3295/2592*zeta(7), (3, 5): 61/19440*pi^5*zeta(3) + 1/1944*pi^4*cl_nielsen(3, 1, 1/3*pi) + 191/5832*pi^3*zeta(5) - 1/9*pi^2*cl_nielsen(5, 1, 1/3*pi) - 3163/7776*pi*zeta(7) + cl_nielsen(5, 3, 1/3*pi)} mgvreductions = {(6, 4): 174605509/4714094246400*pi^10 + 53/12960*pi^4*zeta(3)^2 + 1/36*pi^2*zeta(3)*zeta(5) - 1/18*pi^2*gl_nielsen(6, 2, 1/3*pi) - pi*gl_nielsen(7, 2, 1/3*pi) - 5/2*zeta(3)*zeta(7) - 5/4*zeta(5)^2 + 7/2*gl_nielsen(8, 2, 1/3*pi), (5, 4): -645259/4081466880*pi^9 + 7/216*pi^3*zeta(3)^2 - 1/18*pi^2*gl_nielsen(5, 2, 1/3*pi) + 5/6*pi*zeta(3)*zeta(5) + 5/6*pi*gl_nielsen(6, 2, 1/3*pi) + 5*gl_nielsen(7, 2, 1/3*pi), (2, 6): 5491/37791360*pi^8 + 1/162*pi^3*gl_nielsen(3, 2, 1/3*pi) - 1/108*pi^2*zeta(3)^2 - 1/3*pi*gl_nielsen(5, 2, 1/3*pi) - zeta(3)*zeta(5) - gl_nielsen(6, 2, 1/3*pi), (4, 6): 3291377/30611001600*pi^10 - 79/116640*pi^4*zeta(3)^2 + 1/81*pi^3*gl_nielsen(5, 2, 1/3*pi) - 1/36*pi^2*zeta(3)*zeta(5) - 5/3*pi*gl_nielsen(7, 2, 1/3*pi) - 11/2*zeta(3)*zeta(7) - 11/4*zeta(5)^2 - 7/2*gl_nielsen(8, 2, 1/3*pi), (2, 8): 43809937/2357047123200*pi^10 - 1/29160*pi^5*gl_nielsen(3, 2, 1/3*pi) + 1/11664*pi^4*zeta(3)^2 + 1/162*pi^3*gl_nielsen(5, 2, 1/3*pi) + 1/18*pi^2*gl_nielsen(6, 2, 1/3*pi) - 1/3*pi*gl_nielsen(7, 2, 1/3*pi) - zeta(3)*zeta(7) - 1/2*zeta(5)^2 - gl_nielsen(8, 2, 1/3*pi), (4, 4): 221947/881798400*pi^8 + 1/108*pi^2*zeta(3)^2 - 2/3*pi*gl_nielsen(5, 2, 1/3*pi) - 2*zeta(3)*zeta(5), (3, 6): -2844967/142851340800*pi^9 + 1/1944*pi^4*gl_nielsen(3, 2, 1/3*pi) - 1/648*pi^3*zeta(3)^2 - 1/9*pi^2*gl_nielsen(5, 2, 1/3*pi) + 1/6*pi*zeta(3)*zeta(5) - 5/6*pi*gl_nielsen(6, 2, 1/3*pi) + 3*gl_nielsen(7, 2, 1/3*pi), (2, 2): -23/19440*pi^4, (4, 2): 209/918540*pi^6 - 1/6*zeta(3)^2, (3, 4): -5533/22044960*pi^7 - 1/18*pi^2*gl_nielsen(3, 2, 1/3*pi) + 1/6*pi*zeta(3)^2 + 2*gl_nielsen(5, 2, 1/3*pi), (2, 4): 4069/7348320*pi^6 - 1/3*pi*gl_nielsen(3, 2, 1/3*pi) - 1/3*zeta(3)^2}