puzzle contents
raw puzzle

Original Problem

Find the number of positive integral solutions for the equations

\[\frac{1}{x} + \frac{1}{y} = \frac{1}{N!}\]

Input Format 
An integer N 

Output Format 
The number of positive integral solutions for the above equation modulo 1000007


\(1\leq N\leq 10^6\)


For small integers the problem could be solved with a brute-force solution but since \(N\) can be quite large, \(N!\) would be much larger so we need to find another way. Lets substitute \(n := N!\) and rewrite the original equation such that

\[\begin{array}{rrl} & \frac{1}{x} + \frac{1}{y} &= \frac{1}{n}\\ \Leftrightarrow & \frac{x+y}{xy} &= \frac{1}{n}\\ \Leftrightarrow & xy - xn - yn &= 0\\ \Leftrightarrow & n^2 + xy - xn - yn &= n^2\\ \Leftrightarrow & (x-n)(y-n) &= n^2\\ \end{array}\]

When we now look at the prime factorization of the positive integer \(n^2\) with \(k\) factors \(p_i\) and their exponent \(e_i\) we see

\[n^2 = \prod\limits_{i=1}^k \left(p_i^{e_i}\right)^2 = \prod\limits_{i=1}^k p_i^{2e_i}\]

It follows that the number of integer solutions for the stated equation \(\frac{1}{x} + \frac{1}{y} = \frac{1}{n}\) is given by the number of divisors of \(n^2\).

\[\delta(n^2) = \prod\limits_{i=1}^k (1+2e_i) \]

In order to come back to the original problem, we need to find the exponents of the prime factorization of \(N!\) instead of \(N\). Intuitively, since we multiply \(1\cdot 2\cdot\dots\cdot N\) for \(N!\), any factor has a prime factorization less than or equal to \(N\) and as such also the product of the individual factors must have all prime factors \(\leq N\). From this assertion we can also see that the unique prime factors of \(N!\) must be all primes below \(N\). We can use this and divide each prime factor up to \(N\) to count the exponent \(e_i\). When we use Ruby with its built in eratosthenes prime generator, the code is pretty clear:

require 'prime'

M = 1000007
N = gets.to_i

puts Prime::EratosthenesGenerator
    .take_while {|i| i <= N}
    .inject(1) { |prod, p|

  x = N
  e = 0
  loop do
    x/= p
    break if x == 0
    e+= x

  prod * (2 * e + 1)
} % M