投稿者: koyahata

Project Euler 46

久しぶりに解きました。
あんまり良いプログラムではないけど…

%problem 46
:-lib(ic).
solve:-
	PrimeMax is 100000,
	assert(max(PrimeMax)),
	make_list(PrimeMax,Lst),
	eratosthenes_sieve(Lst,PrimeLst), % エラトステネスのふるいによる素数列挙

	Prime1#::PrimeLst,
	Prime2#::PrimeLst,
	A#=Prime1*Prime2,
	A#=_*2+1, % Aは奇数
	indomain(A,min),
	(
		Prime#::PrimeLst,
		B#>0,
		(A-Prime)/2 #= B*B,
		
		labeling([B])->
		(
%			printf(output,"%d = %d + 2*%d^2\n",[A,Prime,B]),
%			flush(output)
			true
		);
		(
			labeling([Prime1,Prime2]),
			printf(output,"answer %d = %d * %d\n",[A,Prime1,Prime2]),
			flush(output),
			abort
		)
	),
	fail.

mk_lst([],_,[]):-!.
mk_lst([First|Rest],Th,[Rest]):-
	First>Th,
	mk_lst(Rest,Th,Rest).
mk_lst([First|Rest],Th,[First|Rest]):-
	mk_lst(Rest,Th,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).

実行結果(回答伏せます) solve. answer XXXX = XXXX * XXXX Aborting execution ...

2021年1月17日

CHRをGPUで実行させる論文を読みました。

Parallel Execution of Constraint Handling Rules on a Graphical Processing Unit
という論文を読みました。
リンク
リンク切れている場合は論文名で検索してください。

Constraint Handling Rulesの開発者のThom Fr¨uhwirth さん達が書いた論文で、CHRをGPU上で実行させる方法がCUDAのコード込みで書いてあります。

Constraint Handling Rulesに関して自分がわかっているレベル(≒あまりわかっていない)で簡単に説明すると、CHRは基本的には記述された一連の制約を論理的に同等かつよりシンプルなものにしたり、何かに変換（自作の制約で書いた記述 → 組み込み述語のみ使用した記述に変換など）したりします。

例えばECLiPSeではPropiaとは別にCHRを使用しても自作の制約が定義出来ます。

まず以下の3種類のいずれかの構文でルールを記述しておきます。(これが自作制約となります)
CHRには以下の3つしかルールがありません。
・simplification
　　Head1,Head2 <=> Guard | Body
・propagation
　　Head1,Head2 => Guard | Body
・simpagation
　　Head1¥Head2 <=> Guard | Body
（HeadもGuardもBodyも普通のPrologのCallable Termで書きます。
GuardとBodyの部分は複数の述語を書いて良い。Headの部分は最大2つの述語のみ（ECLiPSeの場合。3つ以上が不可なのは以下に述べる述語のピックアップの際の効率の問題らしいです）→すみません訂正します。CHRの昔のライブラリでは2つのみでしたが、新しいECHライブラリでは3つ以上可能でした）

動作の仕組みなのですが、上記で定義した制約を含む述語群を連言（=AND。カンマ区切りの普通のPrologの表現）で記述すると、CHRの機構はそれらの述語群の中から2つピックアップし、
・上記のHead1,Head2にパターンマッチ
・GuardがTrueとなる
を満たすものを見つけると、以下のルールに従ってピックアップした2つの述語を「書き替えて」しまいます（これをfire[発火]と呼ぶ）そしてその際にBodyの実行可能な述語は実行してしまいます。
また、Guardの部分は常時チェックされる番犬のような動作をしていて、例えばguardにnonvar(A)みたいな記述をしていた際、別の処理でAに何か値を設定した場合、その瞬間にGuardのチェックが通り対応するBodyが発火されます（ECLiPSeのsuspend機構を使って実装されている）
simplificationの場合はBodyに書き換え
propagationの場合はHead1,Head2そのままでBodyの述語を追加
simpagationはHead2のみBodyに書き換え(Head1は残す)

上記の変換を適用できるルールがなくなるまで繰り返す（2つ述語ピックアップ→Head1,Head2にマッチングするかチェック→GuardがTrueかチェック→ルールに従いHead1,Head2を書き換え）
どうも書き換え中のルール群をconstraint storeと呼んでいるようです。

