The modified generalized Bessel function (MGBF) is defined in Ref. 1. For indices 1 and 2 and order zero, (which I will write as $I(u,v;t)$) it is expressed as the following series. $$I_0^{12}(u, v; t) = I(u, v; t) = \sum_{k=-\infty}^{\infty}t^kI_{2k}(u)I_k(v),$$ where $I_\nu(x)$ is the 1D modified Bessel function of the first kind of order $\nu$. The function also has the integral representation,
$$I(u, v; e^{i\phi}) = \frac{1}{2\pi}\int_{-\pi}^\pi \exp\left[u\cdot\cos(\theta) + v\cdot\cos(2\theta+\phi)\right]\mathrm{d}\theta$$
I need to numerically evaluate the function for u and v in [0, 10000] or so and for $t = e^{i\phi}$ with $\phi$ in [0, $\pi$].
Due to the exponential-like behavior, this means pulling out some of the asymptotics. That is, in most libraries, there are functions for the 1D modified Bessel function that return $e^{-|x|}I_\nu(x)$ to avoid the problems of floating point overflow when multiplying by a prefactor with $e^{-x}$ behavior. I feel like there should be a way to do the same for the MGBF.
When $\phi = 0$ the prefactor can be seen from the series to basically the same as the 1D case and it $e^{u + v}$. This is also shown in Ref. 1. However, for other parameters and being most extreme when $\phi=\pi$ the function grows at a pretty different rate when both $u$ and $v$ are large. This means that the same overflow issues exist.
Here is a figure of the problem. This is the MGBF multiplied by the prefactor $e^{-u-v}$ evaluated using the integral form. The function is well behaved when $\phi=0$, but otherwise the exponent goes crazy and will cause floating point errors if you evaluate it at larger $u$ and $v$.
My question is, what is the correct asymptotic behavior when $\phi\neq0$ to normalize the function and avoid floating point errors.
[1] Dattoli, G., Torre, A., Lorenzutta, S., Maino, G., & Chiccoli, C. (1991). Theory of generalized bessel functions.-II. Il Nuovo Cimento B, 106(1), 21–51. https://doi.org/10.1007/BF02723125