カテゴリー: PROLOG

問題

最初のほうのはみんな簡単ですね…歯ごたえあるやつまでとばしてやるかな…

プログラム：

% problem 6	
solve:-
	make_list(100,L),
	calc_pow2_sum(L,PowSum),
	calc_sum(L,Sum),
	SumPow is Sum * Sum,
	Ans is SumPow - PowSum ,
	write(Ans).

calc_pow2_sum([],0):-!.
calc_pow2_sum([First|Rest],PowSum):-
	FirstPow is First * First,
	calc_pow2_sum(Rest,PowSumRest),
	PowSum is FirstPow + PowSumRest.
	
calc_sum([],0):-!.
calc_sum([First|Rest],Sum):-
	calc_sum(Rest,RestSum),
	Sum is First + RestSum.
	
% 1 ～ Idx までのリストを作成
make_list(Idx,List):-
	make_list_sub(Idx,List,[]).
	
make_list_sub(0,Ret,Ret):-!.
make_list_sub(Idx,Ret,List):-
	Idx1 is Idx - 1,
	make_list_sub(Idx1,Ret,[Idx | List]).

実行結果：
1 ?- solve.
***********************(解答伏せます)**********************
true.

問題

1～20の最小公倍数を求めているだけ

プログラム：

% problem 5
	
solve:-
	calc_lcm_all([2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20],Lcm),
	write(Lcm).

% すべての数字の最小公倍数を取得
calc_lcm_all([],1):-!.
calc_lcm_all([First|Rest],Lcm):-
	calc_lcm_all(Rest,LcmRest),
	lcm(First,LcmRest,Lcm).

% 最小公倍数
lcm(N, M, Lcm) :-
	Gcd is gcd(N , M),
	Lcm is (N * M) // Gcd, !.

実行結果：

[2] 9 ?- time(solve).
***********************(解答伏せます)**********************
% 79 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
true.

問題

assertのハッシュとselectの順番に頼った力まかせのひどいプログラム。エレガントさのかけらもない。遅いし…

make_pent_list(10000,PentLst)の数字も小さいものから試していって見つかるまでだんだんと増やしてった（ひどい）

答えは合ってた。


:-use_module(library(clpfd)).

main:-
make_pent_list(10000,PentLst),
select(B,PentLst,Rest),
select(A,Rest,_),
A < B, Sum is A+B, p(Sum), Diff is B-A, p(Diff), write([Diff]). make_pent_list(0,[]):-!. make_pent_list(Idx,[First|Rest]):- !, First is Idx * (3*Idx - 1) / 2, assert(p(First)), Idx1 is Idx - 1, make_pent_list(Idx1,Rest). %　実行結果 % %1 ?- time(once(main)). %[5482660] %% 252,016,073 inferences, 295.203 CPU in 296.719 seconds (99% CPU, 853704 Lips) %true.

問題

問題を「3桁の数3つの最大回文積を求めよ」と読み間違っていてずっと解答が合わずに頭を抱えていた。勘違いしていた問題のほうが本題より難しいがこちらの答えは以下のとおり967262769になるようだ。
1 ?- time(solve).
967262769=999*989*979
% 30,373,524 inferences, 16.703 CPU in 16.734 seconds (100% CPU, 1818434 Lips)
true
between述語のデクリメントバージョンdetween_decを作成して対応した。

プログラム：

% problem 4
	
solve:-
	between_dec(999,100,Left3Dig),
	number_chars(Left3Dig,NumChars),
	[D1,D2,D3] = NumChars,
	[D1,D2,D3,D3,D2,D1] = KaibunChars,
	number_chars(Kaibun,KaibunChars),
	
	between_dec(999,100,A),
	Kaibun mod A =:= 0,
	
	B is Kaibun / A,
	
	100 =< B, B =<  999,
	!,
	write(Kaibun),write('='),write(A),write('*'),write(B),nl.

% decrimental version of between
between_dec(Max,Min,Val):-
	integer(Val),
	!,
	Max >= Val,
	Val >= Min.
between_dec(Max,Min,Val):-
	var(Val),
	between_dec_sub(Max,Min,Max,Val).

between_dec_sub(_,Min,Val,_):-
	Val<Min,
	!,
	fail.
between_dec_sub(Max,Min,Val,Val):-
	Min =< Val,Val =< Max.
between_dec_sub(Max,Min,Val,RetVal):-
	Val1 is Val - 1,
	between_dec_sub(Max,Min,Val1,RetVal).
	
	

