# Inverse of difference of two digamma functions

I recently encountered the expression below for which I was interested in solving for $$x$$:

$$$$\psi(x+n+1) - \psi(x+1) =y$$$$

$$\psi$$ is the digamma function, $$n$$ is a positive integer and $$x,y>0$$. Initially, I solved this numerically but I was interested to see if there was the possibility to solve this, either in closed form or as a series. I tried to find some way by using the relationship to harmonic numbers, e.g. rewriting as

$$$$\psi(x+n+1) - \psi(x+1) = \sum_{k=1}^{n} \frac{1}{k+x} = \int_0^1 \frac{t^x(1-t^n)}{1-t}dt$$$$

to see if there was something which I could do, but I didn't notice anything apparent. The only thing I've come up with is that for large values of $$x$$, $$y$$ is small and using $$\psi(x) \sim\ln(x)$$, that $$x \approx \dfrac{n}{e^y - 1}$$.

I'm happy for any suggestions on how to proceed!

• Cases, up to $n=4$, have a closed form inverse, so inverting $\frac1{x+1}+\dots+\frac1{x+n}$ will be like inverting a polynomial. Jan 18 at 13:42
• Yes, you are right @Gary! I changed notation when typing the question. Now it should be fixed. Jan 18 at 14:03
• What is $t$ at the end?
– Gary
Jan 18 at 14:03
• Same goes for that part, I should probably try to stick to my initial notation and not change things :) Jan 18 at 14:04
• If you use the better approximation $\psi(z+1) \approx \log(z+1/2)$ you get $$x\approx \frac{n}{\mathrm{e}^y-1}-\frac{1}{2}.$$
– Gary
Jan 18 at 14:10

