Project Euler 106（1） - inamori’s diary

http://projecteuler.net/index.php?section=problems&id=106

まず、可能なBとCの組合せの数を数えておきましょう。まずBを一つ以上選びます。そして、残った中から一つ以上選びます。

_nC₁(_n-1C₁ + … _n-1C_n-1) + … + _nC_n-1₁C₁
= _nC₁(2^n-1 - 1) + … + _nC_n-1(2¹ - 1)
= _nC₁2^n-1 + … + _nC_n-12¹ - (_nC₁ + … + _nC_n-1)
= 3ⁿ - (2ⁿ + 1) - (2ⁿ - 2)
= 3ⁿ - 2ⁿ⁺¹ + 1

実際にはBとCが入れ替わっただけのペアを重複して数えているので、この半分になります。

(3ⁿ + 1) / 2 - 2ⁿ

さて、例えばA = { a₁, a₂, a₃, a₄ } として、

a₁ < a₂ < a₃ < a₄

とすると、

B = { a₁, a₄ }
C = { a₂, a₃ }

なら、S(B) = S(C)となるかどうかわからないのですが、

B = { a₁, a₃ }
C = { a₂, a₄ }

なら、S(B) < S(C)となるので、調べる必要がありません。
BとCを並べて、（Bが必ずa₁を含むとすると）先頭からインデックスを見て、常にBよりCのほうが小さければ調べる必要がないと判定されます。一つでもそれに反すれば調べる必要があることになります。

また、n = 6のでBとCをあわせて要素数4の場合、6つの中からどの4つを選んでも調べる必要があるBとCの組合せの数は同じです。だから、一つの組合せについて調べて₆C₄倍すればよいです。

from itertools import *
import time

def C(n, m):
    return 1 if m == 0 else C(n, m - 1) * (n - m + 1) / m

def gen_pair_of_sets(n):
    def gen(a, p, q):
        if p == 0:
            yield (), a
        elif q == 0:
            yield a, ()
        else:
            for a1, a2 in gen(a[1:], p - 1, q):
                yield (a[0],) + a1, a2
            for a1, a2 in gen(a[1:], p, q - 1):
                yield a1, (a[0],) + a2
    
    for a1, a2 in gen(tuple(range(1, n)), n / 2 - 1, n / 2):
        yield (0,) + a1, a2

def needs_rule1(a1, a2):
    return not all(a1[k] < a2[k] for k in xrange(len(a1)))

def count_pair_of_sets(n):
    counter = 0
    for m in range(4, n + 1, 2):
        counter += C(n, m) * sum(starmap(needs_rule1,
                                        gen_pair_of_sets(m)))
    return counter

t = time.clock()
N = 20
print count_pair_of_sets(N)
print time.clock() - t