実行結果：

1 ?- time(solve).
回答伏せます
% 420,104 inferences, 0.219 CPU in 0.219 seconds (100% CPU, 1920475 Lips)
true.

問題

解説：
①エラトステネスのふるいによる素数列挙、すべてsosu述語としてassert
②rotate述語で転回した数字をsosu述語で存在するか確認
　すべての転回に関して存在すればretract
　見つかった分カウントをアップ

結構速く動作してると思います（自分の基準では）

% problem 35

% Idx は 9 のみで構成されている必要あり
solve(Idx):-
	retractall(sosu),
	flag(junkan_sosu_cnt, _, 0),
	assert(max(Idx)),
	make_list(Idx,List),
	eratosthenes_sieve(List,List1), % エラトステネスのふるいによる素数列挙
	assert_all(List1),
	rot_chks(List1),
	!,
	flag(junkan_sosu_cnt, Cnt1, Cnt1),nl,
	write(Cnt1).

rot_chks([]):-!.
rot_chks([First|Rest]):-
	rot_chk(First),
	rot_chks(Rest).
rot_chks([_|Rest]):-
	!,
	rot_chks(Rest).
	
rot_chk(Num):-
	!,
	findall(Rot,rotate(Num,Rot),RotLst),
	sort(RotLst,RotLst1),  % delete duplicate numbers like [11,11]
	sosu_assert_chk(RotLst1),
	write(RotLst1),nl,
	length(RotLst1,Len),
	flag(junkan_sosu_cnt, Cnt, Cnt + Len),
	retract_sosu_list(RotLst1).

% リスト内がすべて素数としてassertされているかかェック
sosu_assert_chk([]):-!.
sosu_assert_chk([First|Rest]):-
	sosu(First),
	sosu_assert_chk(Rest).

retract_sosu_list([]):-!.
retract_sosu_list([First|Rest]):-
	retract(sosu(First)),
	retract_sosu_list(Rest).

assert_all([]):-!.
assert_all([First|Rest]):-
	assert(sosu(First)),
	assert_all(Rest).

% 2 ～ Idx までのリストを作成
make_list(Idx,List):-
	make_list_sub(Idx,List,[]).
	
make_list_sub(1,Ret,Ret):-!.
make_list_sub(Idx,Ret,List):-
	Idx1 is Idx - 1,
	make_list_sub(Idx1,Ret,[Idx | List]).
	
% 素数列挙 エラトステネスのふるい
eratosthenes_sieve([],[]):-!.
eratosthenes_sieve([First | Rest],[First | Rest]):-
	max(Max),
	First * First > Max ,
	!.
eratosthenes_sieve([First | Rest], [First | Ret]):-
	erase_baisu(First,Rest,Rest1),
	eratosthenes_sieve(Rest1,Ret).

erase_baisu(_,[],[]):-!.
erase_baisu(Base,[First | Rest],[First | Ret]):-
	First mod Base =\= 0,
	!,
	erase_baisu(Base,Rest,Ret).
erase_baisu(Base,[_ | Rest],Ret):-
	erase_baisu(Base,Rest,Ret).

% Num を回転させた数字が Rot (Choice Pointあり)
% 123 -> 231 , 312
rotate(Num,Rot):-
	number_chars(Num,Chrs),
	length(Chrs,Len),
	rotate_sub(Chrs,Len,Rot).

rotate_sub(_,0,_):-!,fail.
rotate_sub(Chrs,_,RotNum):-
	number_chars(RotNum,Chrs).
rotate_sub([First | Rest],Idx ,Rot ):-
	append(Rest,[First],Rot1),
	Idx1 is Idx -1 ,
	rotate_sub(Rot1,Idx1,Rot).

実行結果：

1 ?- 
% c:/Documents and Settings/Owner/My Documents/Prolog/junkan_sosu.pl compiled 0.02 sec, 25 clauses
1 ?- time(solve(999999)).
***********************(解答伏せます)**********************
% 42,859,931 inferences, 35.094 CPU in 35.156 seconds (100% CPU, 1221298 Lips)
true.

問題

% problem 3

solve:-
	write('600851475143 = '),
	div_mod0(600851475143,2),
	!.

div_mod0(1,_):-
	!,
	write('1').
div_mod0(BaseNum,DivNum):-
	BaseNum mod DivNum =:= 0,
	Divided is BaseNum / DivNum,
	DivNum1 is 2,
	write(DivNum),write(' * '),
	div_mod0(Divided,DivNum1).

