「prolog勉強プロジェクトオイラー161～170」の編集履歴（バックアップ）一覧に戻る
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。
*Problem 161 「トリオミノ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20161
トリオミノで9*12の長方形を埋め尽くすとき何通りの埋め方があるか答える問題。
*解法
左下から右上へ敷き詰めていく。
ビット演算に翻訳して敷き詰め、あとは敷き詰め終わった部分は敷き詰め終わった部分の形（どのようなピースでどのように埋めたかは無視して全体の形だけ）を主キーにして動的計画法でアセプト。
計算に必要なのは３行だけなので9列*3行の２７ビットで管理ができこれはＰｒｏｌｏｇの一番短い整数型28bitより短い。
これで少し高速化を達成できているかもしれない。
Prologコード計算時間 2.112sec。
平凡タイム？
速い人はどうしてるんだろ？
piece(Piece,Row,Col):-
Row=<11,
Col=<8,
Piece is 1+(3<<9).
piece(Piece,Row,Col):-
Row=<11,
1<Col,
Piece is 1+(3<<8).
piece(Piece,Row,Col):-
Row=<11,
Col=<8,
Piece is 3+(1<<9).
piece(Piece,Row,Col):-
Row=<11,
Col=<8,
Piece is 3+(1<<10).
piece(Piece,Row,_):-
Row=<10,
Piece is 1+(1<<9)+(1<<18).
piece(Piece,Row,Col):-
Row=<12,
Col=<7,
Piece is 7.
union_sum([],Y,[Y]):-!.
union_sum([[Set,Perm]|Rest],[Set,Perm1],Result):-
!,
Perm2 is Perm+Perm1,
union_sum(Rest,[Set,Perm2],Result).
union_sum([X|Rest],Y,[Y|Result]):-
union_sum(Rest,X,Result).
one_calc(No,_,Col,No):-0<(No/\(1<<(Col-1))),!.
one_calc(No,Row,Col,No1):-
piece(Piece,Row,Col),
Slide is (Piece<<(Col-1)),
0=:=No /\ Slide,
No1 is No \/ Slide.
calc_next(Memo,Row,Col,[No1,Perm]):-
member([No,Perm],Memo),
one_calc(No,Row,Col,No1).
row_dell([],[]):-!.
row_dell([[No,Perm]|Rest],[[No1,Perm]|Result]):-row_dell(Rest,Result),
No1 is (No>>9).
roop_col(Memo,Row,10):-!,
Row1 is Row+1,
row_dell(Memo,Memo1),
roop_row(Memo1,Row1).
roop_col(Memo,Row,Col):-
Col1 is Col+1,
findall(M,calc_next(Memo,Row,Col,M),Memo1),
msort(Memo1,Memo2),
[Top|Memo3]=Memo2,
union_sum(Memo3,Top,Memo4),
roop_col(Memo4,Row,Col1).
roop_row(Memo,13):-
!,nl,nl,
write(Memo).
roop_row(Memo,Row):-
!,
length(Memo,Len),
write([Row,Len]),nl,
roop_col(Memo,Row,1).
main161:-roop_row([[0,1]],1).
*Problem 162 「16進数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20162
16進法では, 数は以下の16個の数字によって表される
0,1,2,3,4,5,6,7,8,9,A,B,C,D,E,F
16進数の AF は, 10進法での 10x16 + 15 = 175 と等しい.
3桁の16進数 10A, 1A0, A10, A01 には, 0, 1, A の全てが現れている.
10進数で書くときと同様に, 先頭の0は書かないことにする.
0, 1, A がそれぞれ少なくとも1回は現れるような16桁までの16進数はいくつ存在するか?
16進数で答えよ.
(A,B,C,D,E,F は大文字とし, 先頭や末尾の16進数であることを表す記号や先頭の0は許されない. つまり, 1A3F ならOK. 1a3f, 0x1a3f, $1A3F, #1A3F, 0000001A3Fは許されない. )
解法
一桁ずつ伸ばして組み合わせ数を考えます、ある桁数で管理すべきことは。
0が0か1回以上出たの２通り
1が0か1回以上出たの２通り
Aが0か1回以上出たの２通り
で計８通りを管理して
ある桁を計算するとき0が出た１が出たＡがでたそれ以外が出た。
を管理しながら計算していけばよいだけです。
配列のある言語ならBigO(15*8*4)で計算が完了します。
one_calc([ _,N1,Na,Perm],[1,N1,Na ,Perm]).
one_calc([N0, _,Na,Perm],[N0,1,Na ,Perm]).
one_calc([N0,N1, _,Perm],[N0,N1,1 ,Perm]).
one_calc([N0,N1,Na,Perm],[N0,N1,Na,RePerm]):-
RePerm is Perm*13.
next_calc(Memo,Result):-
member(M,Memo),
one_calc(M,Result).
union_sum([],Y,[Y]):-!.
union_sum([[N0,N1,Na,Perm]|Rest],[N0,N1,Na,Perm1],Result):-
!,
Perm2 is Perm+Perm1,
union_sum(Rest,[N0,N1,Na,Perm2],Result).
union_sum([X|Rest],Y,[Y|Result]):-
union_sum(Rest,X,Result).
search(_,16,Ans):-
!,
format('~16r',[Ans]).
search(Memo,Len,Ans):-
findall(M,next_calc(Memo,M),Memo1),
msort(Memo1,Memo2),
[Top|Memo3]=Memo2,
union_sum(Memo3,Top,Memo4),
Len1 is Len+1,
(member([1,1,1,Perm],Memo4)->Ans1 is Ans+Perm;Ans1 is Ans),
search(Memo4,Len1,Ans1).
main162:-
Seed=[[0,0,0,13],[0,1,0,1],[0,0,1,1]],
search(Seed,1,0).
*Problem 164 「どの連続した3桁の和も与えられた数を超えない数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20164
問題
どの連続した3桁の和も9以下のような(9以下は9を含む)20桁の数(先頭の0は認めない)はいくつあるか?
解答
教科書的な動的計画法で何も工夫するところがありません。
簡単すぎる問題なのでサクッとコードを書いてしまいましょう。
Prologは破壊的代入がない言語なので段階を置いて集計する必要があるくらいです。
seed([N1,N2,N3,1]):-
between(1,9,N1),
between(0,9,N2),
between(0,9,N3),
N1+N2+N3=<9.
sum([],0):-!.
sum([[_,_,_,P]|Ps],Result):-sum(Ps,Re),Result is Re+P.
union_sum([],Set,[Set]):-!.
union_sum([[N1,N2,N3,P]|Memo],[N1,N2,N3,P1],Result):-
!,
P2 is P1+P,
union_sum(Memo,[N1,N2,N3,P2],Result).
union_sum([Set|Memo],Set1,[Set1|Result]):-
union_sum(Memo,Set,Result).
calc_next_B(Memo,[N2,N3,N,P]):-
member([_,N2,N3,P],Memo),
between(0,9,N),
N2+N3+N=<9.
calc_next_A(Memo,Memo4):-
findall(Set,calc_next_B(Memo,Set),Memo1),
msort(Memo1,Memo2),
[Top|Memo3]=Memo2,
union_sum(Memo3,Top,Memo4).
calc(20,Memo):-!,sum(Memo,Ans),write(Ans).
calc(N,Memo):-
calc_next_A(Memo,Memo4),
N1 is N+1,
calc(N1,Memo4).
main164:-
findall(Set,seed(Set),Memo),
calc(3,Memo).
*Problem 165 「交点」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20165
全部の線分の交差判定をしました、遅いです。
交点は分数座標で保持することで同じ点かを判断します。
グローバルスタックが足らないのでコードの実行時に拡張してアセプト。
Proloｇにはstd::setがないのでメモリや速度が無駄になるのが難しいところ。
線分の生成規則から何か交点の発生に性質があるのだろうか？
gcd(0, B, G) :- G is abs(B).
gcd(A, B, G) :- A =\= 0, R is B mod A, gcd(R, A, G).
seed(20000,_,Line,[Line]):-!.
seed(N,S,Line,[Line1|Result]):-
0=:= N mod 4,N>0,
!,
S1 is (S*S) mod 50515093,
N1 is N+1,
T is S1 mod 500,
reverse(Line,Line1),
seed(N1,S1,[T],Result).
seed(N,S,Line,Result):-
S1 is (S*S) mod 50515093,
N1 is N+1,
T is S1 mod 500,
seed(N1,S1,[T|Line],Result).
cross(X1,Y1,X2,Y2,Result):-
Result is X1*Y2-X2*Y1.
check([X1,Y1,X2,Y2],[X3,Y3,X4,Y4]):-
SaX is X2-X1,
SaY is Y2-Y1,
SaX3 is X3-X1,
SaY3 is Y3-Y1,
SaX4 is X4-X1,
SaY4 is Y4-Y1,
0>(SaX*SaY3-SaY*SaX3)*(SaX*SaY4-SaY*SaX4).
calc_point([X1,Y1,X2,Y2],[X3,Y3,X4,Y4],[AnsX2,AnsY2,GcdX1,GcdY1]):-
SaX2 is X2-X1,
SaY2 is Y2-Y1,
SaX4 is X4-X3,
SaY4 is Y4-Y3,
cross(SaX4,SaY4,SaX2,SaY2,Down),
SaX31 is X3-X1,
SaY31 is Y3-Y1,
cross(SaX4,SaY4,SaX31,SaY31,Up),
AnsX1 is X1*Down+SaX2*Up,
AnsY1 is Y1*Down+SaY2*Up,
gcd(AnsX1,Down,GcdX),
gcd(AnsY1,Down,GcdY),
!,
AnsX2 is abs(AnsX1//GcdX),
AnsY2 is abs(AnsY1//GcdY),
GcdX1 is abs(Down//GcdX),
GcdY1 is abs(Down//GcdY).
roopB(L1,[L2|_],Ans):-
check(L1,L2),
check(L2,L1),
calc_point(L1,L2,Ans).
roopB(L1,[_|Lines],Result):-
!,
roopB(L1,Lines,Result).
roopA([L1|Lines],Result):-
roopB(L1,Lines,Result).
roopA([_|Lines],Result):-
!,
roopA(Lines,Result).
main165:-X is 512*1024*1024,set_prolog_stack(global,limit(X)),
seed(0,290797,[],Lines),nl,
findall(Ans,roopA(Lines,Ans),Answers),
sort(Answers,Answers1),
length(Answers1,Len),write(Len).
*Problem 166 「十字」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20166
とりあえず答えが見たかったのでC++で実装。
対称性を利用して上２行を計算し終えたら下２行は上２行の１８０度回転として考えて３．３秒を達成。
C++は早いうえに新しいパソコンだからこれは卑怯くさいタイムな気もする。
#include<stdio.h>
#include<map>
#include<iostream>
#include <time.h>
struct SumSet{
int c1,c2,c3,c4,dl,dr;
bool operator<(const SumSet& ss)const {
if(c1!=ss.c1)return c1<ss.c1;
if(c2!=ss.c2)return c2<ss.c2;
if(c3!=ss.c3)return c3<ss.c3;
if(c4!=ss.c4)return c4<ss.c4;
if(dl!=ss.dl)return dl<ss.dl;
if(dr!=ss.dr)return dr<ss.dr;
return false;
}
};
__int64 calc(const int LIMIT){
std::map<SumSet,__int64> memo,nextMemo;
std::map<SumSet,__int64>::iterator it;
SumSet ss,ss2;
__int64 perm,perm2,ans=0;
int temp;
ss.c1=ss.c2=ss.c3=ss.c4=ss.dl=ss.dr=0;
memo[ss]=1;
for(int row=0;row<2;row++){
for(it=memo.begin();it!=memo.end();it++){
ss=(*it).first;
perm=(*it).second;
for(int a=0;a<=9;a++){
temp=ss.c1+a;
if(a>LIMIT)break;
if(temp>LIMIT)continue;
for(int b=0;b<=9;b++){
temp=ss.c2+b;
if(a+b>LIMIT)break;
if(temp>LIMIT)continue;
for(int c=0;c<=9;c++){
temp=ss.c3+c;
if(a+b+c>LIMIT)break;
if(temp>LIMIT)continue;
int d=LIMIT-a-b-c;
if(d<0||9<d)continue;
ss2.c1=ss.c1+a;
ss2.c2=ss.c2+b;
ss2.c3=ss.c3+c;
ss2.c4=ss.c4+d;
if(row==0){
ss2.dl=ss.dl+a;
ss2.dr=ss.dr+d;
}else{
ss2.dl=ss.dl+b;
ss2.dr=ss.dr+c;
}
if(ss2.dl>LIMIT)continue;
if(ss2.dr>LIMIT)continue;
if(nextMemo.find(ss2)==nextMemo.end())nextMemo[ss2]=0;
nextMemo[ss2]+=perm;
}
}
}
}
memo.clear();
memo.insert(nextMemo.begin(),nextMemo.end());
nextMemo.clear();
}
for(it =memo.begin();it!=memo.end();it++){
perm=(*it).second;
ss=(*it).first;
ss2.c4=LIMIT-ss.c1;
ss2.c3=LIMIT-ss.c2;
ss2.c2=LIMIT-ss.c3;
ss2.c1=LIMIT-ss.c4;
ss2.dl=LIMIT-ss.dl;
ss2.dr=LIMIT-ss.dr;
if(memo.find(ss2)==memo.end())continue;
ans+=perm*memo[ss2];
}
std::cout<<LIMIT<<","<<ans<<"\n";
return ans;
}
int main(){
clock_t start,end;
start = clock();
__int64 ans=0;
for(int i=0;i<=36;i++){
ans+=calc(i);
}
end = clock();
std::cout<<"ans"<<" "<<ans<<"\n";
printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC);
}
***問166Prolog版
縦横斜めの合計値を決めて２段目までの組み合わせを求めるところまではC++と同じ。
一つ一つの組み合わせに対し足らない分を縦横斜めで計算してそのセットを保存。
上２段を１８０度回転したものと同じになる組み合わせ同士をcalc_perm述語で照合してかけていけば何とか現実的な速度で解けます。
コード実行時間約１分３０秒。
C++がいかに早いかわかるというものです。
sum([],0):-!.
sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
one_check(Row,N,Limit):-
N =< Limit,
Limit-(3-Row)*9=<N.
naname(0,0,3):-!.
naname(1,1,2):-!.
naname(2,2,1):-!.
naname(3,3,0):-!.
next_one_calc(R,4,Limit,[DL,DR,Perm],Adds,_,[ReDL,ReDR,Perm]):-
!,
naname(R,PL,PR),
nth0(PL,Adds,AddDL),
ReDL is DL+AddDL,
one_check(R,ReDL,Limit),
nth0(PR,Adds,AddDR),
ReDR is DR+AddDR,
one_check(R,ReDR,Limit).
next_one_calc(R,C,Limit,[Col|Cols],Adds,AddSum,[ReCol|Result]):-
!,
between(0,9,Add),
ReCol is Col+Add,
AddSum1 is AddSum+Add,
one_check(R,ReCol,Limit),
one_check(C,AddSum1,Limit),
C1 is C+1,
next_one_calc(R,C1,Limit,Cols,[Add|Adds],AddSum1,Result).
next_calc(R,Limit,Memos,Result):-
member(Memo,Memos),
next_one_calc(R,0,Limit,Memo,[],0,Result).
union_sum([],Y,[Y]):-!.
union_sum([[C1,C2,C3,C4,DL,DR,Perm]|Rest],[C1,C2,C3,C4,DL,DR,Perm1],Result):-
!,
Perm2 is Perm+Perm1,
union_sum(Rest,[C1,C2,C3,C4,DL,DR,Perm2],Result).
union_sum([X|Rest],Y,[Y|Result]):-
!,
union_sum(Rest,X,Result).
calc_perm([],_,Ans,Ans):-!.
calc_perm(_,[],Ans,Ans):-!.
calc_perm([[C1,C2,C3,C4,DL,DR,Perm]|Rest1],
[[C1,C2,C3,C4,DL,DR,Perm2]|Rest2],
Ans,Result):-
!,
Ans1 is Ans+Perm*Perm2,
calc_perm(Rest1,Rest2,Ans1,Result).
calc_perm([[C1,C2,C3,C4,DL,DR,_]|Rest1],
[[C11,C22,C33,C44,DL11,DR11,Perm2]|Rest2],
Ans,Result):-
sort([[C1,C2,C3,C4,DL,DR],[C11,C22,C33,C44,DL11,DR11]],
[[C1,C2,C3,C4,DL,DR]|_]),
!,
calc_perm(Rest1,[[C11,C22,C33,C44,DL11,DR11,Perm2]|Rest2],
Ans,Result).
calc_perm(Memos1,[_|Memos2],Ans,Result):-
!,calc_perm(Memos1,Memos2,Ans,Result).
rev_list(Limit,Memo,[ReC1,ReC2,ReC3,ReC4,ReDL,ReDR,Perm]):-
member([C1,C2,C3,C4,DL,DR,Perm],Memo),
ReC4 is Limit-C1,
ReC3 is Limit-C2,
ReC2 is Limit-C3,
ReC1 is Limit-C4,
ReDL is Limit-DL,
ReDR is Limit-DR.
search(2,Limit,Memos,Ans):-
!,
msort(Memos,Memos1),
findall(M,rev_list(Limit,Memos1,M),RevMemo),
msort(RevMemo,RevMemo1),
calc_perm(Memos,RevMemo1,0,Ans),
write([Limit,Ans]),nl.
search(R,Limit,Memos,Result):-
!,
findall(M,next_calc(R,Limit,Memos,M),Memos1),
msort(Memos1,Memos2),
[Top|Memos3]=Memos2,
union_sum(Memos3,Top,Memos4),
R1 is R+1,
search(R1,Limit,Memos4,Result).
search_w(Ans):-
between(0,36,Limit),
search(0,Limit,[[0,0,0,0,0,0,1]],Ans).
main166:-
findall(Ans,search_w(Ans),Answers),
sum(Answers,AllSum),
write(AllSum).
*Problem 167 「Ulam数列について調べ上げよ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20167%E3%81%88
まずは問題を理解するところから始めてみよう。
何か問題に取り掛かるとき丁寧に列挙すること、丁寧に整理整頓するところから始めるべきだ。
数学の問題だって同じだ。
まずは問題が何を言っているか考えてみる。
∑U(2,2n+1)1000億,2=<n=<10ということは
ΣU(2,5)1000億
ΣU(2,7)1000億
ΣU(2,9)1000億
ΣU(2,11)1000億
ΣU(2,13)1000億
ΣU(2,15)1000億
ΣU(2,17)1000億
ΣU(2,19)1000億
ΣU(2,21)1000億
列挙することで問題が一段シンプルになった。
こういう問題はおそらく小さな値を見て規則性を発見するのが一番素朴な方法だと思う。
まずはそれを試してみよう。
2,5で考えてみる
2,5,7,9,12、、
すぐに気付くのはどの数字も2a+5b(a,bは0以上の自然数)で表現できるということだ。
小さい数から考え、2a+5bで表現できかつ何らかの条件と問題を言い換えてみればよさそうである。
互いに粗だからどのｎもn=2a+5bはaとｂはｎが決まるとa,bは一通りに決まる。
ただ1000億というのが多いし何らかの条件というのも同定するのが難しい。
線形でいっても実行時間がかかるのは明白、ひと工夫必要に思える。
このへんただいま試行錯誤中。
＊Problem 168 「数の循環」 †
自然数 142857 を考える. 最後の桁 (7) を一番前に持っていく, すなわち, 右に循環させると 714285 を得る.
714285 = 5 × 142857 が確認できる.
これは142857の珍しい性質を示している. つまり, 右に循環させた数の約数になっている.
10 < n < 10^100 の n について, この性質をもつ数の総和の下 5 桁を答えよ.
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20168
瞬きする間に答えが出るところまで到達したので満足。
右に循環させた数をNとするとありうる答えはNとｎは桁数が同じなのでN=an aは1～9までの数。
aと下一桁の二つが決まるとほかの桁はすべて一意に決まると気づけばこの問題は簡単。
sum([],0):-!.
sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
calc_next(Memo,KetaP,Mult,Result):-
member([FNum,LNum,Num,Up,Count,Ans],Memo),
Temp is Num*Mult+Up,
NextUp is Temp//10,
NextNum is Temp mod 10,
(KetaP=<6 ->
Ans1 is Ans+Count*NextNum*10^(KetaP-1);
Ans1 is Ans),
Result=[FNum,LNum,NextNum,NextUp,Count,Ans1].
last_calc(Memo,Mult,_,Ans):-
member([FNum,_,Num,Up,_,Ans],Memo),
Num>0,
FNum=:=Num*Mult+Up.
search(Limit,Limit,_,_,_):-!,fail.
search(Limit,KetaP,Mult,Memo,Result):-
KetaP>2,
findall(E,last_calc(Memo,Mult,Limit,E),Es),
sum(Es,Result).
search(Limit,KetaP,Mult,Memo,Result):-
!,
findall(M,calc_next(Memo,KetaP,Mult,M),Memo1),
KetaP1 is KetaP+1,
search(Limit,KetaP1,Mult,Memo1,Result).
seed(Seed,Mult):-
between(1,9,FNum),
LNum is (FNum*Mult) mod 10,
Seed =[FNum,LNum,FNum,0,1,FNum].
calc(Result):-
between(102,102,Limit),
between(1,9,Mult),
findall(Seed,seed(Seed,Mult),Seeds),
search(Limit,2,Mult,Seeds,Result).
main168:-findall(Sum,calc(Sum),Sums),
sum(Sums,Ans),
Ans1 is Ans mod 100000,
write(Ans1).
*Problem 169 「ある数を2のべき乗の和で表せる方法の数を探し当てよ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20169
簡単すぎる問題、何も考えない動的計画法で解けます。
% 答えは１２桁
to_bit2(0,[]):-!.
to_bit2(N,[B|Result]):-
N1 is N//2,
B is N mod 2,
to_bit2(N1,Result).
check(0,0,0).
check(0,0,1).
check(0,1,1).
check(1,0,0).
check(1,1,0).
check(1,1,1).
union_sum([],Set,[Set]):-!.
union_sum([[Add,P]|Memo],[Add,P1],Result):-
!,
P2 is P1+P,
union_sum(Memo,[Add,P2],Result).
union_sum([Set|Memo],Set1,[Set1|Result]):-
union_sum(Memo,Set,Result).
calc_nextB(Memo,Bit,[NextAdd,P]):-
member([Add,P],Memo),
check(Bit,Add,NextAdd).
calc_nextA(Memo,Bit,Memo4):-
findall(Set,calc_nextB(Memo,Bit,Set),Memo1),
msort(Memo1,Memo2),
[Top|Memo3]=Memo2,
union_sum(Memo3,Top,Memo4).
calc(Memo,[1]):-!,member([0,P0],Memo),member([1,P1],Memo),
Ans is P0+P1,
write(Ans).
calc(Memo,[Bit|Bits]):-
!,
calc_nextA(Memo,Bit,Memo4),
write([Bit,Memo4]),nl,
calc(Memo4,Bits).
main169:-
N is 10^25,
to_bit2(N,Bits),
calc([[0,1]],Bits).
