prolog勉強プロジェクトオイラー71～80

プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。

と書きましたが、下記リンク先サイトの方が私のコードよりずっと賢いのでそちらを参考にした方がよいです。
http://d.hatena.ne.jp/nineties/20110119

プロジェクトオイラーについては私のサイトよりこの人のサイトの方がずっと賢いです。
問い75はブロコット木を使えば早くなるとか、76や77は分割数と見るより、1～99までのコインを組み合わせたり素数を組み合わせる動的計画法とみるなど。
読めば目からうろこな発想でかかれています。
物事や問題をシンプルな話に変換できる能力、リンク先の人は私から見れば尊敬に値します。
プロジェクトオイラーの問題は問題番号が大きいほど難しくなるのですが、頭の悪い私でも問題番号100以下の簡単な問題なら解けます。
しかし同じ簡単な問題を解くのでも解き方に頭の良さの差が大きく出るのだと実感しました。

問い71
http://projecteuler.net/problem=71
ファレイ数列の定義そのままです。
ファレイ数列が隙間を埋めていくというのは分かりますが、出てくる数字が既約というのは何だか不思議。
整数論は初歩の話でも面白いことがおおいですね、私の性根はこういう初等的な数学の話を知ることが好きなおとなしい人間だったりします。

farey(U0,D0,_,D1):-
	T is D0+D1,
	T>=1000000,
	!,
	write([U0,D0]).
farey(U0,D0,U1,D1):-
	U2 is U0+U1,
	D2 is D0+D1,
	farey(U2,D2,U1,D1).

main71:-
	farey(2,5,3,7).

問い72　Counting fractions

この問題は一番素朴な方法は試し割でPhi関数を実装、これだと遅いですがメモリ使用量は抑えられますしまあいいかなと。
試し割は無駄な処理が入りますが試し割の無駄を省こうとすると素数リストが必要となり素数リスト関係の計算コストやメモリ使用量が上がります。
整数論に通じてない素朴な身としては試し割が一番シンプルでいいという結論になってしまいました。
なにやら高速に計算する関数があるようですが原理がよくわからないので飛ばしています。

is_prime(N):-
between(2,N,N1),
N2 is N1*N1,
(N<N2->!,true;
0=:=N mod N1,!,fail).
prime_list(N):-between(2,1100,N),is_prime(N).

calc(1000000,Ans,_):-!,write(Ans).
calc(N,Ans,Ps):-N1 is N+1,
	calc_phi(N,Phi,Ps),
	Ans1 is Ans + Phi,
	calc(N1,Ans1,Ps).

main72:-findall(N,prime_list(N),Ps),
        calc(1,0,Ps).

div(N,P,Result):-0=:=N mod P,!, N1 is N//P,div(N1,P,Result).
div(N,_,N):-!.


calc_phi(N,Result,[P|_]):-
	P2 is P*P,
 	N<P2,
	!,
	(N>1 ->
	Result is N-1;
	Result is 1).