div_mod0(BaseNum,DivNum):-
	DivNum1 is DivNum + 1,
	div_mod0(BaseNum,DivNum1).

実行結果：

1 ?- time(solve).
***********************(解答伏せます)**********************
% 18,489 inferences, 0.016 CPU in 0.016 seconds (100% CPU, 1183296 Lips)
true.

問題

最初の問題は簡単だなあ…作業感出てきた…早くもう少し難しくなってほしい

% problem 11

:- dynamic
	gvar/2.
	
gvar(max,0).

set_gvar(Name,X):-
	nonvar(Name),retract(gvar(Name,_)),!,asserta(gvar(Name,X)).
set_gvar(Name,X):-
	nonvar(Name),asserta(gvar(Name,X)).

solve:-
	set_gvar(squares,
		masu(
			hline(08,02,22,97,38,15,00,40,00,75,04,05,07,78,52,12,50,77,91,08),
			hline(49,49,99,40,17,81,18,57,60,87,17,40,98,43,69,48,04,56,62,00),
			hline(81,49,31,73,55,79,14,29,93,71,40,67,53,88,30,03,49,13,36,65),
			hline(52,70,95,23,04,60,11,42,69,24,68,56,01,32,56,71,37,02,36,91),
			hline(22,31,16,71,51,67,63,89,41,92,36,54,22,40,40,28,66,33,13,80),
			hline(24,47,32,60,99,03,45,02,44,75,33,53,78,36,84,20,35,17,12,50),
			hline(32,98,81,28,64,23,67,10,26,38,40,67,59,54,70,66,18,38,64,70),
			hline(67,26,20,68,02,62,12,20,95,63,94,39,63,08,40,91,66,49,94,21),
			hline(24,55,58,05,66,73,99,26,97,17,78,78,96,83,14,88,34,89,63,72),
			hline(21,36,23,09,75,00,76,44,20,45,35,14,00,61,33,97,34,31,33,95),
			hline(78,17,53,28,22,75,31,67,15,94,03,80,04,62,16,14,09,53,56,92),
			hline(16,39,05,42,96,35,31,47,55,58,88,24,00,17,54,24,36,29,85,57),
			hline(86,56,00,48,35,71,89,07,05,44,44,37,44,60,21,58,51,54,17,58),
			hline(19,80,81,68,05,94,47,69,28,73,92,13,86,52,17,77,04,89,55,40),
			hline(04,52,08,83,97,35,99,16,07,97,57,32,16,26,26,79,33,27,98,66),
			hline(88,36,68,87,57,62,20,72,03,46,33,67,46,55,12,32,63,93,53,69),
			hline(04,42,16,73,38,25,39,11,24,94,72,18,08,46,29,32,40,62,76,36),
			hline(20,69,36,41,72,30,23,88,34,62,99,69,82,67,59,85,74,04,36,16),
			hline(20,73,35,29,78,31,90,01,74,31,49,71,48,86,81,16,23,57,05,54),
			hline(01,70,54,71,83,51,54,69,16,92,33,48,61,43,52,01,89,19,67,48)
		)
	),
	set_gvar(max,0),
	\+vartical,
	\+horizontal,
	\+diagonal_topleft2bottomright,
	\+diagonal_topright2bottomleft,
	!,
	gvar(max,Max),
	write(Max).
	
vartical:-
	gvar(squares,Masu),
	between(1,20,H),
	between(1,17,V),
	V1 is V+1,
	V2 is V+2,
	V3 is V+3,
	arg(V,Masu,HLine),
	arg(V1,Masu,HLine1),
	arg(V2,Masu,HLine2),
	arg(V3,Masu,HLine3),
	arg(H,HLine,Val),
	arg(H,HLine1,Val1),
	arg(H,HLine2,Val2),
	arg(H,HLine3,Val3),
	Mul is Val * Val1 * Val2 * Val3,
	gvar(max,Max),
		(Mul > Max -> set_gvar(max,Mul);true),
	fail.

horizontal:-
	gvar(squares,Masu),
	between(1,17,H),
	between(1,20,V),
	H1 is H+1,
	H2 is H+2,
	H3 is H+3,
	arg(V,Masu,HLine),
	arg(H,HLine,Val),
	arg(H1,HLine,Val1),
	arg(H2,HLine,Val2),
	arg(H3,HLine,Val3),
	Mul is Val * Val1 * Val2 * Val3,
	gvar(max,Max),
	(Mul > Max -> set_gvar(max,Mul);true),
	fail.
	

