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

prolog勉強プロジェクトオイラー61~70 - (2013/08/13 (火) 14:37:33) の編集履歴(バックアップ)


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

問い61

この問題は探索でといても十分間に合いますが。
スタック操作だけで動的計画法で解くというチャレンジをしてみました。
どの部分の処理も線形時間か n log(n)で終わるようになっています。
appendを使ってるところだけが効率が悪いくらいです。


calcP3(N1):-between(1,200,N),N1 is (N*(N+1))//2,999<N1,N1<10000.
calcP4(N1):-between(1,100,N),N1 is N*N,999<N1,N1<10000.
calcP5(N1):-between(1,100,N),N1 is (N*(3*N-1))//2,999<N1,N1<10000.
calcP6(N1):-between(1,100,N),N1 is (N*(2*N-1)),999<N1,N1<10000.
calcP7(N1):-between(1,100,N),N1 is (N*(5*N-3))//2,999<N1,N1<10000.
calcP8(N1):-between(1,100,N),N1 is (N*(3*N-2)),999<N1,N1<10000.



mysweep([],[]):-!.
mysweep([[NNum,FNum,Sum]],[]):-integer(NNum),integer(FNum),integer(Sum),!.
mysweep([[NNum,FNum,_],[NNum,FNum,Sum]|Rest],Result):-!,
	mysweep([[NNum,FNum,Sum]|Rest],Result).
mysweep([_,Set|Rest],[Set|Result]):-
	!,
	mysweep([Set|Rest],Result).


calc2([NN,FN,Sum],[[NN,NB]|_],[NB,FN,Sum1]):-
	Sum1 is Sum+NN*100+NB.
calc2(Set,[_|Rest],Result):-
calc2(Set,Rest,Result). 

calc([],_,Ans,Perm):-!,search(Ans,Perm).
calc(_,[],Ans,Perm):-!,search(Ans,Perm).
calc([[NN,FN,Sum]|Rest],[[NN,NB]|Rest1],Ans,Perm):-
	!,
	findall(Re1,calc2([NN,FN,Sum],[[NN,NB]|Rest1],Re1),List),
	append(Ans,List,Ans1),
	calc(Rest,[[NN,NB]|Rest1],Ans1,Perm).
calc([[NN,_,_]|Rest],[[NF,NB]|Rest1],Ans,Perm):-
	NN<NF,!,
	calc(Rest,[[NF,NB]|Rest1],Ans,Perm).
calc(List,[_|Rest],Ans,Perm):-
	!,
	calc(List,Rest,Ans,Perm).

check([[A,A,Ans]|_],Ans):-!.
check([_|Rest],Ans):-
	check(Rest,Ans).

search([],[]):-!,fail.
search(Ans,[]):-check(Ans,Ans1),write([ansis,Ans1]).
search(Ans,Perm):-
	sort(Ans,Ans1),
	[Top|_]=Ans1,
	mysweep(Ans1,Re1),
	Ans2=[Top|Re1],
	select(L,Perm,Rest1),
	calc(Ans2,L,[],Rest1). 

to_date3([],[]):-!.
to_date3([Y|Rest],[[N2,N1,Y]|Result]):-
	N1 is Y // 100,
	N2 is Y mod 100,
	to_date3(Rest,Result).

to_date2([],[]):-!.
to_date2([Y|Rest],[[N2,N1]|Result]):-
	N2 is Y // 100,
	N1 is Y mod 100,
	to_date2(Rest,Result).

to_date([],[]):-!.
to_date(List,Re2):-
	member(L,List),
	to_date2(L,Re1),
	sort(Re1,Re2).

main61:-
	findall(Y,calcP3(Y),L3),
 	findall(Y,calcP4(Y),L4),
	findall(Y,calcP5(Y),L5),
	findall(Y,calcP6(Y),L6),
	findall(Y,calcP7(Y),L7),
	findall(Y,calcP8(Y),L8),
	findall(Re,to_date([L3,L4,L5,L6,L7],Re),Perm),
	to_date3(L8,L88),
	search(L88,Perm).




問い62

Nを小さい数字から始めN^3の各桁の数字の出現数でカウントしていき大きな数字へ続けます。
条件を満たせばそれを暫定解として探索していき暫定解*10<N^3となればそれより小さい答えはありません。
c++のstd::mapにあたる機能がピュアPrologにはないので線形探索をしておりコード実行速度は遅いです。
std::mapに当たる機能があったら非常に高速に解ける問題です。



mycount(0,Temp,Result):-sort(Temp,Result),!.
mycount(N,Temp,Result):-
	N1 is N//10,
	M is N mod 10,
	(select([M,C],Temp,Rest) ->
	C1 is C+1,
	mycount(N1,[[M,C1]|Rest],Result);
	mycount(N1,[[M,1]|Temp],Result)).

search(N,Ans,_):-
	0<Ans,
	T is Ans*10,
	N3 is N*N*N,
	T<N3,
	!,
	write(Ans).


