「prolog勉強プロジェクトオイラー111~120」の編集履歴(バックアップ)一覧に戻る

prolog勉強プロジェクトオイラー111~120 - (2013/11/13 (水) 17:57:06) のソース

プロジェクトオイラーという数学問題の掲載されたサイトの問題を堀江伸一さんがPrologで解くページ。
Prolog言語の特徴。
配列はない、stdコンテナにあたるものはもちろんない。
何があるかっていうとリストしかない。
凄く不便。
しかも破壊的代入がない、そのため動的計画法と物凄く相性が悪い、破壊的代入を前提としたデータ構造とも相性が悪いというか現実的な計算量で実装できない物が多数。
基本ローカル変数しかないから変数が増える問題では関数(Prologでは厳密に述語といい関数とわけて考える)の引数が凄いことになる。
引数を減らすと、解集合を巨大なリストで得てそこから一番いい解を探すような富豪プログラムになる。
利点はバックトラックがあるくらい、このバックトラックという機能は論理的処理には向いており、このためマサチューセッツ工科大学でも論理枠でPrologの講義が一応行われている。
Prologは論理的処理は得意だが、計算処理は苦手、よってプロジェクトオイラーのサイトの問題を解くには向いてない。
Prologという不便な言語でプロジェクトオイラーに挑戦してるのは多分私くらい。
何故Prologでやるかっていうかとこの言語、記述してると妙に楽しいから、それだけ、プロジェクトオイラーは趣味のサイトだからそれでもいいとおもってる。
C++が一番好きだけど書いてて楽しいのはProlog。
Javaは言語使用がヘビーで中々勉強する気にならない。
ちなみに私は整数論の教育を受けたことがないので整数論を中心としたプロジェクトオイラーのサイトでは平凡なコードです。
多様な問題がある会津大学オンラインジャッジ(AOJ)という競技プログラムの過去問を扱ったサイトでは私。
たまにですがコード実行速度1位を取ったり取った後で一月後位に1位を抜かれたりしています。


*Problem 111 「重複桁を持つ素数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20111
この問題例えばMの値が正しいかは置いといて
M(n,1)=2,M(n,2)=5でXは1を2個、2を5個持ち他に何個かの桁がそれ以外の数字で構成された素数があったとして、これが素数なら
1と2両方で足すのかという疑問があったが、計算してみると
M(10,d)はdが6以下になることがないみたいなのでこの疑問は解消されないままアセプト。
解けたものの分かりにくい問題だ。
test述語で駆動。
出てきた答えをExcelで集計してアセプト。

 not_prime(N):-
 	S is floor(sqrt(N)),
 	between(2,S,M),
 	0=:=N mod M,
 	!,
 	true.
 is_prime(N):-not(not_prime(N)).
 
 create_perm(-1,_,_,_,_,_):-!,fail.
 create_perm(0,_,_,_,Num,Num):-
 	is_prime(Num),!.
 
 create_perm(Rest,Nums,D,Rest,Num,Result):-
 	!,
 	Num1 is Num*10+D,
 	Num1>0,
 	Rest1 is Rest-1,
 	create_perm(Rest1,Nums,D,Rest1,Num1,Result).
 create_perm(Rest,Nums,D,M,Num,Result):-
 	member(N,Nums),
 	(N=:=D-> M1 is M-1 ; M1 is M),
 	M1> -1,
 	Num1 is Num*10+N,
 	Num1>0,
 	Rest1 is Rest-1,
 	create_perm(Rest1,Nums,D,M1,Num1,Result).
 
 
 search_m(Nums,D,Ans):-
 	Keta is 10,
 	between(1,9,T),
 	M is Keta-T,
  	findall(Num,create_perm(Keta,Nums,D,M,0,Num),Ans),
 	length(Ans,Len),
 	Len>0,
 	%write([len,M]),
 	!.
 
 sum([],0):-!.
 sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
 
 test:-Nums=[9,8,7,6,5,4,3,2,1,0],
 	between(0,9,D),
 	search_m(Nums,D,Ans),
 	sum(Ans,Ans1),
 	write(Ans1),nl,fail.
 	%write(Ans),nl.


