Number of solutions for $\frac{1}{X} + \frac{1}{Y} = \frac{1}{N!}$ where $1 \leq N \leq 10^6$

Note: this is a programming challenge at this site

For this equation $$\frac{1}{X} + \frac{1}{Y} = \frac{1}{N!}\quad ( N \text{ factorial} ),$$ find the number of positive integral solutions for $(X,Y)$ ?

Note that : $1 \leq N \leq 10^6$.


Equivalently, we want to solve $N!(X+Y)=XY$, which can be rewritten as $$(X-N!)(Y-N!)=(N!)^2.$$

Note that $X\gt N!$ and $Y\gt N!$. For if, for example, $X\lt N!$, then since $X$ and $Y$ are positive, we have $\dfrac{1}{X}+\dfrac{1}{Y} \ge \dfrac{1}{N!}+\dfrac{1}{Y}\gt \dfrac{1}{N!}$.

Now let $s=X-N!$ and $t=Y-N!$. To count the solutions of our original equation, we count the number of solutions of $st=(N!)^2$.

To do this, note that $s$ ranges over the divisors of $(N!)^2$, and that once we know $s$, we know $t$, and therefore $X$ and $Y$.

Thus the number of solutions of our equation is the number $d((N!)^2)$ of (positive) divisors of $(N!)^2$.

That still leaves a great deal of work to do. There is a simple formula for the number of divisors of $n$, once we know the structure of the prime power factorization of $n$. However, there is no simple expression for the structure of the prime power factorization of $(N!)^2$.

However, the problem is not computationally hopeless. For any prime $p$, the highest power of $p$ that divides $N!$ is $$\lfloor N/p\rfloor+\lfloor N/p^2\rfloor+\lfloor N/p^3\rfloor+ \lfloor N/p^4\rfloor+\cdots,$$ and for $(N!)^2$ we just double. The sum above is a finite sum, indeed a rather short sum, since the terms $p^i$ in the denominator grow exponentially, and soon surpass $N$.

Remark: For completeness, so that you can begin to build a program, suppose that $n$ has prime power factorization $$p_1^{a_1}p_2^{a_2}\cdots p_k^{a_k},$$ where the $p_i$ are distinct primes. Then the number of positive divisors of $n$ is $$(a_1+1)(a_2+1)\cdots(a_k+1).$$