包除原理（Inclusion-exclusion principle）でProject Euler Problem 1を解く（その１）

Last updated at 2024-02-26Posted at 2021-09-24

本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題の要約：1000未満の3か5の倍数となる自然数の和を求めよ

これは一見単純な問題で以下のようなプログラムで解けると思いますが。さすがにProject Eulerの最初の問題だけあってよく考えると深い示唆が含まれています。

n = 1000
s = 0
for i in range(1,n):
  if (i % 3 == 0) or (i % 5 == 0):
    s += i
print(s)

n = 1000
# addm: n未満のkの倍数の集合を返す
def setm(n, k):
  s = set()
  for i in range(k,n,k):
    s.add(i)
  return s

print(sum(setm(n, 3) | setm(n, 5)))

いずれも直感的に分かりやすいと思いますが、計算量が$O(n)$なので、nが大きくなると時間がかかることになります。

そこで包除原理（Inclusion-exclusion principle）の登場となります。この問題のように集合がA={3の倍数}とA={5の倍数}の2つであれば、$|A|$を集合$A$に含まれる数の和とすると以下のように表せます。

$|A\cup B|=|A|+|B|-|A\cap B|$

n未満のkの倍数の和を返す関数msumを作ってこれをプログラムにすると、以下のように。

n = 1000
# n未満のkの倍数の和を返す（等差数列の和の公式を使用）
def msum(n, k):
  c = (n-1) // k
  return c*(c+1)*k//2

print(msum(n,3)+msum(n,5)-msum(n,3*5))

ポイントは和の計算に等差数列の公式を使っているのでnが大きくなっても計算量は変わらない（$O(1)$）になることです。Project Eulerの問題はこれからnがとてつもなく大きいもの出てきますから。このような考え方が必要になってきます。

次回のその２では包除原理をもう少し詳しく見てみたいと思います。

以下に各々の解の計算時間を示します。nが大きくなるとその差は歴然となります。

*1 Google Colab RAM不足のため