diagonal_topleft2bottomright:-
	gvar(squares,Masu),
	between(1,17,H),
	between(1,17,V),
	H1 is H+1,
	H2 is H+2,
	H3 is H+3,
	V1 is V+1,
	V2 is V+2,
	V3 is V+3,
	arg(V,Masu,HLine),
	arg(V1,Masu,HLine1),
	arg(V2,Masu,HLine2),
	arg(V3,Masu,HLine3),
	arg(H,HLine,Val),
	arg(H1,HLine1,Val1),
	arg(H2,HLine2,Val2),
	arg(H3,HLine3,Val3),
	Mul is Val * Val1 * Val2 * Val3,
	gvar(max,Max),
	(Mul > Max -> set_gvar(max,Mul);true),
	fail.

diagonal_topright2bottomleft:-
	gvar(squares,Masu),
	between(4,20,H),
	between(4,20,V),
	H1 is H-1,
	H2 is H-2,
	H3 is H-3,
	V1 is V+1,
	V2 is V+2,
	V3 is V+3,
	arg(V,Masu,HLine),
	arg(V1,Masu,HLine1),
	arg(V2,Masu,HLine2),
	arg(V3,Masu,HLine3),
	arg(H,HLine,Val),
	arg(H1,HLine1,Val1),
	arg(H2,HLine2,Val2),
	arg(H3,HLine3,Val3),
	Mul is Val * Val1 * Val2 * Val3,
	gvar(max,Max),
	(Mul > Max -> set_gvar(max,Mul);true),
	fail.