*Problem 112 「はずみ数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20112
とりあえず小さい方から全検証してみた。
割り算は遅いという太古の神話を思い出し、割り算を減らすために
数字を最初からリストにしたら少し早くなるかなと思って計算してみたけれど思ったより早くならなかった。
うーん逆に遅くなった?


 add([],0,[]):-!.
 add([],1,[1]):-!.
 add(Xs,0,Xs):-!.
 add([X|Xs],1,[X1|Result]):-
 	Xt is (X+1), 
 	(Xt=:=10->add(Xs,1,Result),X1=0;add(Xs,0,Result),X1 is Xt).
 check_up([X,X1]):-X1=<X,!.
 check_up([X,X1|Xs]):-X1=<X,check_up([X1|Xs]).
 check_down([X,X1]):-X=<X1,!.
 check_down([X,X1|Xs]):-X=<X1,check_down([X1|Xs]).
 check(Xs):-check_up(Xs),!.
 check(Xs):-check_down(Xs),!.
 
 calc(N,Ca,Cb,_):-Ca*99=:=Cb,N1 is N-1,write([N1,Ca,Cb]),!.
 calc(N,Ca,Cb,Xs):-
 	check(Xs),
  	!,
 	N1 is N+1,
 	add(Xs,1,Xs1),
 	!,
 	Ca1 is Ca+1,
 	calc(N1,Ca1,Cb,Xs1).
 calc(N,Ca,Cb,Xs):-
 	N1 is N+1,
 	add(Xs,1,Xs1),
 	Cb1 is Cb+1,
 	!,
 	calc(N1,Ca,Cb1,Xs1).
 
 main:-calc(100,99,0,[0,0,1]).





*Problem 113 「非はずみ数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20113
10^100以下のはずみ数出ない数の割合を求める問題。
答えは14桁。


 sum([],0):-!.
 sum([X|Xs],Sum1):-sum(Xs,Sum),Sum1 is Sum+X.
 
 next_calc([],0,[]).
 next_calc([Count|Counts],Sum1,[Sum1|Result]):-
 	next_calc(Counts,Sum,Result),
 	Sum1 is Sum+Count.
 calc_up(100,_,Sum,Sum):-!.
 calc_up(N,Counts,Sum,Result):-
 	next_calc(Counts,_,Counts1),
 	sum(Counts1,Add),
 	Sum1 is Sum+Add-1,
 	N1 is N+1,
 	!,
 	calc_up(N1,Counts1,Sum1,Result).
 calc_down(100,_,Sum,Sum):-!.
 calc_down(N,Counts,Sum,Result):-
 	next_calc(Counts,_,Counts1),
 	[_|Counts2]=Counts1,
 	sum(Counts2,Add),
 	Sum1 is Sum+Add-1,
 	N1 is N+1,
  	!,
 	calc_down(N1,Counts1,Sum1,Result).
 
 main113:-Counts=[10,9,8,7,6,5,4,3,2,1],
 	calc_up(2,Counts,0,AnsUp),
 	calc_down(2,Counts,0,AnsDown),
 	Ans is AnsUp+AnsDown+99-8*98,
 	write([AnsUp,AnsDown,Ans]).


*Problem 114 「ブロックの組み合わせ方の数え上げ その1」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20114
長さ 7 ユニットからなる 1 列上に, 最低 3 ユニットの長さを持つ赤ブロックが置かれている. ただしどの赤ブロック同士も, 少なくとも 1 ユニットの黒い正方形が間にある(赤ブロックは長さが異なってもよい). これを敷き詰める方法は, ちょうど 17 通りある.
長さ50では何通りの敷き詰め方があるか?
 
 seed(1):-
 	between(1,52,_).
 
 next_calc(_,_,[],[]):-!.
 next_calc(N,Count,[C1|Counts],[C1|Result]):-
 	N<4,
 	!,
 	N1 is N+1,
 	next_calc(N1,Count,Counts,Result).
 next_calc(N,Count,[C1|Counts],[C2|Result]):-
 	C2 is C1+(N-3)*Count,
 	N1 is N+1,
 	next_calc(N1,Count,Counts,Result).
 
 calc([_,_,_,Ans]):-!,write(Ans).
 calc([Count|Counts]):-
 	next_calc(1,Count,Counts,Counts1),
 	calc(Counts1).
 
 main114:-
 	findall(N,seed(N),Count),
 	calc(Count).