search(N,Ans,Counts):-
	N3 is N*N*N,
	N1 is N+1,
	mycount(N3,[],Re1),
	select([Re1,C,Min],Counts,Rest),!,
	(4 =< C ->
	((Ans>Min;Ans=:=0) -> search(N1,Min,Counts);
	search(N1,Ans,Counts));
	C1 is C+1,
	search(N1,Ans,[[Re1,C1,Min]|Rest])
	).
search(N,Ans,Counts):-
	N3 is N*N*N,
	N1 is N+1,
	mycount(N3,[],Re1),
	search(N1,Ans,[[Re1,1,N3]|Counts]).





問い63

この問題は10より大きな数を考えると10^2>2桁なのでそれより大きな数の乗数は考慮する必要はありません。
1~9まで検証すれば十分です。
Logを使った方がすっきりしますが桁数を簡単に数えてときました。

keta_size(0,0):-!.
keta_size(N,Result):-
	N1 is N//10,
	keta_size(N1,Re),
	Result is Re+1.

search2(A,B,N):-N is A^B,keta_size(N,L),B>L,!,fail.
search2(A,B,N):-N is A^B,keta_size(N,L),B=:=L.
search2(A,B,N):-B1 is B+1,search2(A,B1,N).

search(N):-
	between(1,9,A),
	search2(A,1,N).
main63:-findall(N,search(N),List),sort(List,List1),length(List1,Len),write(Len).




問い64

10000以下の平方数でない数を平方根を連分数展開で計算する。
これを連分数展開した時、連分数展開が奇数周期になる物の個数を答えよ。


gcd(0, B, B) :-!.
gcd(A, B, G) :- A \== 0, R is B mod A, gcd(R, A, G).

calc(N,NSqrt,U,DR,Deep,List,Result):-
	D is N-DR*DR,
	gcd(D,U,G),
	D1 is D//G,
	U1 is U//G,
	L1 is floor((NSqrt+DR)*U1/D1),
	L2 is L1*D1,
	UR2 is L2-DR,
	Deep1 is Deep+1,
	(member([L1,UR2,Deep2,D1],List) ->
	!,Result is (Deep1-Deep2) mod 2
	;calc(N,NSqrt,D1,UR2,Deep1,[[L1,UR2,Deep1,D1]|List],Result)).

search(10001,Ans):-!,write(Ans).
search(N,Ans):-
	sqrt(N,N1),
	N=:= floor(N1)*floor(N1),
	!,
	N2 is N+1,
	search(N2,Ans).
search(N,Ans):-
	sqrt(N,NSqrt),
	FloorSQ is floor(NSqrt),
	calc(N,NSqrt,1,FloorSQ,0,[],Result),
	Ans1 is Ans+Result,
	N1 is N+1,
	search(N1,Ans1).



問い65

素朴に定義通り求めてみました。


gcd(0, B, B) :-!.
gcd(A, B, G) :- A \== 0, R is B mod A, gcd(R, A, G).
add(U,D,U1,D1,ReU,ReD):-
	U2 is U*D1+U1*D,
	D2 is D*D1,
	gcd(U2,D2,G),
	ReU is U2//G,
	ReD is D2//G.


calc([X],U,D):-
	integer(X),
	U is 1,
	D is X.
calc([X|Rest],U1,D1):-
	calc(Rest,U,D),
	add(U,D,X,1,D1,U1).

calc_list(M):-
	between(0,98,N),
	(1=:=N mod 3->
	M is N //3*2+2;
	M is 1).
keta_sum(0,0):-!.
keta_sum(N,Result):-
	N1 is N//10,
	keta_sum(N1,Re),
	Result is Re+N mod 10.

main65:-
	findall(M,calc_list(M),List),calc(List,U,D),add(U,D,2,1,U1,_),
	keta_sum(U1,Ans),write(Ans).





問い66 http://projecteuler.net/problem=66

この問題は色々式変形したり漸化式が立たないか考えたけど自力では解けなかった問題。
ネットに転がっていたjavascriptコードをPrologコードに翻訳して100%カンニングして解いた。
驚くほど高速に答えが出てすこし驚いた。
この手法は連分数展開というものらしい。
大雑把な考え方までは理解したけれど、細かい計算は理解してない状態。
D=(x^2-1)/y^2はyが大きくなると
x^2/y^2と見分けがつかなくなってくる。
√D=x/yを満たす分数を求める問題で近似するという発想らしい。
聞けばなるほどの着想だが厳密解といいう考え方にこだわってると中々出てこない発想だな。
一つ勉強になった。






calc(_,_,_,_,E,F,K,N,_,P2,_,_,X,Result):-
	B1 is (X-E*E)/F,
	C1 is floor((K+E)/B1),
	A1 is C1,
	A1=:=2*K,0=:=N mod 2,!,Result is P2.
 
calc(_,_,_,_,E,F,K,N,P1,P2,Q1,Q2,X,Result):-
	B1 is (X-E*E)/F,
	C1 is floor((K+E)/B1),
	D1 is C1*B1-E,
	A1 is C1,
	E1 is D1,
	F1 is B1,
	P3 is A1*P2+P1,
	Q3 is A1*Q2+Q1,
	N1 is N+1,
	calc(A1,B1,C1,D1,E1,F1,K,N1,P2,P3,Q2,Q3,X,Result).


