What is the Cumulative Distribution Function of the following random variable? Suppose that we have $2n$ iid random variables $X_1,…,X_n,Y_1,…,Y_n$ where $n$ is a large number. 
I want to find $P((k∑_iX_iY_i+(∑_iX_i)(∑_jY_j))<c)$ for any integer c.
Since $n$ is a large number and all the random variables are $iid$, using central limit theorem, we can say that $k∑_iX_iY_i$, $(∑_iX_i)$ and $(∑_jY_j)$ are approximately normal random variables and $(∑_iX_i)$$(∑_jY_j)$ is the product of two normal random variables which would have Normal Product Distribution.
So $k∑_iX_iY_i+(∑_iX_i)(∑_jY_j)$ is the sum of one normal and one normal product random variable which are dependent.
Now the question is how can we find $P((k∑_iX_iY_i+(∑_iX_i)(∑_jY_j)) \le c)$ for any integer c? 
 A: @May: if they are independent (as in the first i of iid) then $\sum_i \sum_j X_i Y_j =\sum_i X_i \times \sum_j  Y_j$. 
For example if they are normally distributed then you have the product of two normal random variables, which is in general not normally distributed and this extends to most other distributions, especially if the sums are over a large number of cases, though if the means are non zero and the standard deviations small compared with the means, this may be difficult to spot.
A: $$Z = \sum_{i=1}^n \sum_{j=1}^n X_i Y_j = \left(\sum_{i=1}^n X_i\right)\left(\sum_{j=1}^n Y_j\right)$$
If $n$ is large, $S_X = \sum_i X_i$ and $S_Y = \sum_j Y_j$ are approximately normal. They have means $n\mu$ and standard deviations $\sqrt{n} \sigma$ where each $X_i$ and $Y_j$ have mean $\mu$ and standard deviation $\sigma$.  Of course they are independent.  Thus
$E[Z] = E[S_X] E[S_Y] = n^2 \mu^2$ and $E[Z^2] = E[S_X^2] E[S_Y^2] = (n^2 \mu^2 + n \sigma^2)^2$, so the variance of $Z$ is $\text{Var}(Z) = E[Z^2] - E[Z]^2 = n^2 \sigma^4 + 2 n^3 \sigma^2 \mu^2$.  
The moment generating function of the product of independent normal random variables with means $n\mu$ and standard deviations $n \sqrt{\sigma}$ has, according to Maple, moment generating function
$$ M_Z(t) = E[e^{tZ}] = \frac{1}{\sqrt{1 - n^2 \sigma^4 t^2}} \exp\left(\frac{n^2 \mu^2 t}{1 - n \sigma^2 t}\right)$$
for  $t < 1/(n \sigma^2)$.
EDIT: If $\mu \ne 0$, it would be better to separate out the effect of the mean.  So let $X_i = \mu + \sigma U_i$ and $Y_i = \mu + \sigma V_i$, where $U_i$ and $V_i$ have mean $0$ and standard deviation $1$.  Then $$Z  = n^2 \mu^2 + n \mu \sigma \sum_{i=1}^n (U_i + V_i) + \sigma^2 \sum_{i=1}^n \sum_{j=1}^n U_i V_j$$
Now $n \mu \sigma \sum_{i=1}^n (U_i + V_i)$ is approximately normal with mean $0$ and 
standard deviation $\sqrt{2} n^{3/2} \mu \sigma$, while $\sigma^2 \sum_{i=1}^n \sum_{j=1}^n U_i V_j$ has mean $0$ and standard deviation $n \sigma^2$.  For large
$n$ this term is negligible compared to the $n^{3/2}$ term.  So a good
approximation to the distribution of $Z$ is normal with mean $n^2 \mu^2$ and standard deviation $\sqrt{2} n^{3/2} \mu \sigma$.
You asked about $ (k−1) \sum_i X_i Y_i+ Z$: call this $(k-1) T + Z$.  If we separate out the effect of the mean, 
$$T = n \mu^2 +  \mu \sigma \sum_{i=1}^n (U_i + V_i) + \sigma^2\sum_{i=1}^n U_i V_i$$
where $\mu \sigma \sum_{i=1}^n (U_i + V_i)$ has mean $0$ and standard deviation $\sqrt{2n} \mu \sigma$ and $\sigma^2 \sum_{i=1}^n U_i V_i$ has mean $0$ and standard deviation $\sqrt{n} \sigma^2$.  Again, these terms are negligible compared to the $n^2$ and $n^{3/2}$ terms.  
