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

prolog勉強プロジェクトオイラー31~40 - (2013/12/15 (日) 06:31:42) のソース

+プロジェクトオイラーの問題をPrologで解くページ。
作者 堀江 伸一
近況 問い36のコードの無駄を修正。
こんなネットの裏の裏の一番端っこの誰からも注目されないWikiに何を書いたって世の中に何の影響も与えないので結構好きかって書いてます。


*Problem 31 「硬貨の和」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2031
イギリスでは硬貨はポンド£とペンスpがあり,一般的に流通している硬貨は以下の8種類である.
1p, 2p, 5p, 10p, 20p, 50p, £1 (100p) and £2 (200p).
以下の方法で£2を作ることが可能である.
1×£1 + 1×50p + 2×20p + 1×5p + 1×2p + 3×1p
これらの硬貨を使って£2を作る方法は何通りあるか?

解法
Prologには破壊的代入がないので動的計画法が素直に使えません。
なのでちょっと苦労しましたがその結果まあまあの速度となりました。
今ならもうちょっと賢い方法で実装しますね


肝心のプログラムは以下の通り。



Hide Code

 fast_reverse(List1,List2):-fr(List1,[],List2).
 fr([Head|Tail],SoFar,Result):-fr(Tail,[Head|SoFar],Result).
 fr([],SoFar,SoFar).
 
 drop_w(Xs,Ys,Result):-length(Xs,Len1),length(Ys,Len2),Len3 is Len1-Len2,
    drop(Xs,Len3,Result).
 drop(Xs, 0, Xs).
 drop([], N, []) :- N > 0.
 drop([_ | Xs], N, Ys) :-
    N > 0, N1 is N - 1, drop(Xs, N1, Ys).
 
 
 add([],[],[]):-!.
 add([Z1|Zs1],[Y|Ys],[Z2|Sum]):-Z2 is Z1+Y,add(Zs1,Ys,Sum).
 
 setY([],Zs1,Ys1,Ans,_,_,Coins):-
        drop_w(Zs1,Ys1,Zs2),
        add(Zs2,Ys1,Zs3),
        !,
        NextZ=[Zs3|Ans],
        flatten(NextZ,NextZ1),
        fast_reverse(NextZ1,NextZ2),
        (Coins==[]->write(NextZ2);
        [Coin|Coins1]=Coins,
        first_set(NextZ2,[],Zs4,Coin,NextZ3),
        setY(NextZ3,Zs4,[],Zs4,0,Coin,Coins1)).
 
 setY(Xs,Zs1,Ys1,Ans,N,N,Coins):-!,
        add(Zs1,Ys1,NextZs),
        setY(Xs,NextZs,[],[NextZs|Ans],0,N,Coins).
 setY([X|Xs],Zs1,Ys1,Ans,Count,N,Coins):-!,
        Count1 is Count+1,
        setY(Xs,Zs1,[X|Ys1],Ans,Count1,N,Coins).
 
 first_set(List,Zs,Zs,0,List):-!.
 first_set([X|List],Zs,Result,N,Result2):-
        N1 is N-1,
        first_set(List,[X|Zs],Result,N1,Result2).
  
 main_31_3:-
        [C|Coins]=[1,2,5,10,20,50,100,200],
        findall(0,between(0,199,_),List),
        List1 = [1|List],
        first_set(List1,[],Zs,C,List2),
        setY(List2,Zs,[],Zs,0,C,Coins).




*Problem 32 「パンデジタル積」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2032
すべての桁に 1 から n が一度だけ使われている数をn桁の数がパンデジタル (pandigital) であるということにしよう: 例えば5桁の数 15234 は1から5のパンデジタルである.
7254 は面白い性質を持っている. 39 × 186 = 7254 と書け, 掛けられる数, 掛ける数, 積が1から9のパンデジタルとなる
掛けられる数/掛ける数/積が1から9のパンデジタルとなるような積の総和を求めよ.
HINT: いくつかの積は, 1通り以上の掛けられる数/掛ける数/積の組み合わせを持つが1回だけ数え上げよ.