calc_phi(N,Result,[P|Ps]):-!,
	div(N,P,N1),
	calc_phi(N1,Re,Ps),
	(N>N1->
	Result is Re*(N//(N1*P))*(P-1);
	Result is Re).

問い73

遅いけどファレイ数列でとりあえず解。
もちょっとましな方法はありそう。

calc(D0,D1,_):-
	D2 is D0+D1,
	D2>12000,!,fail.
calc(D0,D1,D2):-
	D2 is D0+D1.
calc(D0,D1,D3):-
	D2 is D0+D1,
	calc(D0,D2,D3).
calc(D0,D1,D3):-
	D2 is D0+D1,
	calc(D2,D1,D3).
test:-findall(N,calc(3,2,N),List),length(List,Len),write(Len).

問い74

http://projecteuler.net/problem=74
この問題はリスト処理だけで高速に解くのはひと工夫必要でした。
一回目の変換で123456も654321も152346も数字を並び変えたものは全部同じ値になるということで纏めて計算します。
これを利用して計算量を下げました。
纏めをもっと上手にやればもっと計算量が下がりますが、そのためにはコードが膨らむのでこのへんで手打ちとしました。

kaizyo(0,1):-!.
kaizyo(1,1):-!.
kaizyo(2,2):-!.
kaizyo(3,6):-!.
kaizyo(4,24):-!.
kaizyo(5,120):-!.
kaizyo(6,720):-!.
kaizyo(7,5040):-!.
kaizyo(8,40320):-!.
kaizyo(9,362880):-!.

calc(0,0):-!.
calc(N,Result):-M is N mod 10,N1 is N//10,calc(N1,Re),
	kaizyo(M,Re1),Result is Re+Re1.

calc2([],1):-!.
calc2([P|Rest],Result):-
	calc2(Rest,Re),
	kaizyo(P,P1),
	Result is Re*P1.

sum([],Result,Result):-!.
sum([X|Rest],Sum,Result):-Sum1 is Sum+X,sum(Rest,Sum1,Result).

perm_calc(Deep,Perm,Result):-
	Deep1 is Deep-1,
	kaizyo(Deep1,P1),
	reverse(Perm,RevPerm),
	[_|Perm1]=RevPerm,
       select(P2,Perm1,Perm2),
	P3 is P2-1,
	calc2([P3|Perm2],Div1),
	Result is P1//Div1.

perm_calc_w(Deep,Perm,Result):-
	findall(Re,perm_calc(Deep,Perm,Re),List),
	sum(List,0,Result).

perm(_,_,_,_,7,_):-!,fail.
perm(N,Down,K,Perm,Deep,[N,Perm1]):-
	Down>0,
	perm_calc_w(Deep,[K|Perm],Perm1).

perm(N,Down,K,Perm,Deep,Result):-
	K1 is K+1,
	Deep1 is Deep+1,
	kaizyo(Down,P1),
	N1 is N+P1,
	perm(N1,Down,K1,Perm,Deep1,Result).
perm(N,Down,K,Perm,Deep,Result):-
	Deep1 is Deep+1,
	Down1 is Down+1,
	between(Down1,9,Down2),
	kaizyo(Down2,P1),
	N1 is N+P1,
	perm(N1,Down2,1,[K|Perm],Deep1,Result).
test(Result):-perm(0,0,0,[],0,Result).

calc3(_,_,58):-!.
calc3(N,List,Deep):-
	Deep1 is Deep+1,
	calc(N,N1),
	not(member(N1,List)),
	calc3(N1,[N1|List],Deep1).

search(List,P):-
	member([N,P],List),
	calc3(N,[N],0).


main74:-findall(Re,test(Re),List),
	findall(P,search(List,P),Ans),
	sum(Ans,0,AnsN),write(AnsN).

独り言
http://share-wis.com/
この教育系サイトに誰かスポンサーやパトロンや協力者が増えたらいいなと思う、社会的責任投資の一種としてそんなに悪くないと思う。
勉強サイトだけど評価がゲームに近く、一動画が短くサクサク見れて見た動画の数だけプチ勉強した気分になれるというのがいい。

独り言２
ここら辺までのチュートリアル問題と違い問75あたりから難易度がすこしずつ上がり始める。
流石に基本リスト処理しかないPrologでは計算量的に厳しいものが増えてきた。
配列とstd::mapと優先順位付きキューがPrologにもあればまだ楽なんだけどな。

問い75　Singular integer right triangles

http://projecteuler.net/problem=75
この問題は配列があれば高校生でも簡単に解けます。
ListしかないPrologでは一工夫が必要でした。
色々試行錯誤した結果、ListでC++のstd::mapもどきを再現して解きました。
メモリ使用量やCPU使用率を抑えて30秒くらいで答えが出ます。
赤黒木を再現するほどの気力はなかったので375000を真中において左右を生成する方法でstd::mapもどきを生成しています。
赤黒ならもう数倍速くなったかもしれません。

check(N,R,R):-0=:=N mod R,R>5.
check(N,R,D):-0=:=N mod R,D is N // R,D>R,D>5.

farey(LU,LD,RU,RD,_):-
	M is LD+RD,
	N is LU+RU,
	Len is M*(M+N),
	750000<Len,
	!,
	fail.

farey(LU,LD,RU,RD,Result):-
	M is LD+RD,
	N is LU+RU,
	0<N,
	N<M,
	1=:=(M-N) mod 2,
	Result is M*(M+N).
farey(LU,LD,RU,RD,Result):-
	MD is LD+RD,
	MU is LU+RU,
	farey(LU,LD,MU,MD,Result).
farey(LU,LD,RU,RD,Result):-
	MD is LD+RD,
	MU is LU+RU,
	farey(MU,MD,RU,RD,Result).

tree(L,R,[]):-
	Sa is R-L,
	Sa<2,
	!.
tree(L,R,[L1,[M,0],R1]):-
	!,
	M is (L+R)//2,
	tree(L,M,L1),
	tree(M,R,R1).

tree_add(N,[L,[N,C],R],[L,[N,C1],R]):-
	C1 is C+1,!.
tree_add(N,[L,[N1,C],R],[L,[N1,C],R1]):-
	N1<N,
	!,
	tree_add(N,R,R1).
tree_add(N,[L,[N1,C],R],[L1,[N1,C],R]):-
	N<N1,
	tree_add(N,L,L1).

tree_search(N,[_,[N,C],_],C):-!.
tree_search(N,[L,[N1,_],_],Result):-
	N<N1,
	!,
	tree_search(N,L,Result).
tree_search(N,[_,[N1,_],R],Result):-
	N1<N,
	!,
	tree_search(N,R,Result).


tree_add_w(L,U,Tree,Tree):-
	T is L*U,
	750000<T,
	!.
tree_add_w(L,U,Tree,Result):-
	!,
	U1 is U+1,
	T is L*U,
	tree_add(T,Tree,Tree1),
	tree_add_w(L,U1,Tree1,Result).

tree_count([],0):-!.
tree_count([L,[_,C],R],Result):-
	!,
	(C=:=1 -> MC is 1;MC is 0),
	tree_count(L,LC),
	tree_count(R,RC),
	Result is LC+RC+MC.


add_farey_to_tree([],Tree):-!,tree_count(Tree,Result),write(Result).
add_farey_to_tree([L|Lens],Tree):-
	!,
	tree_add_w(L,1,Tree,Tree1),
	add_farey_to_tree(Lens,Tree1).

sum([],Sum,Sum):-!.
sum([X|Rest],Sum,Result):-Sum1 is Sum+X,sum(Rest,Sum1,Result).

問い76

100の分割数-1の数を求めよという問題。
Wikiの定義通り実装。

calc([],_,_,_,0):-!.
calc([X|Rest],K,Deep,Deep,Result):-
	!,
	(K<0->K1 is -K+1;K1 is -K),
 	M is (K1*(3*K1-1))//2,
	Deep1 is Deep+1,
	calc(Rest,K1,Deep1,M,Re),
	(1=:= K mod 2->Add is X;Add is -X),
	Result is Re+Add.
calc([_|Rest],K,Deep,M,Result):-
	!,
	Deep1 is Deep+1,
	calc(Rest,K,Deep1,M,Result).
calc1(List):-
	calc(List,1,1,1,Re),
	length(List,Len),
	write([Len,Re]),nl,
	Ans is Re-1,
	(Len=:=100->!,write([ans,Ans]);calc1([Re|List])).

問い77

http://projecteuler.net/problem=77
ちょっとよんで即座に速度の出る方法を思いつかなかったので他の人の考えを検索。
素数和で素数を小さい順から並べた時末尾の素数と同じより大きな数字を足す動的計画法でアセプト。
という一行を読んで発想法を3秒で理解。
なるほどなるほどとPrologで実装したらすこしめんどくさかった。
std::mapがないから、ソートして集計し直す処理でそれを代用するあたりがめんどくさい。

is_prime(N):-
	between(2,N,N1),
	M is N1*N1,
	(N<M ->!,true;
	0=:=N mod N1,!,fail).

sums([],[N,LP,C],[[N,LP,C]]):-!.
sums([[N,LP,C]|Rest],[N,LP,C1],Result):-
	!,
	C2 is C+C1,
	sums(Rest,[N,LP,C2],Result).
sums([X|Rest],[N1,LP,C1],[[N1,LP,C1]|Result]):-
	!,
	sums(Rest,X,Result).
n_sum(_,[],Result,Result):-!.
n_sum(N,[[N,_,C]|Rest],Sum,Result):-
	!,
	Sum1 is Sum+C,
	n_sum(N,Rest,Sum1,Result).
n_sum(N,[_|Rest],Sum,Result):-
	n_sum(N,Rest,Sum,Result).

calc3(_,[],_,[]):-!.
calc3(_,_,[],[]):-!.
calc3(N,[[N1,LP,C]|List],[P|Ps],[[N,P,C]|Result]):-
	LP=<P,
	N=:=N1+P,
	!,
	calc3(N,List,[P|Ps],Result).
calc3(N,[[N1,_,_]|List],[P|Ps],Result):-
	T is N1+P,
	N>T,
	!,
	calc3(N,List,[P|Ps],Result).

calc3(N,List,[_|Ps],Result):-
	!,
 	calc3(N,List,Ps,Result).
calc2(N,Ps,List):-
	calc3(N,List,Ps,List1),
	append(List,List1,List2),
	msort(List2,List3),
	[Top|List4]=List3,
	sums(List4,Top,List5),
	n_sum(N,List5,0,Perm),
 	write([N,Perm]),nl,
	N<72,
	N1 is N+1,
	calc(N1,Ps,List5).


calc(N,Ps,List):-
	is_prime(N),
	!,
	calc2(N,[N|Ps],List).
calc(N,Ps,List):-
	!,
	calc2(N,Ps,List).

main77:-calc(2,[],[[0,0,1]]).

問い78　Coin partitions

http://projecteuler.net/problem=78
Wikiの公式通り実装しましたが配列と違いリスト処理なので速度が出ません。
p(n+1)とp(n)は足すものが一つずれてるだけなんだからこれを利用して計算量を下げる方法がありそうなのですが。
それを考慮した実装は少しめんどくさそうです。

calc([],_,_,_,0):-!.
calc([X|Rest],K,Deep,Deep,Result):-
!,
 	(K<0->K1 is -K+1;K1 is -K),
	M is (K1*(3*K1-1))//2,
	Deep1 is Deep+1,
	calc(Rest,K1,Deep1,M,Re),
	(1=:= K mod 2->Add is X;Add is 1000000-X),
	Result is (Re+Add) mod 1000000.
calc([_|Rest],K,Deep,M,Result):-
	!,
	Deep1 is Deep+1,
	calc(Rest,K,Deep1,M,Result).
calc1(List):-
	calc(List,1,1,1,Re),
	length(List,Len),
	!,
	write([Len]),
	(0=:=Re mod 1000000 ->!,write([ans,Len]);calc1([Re|List])).

問い79　Passcode derivation

http://projecteuler.net/problem=79
いくつかの桁に重複数字が出たらこの問題難問ですが、重複数字がないので単純なトポロジカルソートで間に合います。

set(A,B,_,A,B).
set(A,_,C,A,C).
set(_,B,C,B,C).
set_w(A,B,C,A2,B2):-set(A,B,C,A1,B1),A2 is A1-48,B2 is B1-48.

read79([A1,B1]):-seen,
     see('e79.txt'),
     repeat,
     get0(A),
     get0(B),
     get0(C),
      get0(D),
      set_w(A,B,C,A1,B1),
      (D== -1 ->!,fail;true).

searchK(List,Ans,Bads,K):-
	between(0,9,K),
	not(member([_,K],List)),
	not(member(K,Ans)),
	not(member(K,Bads)),
	!.

search2(List,Ans,Bads,K):-
	select([K,_],List,List1),
	!,
	search2(List1,Ans,Bads,K).
search2(List,Ans,Bads,_):-
	!,
	search(List,Ans,Bads).

write_ans(Ans):-reverse(Ans,Ans1),nl,nl,write([ans,Ans1]).
search([],Ans,Bads):-!,between(0,9,N),
	not(member(N,Ans)),
	not(member(N,Bads)),
	!,
	write_ans([N|Ans]).

search(List,Ans,Bads):-
	!,
	searchK(List,Ans,Bads,K),
	search2(List,[K|Ans],Bads,K).

bad_nums(List,N):-
	between(0,9,N),
	not(member([N,_],List)),
	not(member([_,N],List)).

main79:-findall(E,read79(E),List),
	findall(BN,bad_nums(List,BN),Bads),
	write(Bads),
	search(List,[],Bads).

問い80

出題側としてはこの問題を通じて連分数展開の意味を学習してほしいのだと思いますが。
連分数展開がめんどくさかったので、ｎ*10＾200の平方根の整数部分だけを取得する方法で解きました（＾＞＾）ふふーん。
速度は出ませんが気にするほどではありません。
そこ手抜きと言わない。
まあそのうち連分数展開のコードも書くと思います。

sqrt_search(Down,Up,N,M):-
	M is (Down+Up)//2,
	T1 is M*M,
	T2 is (M+1)*(M+1),
	T1 =<N,
	N <T2,
	!.
sqrt_search(Down,Up,N,Result):-
	M is (Down+Up)//2,
	T1 is M*M,
	T1<N,
	!,
	sqrt_search(M,Up,N,Result).
 sqrt_search(Down,Up,N,Result):-
	M is (Down+Up)//2,
	sqrt_search(Down,M,N,Result).

to_bignum(Re3):-
	between(1,100,N),
	sqrt(N,N1),
	floor(N1,N2),
	not(N2*N2=:=N),
	Re is N*10^200,
	Down is 10^100,
	Up is N*N*10^100,
sqrt_search(Down,Up,Re,Re1),
 	number_codes(Re1,Re2),
	sum(Re2,0,0,Re3).

sum(_,100,Sum,Sum):-!.
sum([X|Rest],Deep,Sum,Result):-
 	Sum1 is Sum + X-48,
	Deep1 is Deep+1,
	sum(Rest,Deep1,Sum1,Result).

nsum([],Sum,Sum):-!.
nsum([X|Rest],Sum,Result):-
	Sum1 is Sum+X,
 	nsum(Rest,Sum1,Result).

test:-findall(X,to_bignum(X),List),
	nsum(List,0,Ans),write(Ans).

タグ：

+ タグ編集

「prolog勉強プロジェクトオイラー71～80」をウィキ内検索

最終更新：2013年09月02日 07:58

c21coterie @ ウィキ

更新履歴

メニュー

リンク

他のサービス