Find probability that product of two random integers uniformly distributed from $1$ to $n$ is divisible by $n$. I suppose that $$P(n)=\frac{\prod_{i=1}^k (2p_i-1)}{n^2}$$ as $n=\prod_{i=1}^kp_i$ ($p_1$ to $p_k$ are different prime numbers). 
It works with 1,2 or 3 primes, but I am not able to prove it.
 A: Assume that the number of distinct prime divisors of $n$ is $k$, so that:
$$
n=\prod_{i=1}^k p^{r_i}_i,
$$
the overall number of distinct divisors being $\prod_{i=1}^k(r_i+1)$.
Choose some divisor $d$ of $n$ and let $M_d$ be the set of numbers in the range from $1$ to $n$, having $d$ as greatest common divisor with $n$: $m\in M_d: \{1\le m \le n,\ \gcd(m,n)=d\}$. The set will contain $\varphi(n/d)$ numbers, $\varphi(n)$ being the Euler's totient function. Repeating the procedure for all distinct divisors of $n$ one ends up with partition of the set $(1..n)$ into disjoint union of $M_{d_i}$. Observe that for every $m\in M_d$ there are $d$ ways to choose the second number (a multiple of $n/d$) such that the product of both is divisible by $n$. 
Thus, the overall number of "good" pairs is:
$$
N(n)=\sum_{d|n}\varphi\left(\frac{n}{d}\right)d.
$$
$N(n)$ being the Dirichlet convolution of two multiplicative functions $\phi(n)$ and $\text{id}(n)$ is itself a multiplicative function:
$$N\left(\prod_i p_i^{r_i}\right)=\prod_{i}N(p_i^{r_i}),$$
with
$$
N(p^{r})=p^r+\sum_{i=1}^{r} p^{i-1}(p-1)p^{r-i}
=p^r\left[1+\left(1-\frac{1}{p}\right)r\right].
$$
As the overall number of pairs is $n^2$ the probability in question is:
$$P(n)=\frac{N(n)}{n^2}=\frac{\prod_i\left[1+\left(1-\frac{1}{p_i}\right)r_i\right]}{n}.
$$
In your case all $r_i=1$, so that the above expression reduces to:
$$
P(n)=\frac{\prod_i\left(2-\frac{1}{p_i}\right)}{n}=\frac{\prod_i(2p_i-1)}{n^2},
$$
as claimed.
A: Assuming $n$ is square-free, it holds. The function $P$ is multiplicative by the Chinese remainder theorem, and the pairs $(n, m)$ in $1, \dots, p$ with $p| nm$ are exactly the $2p - 1$ ones with $p = n$ or $p = m$.
