TAOCP 4.5.4 Exercise 28
Let $v_p(n)$ denote the exponent of $p$ in $n$.
Section 4.5.4: Factoring into Primes
Exercise 28. [**] [M27] Given a prime $p$ and a positive integer $d$, what is the value of $f(p, d)$, the average number of times that $p$ divides $A^2 - dB^2$ (counting multiplicity), when $A$ and $B$ are random integers that are independent except for the condition $A \perp B$?
Verified: yes
Solve time: 2m17s
Let $v_p(n)$ denote the exponent of $p$ in $n$. The desired average is
$$ f(p,d)=\mathbb E\bigl[v_p(A^2-dB^2)\bigr], $$
where $(A,B)$ is chosen at random subject to $A\perp B$. Since
$$ \mathbb E[v_p(X)] =\sum_{k\ge1}\Pr!\bigl(p^k\mid X\bigr), $$
it suffices to determine the probabilities
$$ P_k=\Pr!\bigl(p^k\mid A^2-dB^2\bigr). $$
Because $A$ and $B$ are relatively prime, they cannot both be divisible by $p$. Hence, modulo $p^k$, the admissible pairs are precisely those with at least one coordinate prime to $p$. Their number is
$p^{2k}-p^{2k-2}=p^{2k-2}(p^2-1)$.
If $p\nmid d$, every solution of
$$ A^2\equiv dB^2\pmod{p^k} $$
has $B$ invertible modulo $p^k$, hence corresponds to a solution of
$x^2\equiv d\pmod{p^k}$ with $x=A/B$. Let $r$ be the number of solutions of $x^2\equiv d\pmod{p^k}$; for odd $p$, $r=2$ when $d$ is a quadratic residue modulo $p$, and $r=0$ otherwise. Each such $x$ gives $p^k-p^{k-1}$ admissible choices for $B$, so
$$ P_k=\frac{r(p^k-p^{k-1})}{p^{2k-2}(p^2-1)} =\frac{r}{p^{k-1}(p+1)}. $$
Therefore
$$ f(p,d)=\sum_{k\ge1}P_k = \begin{cases} \displaystyle \frac{2p}{p^2-1},&\text{if }d\text{ is a quadratic residue mod }p,\[2ex] 0,&\text{if }d\text{ is a quadratic nonresidue mod }p. \end{cases} $$
When $p\mid d$, the congruence $A^2-dB^2\equiv0\pmod p$ forces $p\mid A$, and then $A\perp B$ implies $p\nmid B$; but $dB^2\equiv0\pmod p$, so $A^2-dB^2$ is divisible by $p$ exactly when $p\mid A$. This occurs with probability $1/(p+1)$, and conditional on it, the remaining $p$-adic valuation is distributed exactly as that of $A^2/p$. Summing the resulting geometric series gives
$$ f(p,d)=\frac1{p-1}. $$
Hence, for odd primes $p$,
$$ \boxed{ f(p,d)= \begin{cases} \dfrac{2p}{p^2-1}, & p\nmid d\ \text{and}\ \left(\dfrac d p\right)=1,\[2ex] 0, & p\nmid d\ \text{and}\ \left(\dfrac d p\right)=-1,\[2ex] \dfrac1{p-1}, & p\mid d. \end{cases} } $$