ちなみに上記一連の処理中にはバックトラックは発生しません。
CHRの動作を「コミッテドチョイス」と言うらしいのですが、CHRの機構がどの述語をピックアップするか勝手に決めるのでプログラマはコントロールできない（選択を委託する）という意味かもしれません。
また、Headの部分は２つじゃなくて1つでも良いです。

CHRは基本的には上記の書き換えを行うのみなので、用途は特に制約の表現のみにとどまらず、サンプルを見ると
　最大公約数の計算
　ガウスの消去法
　高速フーリエ変換
　マージソート
など、色々な応用があります。

めちゃめちゃ端折ってますが上記がCHRの基本動作です。

記述した制約群をどのような組み合わせでピックアップ・変換したとしても必ず最終的に生成される述語群が同じになるようなとき、Confluence（合流）性を持つと言うそうです。

Confluence性はどのようにすれば担保されるかなのですが、以下の論文に書いてありますが自分はまだ理解できていない（論理学の用語がたくさん使用されてて難しいため）どうも冗長なルールを追加することにより合流性が担保されるようです。

Theory and practice of constraint handling rules
リンク切れしている場合は論文名で検索してください。

Confluence性を満たす場合CHRが並列処理に非常に向いているとのことで、理屈では確かにそうだなと自分も納得しました。

前置きが長くなりましたが、上記を踏まえたうえでGPUでCHRを使った上記の論文を読んだわけです（GPUは数千コア積んでいるものもあり単純な計算を同時に行うのに非常に有力）
論文の内容はCUDAコードが実際に記載されていて実装イメージがしやすいものでした。
しかし、ルールの述語自体をCUDAコンパイル用のコードで書いているため、あまり汎用性がないのではないという気がしました。
使う側は当然任意の述語を使いたいのですが、論文のやり方ではそのたびに新たなプログラムをCUDA用に書いてコンパイルしなくてはならないのでは？と感じました。
あまりGPUのことわかってないですが、if文を多用するプログラムは遅くなるなど、性能を100%出すためにいろいろ条件があるようです。
現状ではCHRをGPUの速度を100%出し切って実行するのは難しいのかな、と感じました。
ただ理論的にCHRのいろいろな性質が証明されているので、ハードウェア側でCHRに向いた実装が実現されたら、かなり望みがあるような気がします。

余談ですが「CHRには3つしかルールがない」と書きましたが、
・simpagation
　　Head1¥Head2 <=> Guard | Body
の、
・simplification はこれの Head1 がないバージョン（Head2がBodyに書き換えられる)
・propagation はこれの Head2 がないバージョン(Head1そのままでBodyが追加される)
となるので実質はsimpagationだけで良いみたいな文章をどこかで読みました。

2021年1月2日
labelingに関して

制約論理プログラミングで必ず出てくるlabelingに関して説明しようと思います。
だだーっとマニュアル見ないでいい加減に書くので多分嘘も書いてますが大体で読んでください。なるべくプログラム自体も書かない記述にします（めんどくさいので）

概要的にはECLiPSe CLP のtutorialに書いてあることを抜粋して書くので英語読める人はそちらを読んでもらう方が絶対に良いです。とても分かりやすいです。
ECLiPSeのtutorial

制約論理プログラミングでは基本的なプログラムの構造は以下のような形になってます。
１．制約の記述
２．解のサーチ
３．解の表示

labelingというのはこれの2．解のサーチに相当します。

labeling述語には制約変数のリストを渡します。
例えば、（便宜的に整数のみ扱います。というか整数しか自分は扱い方が良くわかりません。）
変数Aが１から３
変数Bが１から３
変数Cが１から３
といった制約をつけた後、labeling([A,B,C])というような呼び出しをすると、A=1,B=1,C=1という解が設定されます。
Prologをやっている人はピンとくると思いますが、labelingの呼び出しでバックトラックが発生し、failすると次のA=2,B=1,C=1という別解が設定されます。バックトラックするごとに制約を満たす解が順番に設定されます。

ちょっと余談で、ECLiPSe CLPの資料を読んでるとことあるごとに口酸っぱく書いてあるのですが
「１．の制約の記述の部分は、バックトラックを発生させてはいけない」とのことです。制約の部分でバックトラックを発生させる作りにすると正確な解が得られないとのことです。
バックトラックの有無がPrologの一般的な述語と制約論理プログラミングの制約の大きな違いだと思います。
ですので、ECLiPSe CLPにもSWI-PrologのCLPFDにもたくさん制約用の述語が定義されてますが、そのどれもがバックトラックをしない述語になってるはずです。

