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

prolog勉強プロジェクトオイラー101~110 - (2013/10/16 (水) 09:35:54) のソース

プロジェクトオイラーの問題をPrologで解く堀江伸一さんのページ。
学がないのでリンク先みたいなサイトを見るたびうらやましくなるような人間です。
http://d.hatena.ne.jp/Lost_dog/touch/searchdiary?word=*%5BHaskell%5D



*問い101「最適多項式」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20101
ニュートン補間を題材にした練習問題。
計算法や原理は分かったものの、Prologでどう実装したらよいかわからなかったのでとりあえずわかりやすさ重視でコードを実装。






 f(N,Re):-Re is 1-N+N^2-N^3+N^4-N^5+N^6-N^7+N^8-N^9+N^10.
 f2(N,Re):-Re is N^3.
 
 g(X,[X1],Result):-Result is X-X1,!.
 g(X,[X1|Xs],Result):-
 	!,
 	g(X,Xs,Re),
 	Result is Re*(X-X1).
 newton_calc(_,_,[C],C):-!.
 newton_calc(X,[_|Xs],[C|Cs],Result):-
 	newton_calc(X,Xs,Cs,ReY),
 	g(X,Xs,MultX),
 	Result is ReY+C*MultX.
 
 
 newton(_,[Y],[Y]):-!.
 newton([X|Xs],[Y|Ys],[C|Cs]):-
 	!,
 	newton(Xs,Ys,Cs),
 	newton_calc(X,Xs,Cs,ReY),
 	g(X,Xs,MultX),
 	C is (Y-ReY)/MultX.
 
 y_ calc(_,[_],[C],C):-!.
 y_calc(X,[_|Xs],[C|Cs],ResultY):-
 	y_calc(X,Xs,Cs,Re1),
 	g(X,Xs,MultX),
  	ResultY is Re1+MultX*C.
 
 
 
 
 calc(_,10,_,Sum):-write(Sum),!.
 calc(Ys,X,Xs,Sum):-
 	X1 is X+1,
  	f(X1,Y),
 	Ys1=[Y|Ys],
 	Xs1=[X1|Xs],
 	X2 is X1+1,
 	newton(Xs1,Ys1,Cs),
  	y_calc(X2,Xs1,Cs,Add),
 	Sum1 is Sum+Add,
 	write(Sum1),nl,
 	calc(Ys1,X1,Xs1,Sum1).
 
 main101:-calc([],0,[],0).



*問い102 Triangle containment
http://projecteuler.net/problem=102
テキストファイルが読みにくかったのでテキストの内容をエディタでPrologのリスト形式に加工してから解。
解き方は自由のはずなのでこれはルール違反ではないはず。
原点Oが三角形ABCの内部にあるとはOAとOBの内積、OBとOCの内積、OCとOAの内積が全部同じ符号であるということ。
辺上にある場合の処理が問題に書いてないが、その場合も三角形内部にあるものと判定して計算でも合格できました。



 calc([],Ans):-!,write(Ans).
 calc([X1,Y1,X2,Y2,X3,Y3|Rest],Ans):-
 	I1 is X1*Y2-X2*Y1,
 	I2 is X2*Y3-X3*Y2,
 	I3 is X3*Y1-X1*Y3,
  	((0=<I1,0=<I2,0=<I3);(I1=<0,I2=<0,I3=<0)),
 	!,
 	Ans1 is Ans+1,
 	calc(Rest,Ans1).
 calc([_,_,_,_,_,_|Rest],Ans):-
 	calc(Rest,Ans).
 
 main102:-
 	seen,
 	see('e102.txt'),
 	read(X),
 	seen,
 	calc(X,0).




*問い104 「両端がパンデジタルなフィボナッチ数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20104
先頭桁9桁に影響を与えるのは大きな桁のほうだけなので先頭桁を計算するのに必要な桁数だけ先頭を計算すればC++のint64でも計算が出来ます。
Prologではそれをする意味があるのか考えてもよくわからなかったのですが、実装がわかりやすいという理由でこの必要なケタだけ計算を採用しました。


 cut_div(N,R,N,1):-Up is 10^R,N<Up,!.
 cut_div(N,R,Result,D1):-N1 is N//10,cut_div(N1,R,Result,D),D1 is D*10.
 check(0,[]):-!.
 check(N,Nums):-
 	N1 is N mod 10,
 	select(N1,Nums,Nums1),
 	N2 is N //10,
 	check(N2,Nums1).
 
 fina(N,Top,_,Tail,_):-
 	cut_div(Top,9,Top1,_),
 	check(Top1,[1,2,3,4,5,6,7,8,9]),
 	check(Tail,[1,2,3,4,5,6,7,8,9]),
 	!,write([N,Top1,Tail]).
 fina(N,Top,Top1,Tail,Tail1):-
 	Top2 is Top+Top1,
 	cut_div(Top2,18,TopB,D),
         TopA is Top1//D,
 	Tail2 is (Tail+Tail1) mod 10^9,
 	N1 is N+1,
 	!,
 	fina(N1,TopA,TopB,Tail1,Tail2).



