 Solution

Introduce $t=\mu y/(2\lambda )\ge 0$. Solving the quadratic and using its product of roots givesBoth roots are positive, with $x_{-}\le \mu \le x_{+}$; equality occurs only for $y=0$. For computation, the smaller root has the equivalent cancellation-free formThis is stable root evaluation for inverse Gaussian sampling.

Define $h(x)=\lambda (x-\mu {)}^2/({\mu }^{2}x)$. It decreases from infinity to zero on $(0,\mu )$ and increases from zero to infinity on $(\mu ,\infty )$, sinceThe two roots are precisely its inverse branches. This supplies the Jacobians needed for the reciprocal-root inverse Gaussian sampler.