Probability that product is a perfect square 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\}$?

 A: The main term certainly exceeds $1/n$ because $0$ appears $2n+1$ times . . .
I assume, then, that you meant to choose $(x,y)$ uniformly from $[1,n]^2$,
not $[0,n]^2$.  But it's still more than $1/n$: you get
$n$ solutions from just $x=y$, and another $n - O(n^{1/2})$ from 
$(x,y)=(a^2,b^2)$, and these have an overlap of only $\lfloor n^{1/2} \rfloor$
so we get a probability of at least $2/n - O(n^{-3/2})$.
In fact the main term is $(C \log n) / n$ where $C = 1/\zeta(2) = 6/\pi^2$.
Indeed $x,y$ such that $xy$ is a square are parametrized by
$(ma^2,mb^2)$ with $\gcd(a,b)=1$.  Given $m \leq n$, there are
$\lfloor(n/m)^{1/2}\rfloor$ choices for each of $a$ and $b$,
and they are relatively prime with probability about $1/\zeta(2) = C$
as long as $m = o(n)$.  Hence the count is asymptotic to
$C \sum_{m=1}^n n/m \sim C n \log n$, and the probability is
asymptotic to $(C\log n) / n$ as claimed.
A: Using the method of contour integration, one can get
$$P(N)=\frac{\log N}{\zeta(2)N}+\frac{A}{N}+O(N^{-23/22+o(1)}),$$
where $P(N)$ is our probability and $A=\frac{3\gamma}{\zeta(2)}-1-\frac{2\zeta'(2)}{\zeta(2)^2}-\frac{1}{\zeta(2)}$. To prove this formula, lets note that if
$$N^2P(N)-(N-1)^2P(N-1)=a(N)$$
then $a(N)-1$ is equal to the number of pairs of the form $(x,N)$ or $(N,x)$ with $xN=y^2$ for some $y$, $x<N$. One can easily see that if $D(N)^2$ is a maximal square divisor of $N$ then $a(N)-1=2D(N)-2$. This is because $N=MD(N)^2$ with squarefree $M$ and if $xN=y^2$, then we should have $x=Mz^2$ for some integer $z$. There are exactly $D(N)-1$ admissible values for $z$, hence there are $2(D(N)-1)$ pairs of the form $(x,N)$ or $(N,x)$ with $xN=y^2$ and $x<N$. So we get
$$a(N)=2D(N)-1.$$
Now let us consider the generating Dirichlet series of the function $a(N)$:
$$f(s)=\sum_{n\in \mathbb N} \frac{a(n)}{n^s}=2\sum_{n \in \mathbb N} D(n)/n^s-\zeta(s).$$
The function $D(n)$ is multiplicative and for any prime $p$ we have $D(p^k)=p^{[k/2]}$. Therefore,
$$\sum_{n \in \mathbb N} D(n)n^{-s}=\prod_p (1+p^{-s}+p^{1-2s}+p^{1-3s}+\ldots)=\prod_p(1+p^{-s})(1+p^{1-s}+p^{2-2s}+\ldots)=\frac{\zeta(s)\zeta(2s-1)}{\zeta(2s)}.$$
Next, using the truncated version of Perron's formula and the fact that $D(n) \leq \sqrt n$, we obtain for any $b>1$ (and fixed) and any $T$ the following:
$$\sum_{n \leq N} D(n)=\frac{1}{2\pi i}\int_{b-iT}^{b+iT} N^s\frac{\zeta(s)\zeta(2s-1)}{s\zeta(2s)}ds+O\left(\frac{N^{3/2}\log N}{T}\right).$$
Moving the contour of integration to the line $\mathrm{Re}\,s=1/2+\varepsilon$ for any $\varepsilon>0$ and using the estimates $\zeta(s) \ll |s|^{1/6}$, $\zeta(2s-1) \ll |s|^{1/2}$ and $\zeta(s)^{-1} \ll 1$ one can prove that
$$\sum_{n \leq N} D(n)=\mathrm{Res}_{s=1}N^s\frac{\zeta(s)\zeta(2s-1)}{s\zeta(2s)}+O(N^\varepsilon(\sqrt N T^{5/6}+NT^{-1/3}+N^{3/2}T^{-1})).$$
Choosing $T=N^{6/11}$, we finally obtain
$$\sum_{n \leq N} D(n)=\mathrm{Res}_{s=1}N^s\frac{\zeta(s)\zeta(2s-1)}{s\zeta(2s)}+O(N^{21/22+\varepsilon}).$$
It remains to compute the residue.
For $s \to 1$ we have
$$\frac{N^s}{s}=N+(N\log N-N)(s-1)+O((s-1)^2),$$
$$\zeta(s)=\frac{1}{s-1}+\gamma+O(s-1),$$
$$\zeta(2s-1)=\frac{1}{2(s-1)}+\gamma+O(s-1)$$
and
$$\frac{1}{\zeta(2s)}=\frac{1}{\zeta(2)}-\frac{2\zeta'(2)}{\zeta(2)^2}(s-1)+O((s-1)^2).$$
Multipying this gives the result
$$\frac{N}{2\zeta(2)^2(s-1)^2}+\frac{N\log N/2\zeta(2)+N(A+1)/2}{s-1}+O(1),$$
so
$$\sum_{n \leq N} D(n)=\frac{N\log N}{2\zeta(2)}+N(A+1)/2+O(N^{21/22+o(1)})$$
and
$$N^2P(N)=\sum_{n \leq N} a(N)=2\sum_{n \leq N} D(n)-N=\frac{N\log N}{\zeta (2)}+AN+O(N^{21/22+o(1)})$$
which is equivalent to the stated result.
Analogous estimates for $k$-th powers with $k>1$ will give us something like
$$\frac{2\zeta(k-1)}{\zeta(k)N}-\frac{1}{N}+O(N^{-1-e_k})$$
with $e_k>0$, as the corresponding generating function is $2\frac{\zeta(s)\zeta(ks-1)}{\zeta(ks)}-\zeta(s)$. Also, I think that the error term could be significantly improved under the assumption of the Riemann Hypothesis.
A: 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
and the Gauss
circle problem.