$$\def\F{\operatorname F}\def\H{\operatorname H}$$Derivation: After different Lagrange reversion setups, this one had the simplest derivation without needing to factor. Additionally, less series are possible with Tschirnhausen transformations, but applying them to a degree $$m$$ polynomial in general is complicated. The harmonic number H$$_m$$ set up includes $$(0,0)$$: $$e^{\H_m-\sum\limits_{k=1}^m\frac1{y+k}}-1=x\implies y=x+\sum_{n=1}^\infty \frac1{n!}\frac{d^{n-1}}{dx^{n-1}}\left(x+1-e^{H_m}\prod_{k=1}^me^{-\frac1{x+k}}\right)^n$$ shown here, so here is a solution directly without reducing the equation. Using the binomial theorem on the derivative: $$\frac{d^{n-1}}{dx^{n-1}}\left(x+1-e^{\H_m}\prod_{k=1}^me^{-\frac1{x+k}}\right)^n=\sum_{j=0}^n\binom nj\left(-e^{\H_m}\right)^j\frac{d^{n-1}}{dx^{n-1}}\left((x+1)^{n-j} \prod_{k=1}^me^{-\frac j{x+k}}\right)$$ Applying the general Leibniz rule and factorial power $$a^{(b)}$$: $$\frac{d^{n-1}}{dx^{n-1}}\left((x+1)^{n-j} \prod_{k=1}^me^{-\frac j{x+k}}\right)=\sum_{v=0}^{n-1} \binom{n-1}v\frac{d^{n-1-v}}{dx^{n-1-v}}((x+1)^{n-j}) \frac{d^v}{dx^v}\left(\prod_{k=1}^me^{-\frac j{x+k}}\right)= \sum_{n_0=0}^{n-1} \binom{n-1}{n_0} (x+1)^{n_0-j+1}(n-j)^{(n-1-n_0)}\frac{d^{n_0}}{dx^{n_0}}\left(\prod_{k=1}^me^{-\frac j{x+k}}\right)$$ Also, we use the generalized product rule: $$\frac{d^{n_0}}{dx^{n_0}}\left(\prod_{k=1}^me^{-\frac j{x+k}}\right)= \sum_{n_1=0}^{n_0}\sum_{n_2=0}^{n_1}\dots\sum_{n_{m-1}=0}^{n_{m-2}}\frac{{n_0}!f_1^{(n_0-n_1)}f_m^{(n_{m-1})}}{(n_0-n_1)!n_{m-1}!}\prod_{k=2}^{m-1}\frac{f_k^{(n_{k-1}-n_k)}}{(n_{k-1}-n_k)!}$$ A formula with the confluent hypergeometric $$_1\F_1(a;b;z)$$ function and Kronecker $$\delta_x,\delta_{x,y}$$: $$\frac{d^r}{dx^r}e^{-\frac j{x+k}}=\sum_{u=0}^\infty\frac{(-j)^u}{u!}\frac{d^r}{dx^r}(x+k)^{-u}=\sum_{u=0}^\infty \frac{(-j)^u (-u)^{(r)}(x+k)^{-r-u}}{u!}=(x+k)^{-r}0^{(r)}-j(x+k)^{-r-1}(-1)^{(r)}\,_1\text F_1\left(r+1;2;-\frac j{x+k}\right)=\delta_r-j(-1)^r r!(x+k)^{-r-1}\,_1\text F_1\left(r+1;2;-\frac j{x+k}\right)$$ Simple Verification: Unfortunately, Wolfram Alpha cannot evaluate the full sums, but does for a partly written sum. Let’s verify that the above steps work for $$m=2$$. Testing $$n=3,x=0.3$$ gives matching series coefficients. Thereafter, testing $$\frac{d^v}{dx^v}\left(e^{-\frac j{x+1}-\frac j{x+2}}\right),j=4,v=3,x=0.3$$ gives the same result. Similarly when $$m=3$$, we get matching series coefficients. Testing $$\frac{d^v}{dx^v}\left(e^{-\frac j{x+1}-\frac j{x+2}-\frac j{x+3}}\right),v=4,j=5,x=0.3$$ gives the same result. Wolfram Alpha cannot evaluate the derivatives as confluent hypergeometric functions, but the derived $$\frac{d^r}{dx^r} e^{-\frac j{x+k}}$$ formula is true presents no convergence issues. Other values of $$v,j,x$$, including $$v,j=0$$, work when $$m=2,3$$ and combining verification links give evidence for the Lagrange reversion result working when $$m>3$$ Formula: Therefore we should get: $$\sum_{k=1}^m\frac1{y+k}=x\implies y=\H_m-\ln\left(1+\sum_{n=1}^\infty\sum_{j=0}^n\binom nj\left(-e^{\H_m}\right)^j \sum_{n_0=0}^{n-1}(x+1)^{n_0-j+1}(n-1)^{(n_0)}(n-j)^{(n-1-n_0)} \sum_{n_1=0}^{n_0}\sum_{n_2=0}^{n_1}\dots\sum_{n_{m-1}=0}^{n_{m-2}}\left(\frac{\delta_{n_0,n_1}} {(n_0-n_1)!}-j(-1)^{n_0-n_1}(x+1)^{n_1-n_0-1}\,_1\text F_1\left(n_0-n_1+1;2;-\frac j{x+1}\right)\right) \left(\frac{\delta_{n_{m-1}}} {n_{m-1}!}-j(-1)^{n_{m-1}}(x+m)^{-n_{m-1}-1}\,_1\text F_1\left(n_{m-1}+1;2;-\frac j{x+m}\right)\right)\prod_{k=2}^{m-1}\left(\frac{\delta_{n_{k-1},n_k}} {(n_{k-1}-n_k)!}-j(-1)^{n_{k-1}-n_k}(x+k)^{n_k-n_{k-1}-1}\,_1\text F_1\left(n_{k-1}-n_k +1;2;-\frac j{x+k}\right)\right) \right)= \boxed{\H_m-\ln\left(1+\sum_{n=1}^\infty\sum_{j=0}^n\sum_{n_0=0}^{n-1}\frac{\Gamma(n)n! (x+1)^{n_0-j+1} \left(-e^{\H_m}\right)^j}{(v-j+1)!\Gamma(n-v)j!}\sum_{0\le n_{m-1}\le\dots\le n_0}\left(\frac{\delta_{n_{m-1}}} {n_{m-1}!}-j(-1)^{n_{m-1}}(x+m)^{-n_{m-1}-1}\,_1\text F_1\left(n_{m-1}+1;2;-\frac j{x+m}\right)\right)\prod_{k=1}^{m-1}\left(\frac{\delta_{n_{k-1},n_k}} {(n_{k-1}-n_k)!}-j(-1)^{n_{k-1}-n_k}(x+k)^{n_k-n_{k-1}-1}\,_1\text F_1\left(n_{k-1}-n_k +1;2;-\frac j{x+k}\right)\right) \right)}$$ Closed form?: Since solving $$\sum\limits_{k=1}^m\frac1{y+k}=x$$ is solving a polynomial, you use the Riemann theta function, but this is too advanced for some.