*Problem 115 「ブロックの組み合わせ方の数え上げ その2」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20115
リスト操作しかないのが逆転の発想を与えてくれた問題。
ユニットの長さを伸ばす操作を手前に一つ追加する操作だと考えればすっきりコードが書ける。
SF作家であり医療研究機関のプログラマでもあったグレッグイーガンが作中、セルオートマトンでセルの計算を行うにはセルオートマトンのサイズが巨大になるほど、近傍セルへのアクセス時間が延びるという記述があった。
多分彼は配列型言語ではなく、リスト操作型言語のプログラマだったんだろう。
配列ならアクセス速度は1のはずだからである。
まあそんなのはどうでもいいことなのでコード記述。
計算量はもう少し下げられる可能性はあるがまあ十分高速なのでいいかな。


 next_calc(_,[],0):-!.
 next_calc(N,[Count|Counts],Sum1):-
 	N=<0,
 	!,
 	next_calc(N,Counts,Sum),
 	Sum1 is Sum+Count.
 next_calc(N,[_|Counts],Sum):-
  	!,
 	N1 is N-1,
 	next_calc(N1,Counts,Sum).
 seed(1):-
 	between(1,51,_).
 calc(N,Sum,_):-1000000=<Sum,!,N1 is N-3,write(N1).
 calc(N,Sum,Counts):-
 	N1 is N+1,
  	next_calc(50,Counts,Sum1),
 	Sum2 is Sum+Sum1,
 	calc(N1,Sum2,[Sum2|Counts]).
 
 main115:-
 	findall(N,seed(N),Counts),
 	calc(52,1,Counts).





*Problem 116 「赤タイル, 緑タイル, あるいは青タイル」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20116
%風をひいて頭がうまく動かない カユ ウマ
とりあえず短いコードなので手っ取り早く書いてみた。


 seed(1):-
 	between(1,51,_).
 next_calc(_,_,[],[]):-!.
 next_calc(N,Add,[Count|Counts],[Count1|Result]):-
 	N=<1,
 	!,
 	Count1 is Count+Add,
 	next_calc(N,Add,Counts,Result).
 next_calc(N,Add,[Count|Counts],[Count|Result]):-
 	N1 is N-1,
 	next_calc(N1,Add,Counts,Result).
 
 calc(Len,Counts,Result):-length(Counts,Len),!,nth1(Len,Counts,Result).
 calc(Len,[Count|Counts],Result):-
 	next_calc(Len,Count,Counts,Counts1),
  	!,
 	calc(Len,Counts1,Result).
 
 main116:-
 	findall(Count,seed(Count),Counts),
  	calc(2,Counts,Re1),
 	calc(3,Counts,Re2),
 	calc(4,Counts,Re3),
 	Ans is Re1+Re2+Re3-3,
 	write(Ans).



*problem 117 「赤タイル, 緑タイル, そして青タイル」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20117
これも簡単。

 seed(1):-
 	between(1,51,_).
 next_calc(_,_,[],[]):-!.
 next_calc(N,Add,[Count|Counts],[Count1|Result]):-
  	N=<1,
 	!,
 	Count1 is Count+Add,
 	next_calc(N,Add,Counts,Result).
 next_calc(N,Add,[Count|Counts],[Count|Result]):-
 	N1 is N-1,
 	next_calc(N1,Add,Counts,Result).
 
 calc([_,Ans]):-!,write(Ans).
 calc([Count|Counts]):-
 	next_calc(2,Count,Counts,Counts1),
  	next_calc(3,Count,Counts1,Counts2),
 	next_calc(4,Count,Counts2,Counts3),
 	calc(Counts3).
 
 main117:-findall(N,seed(N),Counts),
 	calc(Counts).





