「prolog勉強プロジェクトオイラー161~170」の編集履歴(バックアップ)一覧はこちら
追加された行は緑色になります。
削除された行は赤色になります。
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。
*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
全部の線分の交差判定をしました、遅いです。
交点は分数座標で保持することで同じ点かを判断します。
グローバルスタックが足らないのでコードの実行時に拡張してアセプト。
Prologには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 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もn=2a+5bはaとbはnが決まると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は桁数が同じなので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
簡単すぎる問題、何も考えない動的計画法で解けます。
% 答えは12桁
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).
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。
*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
全部の線分の交差判定をしました、遅いです。
交点は分数座標で保持することで同じ点かを判断します。
グローバルスタックが足らないのでコードの実行時に拡張してアセプト。
Prologには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 「十字」 †
とりあえず答えが見たかったのでC++で実装。
対称性を利用して上2行を計算し終えたら下2行はその反転として考えて3.3秒を達成。
C++は早いうえに新しいパソコンだからこれは卑怯くさいタイムな気もする。
今のところstd::mapをPrologに移植する手段を思いつかない。
意外と全探索で何とかなったりして?
#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);
}
*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もn=2a+5bはaとbはnが決まると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は桁数が同じなので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
簡単すぎる問題、何も考えない動的計画法で解けます。
% 答えは12桁
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).