In this paper, we derive the joint (amplitude, phase) distribution of the product of two independent non-zero-mean Complex Gaussian random variables. We call this new distribution the complex Double Gaussian distribution. This probability distribution function (PDF) is a doubly infinite summation over modified Bessel functions of the first and second kind. We analyze the behavior of this sum and show that the number of terms needed for accuracy is dependent upon the Rician k-factors of the two input variables. We derive an upper bound on the truncation error and use this to present an adaptive computational approach that selects the minimum number of terms required for accuracy. We also present the PDF for the special case where either one or both of the input complex Gaussian random variables is zero-mean. We demonstrate the relevance of our results by deriving the optimal Neyman-Pearson detector for a time reversal detection scheme and computing the receiver operating characteristics through Monte Carlo simulations, and by computing the symbol error probability (SEP) for a single-channel M-ary phase-shift-keying (M-PSK) communication system.