prolog勉強プロジェクトオイラー201～210 - c21coterie @ ウィキ

プロジェクトオイラーの問題を堀江伸一子と私がProlog言語で解くページ。
最近スランプ気味。

Problem 201 「唯一の和を持つ部分集合」 †

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20201
数の集合Aについて, sum(A)でAの要素の和を表す. 集合B = {1,3,6,8,10,11}を考えよう. Bの3要素の部分集合は20個あり, それぞれ和は以下になる.

sum({1,3,6}) = 10,
sum({1,3,8}) = 12,
sum({1,3,10}) = 14,
sum({1,3,11}) = 15,
sum({1,6,8}) = 15,
sum({1,6,10}) = 17,
sum({1,6,11}) = 18,
sum({1,8,10}) = 19,
sum({1,8,11}) = 20,
sum({1,10,11}) = 22,
sum({3,6,8}) = 17,
sum({3,6,10}) = 19,
sum({3,6,11}) = 20,
sum({3,8,10}) = 21,
sum({3,8,11}) = 22,
sum({3,10,11}) = 24,
sum({6,8,10}) = 24,
sum({6,8,11}) = 25,
sum({6,10,11}) = 27,
sum({8,10,11}) = 29.
これらの和は1つしか現れない場合もそうでない場合もある. 集合Aについて, U(A,k)で, Aのk要素の集合全体について和を取ったときに1回のみ現れる和の集合を表す. 上の例をとると, U(B,3) = {10,12,14,18,21,25,27,29}であり, sum(U(B,3)) = 156となる.
今, 100個の要素を持つ集合 S = {1^2, 2^2, ..., 100^2}を考える. Sの50要素の部分集合は100891344545564193334812497256個ある.
50要素の部分集合の和の中で1回のみ現れる和の集合の総和を求めよ. すなわち, sum(U(S,50))を求めよ.

解法
正しい解法は数式的なものや判別関数的なものかもしれませんがとりあえず、その数字をｒ個の数で作れたとき。
0.1.2以上はたくさんと数えて２ビットごとに分割して数えました。
[Num,Counts]のNumは作れた数字。
Counts変数の0～１ビット目は０個の数で数Numを作れた。
２～３ビット目は１個の数で作れたNumの種類数、２r～２r+1ビット目はr個の数でNumを作れた種類数を０個、１個、２個（３個以上を含む）以上の３種類に分けて記録してメモ化計算をしていきます。
あとはこれを利用して組み合わせを検証していくだけです。
こんなアクロバティックな解法私だけかもしれません。
普通の人は数学的知識を使って解くような気がします。
コード実行時間59.502sec。
C++なら１秒切れます。

maskA(1690200800304305868662270940501):-!.
maskB(5070602400912917605986812821503):-!.

count_union(Counts1,Counts2,Result):-
	maskA(MaskA),
	maskB(MaskB),
	Result is
	MaskB /\ (((MaskA /\ Counts1)+(MaskA /\ Counts2))
	\/ (Counts1\/Counts2)).

count_union_all([],Y,[Y]):-!.
count_union_all([[Num,Counts]|Rest],[Num,Counts1],Result):-
	!,
	count_union(Counts,Counts1,Counts2),
	count_union_all(Rest,[Num,Counts2],Result).
count_union_all([X|Rest],Y,[Y|Result]):-
	!,
	count_union_all(Rest,X,Result).

one_calc(_,X,X).
one_calc(N,[Num,Counts],[NextNum,NextCounts]):-
	maskB(MaskB),
	NextNum is Num+N^2,
	NextCounts is (Counts*4)/\MaskB.

calc_next(N,Memos,Result):-
	member(M,Memos),
	one_calc(N,M,Result).

check_sum([],0):-!.
check_sum([[Num,Counts]|Rest],Result):-
	1=:=((Counts>>100) /\ 3),
	!,
	check_sum(Rest,Re),Result is Re+Num.
check_sum([_|Rest],Result):-
	!,
	check_sum(Rest,Result).

search(101,Memos):-!,check_sum(Memos,Ans),write(Ans).
search(N,Memos):-
	findall(M,calc_next(N,Memos,M),Memos1),
	msort(Memos1,Memos2),
	[Top|Memos3]=Memos2,
	count_union_all(Memos3,Top,Memos4),
	N1 is N+1,
	length(Memos4,Len),
	write(Len),nl,
	search(N1,Memos4).
main201:-search(1,[[0,1]]).

問201のC++版

原理は同じだが、速度は0.06秒と段違い。
やっぱC++は早い。

#include<stdio.h>
#include<math.h>
#include<iostream>
#include<time.h>

const int LIMIT=340000;
__int64 count2[LIMIT]={0,},count1[LIMIT]={0,};

int main(){	
	__int64 mask=pow(2,51),up=0,add,ans=0;
 	mask--;
count1[0]=1;
 	clock_t start,end;
 	start = clock();
	
	for(int i=1;i<=100;i++){
		add=i*i;
		for(int k=up;k>=0;k--){
			int p=k+add;
 			count2[p]=mask&(count2[p] | (count2[k]<<1) | (count1[p] & (count1[k]<<1)));
			count1[p]=mask&(count1[p]  ^ (count1[k]<<1));
		}
		up+=add;
	}
 	for(int i=1;i<=up;i++){
		if((((count2[i]>>50)&1)==0) && (((count1[i]>>50)&1)==1))ans+=i;
	}
	end= clock();
	std::cout<<ans;
	 printf("\n%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC);
}

Problem 205 「サイコロゲーム」 †

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20205
ピーターは4面のサイコロを9つ持っている. サイコロの各面には1, 2, 3, 4と書いてある. コリンは6面のサイコロを6つ持っている. サイコロの各面には1, 2, 3, 4, 5, 6と書いてある.
ピーターとコリンはサイコロを投じ, 出た目の合計を比べる. 合計が多い方が勝ちである. もし出た目の合計が等しければ勝負は引き分けになる.
ピーターがコリンに勝つ確率はいくつだろうか? 10進7桁にroundし, 0.abcdefgという形で回答欄に入力せよ.

解法
今日は完全にスランプなので一番簡単な問題を探して解。
なぜこんな簡単な問題が２００番台にあるのだろう、５０番台位の問題に思える？
とりあえずC++で書いて後でPrologに移植することにする。
サイコロの和の組み合わせ数を求めて、あとはピーターがコリンに勝つ確率は
Σコリンが出した合計値＜ピーターが出したある合計値
を全体の組み合わせで割れば答えが出ます。

計算量は9*32*4+6*30*6+37+37+pow関数２回で2500回+pow２回くらい。

#include<stdio.h>
#include<math.h>

int main(){
	double xai4[38]={0,};
	double xai6[38]={0,};

	
	xai4[0]=1;
	xai6[0]=1;
	for(int j=0;j<9;j++){
		for(int i=32;i>=0;i--){
			for(int k=1;k<=4;k++)xai4[i+k]+=xai4[i];
			xai4[i]=0;
		}
	}
	for(int j=0;j<6;j++){
 		for(int i=30;i>=0;i--){
			for(int k=1;k<=6;k++)xai6[i+k]+=xai6[i];
			xai6[i]=0;
		}
	}

	for(int i=2;i<37;i++){
		xai6[i]+=xai6[i-1];
 	}
	double perm=0;
	for(int i=2;i<=36;i++){
		perm+=xai4[i]*xai6[i-1];
	}
	
	printf("%.7f",perm/(pow(4,9)*pow(6,6)));
}