*問い105「特殊和集合:検査」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20105
とりあえずいい方法は実装がめんどくさそうだったので手軽に全探索してみた。
f(全探索,馬鹿).
test:-f(全探索,X),write(X).
そんな感じ?
まあデータサイズが小さいから答えはすぐに出たけれど、もう少し高速化したいところ。
それと例のごとくデータファイルをProlog形式に加工しております。


 sum([],0):-!.
 sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
 
 subset_sum(0,As,0,As):-!.
 subset_sum(N,[A|As],BSum,[A|ReAs]):-
 	subset_sum(N,As,BSum,ReAs).
 subset_sum(N,[B|As],BSum1,ReAs):-
 	N1 is N-1,
 	!,
 	subset_sum(N1,As,BSum,ReAs),
 	BSum1 is BSum+B.
 
 is_bad2(_,_,BSum,BSum):-!.
 is_bad2(N1,N2,BSum,CSum):-N1<N2,CSum=<BSum,!.
 is_bad2(_,_,_,_):-fail.
 
 is_bad(Data):-
 	length(Data,Len),
 	HerfLen is Len//2,
 	between(1,HerfLen,N1),
 	subset_sum(N1,Data,BSum,Data2),
  	RestLen is Len-N1,
 	between(N1,RestLen,N2),
 	subset_sum(N2,Data2,CSum,_),
 	is_bad2(N1,N2,BSum,CSum),
 	!,
 	true.
 
 
 check([],Ans):-write(Ans).
 check([Data|Datas],Ans):-
 	write(Data),nl,
 	not(is_bad(Data)),
 	write(ok),
 	!,
 	sum(Data,Add),
 	Ans1 is Ans+Add,
 	check(Datas,Ans1).
 check([_|Datas],Ans):-!,check(Datas,Ans).
 
 main105:-
 	seen,
 	see('e105.txt'),
 	read(Datas),
 	check(Datas,0).



*問い106 「特殊和集合:メタ検査」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20106
特殊和集合でn=12の時検査すべき部分和集合の数を答える問題。
特殊和集合の数列に番号をつけて手前からm個とm個を拾いbi,ciとし、その数列が
全てのbiとciでbi<ciかbi>ciを満たしているならそれはチェックする必要がない。
それ以外はすべてチェックする必要がある。
コードを実行して出てきた答えを2で割ったのが正式な答え。

 check_up([],[]):-!.
 check_up([B|Bs],[C|Cs]):-
 	B<C,check_up(Bs,Cs).
 
 set1(0,As,[],As):-!.
 set1(N,[A|As],Bs,[A|ReAs]):-
 	set1(N,As,Bs,ReAs).
 set1(N,[B|As],[B|Bs],ReAs):-
 	N1 is N-1,
 	set1(N1,As,Bs,ReAs).
 test(1):-
 	between(2,6,N),
 	set1(N,[1,2,3,4,5,6,7,8,9,10,11,12],Bs,As1),
 	set1(N,As1,Cs,_),
        not(check_up(Bs,Cs)),
 	not(check_up(Cs,Bs)).
 main106:-findall(C,test(C),Cs),
 	length(Cs,Len),
 	write(Len).



*Problem 107 「最短ネットワーク」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20107
この問題 破壊的代入なしでコードを書くのはかなりめんどくさい。
C++なら3重ループで済むところを素直に直訳すると三種類の再帰の組み合わせに翻訳される。
なので一応楽チンに記述できる形でかいてみた。
グラフのサイズが小さいのでこれでも問題なし。


 roop(Points,Cons,[P,Cost]):-
 	member(Point,Points),
 	nth0(Point,Cons,Con),
 	between(0,39,P),
 	nth0(P,Con,Cost),
 	integer(Cost),
 	not(member(P,Points)).
 
 
 min_point(_,Cost1,Point2,Cost2,Point2,Cost2):-Cost2<Cost1,!.
 min_point(Point1,Cost1,_,_,Point1,Cost1):-!. 
 
 select_min([PointCost],Point,Cost):-!,[Point,Cost]=PointCost.
 select_min([[Point,Cost]|PointCost],Point2,Cost2):-
 	!,
 	select_min(PointCost,Point1,Cost1),
 	min_point(Point,Cost,Point1,Cost1,Point2,Cost2).
 
 
 search(_,[],_,Ans):-!,write(Ans).
 search(Points,RestPoints,Cons,Ans):-
 	findall(PointCost,roop(Points,Cons,PointCost),PointCosts),
 	select_min(PointCosts,Point,Cost),
 	Ans1 is Ans-Cost,
 	select(Point,RestPoints,RestPoints1),
 	search([Point|Points],RestPoints1,Cons,Ans1).
 sum([],0):-!.
 sum([X|Xs],Sum1):-integer(X),!,sum(Xs,Sum),Sum1 is Sum+X.
 sum([_|Xs],Sum):-sum(Xs,Sum).
 
 all_sum(Cons,Sum):-
 	flatten(Cons,Cons1),
 	sum(Cons1,Sum).
 
 main107:-seen,
 	see('e107.txt'),
 	read(Cons),
         findall(N,between(1,39,N),RestPoints),
 	all_sum(Cons,Ans),
 	Ans1 is Ans//2,
 	search([0],RestPoints,Cons,Ans1).



