約数の総和は約数関数（Wikipedia）によると以下のように$\sigma$で表されます。
\begin{align}
&\sigma(n) = \sum_{d|n}d \\
&nの小さい順にならべると \\
&1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 12, 28, 14, 24, 24 … \\
&となります
\end{align}
ここで例題を考えます。
【例題１】$N = 10^3$以下の正の整数でその約数の総和が67の倍数になっているものの和を求めよ
これを素直に求めると以下の4つの整数が求まります。
import sympy as sp
def sigma(n):
return sum(sp.divisors(n))
m, N =67, 10**3
ans = 0
for n in range(1,N+1):
s = sigma(n)
if s % m == 0:
print(f"n = {n}: sum({sp.divisors(n)}) = {s} = {m} x {s//m}")
ans += n
print(f"Answer = {ans}")
# n = 401: sum([1, 401]) = 402 = 67 x 6
# n = 802: sum([1, 2, 401, 802]) = 1206 = 67 x 18
# n = 841: sum([1, 29, 841]) = 871 = 67 x 13
# n = 937: sum([1, 937]) = 938 = 67 x 14
Answer = 2981
ただこれではすべての約数を求めているのでNが大きくなると効率が悪くなります。そこで約数の和の公式を使って和を積の形に変形します。
\begin{align}
& n = p_1^{a_1}p_2^{a_2}...p_r^{a_r} \\
&と表されるとき \\
&\sigma(n) = \sum_{d|n}d \\
&=(1+p_1+p_1^2+...+p_1^{a_1})(1+p_2+p_1^2+...+p_2^{a_2})...(1+p_r+p_r^2+...+p_r^{a_r})\\
&= \prod_{i=1}^{r}\sum_{j=0}^{a_j}p_{i}^{j}
\end{align}
積の形になれば以下の順序で67の倍数の約数の和を探すことができます。
まず$(1+p+p^2+...+p^{a})$は何度も出てくるので関数$s(p,a)$を定義します
\begin{align}
s(p,1) &= (1+p) \\
s(p,2) &= (1+p+p^2) \\
&: \\
s(p,a) &= (1+p+p^2, ..., p^a)
\end{align}
- $s(p,a)$が67の倍数になる$(p,a)$を求める。すると$p^{a}$の倍数の$\sigma$は67の倍数になるのでその総和を求める
- ただし$p^{a+1}$を因数にもつとは67の倍数ではなくなるので$p^{a+1}$ の倍数を除く
-
- で求めた$p^a$を2つ素因数に持つ場合は重複して数えているのでその倍数をやはり除く（包除原理なので3つ素因数に持つ数がある場合は、さらに加える必要がある）
こう書くと分かりづらいので例題の場合で示します。
1. s(p,a)が67の倍数になる(p,a)を求める
m, N =67, 10**3
for a in range(1,3+1):
p = 2
while p <= N:
if p**a > N: break
x = sum([p**k for k in range(a+1)])
if x % m == 0:
print(f"s({p},{a}) = {x} = {m}x{x//67}, {p}^{a} = {p**a}")
p = sp.nextprime(p)
# s(401,1) = 402 = 67x6, 401^1 = 401
# s(937,1) = 938 = 67x14, 937^1 = 937
# s(29,2) = 871 = 67x13, 29^2 = 841
このプログラムの結果以下のように3つの$p$が求まります。その倍数を含む4つの数$401, 802, 937, 841$が求まります。
| $s(p,a)$ | $n=kp^a$ |
|---|---|
| $s(401,1) = 402=67 \times 6$ | $401, 802$ |
| $s(937,1) =938=67*14$ | $937$ |
| $s(29,2) =871=67*13$ | $29^2=841$ |
$N=10^3$ではこれですべて求まりますが$N=10^6$になると２と３のケースが出てきますので【例題２】を考えます
【例題２】$N = 10^6$以下の正の整数でその約数の総和が67の倍数になっているものの和を求めよ
2. p^(a+1)の倍数を除く
$n = 29^3 = 24389$の場合には$s(29,3)=25260=377 \times 67 + 1$となり67の倍数にならなくなります。従って16411の倍数を除く必要があります。
3. s(p,a)が67の倍数のp^aを2つ素因数に持つ数を除く
$n = 401 \times 841 = 337241$ のような$p$を2つ持つ場合は重複しているのでこれを除きます。（包除原理）
さらにこれが３つある場合にはまた足す必要がありますが。最小の$p=401$でも$p^3=401^3=64481201>10^6$なので、今回は考えなくて良さそうです。
これをコードにしたものが以下になります。(2)のケースが４個、(3)のケースが８個あることが分かります。
m, N =67, 10**6
def s(p,a):
return sum([p**k for k in range(a+1)])
def multsum(x, N):
n = N // x
return x*(1 + n)*n//2
pl, ans = [], 0
print("--- (2). Subtract sum of multiples of p^(a+1) ---")
for a in range(1,3+1):
p = 2
while p**a <= N:
s1 = s(p,a)
if s1 % m == 0:
pl.append(p**a)
m1, m2 = multsum(p**a,N), multsum(p**(a+1),N)
ans += m1 # (1). Add sum of multiples of p^a
ans -= m2 # (2). Subtract sum of multiples of p^(a+1)
if m2 > 0:
print(f"s({p},{a}) = {s1}, p^a = {p**a}, p^(a+1)={p**(a+1)}, m1 = {m1}, m2 = {m2}")
p = sp.nextprime(p)
print(f"--- (3). Subtract sum of p_i x p_j multiples ---")
for i in range(len(pl)):
for j in range(i+1,len(pl)):
if pl[i]*pl[j] <= N:
print("+++",pl[i],pl[j], N // (pl[i]*pl[j]))
ans -= multsum(pl[i]*pl[j], N)
print(f"Answer: {ans}")
# --- (2). Subtract sum of multiples of p^(a+1) ---
# s(401,1) = 402, p^a = 401, p^(a+1)=160801, m1 = 1246617171, m2 = 3376821
# s(937,1) = 938, p^a = 937, p^(a+1)=877969, m1 = 533881986, m2 = 877969
# s(29,2) = 871, p^a = 841, p^(a+1)=24389, m1 = 594969655, m2 = 20998929
# s(37,2) = 1407, p^a = 1369, p^(a+1)=50653, m1 = 365269735, m2 = 9624070
# --- (3). Subtract sum of p_i x p_j multiples ---
# +++ 401 937 2
# +++ 401 1607 1
# +++ 401 1741 1
# +++ 401 2143 1
# +++ 401 2411 1
# +++ 401 841 2
# +++ 401 1369 1
# +++ 937 841 1
#Answer: 7740907171
（開発環境：Google Colab）
この考え方はProject Euler Problem 565: Divisibility of Sum of Divisorsを解くのに役に立ちます