Project Euler Problem 561

Let S(n) be the number of pairs (a,b) of distinct divisors of n such that a divides b.

Project Euler Problem 561

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