I wish to work with the canonical form of the polygamma function when it is being evaluated. Namely, $$\textrm{polygamma}(m,n)=(-1)^{m+1}m!\zeta(m+1,n)$$, with $\zeta(m+1,n)$ the Hurwitz Zeta function. This is in turn related to the Riemann Zeta function by $$\zeta(m,n+1) = \zeta(m) - \textrm{harmonic}(n,m)$$. Therefore, we should have the relationship $$\textrm{polygamma}(m,n)=(-1)^{m+1}\Gamma(m+1)*[\zeta(m+1)-\textrm{harmonic}(n-1,m+1)]$$.
In SymEngine, I would wish to see:
import symengine as se
m,n = se.symbols('m,n')
f = se.polygamma(m,n)
print(f.subs({m:3,n:15}))
>> 6*(-18249859383918836502097/16863445484161436160000 + (1/90)*pi**4)
print(f.subs({m:4,n:15}))
>>6301283671733562325696376933/253204633944683963942400000−24*se.zeta(5,1)
However, for all even values of m, we simply return the polygamma function itself:
print(f.subs({m:4,n:15}))
>>polygamma(4,15)
One way to work around this is to use SymPy capabilities:
print(f.subs({m:4,n:15})._sympy_().rewrite(sp.harmonic))
>>6301283671733562325696376933/253204633944683963942400000−24𝜁(5)
However, SymEngine natively has a function to rewrite this as a Zeta value (see rewrite_as_zeta()). I was looking at the SymEngine.py constructor and was wondering if it would be possible to add this either as an option, or a class method for PolyGamma. I haven't worked much with conversions between already-written C code and Python, and was wondering how feasible this would be.
Alternatively, I would be interested in using the .replace() or .xreplace() functionality, but I have been getting errors when I try to use a lambda function to replace every component of the expression appropriately. Moreover, I'm not sure that this would work, as it seems like although we call the C code when creating the polygamma function, we don't have access to the C code for the harmonic numbers found here.
Is there any straightforward way to get this form using SymEngine by itself?
Thank you!