*Problem 109 「ダーツ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20109
ダブルで終わるという条件を満たすために得点を最後から計算していきます。
後は探索空間が狭いので全探索でフィニッシュです。
動的計画法やメモ化計算を使うまでもありません。

 ps_sdt(0,0):-!.
 ps_sdt(25,P):-member(P,[25,50]).
 ps_sdt(P,P1):-0<P,P<21,member(M,[1,2,3]),P1 is M*P.
 ps_d(P,P1):-((0<P,P<21);P=:=25),P1 is P*2.
 
 points_create(P):-
 	between(1,25,N),
 	ps_sdt(N,P).
 
 search(_,[],_,_,_):-!,fail.
 search(_,_,Score,_,_):-
 	100=<Score,!,fail.
 search(4,_,_,_,_):-
 	!,fail.
 search(_,_,Score,Score,1).
 search(Turn,[Point|Points],Score,Result,_):-
 	Turn1 is Turn+1,
 	Score1 is Score+Point,
 	search(Turn1,[Point|Points],Score1,Result,1).
 search(Turn,[Point|Points],Score,Result,_):-
 	Check is Score+Point,
 	Check<100,
 	search(Turn,Points,Score,Result,0).
 
 search_top(Points,LastScore):-
 	between(1,25,N),
 	ps_d(N,Score),
 	search(1,Points,Score,LastScore,1).
 
 main109:-
 	findall(P,points_create(P),Points),
 	msort(Points,Points1),
 	findall(LastScore,search_top(Points1,LastScore),LastScores),
 	length(LastScores,Ans),write(Ans).



*問い108 110「ディオファントス逆数 その2」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20110
Piを素数、Aiを自然数とし。
n=P1^A1*P2*A2*,,*Pm*Am でPiの値にかかわらずAiの値と個数だけから答えが決まる。
ここまでは小さな値での規則性からみてわかったものの、なぜそうなるのかよくわからなかったのでYahoo知恵袋で質問
 n(x+y)=xy 
 ∴(x-n)(y-n)=n^2
 (x-n,y-n)=(1,n^2),(a1,n^2/a1),....,(n^2,1)
というヒントを教えてもらう。
このヒントからnが決まると
int((2A1+1)(2A2+1),,,(2Am+1)/2)+1個になる。
と理解。
とりあえず速度より答えを出てきたらいいかな程度でコードを記述。
遅かったらコード高速化だったけど手元のパソコンで5秒で答えが出たのでまあいいかともおもう。

 multi([],1):-!.
 multi([X|Xs],Result):-multi(Xs,Re),Result is Re*X.
 
 search(_,[],_,_,_,_):-!,fail.
 search(Up,_,Num,_,_,_):-
 	Up=<Num,!,fail.
 search(_,_,Num,Multi,A,Num):-
 	Multi1 is (Multi*(2*A+1))//2+1,
 	4000000=<Multi1,!.
 search(Up,[P|Ps],Num,Multi,A,Result):-
 	Num1 is Num*P,
 	A1 is A+1,
 	search(Up,[P|Ps],Num1,Multi,A1,Result).
 search(Up,[_,P|Ps],Num,Multi,A,Result):-
 	Multi1 is Multi*(2*A+1),
 	Num1 is Num*P,
 	search(Up,[P|Ps],Num1,Multi1,1,Result).
 
 main110:-
 	Primes=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47],
 	multi(Primes,Up),
 	findall(Re,search(Up,Primes,2,1,1,Re),Ans),
 	length(Ans,Len),write(Len),nl,
 	sort(Ans,Ans1),
 	[Ans2|_]=Ans1,write(Ans2).