• If one wants to reduce the polynomial equation, then it will give a series expansion with less series than above. However, it is tedious and new to me for a general $m$ degree polynomial from solving $\sum\limits_{n=1}^m\frac1{y+m}=x$. There are $m+2$ sums in this solution and with these transformations, the Lagrange reversion solution could give $m-4$ sums. Feb 25 at 3:24
• Thank you @TymaGaidash! I’ve read through your calculations a couple of times and seem to understund the steps. The answer is quite a brute though! I’ll try to repeat your calculations so I understand it in detail. Regarding the theta function, I’ve only encountered it once briefly, but maybe that’s something to look at one day :) Feb 26 at 5:39

$$\psi(x+n+1) - \psi(x+1) =y$$

For large values of $$x$$, the lhs write $$\frac{n}{x}-\frac{n (n+1)}{2 x^2}+\frac{n (n+1) (2 n+1)}{6 x^3}-\frac{n^2 (n+1)^2}{4 x^4}+O\left(\frac{1}{x^5}\right)$$ Using series reversion $$x_{(1)}=\frac{n}{y}-\frac{n+1}{2}+\frac{\left(n^2-1\right)}{12 n}y+O\left(y^3\right)$$

Comparing with your approximation $$x_{(2)}= \dfrac{n}{e^y - 1}=\frac{n}{y}-\frac{n}{2}+\frac{n }{12}y+O\left(y^3\right)$$ It is really good since $$x_{(2)}-x_{(1)}=-\frac{1}{2}-\frac{y}{12 n}+O\left(y^3\right)$$

Using the above for $$n=1000$$ and $$y=1$$, the exact solution is $$581.477$$ while $$x_{(2)}=582.833$$ and $$\frac{1000}{e-1}=581.977$$.

Edit

We could do better using more terms in the series expansion, then series reversion and transformation of the result into a $$[n,n+1]$$ Padé approximant.

This would give $$x_{(3)}=\frac{60 (n-1) n^3 -24 n^2 \left(n^2-4\right)y+3 n \left(n^3-3 n^2-9 n-13\right)y^2} {60 (n-1) n^2 y+6 n \left(n^2+11\right)y^2+(n-1) \left(n^2+11\right)y^3 }$$

For the worked example, this gives $$x_{(3)}=581.590$$.

• Thank you! I changed my estimate using Garys suggestion in order to improve it, but I have to admit that the solution in your edit works very well for a large range of values of $y$. Jan 19 at 8:08
• @AxelT. Glad to know ! In fact we could extend the Padé approximant to any higher order. Cheers :-) Jan 19 at 8:32
• I've seen that you've put the Padé approximant to good use previously with good success, so I guess I'd better learn it. I reccon that I would be able to use the same method for small values of $x$ as well. Jan 19 at 8:42