これもまた余談ですがECLiPSeには自分が任意に定義したPrologの述語を制約に変換することが出来るPropiaという仕組みがあり、とても便利です。これを使うとバックトラックを生成させてしまう術語がバックトラックを生成しない制約に変わります。

話は戻って、labelingなのですが、既定では与えられた変数群の解を生成する際、
・リストの先頭の変数から順番に値を確定
・値はドメイン（取りうる範囲）の小さいほうから大きいほうへ
決定しています。

ただ、これが一番良いやり方かどうかはわかりません。
SWI-PrologのclpfdでもECLiPSe CLPでも「どの変数の値を優先的に決定するか」「値をどのように設定するか」の2つの観点でいろいろなオプションが用意されてます。
例えば以下の観点のオプションなどあります。
・一番ドメインの範囲（取りうる値の範囲）が小さい変数から決定する
・値を中央から外側に向かうように設定する
・値をランダムに設定する

SWI-PrologのCLPFDではここまでの話でlabelingに関しての説明は完了なのですが、ECLiPSe CLPはこのlabelingの部分を完全にプログラム制御できます。tutorialのTree Search Methodの項を見てください。

まずlabelingを使用せずに、1つの制約変数を指定して、
indomain(変数)という書き方ができます。
これは変数の取りうる値を小さいほうから順に設定してfailするたびバックトラックをさせるという、labelingの1変数バージョンのような述語です。そしてindomainの2番目の引数にオプションを渡すと、いろいろな順番で値を決めることが出来ます（大きい順、ランダム、真ん中から、etc…）

また、deleteという紛らわしい名前の述語があるのですが、これにfirst_failなどのオプションを指定して制約変数のリストを渡すと、指定したオプションに応じたやり方で1つ変数を取り出し、それ以外の変数をRestに返してくれます。
first_failの場合はリスト内の一番ドメインが少ない変数を選択し、残りをRestで返します。
（「次に値を決める変数の選択」を行ってくれる述語）

ECLiPSeではこのdeleteとindomainを駆使して、細やかに変数のラベリングをコントロールできます。ECLiPSeのtutorialには、これらを使用してN-Queenの問題をいろいろなラベリング方法で解き、発生したバックトラックの回数を計測するという興味深い記事が出てますので是非読んでみてください。

あとsearch/6という似たようなことする術語もあるのですがまだよくわかってません。
あとECLiPSeでは整数以外の値が扱えるのですが、そのときにindomainの実数バージョンのようなlocateという術語があるのですが、これもまだ詳細に使いこなせてません。

これ読んでる方もぜひ一緒に勉強しましょう（そして教えて下さい）

2020年12月4日
DOMINO PUZZLEソルバーの解説

ここ最近、勉強で ECLiPSe CLPのExampleページのドミノタイル問題ソルバーのソースを解析していました。大体理解したので解説したいと思います。

問題の解説：
0-0、0-1、0-2…..6-6 の全部で28枚のドミノタイルがある。

このタイルが横８×縦７の長方形の中にランダムに敷き詰められている。
縦、横、左右逆、上下逆のどの置き方でも良い。

どのように置いたかは不明で、並べた結果全体の数字が以下のようになっていた場合、タイルがどのように配置されていたかを解け。
3 1 2 6 6 1 2 2
3 4 1 5 3 0 3 6
5 6 6 1 2 4 5 0
5 6 4 1 3 3 0 0
6 1 0 6 3 2 4 0
4 1 5 2 4 3 5 5
4 1 0 2 4 5 2 0

プログラムの解説：

プログラムの概要を説明します。
①述語one_by_two_tiling/4でLeftとUpというマス同士の接続状態を表す配列を作成、制約を設定します。

Leftの配列は以下の色のついた部分の接続状態を表す2次元配列です。
1つのドミノタイルとしてつながっていれば１、つながっていなければ０になります。

Upは以下の色のついた部分の接続状態を表す2次元配列です。

制約のつけ方なのですが、あるマスに注目し上下左右の接続状態を考えたとき、必ず上下左右のいずれか1つの方向のみ接続されていることになるので、そのような制約を設定します（上下左右の変数の和が１）

②２つのdo文で、隣接した2マスに注目しながら、「2つの数字-対応する接続変数」の組み合わせをどんどんリストに集めます。このとき、タイルは数字が小-大の順番になるようにして保存します。
その部分にタイルが配置されていた場合は、該当する接続状態の変数が1になるというわけです。

