As we know, the joint distribution of standard Brownian motion $B_t$ and its running maximum $M_t=\sup_{s\in[0,t]} B_s$ is given by $$P(B_t\in da,M_t\in db)=\frac{2(2b-a)}{\sqrt{2\pi t^3}}e^{-\frac{(2b-a)^2}{2t}}da\,db$$ What is the joint distribution of $B_t$ and $M_t(M_t-B_t)$, i.e., how to compute $P(B_t\in da,M_t(M_t-B_t)\in dc)$?
I have tried some formal way like solving $M_t$ by letting $B_t=da$ in $M_t(M_t-B_t)=dc$, then I get $M_t=(da+\sqrt{(da)^2+4dc})/2$. I let $db=(da+\sqrt{(da)^2+4dc})/2$ in the given formula, and I cannot proceed.