*Problem 118 「パンデジタル素数集合」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20118
1 から 9 の全ての数字を使い, 自由につなげることで 10 進数の数字を作り, 複数の集合を作ることができる. 集合 {2,5,47,89,631} は面白いことに全ての要素が素数である.
1 から 9 の数字をちょうど 1 個ずつ含み, 素数の要素しか含まない集合はいくつあるか?

うーん?
全探索してみましたが速度が出ません。
もう少し賢い方法があるのでしょうか?

 not_prime(N):-N<2,!.
 not_prime(N):-
 	S is floor(sqrt(N)),
 	between(2,S,M),
 	0=:=N mod M,
 	!,
 	true.
 is_prime(N):-not(not_prime(N)).
 
 
 search_list(N,[],P,1):-!,
 	(P<N;N=:=0),
 	(N=:=0;is_prime(N)).
 search_list(N,Nums,P,Result):-
 	P<N,
 	is_prime(N),
 	search_list(0,Nums,N,Result).
 search_list(N,Nums,P,Result):-
 	select(N1,Nums,Nums1),
 	N2 is N*10+N1,
 	search_list(N2,Nums1,P,Result).
 
 main118:-findall(Add,search_list(0,[1,2,3,4,5,6,7,8,9],0,Add),Ans),
 	length(Ans,Len),write(Len).





*Problem 119 「数字和のべき乗」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20119
A^Bの桁数和がAになる数を探せばよい、それだけの問題。
適当に桁数上限を決め手抜きアセプト。
本当なら、最小値を確定しながら逐次的に答えを伸ばすべきなのだがそれをしないでも答えは出るのでいいかと思う。
それをするとコードがめんどくさいことになる。
それならいっそ大きな範囲を決めていっきに答えた方が楽。

 keta_sum(0,0):-!.
 keta_sum(N,Sum1):-N1 is N//10,keta_sum(N1,Sum),
 	Sum1 is Sum+(N mod 10).
 
 calc2(N,A,B):-
 	between(2,100,B),
 	N is A^B,
 	(N>10^20 ->(!,fail);true),
 	keta_sum(N,Sum),
 	Sum=:=A.
 
 calc([N,A,B]):-
 	between(2,180,A),
 	calc2(N,A,B).
 print_ans(31,_):-!.
 print_ans(N,[Ans|Ansers]):-write([N,Ans]),nl,
 	N1 is N+1,
 	print_ans(N1,Ansers).
 
 main119:-findall(Ans,calc(Ans),Ansers),
 	msort(Ansers,Ansers1),
 	print_ans(1,Ansers1).





*Problem 120 「自乗で割った余り」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20120
モード演算ですから計算結果は周期になります。
その時を捕まえれば一応全探索で答えが出ましたがそのコードは素朴すぎて却下。
もうちょっと頭のいい方法がありそうなので調べました。

http://www.mathblog.dk/project-euler-120-maximum-remainder/
とりあえず英語が読めないのですがこういうサイトくらいの解説は私のサイトも行うべきなのかもしれません。
リンク先の方法は2項定理で数式が打ち消し合ってその余りを取れば2anか2が残るところまでは理解。
ここで式は
2an=r mod a^2
となるので余りはaの倍数でしかありえない。
n=(a-1)/2の時余りは最大となるからこれが答えとなるようです。
というわけで今回はカンニングして解きました。
コードもシンプルなものです。
いかにもプロジェクトオイラーらしい答えですね。
これくらいなら良く考えれば自分でも思いつけたかも。


 search(Result):-
 	between(3,1000,A),
 	Result is 2*A*((A-1)//2).
 
 sum([],0):-!.
 sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
 
 main120:-findall(Ans,search(Ans),Ansers),
 	sum(Ansers,Ans2),write(Ans2).









*事実
会津大学オンラインジャッジでコード実行速度1位を取ったりすることもある私ですが。
兵庫県加古川市加古川町南備後在住の森元さんや藤村さんに言わせれば「AOJで一位を取るくらい馬鹿でもできる、ぷぷぷ笑い物」
だそうです、私より成績の悪い人も全員下等な見下された存在に分類されるそうです。
彼らは頭がいいんでしょうねきっと