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

集合が３つの時の包除原理

前回の投稿では、集合が２つの時の包除原理（Inclusion-exclusion principle）を使って問題を解きましたが集合が3つになるとちょっと複雑になります。Wikiのリンクから3つの場合は以下のようになります。

\begin{align}
|A\cup B \cup C|&=|A|+|B|+|C| \\
&-|A\cap B|-|B\cap C|-|C\cap A| \\
&+|A\cap B\cap C|
\end{align}

したがって前回の問題に７の倍数を含めるようにすると。包除原理を使ったプログラムは解.４のようになります。（関数msumは前回のものを使う）

問題２ 1000未満の、3か5か7の倍数となる自然数の和を求めよ

解.４ 3集合の包除原理を使ったプログラム

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

これで分かるように偶数の数の積では引いて、奇数の数の積では足せばよいことが分かるので、互いに素な数が[3,5,7,11,13]のようにリストで与えられた場合のプログラムば以下のように書けます。規則性がつかめればプログラムにできるのでそれほど複雑には見えなくなりますね。

解.５ N集合の包除原理を使ったプログラム

import itertools
import numpy as np
n = 1000
pl = [3,5,7,11,13]

mulsum = 0
for p in range(1,len(pl)+1):
  for pp in itertools.combinations(pl,p):
    mulsum += (-1)**(p-1)*msum(n,np.prod(pp))

print(mulsum)