読者です 読者をやめる 読者になる 読者になる

本履歴

購入した古本の履歴と時々プログラミング

プログラミングコンテストの練習

ポリアの組合せ論入門(-組合せ論入門)に面白い問題があったので、ちょっと考えよう。

二項係数の偶奇性を求めるアルゴリズムを考えよ

二項係数はn個の中からk個選ぶ組合せの数n_C_kだったね。n_C_k%2、なんだ簡単じゃんといかないのが面白い。

愚直に

まず検算用にn_C_kをそのまま計算し、%2する方法。Rubyならメモリの許す限り頑張れる。
(n*(n-1)*...*(n-r+1)/k*(k-1)*...*2*1)%2

2の因数の数

次に、上の愚直計算の中で、奇数は解答に全く寄与しないことがわかる。というのも、我々が欲しい結果は2のあまりだけだから、もっと言えば、2の因数の数が分母と分子で同じかどうか、だけ知ることができればOK。
分母のほうが大きくなることはない、というのも二項係数は整数だからね。

パスカルの三角形を睨んで再帰

次に、有名なmod2で作ったパスカルの三角形がまさしく答えだね。この図形はシェルピンスキーのガスケットとも呼ばれフラクタルな構造だ。つまり、再帰が使える!
n=2^k+l
と表されていた場合、r>=2^kならl_C_(r-2^k)だし、それ以外はl_C_rだから再帰が止まるように終了条件を設定すればよさそう。
ちなみにnから2^kを求めるのは、ハッカーのたのしみ―本物のプログラマはいかにして問題を解くか に載っていたn&(n-1)で一番右の1のビットが消えるというのを使った。これを繰り返していって、0になる直前の数が2^k(一つのビットだけ立っている)だから。

ポリア先生の解法

nが偶数でrが奇数ならばn_C_rは偶数、それ以外の時には、[n/2]_C_[r/2]が偶数の時その時に限りn_C_rは偶数

という命題を証明していた。上のパスカルの三角形の方法と同値なのかどうかがわからないなあ。あとで考える。

Rubyでプログラムを組んでみた。

#return parity of combination(n,r)

#straightforward
def combination_parity1(n,r)
	ret=1
	(n-r+1).upto(n){|i|ret*=i}
	2.upto(r){|i|ret/=i}
	ret&1
end

#require "memoize"
#include Memoize

def num_of_factor_2(n)
	ret=0
	while n&1==0
		n/=2
		ret+=1
	end
	ret
end

#memoize(:num_of_factor_2)

#count factor 2
def combination_parity2(n,r)
	factor2=0
	(n-r+1).upto(n){|i|factor2+=num_of_factor_2(i)}
	2.upto(r){|i|factor2-=num_of_factor_2(i)}
	factor2>0 ? 0 :1
end

#return n以下の2の冪乗
def flp2(n)
	ret=n
	while n>0
		ret=n
		n=n&(n-1)
	end
	ret
end

#recursion on (mod 2)-Pascal's triangle
def combination_parity3(n,r)
	return 0 if n<r
	return 1 if (n==r)||(r==0)
	flp=flp2(n)
	return r>=flp ? combination_parity3(n-flp,r-flp) : combination_parity3(n-flp,r)
end

#from Notes on Introductory Combinatorics by G. Polya
def combination_parity4(n,r)
	until r==0
		return 0 if (n&1==0)&&(r&1==1)
		n/=2
		r/=2
	end
	return 1
end


require "benchmark"

M=300

Benchmark.bm do |x|
  x.report(":4"){1.upto(M){|n|0.upto(n){|r|combination_parity4(n,r)}}}
  x.report(":3"){1.upto(M){|n|0.upto(n){|r|combination_parity3(n,r)}}}
  x.report(":2"){1.upto(M){|n|0.upto(n){|r|combination_parity2(n,r)}}}
  x.report(":1"){1.upto(M){|n|0.upto(n){|r|combination_parity1(n,r)}}}
end

M=300での結果は、
user system total real
:4 0.030000 0.000000 0.030000 ( 0.032172)
:3 0.050000 0.000000 0.050000 ( 0.055248)
:2 2.920000 0.000000 2.920000 ( 3.002428)
:1 5.360000 0.010000 5.370000 ( 5.370607)


M=3000でパスカル v.s. ポリアの結果は(1,2はこの時点で既に宇宙が終わるまで終わんない)
user system total real
:4 3.450000 0.000000 3.450000 ( 3.508069)
:3 6.820000 0.010000 6.830000 ( 6.835744)

だいたい2倍差でポリアの勝ち。再帰がオーバーヘッドになっているとか、flp2を求めなくていいとかあるかな。
でもflp2をmemoizeしたら逆に遅くなったw
M=300と合わせると計算量は両者同じか・・・