The probability a given integer in $[0,n]$ is a square is $\frac1{\sqrt n}$. What is the probability that if you take two integers uniformly then their product is square? 


I know the main term is $\frac1n$. I am also looking for correction terms. The difficulty is an average integer has $\omega(\log\log n)$ factors which can go as high as $\frac{\log n}{\log\log n}$ and as low as $2$. Then we are looking at the probability that the prime factors with odd parity in both integers come in pairs. Is there a way to organize the computations and get one precise quantity?



> 
> In general what is the number of $m$-tuples whose product is a $k$ power where each integer $x\_i$ is in $[a\_i,b\_i]$ where $i\in\{1,\dots,m\}$?
> 
> 
> 



Probability that product is a perfect square
For the case $m=k=2$, in which we seek the number $N(n)$ of pairs
$(x,y) \in [1,n]^2$ for which $xy$ is a square, we give
an elementary estimate
$$
N(n) = Cn \log n + An + O(n^{2/3}),
$$
where $C = 1/\zeta(2) = 6/\pi^2$ and
$$
A = \frac{3\gamma-1}{\zeta(2)}
 - \frac{2\zeta'(2)}{\zeta(2)^2}
 - 1
= 0.1377775\ldots \, .
$$
This agrees with the analytic calculation of **Asymptotiac K**
(and is corroborated by numerical computation up to $n = 2^{30}$),
and improves the error term (the contour-integral analysis gave
an error estimate equivalent to $n^{21/22+o(1)}$ rather than $n^{2/3}$).


Recall that $xy$ is a square **iff** $(x,y) = (ma^2,mb^2)$
for some positive integers $m,a,b$, and the representation
can be made unique by requiring either that $a,b$ be coprime
or that $m$ be squarefree, which is the source of the factor $1/\zeta(2)$.
Let $M(n)$, then, be the number of triples $(m,a,b)$ of positive integers
such that $ma^2\leq n$ and $mb^2\leq n$, without the additional coprimality or squarefree condition.
Then Möbius inversion (applied with either condition $\gcd(a,b)=1$
or $\mu(m)^2 = 1$ yields
$$
N(n) = \sum\_{d=1}^{\lfloor\sqrt{n}\rfloor} \mu(d) M(\lfloor n/d^2 \rfloor).
$$
We show:


**Proposition.** *$M(n) = n \log n + Bn + O(n^{2/3})$, where
$B = 3\gamma - 1 - \zeta(2) = -0.913287\ldots$* (and again
$\gamma$ is Euler's constant $0.5772156649\ldots$).


*Proof:* Let $R = \lfloor n^{1/3} \rfloor$. If $m\cdot\max(a,b)^2 \leq n$
then either $m \leq R$ or $\max(a,b) \leq R$.
Given $m$, the number of $(a,b)$ pairs is
$\lfloor \sqrt{n/m} \rfloor^2 = n/m - O(\sqrt{n/m})$;
summing this over $m \leq R$ gives $n H\_R - O(n^{2/3})$
where $H\_R$ is the harmonic sum $\sum\_{m=1}^r 1/m$.
In the other direction, each $k \leq R$ occurs $2k-1$ times as
$\max(a,b)$ for positive integers $a,b$, each of which accounts for
$\lfloor n/k^2 \rfloor$ solutions, for a sum of
$$
\sum\_{k=1}^R (2k-1) \, (n/k^2 - O(1)) = (2 H\_R - \zeta(2)) n - O(n^{2/3})
$$
solutions. For the total count we add this to $n H\_r - O(n^{2/3})$,
and subtract $R^3 = n - O(n^{2/3})$ which is the number of solutions
for which $a,b,m$ are all $\leq R$, finding
$(3 H\_R - \zeta(2) - 1) n + O(n^{2/3})$. The Proposition then follows from
$H\_R = \log R + \gamma + O(1/R) = \frac13 \log n + \gamma + O(n^{-1/3})$. $\Box$


The estimate for $N(n)$ then follows from
$N(n) = \sum\_{d=1}^{\lfloor\sqrt{n}\rfloor} \mu(d) M(\lfloor n/d^2 \rfloor)$;
the $2\zeta'(2) / \zeta(2)^2$ comes from the second term of
$$
\frac{n}{d^2} \log \frac{n}{d^2} = \frac1{d^2} n \log n - 2n \frac{\log d}{d^2}
$$
because $\sum\_{d=1}^\infty \mu(d) \log d \, / \, d^2$
is the derivative at $s=2$ of $-1/\zeta(s)$.


One can probably use bounds on exponential sums to further reduce the
$O(n^{2/3})$ error, both in worst and average case, as is done for the
[Dirichlet divisor problem](https://en.wikipedia.org/wiki/Divisor_summatory_function#Dirichlet.27s_divisor_problem)
and the [Gauss
circle problem](https://en.wikipedia.org/wiki/Gauss_circle_problem).


