How to prove that : $$ \int x^{n-1}W(x)dx = \frac {x^ne^{[-nW(x)]}[-nW(-x)]^{-n}[n\Gamma(n+1, -nW(x)- \Gamma(n+2, -nW(x))]} {n^2} $$ Where $W(x)$ is the Lambert-W function https://en.wikipedia.org/wiki/Lambert_W_function and $\Gamma$ is the incomplete gamma function.
I found this on wolfram alpha https://www.wolframalpha.com/input/?i=int+(x)%5E(n-1)W(x), but have no idea how it got this.
For those who are wondering why I want to know the proof of this integral it is because with this integral I can find the integral of the infinite titration of x as :
$ \int x^{x^{x^{.^{.......}}}} dx = \sum_{n=0}^{\infty}-\frac {(-1)^n( \ln x)^ne^{[-nW(-\ln x)]}[-nW(-\ln x)]^{-n}[n\Gamma(n+1, -nW(-\ln x)- \Gamma(n+2, -nW(-\ln x))]} {(n!)n^2} $
Lol the answer is so long that it doesnt fit on one line.