calc_first([Result,D]):-
	between(2,1000,D),
	sqrt(D,T),
 	K is floor(T),
	K2 is K*K,
	D\==K2,
	A is K,
	E is K,
	F is 1,
	N is 1,
	P1 is 1,
	P2 is K,
	Q1 is 0,
	Q2 is 1,
	calc(A,_,_,_,E,F,K,N,P1,P2,Q1,Q2,D,Result).
main66_2:-
	findall([D,X],calc_first([D,X]),List),
	sort(List,List1),write(List1).






問い67

問題自体は簡単。
ただPrologには破壊的代入がないせいで少しめんどくさいので只今効率的なコードとは何か熟慮中。
とりあえずProlog風味なファイル読み込みとは何かと考えてコーディングをしてみました。
データにゴミや入力ミスが入ってないから練習問題って楽でいいですね。
次の行の計算が少し混乱してるがまあ正しい答えが出てるからいいかな?


one_read(Result,State):-get0(C),((C =:= -1 ; C=:=10;C=:=32)->Result=[],State is C;
			  one_read(Re,State),Result=[C|Re]).
row_read([Re1|Result],ReState):-
	one_read(Re,State),(Re==[]->Re1=[];to_num(Re,0,Re1)),((State=:= -1;State=:=10)->ReState is State,Result=[];
			       row_read(Result,ReState)).
all_read([Result|Re]):-row_read(Result,State),(State=:= -1 ->Re=[];
					 all_read(Re)).
to_num([],Sum,Sum):-!.
to_num([X|Rest],Sum,Result):-
	Sum1 is Sum*10+(X-48),to_num(Rest,Sum1,Result).

calc_row(_,[],Re,Result):-!,reverse(Re,Result).
calc_row([L,R|Rest],[Z|Sums],[RS|Rest1],Result):-
	!,
	L1 is L+Z,
	R1 is R+Z,
	(L1<RS -> R2 is RS;R2 is L1),
	calc_row([R|Rest],Sums,[R1,R2|Rest1],Result).


calc_row_first([L,R|Rest],[Z|Sums],Result):-
	L1 is L+Z,
	R1 is R+Z,
	calc_row([R|Rest],Sums,[R1,L1],Result).

 
calc([[[]]],Sums):-!,write(a),sort(Sums,Sums1),write(Sums1).
calc([Row|Rest],Sums):-
       calc_row_first(Row,Sums,Sums1),
       calc(Rest,Sums1).

main67:-seen,see('e67.txt'),all_read(Re),
	[Sums|Rest]=Re,
	calc(Rest,Sums).





問い68

こういう探索系の問題はPrologは強いですね。
いかにもPrologらしいという感じで処理をかけた気がします。

search([],[],[],Sum,Sum,Ans):-!,write(Ans).
search([],[L|Perm],Nums,Sum,Sum,Ans):-!,
       [Xi|L1]=L,
	[X1|_]=Ans,
	select(Xi,Nums,NumsRest),
	(X1=:=10;X1<Xi),
	search(L1,Perm,NumsRest,Xi,Sum,Ans).
search([Xi|Rest],Perm,Nums,NowSum,Sum,Ans):-
	integer(Xi),
	!,
	Xi<10,
	NowSum1 is NowSum+Xi,
	search(Rest,Perm,Nums,NowSum1,Sum,Ans).
search([Xi|Rest],Perm,Nums,NowSum,Sum,Ans):-
	select(Xi,Nums,RestNum),
	NowSum1 is NowSum+Xi,
	search(Rest,Perm,RestNum,NowSum1,Sum,Ans).

main68:-
	Ans=[X1,X2,X3, X4,X3,X5, X6,X5,X7, X8,X7,X9, X10,X9,X2],
	Perm=[[X4,X3,X5],[X6,X5,X7],[X8,X7,X9],[X10,X9,X2]],
	select(X1,[9,8,7,6,5,4,3,2,1,10],Rest),
	select(X2,Rest,Rest1),
	select(X3,Rest1,Rest2),
	X2<10,
	X3<10,
	Sum is X1+X2+X3,
	search([],Perm,Rest2,Sum,Sum,Ans).




問い69

この問題はわざわざプログラムを書くまでもありません。
φ(n)=n*(1-1/p1)*(1-1/p2),,,*(1-1/pn)なのですから
A=n/φ(n)=1/(1-1/p1)*(1-1/p2),,,*(1-1/pn)です。
Aが最大化されるのは右辺の分母が最小になるとき。
素数を小さい方からかけていき、2,3,5,7,11,13、、、とかけていき100万を超える手前が答えとなります。



問い70

この問題、正答者掲示板を見ましたが2数の素数積が答えとなると英語で書かれていたのですが理由がよくわかりませんでした。
例えば3素因数で条件を満たす数があったとして、このサブセットの積が条件を満たすとは限りません。
2素因数でもp1^n*p2^mでn、mが1より大きい時だけ条件を満たす場合もあり得ます。
私には2素因数だけを検討すれば良い理由がわからなかったので、全部検証することにしました。