『続・とっておきの数学パズル』を力技で解く、またはべき集合の列挙を実装してみた

帰省の暇つぶしに、と衝動買い。

続・とっておきの数学パズル

新鮮でバラエティに富んだ、しかも歯応えのある問題が多く、
頭の体操――かなりハードな――として非常にお買い得な一冊でした。

以下は本書の問3.5「真に平等な分割」という問題。

1から16までの整数を、要素の個数が等しい2つの集合に分ける。
2つの集合は、それぞれの和が等しく、各々の数の2乗の和も等しく、さらに各々の数の3乗の和も等しいように分けたい。
さて、どのように分ければいいだろうか？

これがなかなか解けなかったので、あきらめて力技で解くことにしました。
1〜16の集合を2分割する全てのパターン、要はべき集合について各種の和を求めていくことになります。
べき集合の要素を順に返すイテレータを作成する方針で実装してみました。

from datetime import datetime

def main():

  #問題3.5の各種定義
  def square_sum(s): #数値配列の二乗和を返す
    return sum(map(lambda x: x*x, s))
  def cube_sum(s): #数値配列の三乗和を返す
    return sum(map(lambda x: x*x*x, s))
  numbers_to_divide = range(1,17) #分割対象となる配列([1,2,3,..,16])
  sum1 = sum(numbers_to_divide) / 2
  sum2 = square_sum(numbers_to_divide) / 2
  sum3 = cube_sum(numbers_to_divide) / 2
    
  #配列のべき集合を走査
  
  start = datetime.now() #計算時間測定用

  for s in PowersetProvider(numbers_to_divide):
    if sum(s) == sum1 and square_sum(s) == sum2 and cube_sum(s) == sum3:
      print(s)
         
  end = datetime.now()
  print("Time: " + str(end - start))

#与えられた集合（※）のべき集合の要素を順に返すイテレータ
#（※ただし、入出力ともに集合は配列の形式をとる）
#nextでの出力対象を決定するための、マスキング変数（※）をメンバとして持つ。
#
#※マスキング変数を2進表示した数をb(bn,...,b2,b1)とすると、
#　与えられた配列のi番目を次回出力に...
#　含める:bi桁目=1, 含めない:bi=0 の値を持つ
#　例）
#　　与えられた集合が[3,1,4,5]でマスキング変数が9(2進で1010)の場合
#　　[1,5]を次回のイテレーションで出力する。
class PowersetProvider:
  
  #出力を制御するマスキング変数を初期化する。
  
  def __init__(self, provided_set):
    #マスキング変数を初期化（初回出力は空集合となる）
    self.mask = 0
    self.set = provided_set
    
  def __iter__(self):
    return self

  #現在のマスキング変数に相当する集合を出力し、
  #その後マスキング変数を1加算する。
  def __next__(self):

    #イテレーションが終了していた場合、StopIteration例外を投げる
    if self.mask == 2 ** len(self.set):
      raise StopIteration
    
    rtn = []
    #マスキング変数を2進数と見て、各桁の0/1により出力を制御する
    for i in range(0, len(self.set)):
      if (self.mask >> i & 1) == 1:
        rtn.append(self.set[i])

    self.mask += 1
    return rtn

if __name__ == '__main__': main()

結果は下記の通り。
べき集合の全要素（65536集合）を走査している割には、思ったより早いですね。
（ネタバレになるので、出力部分は＊で伏せています）
再帰でのべき集合列挙も考えられるので、また暇なときにやってみようと思います。

[＊＊＊＊＊＊＊＊＊＊＊＊＊]
[＊＊＊＊＊＊＊＊＊＊＊＊＊]
Time: 0:00:00.738000