ヨコ方向の隣接マスに注目した後は、タテ方向で同様の操作を行います。どちらも同じ配列にまとめます。
このようにまとめると当然 (1,4)-変数というような組み合わせが複数見つかりますが、いくつか見つかった中のどれか1つの組み合わせだけが実際に配置されているタイルの情報だということになります。

③まとめた配列をkeysortでソートします。タイルごとにまとまります。

④group_same_key_values/2で、ソート後の
[…(1,4)-接続状態変数1,(1,4)-接続状態変数2,(1,4)-接続状態変数3,(1,5)-接続状態変数4,(1,5)-接続状態変数5…]
のような状態のリストを
[…(1,4)-[接続状態変数1,接続状態変数2,接続状態変数3],(1,5)-[接続状態変数4,接続状態変数5]…]
のような形式にまとめます。

⑤
[接続状態変数1,接続状態変数2,接続状態変数3]　の和が１
[接続状態変数4,接続状態変数5]　の和が１
というような制約を設定します(foreach)

⑥term_variableで全変数を1つのリストにまとめた後labelingですべての接続状態変数の解を探します。

⑦
prety_print出力結果：

2020年12月1日
四角に切れ（パズル）を解くプログラム

昔いくつか書いたPrologのパズルソルバーのコードがリンク切れになってたので再掲します。

今回は四角に切れ（パズル）のソルバーです。

プログラム
 解説と実行結果

これも数年後、英語のStackOverFlow で以下の ECLiPSe CLP で書かれたプログラムを発見しました。
ECLiPSe の開発メンバーのjschimpfさんのコードですね。
ECLiPSeで書かれた四角に切れソルバー

2020年12月1日
RUSH HOUR(パズル)を解くプログラム

昔いくつか書いたPrologのパズルソルバーのコードがリンク切れになってたので再掲します。

今回はラッシュアワーというパズルを解くプログラムです。

プログラム
 解説と実行結果

2020年12月1日
ペントミノを解くプログラム

昔いくつか書いたPrologのパズルソルバーのコードがリンク切れになってたので再掲します。

今回はペントミノのソルバーです。

プログラム
 解説と実行結果

2020年12月1日

数独を解くプログラム

昔いくつか書いたPrologのパズルソルバーのコードがリンク切れになってたので再掲します（これから4回分続きます）。

まずは数独のソルバーです。

これ書いたときはまだ制約論理プログラミングのことを知らなかったのでPure Prologで書いてます（SWI-Prolog)。

その後SWI-Prologのマニュアルのclpfdライブラリの説明の箇所の数独ソルバーのプログラムを発見して衝撃を受け、制約論理プログラミングの勉強を開始したのであった。

プログラム
 解説と実行結果

以下は僕が衝撃を受けたSWI-Prologのマニュアルに載っている数独ソルバー（僕のプログラムよりエレガントかつ速い）

:- use_module(library(clpfd)).
sudoku(Rows) :-
	length(Rows, 9), maplist(same_length(Rows), Rows),
	append(Rows, Vs), Vs ins 1..9,
	maplist(all_distinct, Rows),
	transpose(Rows, Columns),
	maplist(all_distinct, Columns),
	Rows = [As,Bs,Cs,Ds,Es,Fs,Gs,Hs,Is],
	blocks(As, Bs, Cs),
	blocks(Ds, Es, Fs),
	blocks(Gs, Hs, Is).

blocks([], [], []).
blocks([N1,N2,N3|Ns1], [N4,N5,N6|Ns2], [N7,N8,N9|Ns3]) :-
	all_distinct([N1,N2,N3,N4,N5,N6,N7,N8,N9]),
	blocks(Ns1, Ns2, Ns3).
problem(1, [[_,_,_,_,_,_,_,_,_],
	[_,_,_,_,_,3,_,8,5],
	[_,_,1,_,2,_,_,_,_],
	[_,_,_,5,_,7,_,_,_],
	[_,_,4,_,_,_,1,_,_],
	[_,9,_,_,_,_,_,_,_],
	[5,_,_,_,_,_,_,7,3],
	[_,_,2,_,1,_,_,_,_],
	[_,_,_,_,4,_,_,_,9]]).

実行結果：