実行結果(解答伏せています）：

1 ?- time(solve).
XXXXXXXX
% 12,531 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
true.

問題

プロジェクトオイラーは一度素数列挙ロジックを作るとあちこちで使いまわせる。

% problem 10

solve(Idx):-
	assert(max(Idx)),
	make_list(Idx,List),
	eratosthenes_sieve(List,List1), % エラトステネスのふるいによる素数列挙
	add_all(List1,Sum),
	write(Sum).

add_all([],0):-!.
add_all([First|Rest],Ret):-
	add_all(Rest,RestSum),
	Ret is First + RestSum.

retract_sosu_list([]):-!.
retract_sosu_list([First|Rest]):-
	retract(sosu(First)),
	retract_sosu_list(Rest).

assert_all([]):-!.
assert_all([First|Rest]):-
	assert(sosu(First)),
	assert_all(Rest).

% 2 ～ Idx までのリストを作成
make_list(Idx,List):-
	make_list_sub(Idx,List,[]).
	
make_list_sub(1,Ret,Ret):-!.
make_list_sub(Idx,Ret,List):-
	Idx1 is Idx - 1,
	make_list_sub(Idx1,Ret,[Idx | List]).
	
% 素数列挙 エラトステネスのふるい
eratosthenes_sieve([],[]):-!.
eratosthenes_sieve([First | Rest],[First | Rest]):-
	max(Max),
	First * First > Max ,
	!.
eratosthenes_sieve([First | Rest], [First | Ret]):-
	erase_baisu(First,Rest,Rest1),
	eratosthenes_sieve(Rest1,Ret).

erase_baisu(_,[],[]):-!.
erase_baisu(Base,[First | Rest],[First | Ret]):-
	First mod Base =\= 0,
	!,
	erase_baisu(Base,Rest,Ret).
erase_baisu(Base,[_ | Rest],Ret):-
	erase_baisu(Base,Rest,Ret).

実行結果：

1 ?- time(solve(2000000)).
XXXXXXXXXXX
% 91,150,992 inferences, 76.484 CPU in 77.969 seconds (98% CPU, 1191760 Lips)
true.

問題

CLPFDライブラリ読み込んでるけど使ってなかった…

:-use_module(library(clpfd)).

adding(1,0):-!.
adding(Idx,Sum):-
	(Idx mod 3 =:= 0 ; Idx mod 5 =:= 0),
	!,
	Idx1 is Idx - 1,
	adding(Idx1,Sum1),
	Sum is Sum1 + Idx.

adding(Idx,Sum):-
	!,
	Idx1 is Idx - 1,
	adding(Idx1,Sum).

実行結果
[1] 2 ?- adding(999,R).
***********************(解答伏せます)**********************

2週間ほど前にProject Eulerの存在を知り、面白いのでチャレンジしている。本日の時点では26問解いた。全部Prologで解いている。

ところどころCLPFDを使える問題があり、使っている。

制約論理プログラミングの動作の特徴？みたいなのを感じられる面白い問題があったので、紹介しようと思う。

Probrem 9 特別なピタゴラス数

この問題は、CLPFDだと探索の処理を一切書かずに以下のようにただ定義だけをずらずら書くだけで解ける。

:-use_module(library(clpfd)).

solve(Sum):-
	[A,B,C] ins 1..1000,
	A+B+C#=Sum,
	A*A + B*B #= C*C,
	A #< B,
	B #< C,
	label([A,B,C]),
	write([A,B,C]).

問題の条件をそのまま書いただけのようなプログラムだけど、本当にこれだけで解けます。

ただ、これだと非常に処理時間がかかかる。僕のしょぼいマシンでは、
以下のように150秒かかった。

1 ?- time(solve(1000)).
回答伏せます
% 69,897,996 inferences, 142.641 CPU in 153.188 seconds (93% CPU, 490029 Lips)
true 

ネットで調べるとこのピタゴラス数の法則のようなものがいくつかヒットするのですが、そのひとつが以下の図のようなもので、

ようするにBとCの関係に関してなのですが、長さBの正方形の辺を1つづつ大きくして一回りづつ囲むように配置してゆくといつか長さCの正方形と等しくなるというもので、しかもこの大きくした分の面積はA × Aと等しいというもの。

大きくした分の面積は、N(=1,2,3..)の数列として以下のように表すことができる。
En = 2 × N × B + N × N
そして、A × A = En

前述のプログラムにこの条件を追加して動作させると、答えを出す速度が段違いに速くなる。

:-use_module(library(clpfd)).

solve(Sum):-
	[A,B,C] ins 1..1000,
	A+B+C #= Sum,
	A #< B,
	B #< C,
	A * A + B * B #= C * C,
	A*A #= 2*N*B+N*N,　% 追加した条件
	N in 1..1000,　　　% 追加した条件　
	label([A,B,C]),
	write([A,B,C]).

実行結果：

[4] 6 ?- time(solve(1000)).
回答伏せます
% 14,116,046 inferences, 16.094 CPU in 16.516 seconds (97% CPU, 877114 Lips)
true 

10倍近く速くなった。
多分Aの探索空間が新しい条件の数列で絞り込まれて速くなるのだと思う。

「思う」と書いたのも制約論理プログラミングの特徴といえば特徴で、条件を書くたびに大量の制約伝播用のプログラムが内部で自動的に生成されるため、厳密な動作が非常に追いずらい。ステップ実行などのデバックもしずらい(CLPFDライブラリ内部で自動的に大量のバックトラックが発生している）。

これは悪い面で、動作が見積もれないことが原因となり正確・確実な動作を期待される産業分野でなかなか採用されないかもしれない。ちゃんとした企業ほど異常動作のときなどの原因追及フェーズを徹底的に行っているので、そのときに｢思う」とか｢これ入れたら速くなるかも」とかは通用しないわけです。動作確認のためのログを埋め込むのも結構大変です（ライブラリに直接埋め込まないと駄目だし制約伝播状態を出力した解析困難なログになることが予想される）

ちなみに制約伝播というのは、prolog の attributed_variableの仕組みを使って変数の探索空間が変更された時点でキックされる述語が予約されていて、ドミノ倒しのように次々と他の変数の探索空間をせばめてゆくという手法です。

CLPFDのソースを読んでゆくとわかるのですが A in 1..500 などという探索空間は木構造で表現されており、探索空間の変更はこの木構造を変えることで行っている。

たとえば A の木構造が最初 左の枝が1、右の枝が500 の状態で、この状態から100～200の可能性がないことが判明した時点で、右の500の枝の部分に100..200のノードが新しく追加される感じ？

Prologの変数は基本的に再代入不可なのですが、attributed_variableは再代入可能かつ履歴情報を保持していてバックトラック時に直前の値に戻る（Prologの通常の自由変数はバックトラック時は未設定状態になるだけで履歴情報は持っていない）ようなのでこれを利用しているのでしょう。

CLPFDは非常に簡単に記述できて問題が解けるため、勉強し始めのころはまるで魔法のツールのようだと感じたが、上記のような、問題自身の性格を表す条件を入れないとすぐに処理時間が膨大になってしまうように感じる。

当たり前のことだが、問題に対する洞察・分析も非常に大切ということだろう。