Project Euler Problem 561
Let S(n) be the number of pairs (a,b) of distinct divisors of n such that a divides b.
Solution
Answer: 452480999988235494
Let
$$n=\prod_{j=1}^m p_j^{,r_j}.$$
Every divisor of $n$ corresponds to a vector
$$(e_1,\dots,e_m), \qquad 0\le e_j\le r_j.$$
For two divisors $a,b$, the condition $a\mid b$ is equivalent to
$$e_j(a)\le e_j(b)\quad \text{for all }j.$$
We are asked about
$$S(n)=#{(a,b): a\ne b,\ a\mid b,\ a,b\mid n}.$$
For Project Euler 561,
$$n=(p_m)^k=\prod_{j=1}^m p_j^k,$$
so all exponents are equal to $k$.
Step 1: Closed formula for $S((p_m)^k)$
For one prime coordinate, the number of pairs $(x,y)$ with
$$0\le x\le y\le k$$
is
$$\frac{(k+1)(k+2)}2.$$
Since coordinates are independent, the total number of divisor pairs with $a\mid b$ is
$$\left(\frac{(k+1)(k+2)}2\right)^m.$$
This includes the $a=b$ cases, of which there are $(k+1)^m$. Hence
$$\boxed{ S((p_m)^k)= \left(\frac{(k+1)(k+2)}2\right)^m-(k+1)^m }$$
or equivalently
$$\boxed{ S((p_m)^k) =(k+1)^m\left(\left(\frac{k+2}{2}\right)^m-1\right). }$$
Inline formula visualization:
$S\left((p_m)^k\right)=\left(\frac{(k+1)(k+2)}{2}\right)^m-(k+1)^m$
We need
$$E(m,k)=v_2!\left(S((p_m)^k)\right),$$
where $v_2(x)$ is the exponent of 2 dividing $x$.
Here
$$m=904961,$$
which is odd.
Step 2: Determine $E(m,k)$
We split into parity cases.
Case A: $k=2r+1$ (odd)
Then
$$k+1=2(r+1),$$
so
$$S=(r+1)^m\bigl((2r+3)^m-2^m\bigr).$$
Since $(2r+3)^m$ is odd and $2^m$ is even,
$$(2r+3)^m-2^m$$
is odd.
Therefore
$$E(m,k)=m,v_2(r+1).$$
Because
$$r+1=\frac{k+1}{2},$$
we get
$$\boxed{ E(m,k)=m\left(v_2(k+1)-1\right) \qquad (k\text{ odd}). }$$
Case B: $k=2r$ (even)
Then
$$S=(2r+1)^m\bigl((r+1)^m-1\bigr).$$
The first factor is odd, so
$$E(m,k)=v_2((r+1)^m-1).$$
Since $m$ is odd:
- if $r$ is odd, then $r+1$ is even, so $(r+1)^m-1$ is odd:
$$E=0.$$
- if $r$ is even, then $r+1$ is odd, and LTE gives
$$v_2((r+1)^m-1)=v_2(r)+v_2(m)=v_2(r),$$
because $m$ is odd.
Thus:
$$\boxed{ E(m,k)= \begin{cases} v_2(k)-1, & 4\mid k,\[4pt] 0, & k\equiv2\pmod4. \end{cases} }$$
Step 3: Sum $Q(n)$
We need
$$Q(N)=\sum_{k=1}^N E(904961,k), \qquad N=10^{12}.$$
Separate odd and even contributions.
Odd $k$
Write
$$k=2j-1.$$
Then
$$v_2(k+1)-1=v_2(2j)-1=v_2(j).$$
So the odd contribution is
$$904961\sum_{j\le N/2} v_2(j).$$
Even $k$
Only multiples of 4 contribute. Write
$$k=4j.$$
Then
$$v_2(k)-1=v_2(4j)-1=v_2(j)+1.$$
Hence the even contribution is
$$\sum_{j\le N/4}(v_2(j)+1).$$
Therefore
$$\boxed{ Q(N)=904961,A(N/2)+A(N/4)+N/4 }$$
where
$$A(n)=\sum_{j\le n} v_2(j).$$
But
$$A(n)=v_2(n!)=n-s_2(n),$$
with $s_2(n)$ the binary digit sum.
So
$$Q(N)=904961,(N/2-s_2(N/2)) +(N/4-s_2(N/4)) +N/4.$$
Step 4: Evaluate at $N=10^{12}$
We need:
$$N/2=500000000000, \qquad N/4=250000000000.$$
Binary digit sums:
$$s_2(500000000000)=13, \qquad s_2(250000000000)=13.$$
Thus
$$A(500000000000)=499999999987,$$
$$A(250000000000)=249999999987.$$
Now compute:
$$Q(10^{12}) = 904961\cdot 499999999987 +249999999987 +250000000000.$$
This equals
$$\boxed{452480999988235494}.$$
Python implementation
def v2_sum(n):
# Sum_{k<=n} v2(k) = n - popcount(n)
return n - bin(n).count("1")
def Q(N):
m = 904961
return (
m * v2_sum(N // 2)
+ v2_sum(N // 4)
+ N // 4
)
print(Q(8)) # 2714886
print(Q(10**12)) # 452480999988235494
Code walkthrough
v2_sum(n) uses the classical identity
$$\sum_{k=1}^n v_2(k)=v_2(n!)=n-s_2(n).$$
Then the function Q(N) directly evaluates the derived closed formula:
$$Q(N)=904961,A(N/2)+A(N/4)+N/4.$$
The sample check:
Q(8)
returns
$$2714886,$$
matching the problem statement.
Final verification
- The answer is positive and roughly of size
$$904961 \times 5\times10^{11}\approx 4.5\times10^{17},$$
which matches the computed magnitude.
- The sample $Q(8)$ checks exactly.
- The derivation uses only exact valuation identities and LTE.
Answer: 452480999988235494