解法
うーん?
ヒュースティリックに解いてしまった。
コードはむやみに長いし、自動的に解の探索範囲がもとまらず探索範囲を人手で指定しているし。
これっていいコードなんだろうか?

 sets([_],[_,_,_,_]).
 sets([_,_],[_,_,_]).
 
 selects([],Rest,Rest).
 selects([X|Rest],Nums,Result):-select(X,Nums,NumRest),selects(Rest,NumRest,Result).
 
 toNum([],C1,C1).
 toNum([X|Rest],C1,Result):-
 	C2 is C1*10+X,toNum(Rest,C2,Result).
 
 calc(Z):-
 	sets(Xs,Ys),
  	selects(Xs,[1,2,3,4,5,6,7,8,9],Rest1),
 	toNum(Xs,0,X),
 	selects(Ys,Rest1,Rest2),
 	toNum(Ys,0,Y),
 	Z is X*Y,
 	Z<10000,
 	Z1 is Z mod 10,
 	Z2 is (Z//10) mod 10,
 	Z3 is (Z//100) mod 10,
 	Z4 is (Z//1000) mod 10,
 	selects([Z1,Z2,Z3,Z4],Rest2,_).
 
 sum([],Sum,Sum).
 sum([X|Rest],Sum,Result):-Sum1 is Sum+X,sum(Rest,Sum1,Result).
 main32:-
 	findall(Z,calc(Z),List),sort(List,List1),
 	sum(List1,0,Ans),write(Ans).




*Problem 33 「桁消去分数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2033
49/98は面白い分数である.「分子と分母からそれぞれ9を取り除くと, 49/98 = 4/8 となり, 簡単な形にすることができる」と経験の浅い数学者が誤って思い込んでしまうかもしれないからである. (方法は正しくないが,49/98の場合にはたまたま正しい約分になってしまう.)
我々は 30/50 = 3/5 のようなタイプは自明な例だとする.
このような分数のうち, 1より小さく分子・分母がともに2桁の数になるような自明でないものは, 4個ある.
その4個の分数の積が約分された形で与えられたとき, 分母の値を答えよ.


解法
この問題コードは圧倒的に短くなるらしい。
とりあえず全部検討してみた。
search述語で重複が出なかったのでそのまま答えとした。


 gcd(0, B, B).
 gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G).
 calc(U,D,N,ReU,ReD):-ReU is U*10+N,ReD is D*10+N.
 calc(U,D,N,ReU,ReD):-ReU is U*10+N,ReD is D+N*10.
 calc(U,D,N,ReU,ReD):-ReU is U+N*10,ReD is D*10+N.
 calc(U,D,N,ReU,ReD):-ReU is U+N*10,ReD is D+N*10.
  
 search([U1,D1]):-
 	between(1,8,U),
 	T is U+1,
 	between(T,9,D),
 	between(1,9,N),
 	calc(U,D,N,U1,D1),
 	gcd(U1,D1,G1),
 	gcd(U,D,G),
 	U//G=:=U1//G1,
 	D//G=:=D1//G1.
 calc([],AnsU,AnsD):-gcd(AnsU,AnsD,G),Ans is AnsD//G,write(Ans).
 calc([[U,D]|Rest],AnsU,AnsD):-
 	AnsU1 is U*AnsU,
 	AnsD1 is D*AnsD,
 	calc(Rest,AnsU1,AnsD1).
 main33:-findall([U1,D1],search([U1,D1]),List),write(List),calc(List,1,1).





*Problem 34 「桁の階乗」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2034
145は面白い数である. 1! + 4! + 5! = 1 + 24 + 120 = 145となる.
各桁の数の階乗の和が自分自身と一致するような数の和を求めよ.
注: 1! = 1 と 2! = 2 は総和に含めてはならない.


問い30のコードをそのまま流用。
特に書くことなし。

 check([_|_],0):-!,fail.
 check([],0):-!.
 check(Nums,Sum):-
 	T is Sum mod 10,
 	select(T,Nums,Rest),!,
 	Sum1 is Sum//10,
 	check(Rest,Sum1).
 search(Nums,_,_,Sum,Sum):-
 	check(Nums,Sum).
 search(_,_,7,_,_):-!,fail. 
 
 search(Nums,P,Deep,Sum,Result):-
 	between(P,9,N),
 	nth0(N,[1,1,2,6,24,120,720,5040,40320,362880],N1),
 	Sum1  is Sum+N1,
 	Deep1 is Deep+1,
 	search([N|Nums],N,Deep1,Sum1,Result).
 
 sum([],Sum,Sum):-!.
 sum([X|Rest],Sum,Result):-
 	Sum1 is Sum+X,sum(Rest,Sum1,Result).
 
 main34:-findall(Re,search([],0,0,0,Re),List),
 	sum(List,0,Ans),Ans1 is Ans-3,write(Ans1).





*Problem 35 「巡回素数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2035
197は巡回素数と呼ばれる. 桁を回転させたときに得られる数 197, 971, 719 が全て素数だからである.
100未満には巡回素数が13個ある: 2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, および97である.
100万未満の巡回素数はいくつあるか?


計算量低減テクニックで書いたコードといっても。
どの桁も奇数であるというものだけを探索して検討するようにしただけ。
2は探索では検討せず最後に足している。
テクニックがみについたら無意味に使いたがる、そういうコードかもしれない。
テクニックで遊ぶのも技術習得には必要な気もするのでこれもありかと。
 
 
 
 round_list([],[],_).
 round_list([X|L],R,Result):-round_list(L,[X|R],Result).
 round_list([X|L],R,[X|Result]):-
 	reverse(R,R1),
 	append(L,R1,Result).
 
 not_prime(N):-N<2,!.
 not_prime(N):-
 	sqrt(N,N1),
 	N2 is floor(N1),
 	between(2,N2,N3),
 	N4 is N mod N3,
 	N4=:=0,!.
 
 to_num([],Result,Result):-!.
 to_num([X|List],Num,Result):-Num1 is Num*10+X,to_num(List,Num1,Result).
 
 check_exe(List,Ans):-
 	round_list(List,[],Perm),
 	to_num(Perm,0,Num),
 	(member(Num,Ans)->!,true;true),
 	(not_prime(Num)->
 	true,!;fail).
 
 result_list(List,_,_,List).
 result_list(List,Deep,Ans,Result):-
 	create_perm(List,Deep,Ans,Result).
 
 check_perm(List,Deep,Ans,Result):-
 	not(check_exe(List,Ans)),!,
 	to_num(List,0,Num),
 	result_list(List,Deep,[Num|Ans],Result).
 check_perm(List,Deep,Ans,Result):-
 	create_perm(List,Deep,Ans,Result).
 
 create_perm(_,6,_,_):-!,fail.
 create_perm(List,Deep,Ans,Result):-
 	Deep<6,
 	between(1,9,N),
 	T1 is N mod 2,
 	T1\==0,
 	Deep1 is Deep+1,
 	check_perm([N|List],Deep1,Ans,Result).
 main35:-findall(L,create_perm([],0,[],L),List),write(List),
 	length(List,Len),Len1 is Len+1,write(Len1).



*問い36 Double-base palindromes
http://projecteuler.net/problem=36
Problem 36 「二種類の基数による回文数」 †
585 = 10010010012 (2進) は10進でも2進でも回文数である.
100万未満で10進でも2進でも回文数になるような数の総和を求めよ.
(注: 先頭に0を含めて回文にすることは許されない.)

解法
Prolog的コードで書いたつもり。
前載せてたコードに恥ずかしい部分(実装上の無駄)があったので単純修正。

 

 set([_]).
 set([X1,X1]).
 set([X1,_,X1]).
 set([X2,X1,X1,X2]).
 set([X2,X1,_,X1,X2]).
 set([X2,X1,X0,X0,X1,X2]).
 
 num_perm([]).
 num_perm([X|Rest]):-member(X,[0,1,2,3,4,5,6,7,8,9]),num_perm(Rest).
 
 to_num([],Sum,Sum).
 to_num([X|Rest],Sum,Result):-Sum1 is Sum*10+X,to_num(Rest,Sum1,Result).
 
 to_bit2(0,[]):-!.
 to_bit2(Y,[B|Result]):-Y1 is Y//2,B is Y mod 2,to_bit2(Y1,Result).
 
 calc(X):-
 	set(A),
  	num_perm(A),
 	[Top|_]=A,
 	0<Top,
 	to_num(A,0,X),
 	to_bit2(X,Re),
 	reverse(Re,Re).
 sum([],Sum,Sum).
 sum([X|Rest],Sum,Result):-Sum1 is X+Sum,sum(Rest,Sum1,Result).
 
 main36:-
 	findall(Y,calc(Y),List),write(List),
 	sum(List,0,Ans),write(Ans).




*Problem 37 「切り詰め可能素数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2037
3797は面白い性質を持っている. まずそれ自身が素数であり, 左から右に桁を除いたときに全て素数になっている (3797, 797, 97, 7). 同様に右から左に桁を除いたときも全て素数である (3797, 379, 37, 3).
右から切り詰めても左から切り詰めても素数になるような素数は11個しかない. 総和を求めよ.
注: 2, 3, 5, 7を切り詰め可能な素数とは考えない.

解法
ちょっとコードが冗長かもしれない。
速度はそれなり。
短い方から左へ数字列を伸ばしてそれが条件を満たす限りのばす単純な探索で作成。
右側を消したとき素数でなくなる数でも左側へ無限に数字を加えたらどこかで素数になるものがあるかもしれないので。
素数が11個という情報がなかったらこの問題は恐ろしい難問となるはず。
探索は深さ制限で実行。


 not_prime(N):-N<2,!.
 not_prime(N):-
 	sqrt(N,N1),
 	N2 is floor(N1),
 	between(2,N2,N3),
 	N4 is N mod N3,
 	N4=:=0,!.
 is_prime(N):-not(not_prime(N)).
 
 to_num([],Sum,Sum).
 to_num([X|Rest],Sum,Result):-Sum1 is Sum*10+X,to_num(Rest,Sum1,Result).
 
 
 
 check_exe_R(0):-!.
 check_exe_R(Num):-
 	Num1 is Num//10,
 	is_prime(Num),
 	check_exe_R(Num1).
 
 search(List,List,_):-
 	to_num(List,0,Num),
 	10<Num,
 	check_exe_R(Num).
 
 search(List,Ans,Deep):-
 	Deep>0,
  	Deep1 is Deep-1,
 	member(N,[1,2,3,5,7,9]),
 	to_num([N|List],0,Num),
 	is_prime(Num),
 	search([N|List],Ans,Deep1).
 
 sum([],Sum,Sum).
 sum([X|Rest],Sum,Result):-to_num(X,0,Num),Sum1 is Sum+Num,sum(Rest,Sum1,Result).
 
 main37:-between(2,10,N),
 	findall(Ans,search([],Ans,N),List),length(List,11),
 	!,write(List),sum(List,0,Ans),write(Ans).





*Problem 38 「パンデジタル倍数」 †
192 に 1, 2, 3 を掛けてみよう.
192 × 1 = 192
192 × 2 = 384
192 × 3 = 576
積を連結することで1から9の パンデジタル数 192384576 が得られる. 192384576 を 192 と (1,2,3) の連結積と呼ぶ.
同じようにして, 9 を 1,2,3,4,5 と掛け連結することでパンデジタル数 918273645 が得られる. これは 9 と (1,2,3,4,5) との連結積である.
整数と (1,2,...,n) (n > 1) との連結積として得られる9桁のパンデジタル数の中で最大のものはいくつか?

解法
探索するだけの問題です。
最後に出てきた数字の並びが答えとなります。

 num_add(0,List,List):-!.
 num_add(Num,List,[X|Result]):-
 	X is Num mod 10,
 	Num1 is Num//10,
 	num_add(Num1,List,Result),
 	X>0,
 	not(member(X,Result)).
 
 
 search(_,_,List,_):-
 	length(List,Len),9<Len,!,fail.
 search(_,_,List,List1):-
 	length(List,9),reverse(List,List1),!.
 search(N,P,List,Result):-
 	N1 is N*P,
 	P1 is P+1,
 	num_add(N1,List,NextList),
 	search(N,P1,NextList,Result).
 search_start(Re):-
 	between(1,9999,N),
 	search(N,1,[],Re).
 
 main38:-
 	findall(Re,search_start(Re),List),
 	sort(List,List1),write(List1).






*Problem 39 「整数の直角三角形」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2039
辺の長さが {a,b,c} と整数の3つ組である直角三角形を考え, その周囲の長さを p とする. p = 120のときには3つの解が存在する:
{20,48,52}, {24,45,51}, {30,40,50}
p ≤ 1000 のとき解の数が最大になる p はいくつか?


解法
答えはめんどくさいので全部表示した。
一番最後に表示されているのが答えとなる。
このコードはWikiの
http://ja.wikipedia.org/wiki/%E3%83%94%E3%82%BF%E3%82%B4%E3%83%A9%E3%82%B9%E3%81%AE%E5%AE%9A%E7%90%86
を参考に記述した。



 %1000=>2*M*(M+N)
 %500/M=>M+N
 %500/M-M=>N
 gcd(0, B, G) :- G is abs(B).
 gcd(A, B, G) :- A =\= 0, R is B mod A, gcd(R, A, G).
 
 qsort([], []).
 qsort([X|L], S) :-
        mypartition(L, X, L1, L2),
        qsort(L1, S1),
        qsort(L2, S2),
        append(S1, [X|S2], S).
 
 mypartition([], _, [], []).
 mypartition([Y|L], X, [Y|L1], L2) :-
        Y < X,
        mypartition(L, X, L1, L2).
 mypartition([Y|L], X, L1, [Y|L2]) :-
        Y >= X,
        mypartition(L, X, L1, L2).
 set_MN([M,N]):-
 	between(2,500,M),
 	Up is floor(500/M-M)+1,
 	between(1,Up,N),
 	M>N,
 	1=:=(M-N) mod 2,
 	gcd(M,N,1),
 	L1 is 2*M*(M+N),
 	L1<1001.
 
 calc2(L,  _,[]):-L>1000,!.
 calc2(L,Add,[L|Result]):-
 	L1 is L+Add,
 	calc2(L1,Add,Result).
 
 calc([],[]):-!.
 calc([[M,N]|Rest],[List|Result]):-
 	L is 2*M*(M+N),
 	calc2(L,L,List),
 	calc(Rest,Result).
 
 count([],Count,L,[Count,L]):-!.
 count([L|Rest],Count,L,Result):-!,
 	Count1 is Count+1,
 	count(Rest,Count1,L,Result).
 count([L|Rest],Count,OldL,[[Count,OldL]|Result]):-
 	count(Rest,1,L,Result).
 
 
 main39:-
 	findall([M,N],set_MN([M,N]),List),
 	calc(List,Re1),
 	flatten(Re1,Re2),
 	qsort(Re2,Re3),
 	count(Re3,0,0,Result),write(Result),sort(Result,Ans),write(Ans).






*Problem 40 「チャンパーノウン定数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2040
正の整数を順に連結して得られる以下の10進の無理数を考える:
0.123456789101112131415161718192021...
小数第12位は1である.
dnで小数第n位の数を表す. d1 × d10 × d100 × d1000 × d10000 × d100000 × d1000000 を求めよ.


解法
割り算とユニフィケーションどちらが遅いのだろうか?
まあ値が小さいからあっという間に答えは出るけれど、数式的に解くとどうなんだろう?
と思って調べてみたら高速に解くということは数字列のn番目の要素を高速に計算する関数を作ることだったようです。
ま妥当な話ですね。
とりあえず素朴に計算する例。



 len(N,1):-N<10,!.
 len(N,2):-N<100,!.
 len(N,3):-N<1000,!.
 len(N,4):-N<10000,!.
 len(N,5):-N<100000,!.
 len(N,6):-N<1000000,!.
 len(N,7):-N<10000000,!. 
 
 no(0,N,Result):-!,Result is N mod 10.
 no(D,N,Result):-D1 is D-1,N1 is N//10,no(D1,N1,Result).
 
 calc([],_,_,Ans):-!,write(Ans).
 calc([D|Rest],N,Len,Ans):-
 	len(N,Add),
 	Len1 is Len+Add,
 	D=<Len1,!,
 	D2 is Len1-D,
 	no(D2,N,F),
 	Ans1 is Ans*F,
 	N1 is N+1,
 	calc(Rest,N1,Len1,Ans1).
 calc(Ds,N,Len,Ans):-
 	len(N,Add),
 	Len1 is Len+Add,
 	N1 is N+1,
 	calc(Ds,N1,Len1,Ans).
 main40:-
 	calc([1,10,100,1000,10000,100000,1000000],1,0,1).