?- problem(1, Rows), sudoku(Rows), maplist(portray_clause, Rows).
[9, 8, 7, 6, 5, 4, 3, 2, 1].
[2, 4, 6, 1, 7, 3, 9, 8, 5].
[3, 5, 1, 9, 2, 8, 7, 4, 6].
[1, 2, 8, 5, 3, 7, 6, 9, 4].
[6, 3, 4, 8, 9, 2, 1, 5, 7].
[7, 9, 5, 4, 6, 1, 8, 3, 2].
[5, 1, 9, 2, 8, 6, 4, 7, 3].
[4, 7, 2, 3, 1, 9, 5, 6, 8].
[8, 6, 3, 7, 4, 5, 2, 1, 9].
Rows = [[9, 8, 7, 6, 5, 4, 3, 2|...], ... , [...|...]].

2020年12月1日

Project Euler Problem 62

Problem 62

ECLiPSe CLP で解きました。
ほぼPure Prolog で、制約論理の述語は使用してません。
自分のマシンで約9秒。

プログラム：

:-lib(util).
:-dynamic(cubic_sorted/2).

go:-
	retractall(cubic_sorted(_,_)),
	between(1,1000000000000,I),
		CubicNum is I * I * I,
		digit_sort(CubicNum,Sorted),
		assert(cubic_sorted(Sorted,I)),
		findall(Num,cubic_sorted(Sorted,Num),NumLst),
		length(NumLst,NumLen),
		(
			NumLen>=5->
				sort(NumLst,NumLstSorted),
				[First|_]=NumLstSorted,
				CubeFirst is First * First * First,
				printf(output, "CubicNum %d ,Lst %w", [CubeFirst, NumLst]),nl,
				true
				;
				fail
		)
	.

digit_sort(Num,Sorted):-
	num2lst(Num,Lst),
	msort(Lst,Sorted).

% convert number to list(reversed)
num2lst(0,[]):-!.
num2lst(Num,[Dig|DigRst]):-
	Dig is Num mod 10,
	RestNum is Num // 10,
	num2lst(RestNum,DigRst).

コマンドライン：
?- time(go). Yes (9.22s cpu)

標準出力：
<解答伏せます。>

Success, times: 9.218750s thread, 9.219+0.000s process, 9.224s real

2020年7月12日

Project Euler Problem 59

Problem 59

全体の文字数 – 「何文字がA-Z a-zの範囲だったか」　の数値をコストとして計算し、一番コストが小さくなるものをbb_minで調べました。
あんまり効率よくないプログラムで、1900秒くらいかかりました…遅いな…
あとテキストの文字数が3の倍数というところは前もって調べてて、それをあてにしたプログラムになってます（文字数が3の倍数でないとdecryptでこける）

ソース

:-lib(ic).
:-lib(ic_sets).
:-lib(propia).
:-lib(branch_and_bound).
:-lib(csv).

go(Str,Cost,Sum) :-
	csv:(csv_read("p059_cipher.txt", [InLst], Option)),
	char_code('a',AscA),
	char_code('z',AscZ),
	[Key1,Key2,Key3]#::AscA..AscZ,
	decrypt(InLst,OutLst,[Key1,Key2,Key3]),
	OutLst#::0..127,
	OutLst#::0..127,
	letter_count(OutLst,Cnt),
	length(OutLst,Len),
	Cost #= Len - Cnt, 
	flatten([Key1,Key2,Key3,OutLst],AllVals),
	labeling(AllVals),
	iso_light:(atom_codes(Str,OutLst)),
	mysum(OutLst,Sum).

decrypt([],[],_).
decrypt([N1,N2,N3|NRest],[D1,D2,D3|DRest],[K1,K2,K3]):-
	xor(N1,K1,D1) infers most,
	xor(N2,K2,D2) infers most,
	xor(N3,K3,D3) infers most,
	decrypt(NRest,DRest,[K1,K2,K3]).

letter_count([],0).
letter_count([First|Rest],Cnt):-
	letter_count(Rest,Cnt1),
	#::(First, [65..90, 97..122],Hit),
	Cnt #= Cnt1+Hit.

mysum([],0).
mysum([First|Rest],Sum):-
	mysum(Rest,Sum1),
	Sum is Sum1 + First.

実行結果

?- bb_min(go(A, Cost, Sum), Cost, _223).
A = 復号化した文字列（伏せます）
Cost = 伏せます
Sum = 伏せます
Yes (1893.05s cpu)

2020年6月23日