Given particle undergoing Geometric Brownian Motion, want to find formula for probability that max-min > z after n days

Consider a particle undergoing geometric brownian motion with drift $\mu$ and volatility $\sigma$ e.g. as in here. Let $W_t$ denote this geometric brownian motion with drift at time $t$. I am looking for a formula to calculate:
$$ \mathbb{P}\big(\max_{0 \leq t \leq n} W_t - \min_{0\leq t \leq n} W_t > z\big) $$ The inputs to the formula will be $\mu$, $\sigma$, $z$, and $n$.

The following horrible formula for the joint distribution of max, min and end value of a Brownian motion was copied without guarantees from the Handbook Of Brownian Motion (Borodin/Salminen), 1.15.8, p.271. First, for simplicity, this is only written for $\sigma=1,t=1$, and the more general case comes directly from scaling. If we shorten W as the Brownian Botion at t=1, m as the minimum and M as the maximum over $[0,1]$, then for $a < min(0,z) \le max(0,z) < b$ it holds $$ P(a < m, M < b, W \in dz) = \frac{1}{\sqrt{2\pi}}e^{(\mu z-\mu^2/2)} \cdot \sum_{k =-\infty}^{\infty} \Bigl(e^{-(z+2k(b-a))^2/2} - e^{(z-2a + 2k(b-a))^2/2} \Bigr) dz\; . $$ (Apologies for using z here in a different context.) If one really wants to, one can compute from this an even more horrible formula for the above probability. It is now in principle possible to derive from this a formula for what you want, by finding the density function $p_{m,M,W}$, and using $$ P(e^M-e^m\le r) = \int_{(x,y,z)\ :\ e^x \le e^z \le e^y \le e^x + r} p_{m,M,W}(x,y,z) d(x,y,z)\;, $$ but I shudder at the monster I expect to fall out from this. It might be better to give up and simulate the probability in question, and find some asymptotics.

However, if you would like to proceed with it, I suggest you look not into the Handbook Of Brownian Motion, but rather into this paper, as it is much more readable.