「prolog勉強プロジェクトオイラー41～50」の編集履歴（バックアップ）一覧はこちら
追加された行は緑色になります。
削除された行は赤色になります。
プロジェクトオイラーの問題を堀江伸一がPrologで解くページ。
私は整数論の知識がほぼ0なのでコードは基本素朴な実装しかありません。
Prologプログラムもまだ休んでた時期を除くと２カ月になりません。
というわけで９０年代Prolog黄金期を体験した熟練Prologプログラマから見たらおしかりを受けるようなコードだと思いますがまあ素人はこの程度の発想に到達するのだと思ってお読みいただければと思います。
*問い41
パンデジタル数の問題。
全部試し割。
全ての桁を足したものが３の倍数になると素数にならないので7桁か4桁しか検討する必要がない。
sets([7,6,5,4,3,2,1]).
sets([4,3,2,1]).
nums([],Result,Result):-!.
nums(Nums,Num,Result):-select(X,Nums,NumRest),
Num1 is Num*10+X,
nums(NumRest,Num1,Result).
not_prime(N):-N<2.
not_prime(N):-
sqrt(N,D),
D1 is floor(D),
between(2,D1,R),
0=:=N mod R.
all_get(N):-
sets(A),
nums(A,0,N),
not(not_prime(N)).
main41:-
findall(N,all_get(N),List),sort(List,List1),write(List1).
*問い42
ファイルを読み込み名前を番号に変換して、それが三角数であるかを判定してその個数を数える問題。
コードは難しく見えるがこの問題やってることを考えると小学生みたいな気分になる。
last(-1).
spliter(34).
read_name(Names,Name):-
get0(C),
(spliter(C)->
karayomi_read([Name|Names]);
Name1 is Name+C-"A"+1,
read_name(Names,Name1)).
karayomi_read(Names):-get0(C),karayomi(Names,C).
karayomi(Names,C):-last(C),!,count(Names,0).
karayomi(Names,C):-spliter(C),!,read_name(Names,0).
karayomi(Names,_):-karayomi_read(Names).
isDelta(T,_):-T<0,!,fail.
isDelta(0,_):-!.
isDelta(T,Dell):-T1 is T-Dell,Dell1 is Dell+1,isDelta(T1,Dell1).
count([],Ans):-!,write(Ans).
count([Num|Rest],Ans):-isDelta(Num,1),!,Ans1 is Ans+1,count(Rest,Ans1).
count([_|Rest],Ans):-count(Rest,Ans).
main42:-seen,see('e42.txt'),karayomi_read([]).
*問い43
枝刈しながら求めるだけ。
行列演算を解くふうみな何かがあるのかも？
to_num([],Result,Result):-!.
to_num([D|Rest],N,Result):-N1 is N*10+D,to_num(Rest,N1,Result).
search([D1,D2],_,Ans,[],Result):-Ans1=[D2,D1|Ans],
reverse(Ans1,Ans2),
to_num(Ans2,0,Result).
search([D1,D2],Nums,Ans,[P|Ps],Result):-
select(D3,Nums,Rest),
to_num([D1,D2,D3],0,Num),
0=:=Num mod P,
search([D2,D3],Rest,[D1|Ans],Ps,Result).
search_w(N):-
select(D1,[0,1,2,3,4,5,6,7,8,9],Rest),
0<D1,
select(D2,Rest,Rest1),
search([D1,D2],Rest1,[],[1,2,3,5,7,11,13,17],N).
sum([],0).
sum([X|Rest],Result):-sum(Rest,Re),Result is Re + X.
main43:-
findall(N,search_w(N),List),write(List),sum(List,Ans),write(Ans).
*問い44
http://projecteuler.net/problem=44
条件を見たすPn、Pmを考えた時
m=n-aとし、5角数ｄを考えると
d=n(3n-1)/2-(n-a)(3(n-a)-1)/2
を展開して整理して
① 2d=6an-3a^2-a
②a=(6n-1 (+ or -)sqrt((6n-1)^2-24d))/6
c=6n-1とすると
c^2-24d=b^2
c^2-b^2=24d
(c+b)(c-b)=24d
ここで
1からint(sqrt(24d))の範囲で24dを試し割して余りが出なければ
c-bをその割った数と確定ここから割り算のc+bが反対に出てくるので、ここでｃが確定する。
この確定したｃが条件①②を満たすか？
とやるとまあまあの計算量。
まだ計算量が多いがそれなりの時間で答えが出るのでまあいいかな。
多分条件を満たすc^2-b^2を連分数展開とかで求めるのではないかと思う。
以下効率の悪いコード。
手続き的で何かProlog的でないような気がします。
isPenta(N):-
N>0,
N1 is N*24+1,
sqrt(N1,N2),
NN is floor(N2),
NN =:=N2,
0 =:= (NN + 1) mod 6 .
penta(N,Result):-Result is (N*(3*N-1))//2.
calc(N,A,Result):-
T is 6*N-1+A,
T>0,
0=:=T mod 6,
Result is T //6.
calc(N,A,Result):-
T is 6*N-1-A,
T>0,
0=:=T mod 6,
Result is T //6.
calcN(D,A,Result):-
A>0,
T is 2*D+3*A*A+A,
0=:=T mod (6*A),
Result is T//(6*A).
main44:-
between(1,10000000,Seed),
penta(Seed,D),
D24 is D*24,
sqrt(D24,Up),
Up1 is floor(Up),
between(1,Up1,R),
0 =:= D24 mod R,
L is D24 // R,
0=:=(L+R) mod 2,
B is (L+R)//2,
0=:=(B+1) mod 6,
N is (B+1)//6,
A is B*B-D24,
sqrt(A,A1),
A2 is floor(A1),
A2=:=A1,
calc(N,A2,A3),
calcN(D,A3,N3),
N4 is N3-A3,
N4>0,
L1 is N3*(3*N3-1),
L2 is N4*(3*N4-1),
0=:=L1 mod 2,
0=:=L2 mod 2,
L11 is L1 // 2,
L22 is L2 // 2,
L33 is L11+L22,
isPenta(L33),
Ans is L11-L22,
write([L22,L11,Ans]).
*問い45
5角数と6角数が同じ値になるもののうち下から３つ目の数を答えよ。
penta(N,Result):-Result is (N*(3*N-1))//2.
hexa(N,Result):-Result is N*(2*N-1).
search(P,H):-
penta(P,PP),
hexa(H,HP),
PP=:=HP,!,
write(HP).
search(P,H):-
penta(P,PP),
hexa(H,HP),
PP<HP,!,
P1 is P+1,
search(P1,H).
search(P,H):-
penta(P,PP),
hexa(H,HP),
HP<PP,!,
H1 is H+1,
search(P,H1).
main45:-search(166,144).
*問い46
下から全部試すだけです
notPrime(N):-N<2,!.
notPrime(N):-
sqrt(N,N1),
N2 is floor(N1),
between(2,N2,N3),
0=:=N mod N3,
!.
isPrime(N):-not(notPrime(N)).
check(N):-
sqrt(N,N1),
N2 is floor(N1),
between(1,N2,N3),
N4 is N-N3*N3*2,
N4>0,
isPrime(N4),
!.
search(N):-
(isPrime(N);0=:=N mod 2),
!,
N1 is N+1,
search(N1).
search(N):-
(check(N)->
N1 is N+1,
!,
search(N1);
write(N)).
main46:-search(35).
*問い47
エレガントな解法はリンク先をご覧ください。
http://tsumuji.cocolog-nifty.com/tsumuji/2010/09/project-euler-p.html
ハクセル的発想というものを知らなければまず思いつけない重要な発想です。
こちらは私の書いた素朴そのもののコード。
下から愚直に数えてるだけ。
流石にこれでは恥ずかしいので４種類の素因数からなる合成数を小さい方から生成するルールでも考えてみます。
count_prime(1,_,List,Result):-!,sort(List,List1),length(List1,Result).
count_prime(N,R,List,Result):-
N<R*R,!,
count_prime(1,R,[N|List],Result).
count_prime(N,R,List,Result):-
0=:=N mod R,
!,
N1 is N//R,
count_prime(N1,R,[R|List],Result).
count_prime(N,R,List,Result):-
R1 is R+1,
!,
count_prime(N,R1,List,Result).
search(N,R):-
count_prime(N,2,[],4),
!,
N1 is N+1,
(R>2->
Ans is N-3,
write(Ans);
R1 is R+1,
search(N1,R1)).
search(N,_):-
N1 is N+1,
search(N1,0).
main47:-search(2,0).
*問い48
Σn^n(n=1..1000)の下１０ケタを答えよ
無制限の数整数が標準でサポートされている処理系ならこれは何の問題もない問題。
sum(1001,Sum,Result):-!,Result is Sum mod (10^10).
sum(N,Sum,Result):-Sum1 is Sum+N^N,N1 is N+1,sum(N1,Sum1,Result).
main48:-sum(1,0,Ans),write(Ans).
*問い49
３つの数列が４桁かつ等差数列かつ、どの数字も各桁を並べ替えたものとなる数列は２種類ある。
片方は問題文で提示されているのでもう片方を答えよ。
初項Nを各桁の数字のリストになおして並べ替えたものN1を1つ求める。
N<N1を確認してこれから３つ目の数字を２*N1-Nから求める。
この３つが等差数列になってて全て素数ならそれが答えです。
５桁で試してみてもこのコードはまあまあの速度で動きます。
notPrime(N):-N<2,!.
notPrime(N):-
sqrt(N,N1),
N2 is floor(N1),
between(2,N2,N3),
0=:=N mod N3,
!.
isPrime(N):-not(notPrime(N)).
toNum([],0):-!.
toNum(List,Result):-select(X,List,Rest),toNum(Rest,Num),Result is Num*10+X.
toList(0,[]):-!.
toList(Num,[X|Rest]):-X is Num mod 10,Num1 is Num //10,toList(Num1,Rest).
eqlist([],[]):-!.
eqlist([X|Rest],List):-select(X,List,Rest1),eqlist(Rest,Rest1).
search([N,N1,N2]):-between(1000,9997,N),
isPrime(N),
toList(N,List),
toNum(List,N1),
N<N1,
isPrime(N1),
N2 is 2*N1-N,
N2<10000,
toList(N2,List1),
eqlist(List,List1),
isPrime(N2).
main49:-findall(Set,search(Set),List),sort(List,List1),write(List1).
*問い50
連続した素数の和で表せる１００万以下の素数を考える。
この時、足す素数の個数が最長となる数を求めよ。
2～5万までの間にある素数リストを得る（問題で1000いかに２１項があると分かってるので答えが５万以上の素数から始まることはない）。
先頭を１つ消した素数リストとこれをたし（一つずらして足す）100万以下のもの求めてsearch述語に再代入する。
後はこれを繰り返せば答えとなります。
つまり全部検証してます。
配列があればもうすこし楽かつかなり高速に実装できるのですがピュアPrologには配列がなくリストしかないので処理はかなり遅いです。
配列があればもっと早いのですが。
notPrime(N):-N<2,!.
notPrime(N):-
sqrt(N,N1),
N2 is floor(N1),
between(2,N2,N3),
0=:=N mod N3,
!.
isPrime(N):-not(notPrime(N)).
prime_list(N):-
between(2,50000,N),
isPrime(N).
adds([],_,[]):-!.
adds(_,[],[]):-!.
adds([X|Rest],[Y|Rest1],[Z|Result]):-
Z is X+Y,Z<1000000,!,adds(Rest,Rest1,Result).
adds([_|Rest],[_|Rest1],Result):-!,adds(Rest,Rest1,Result).
prime_search([X|_],X):-isPrime(X),!.
prime_search([_|Rest],Result):-prime_search(Rest,Result).
search(_,[],Ans):-!,write(Ans).
search([],_,Ans):-!,write(Ans).
search(List,[_|List1],Temp):-
adds(List,List1,Next),
!,
(prime_search(Next,Result)->
search(Next,List1,Result);
search(Next,List1,Temp)).
プロジェクトオイラーの問題を堀江伸一がPrologで解くページ。
私は整数論の知識がほぼ0なのでコードは基本素朴な実装しかありません。
Prologプログラムの勉強も休んでた時期を除くと２カ月になりません。
というわけで９０年代Prolog黄金期を体験した熟練Prologプログラマから見たらおしかりを受けるようなコードだと思いますがまあ素人はこの程度の発想に到達するのだと思ってお読みいただければと思います。
*問い41
パンデジタル数の問題。
全部試し割。
全ての桁を足したものが３の倍数になると素数にならないので7桁か4桁しか検討する必要がない。
sets([7,6,5,4,3,2,1]).
sets([4,3,2,1]).
nums([],Result,Result):-!.
nums(Nums,Num,Result):-select(X,Nums,NumRest),
Num1 is Num*10+X,
nums(NumRest,Num1,Result).
not_prime(N):-N<2.
not_prime(N):-
sqrt(N,D),
D1 is floor(D),
between(2,D1,R),
0=:=N mod R.
all_get(N):-
sets(A),
nums(A,0,N),
not(not_prime(N)).
main41:-
findall(N,all_get(N),List),sort(List,List1),write(List1).
*問い42
ファイルを読み込み名前を番号に変換して、それが三角数であるかを判定してその個数を数える問題。
コードは難しく見えるがこの問題やってることを考えると小学生みたいな気分になる。
last(-1).
spliter(34).
read_name(Names,Name):-
get0(C),
(spliter(C)->
karayomi_read([Name|Names]);
Name1 is Name+C-"A"+1,
read_name(Names,Name1)).
karayomi_read(Names):-get0(C),karayomi(Names,C).
karayomi(Names,C):-last(C),!,count(Names,0).
karayomi(Names,C):-spliter(C),!,read_name(Names,0).
karayomi(Names,_):-karayomi_read(Names).
isDelta(T,_):-T<0,!,fail.
isDelta(0,_):-!.
isDelta(T,Dell):-T1 is T-Dell,Dell1 is Dell+1,isDelta(T1,Dell1).
count([],Ans):-!,write(Ans).
count([Num|Rest],Ans):-isDelta(Num,1),!,Ans1 is Ans+1,count(Rest,Ans1).
count([_|Rest],Ans):-count(Rest,Ans).
main42:-seen,see('e42.txt'),karayomi_read([]).
*問い43
枝刈しながら求めるだけ。
行列演算を解くふうみな何かがあるのかも？
to_num([],Result,Result):-!.
to_num([D|Rest],N,Result):-N1 is N*10+D,to_num(Rest,N1,Result).
search([D1,D2],_,Ans,[],Result):-Ans1=[D2,D1|Ans],
reverse(Ans1,Ans2),
to_num(Ans2,0,Result).
search([D1,D2],Nums,Ans,[P|Ps],Result):-
select(D3,Nums,Rest),
to_num([D1,D2,D3],0,Num),
0=:=Num mod P,
search([D2,D3],Rest,[D1|Ans],Ps,Result).
search_w(N):-
select(D1,[0,1,2,3,4,5,6,7,8,9],Rest),
0<D1,
select(D2,Rest,Rest1),
search([D1,D2],Rest1,[],[1,2,3,5,7,11,13,17],N).
sum([],0).
sum([X|Rest],Result):-sum(Rest,Re),Result is Re + X.
main43:-
findall(N,search_w(N),List),write(List),sum(List,Ans),write(Ans).
*問い44
http://projecteuler.net/problem=44
条件を見たすPn、Pmを考えた時
m=n-aとし、5角数ｄを考えると
d=n(3n-1)/2-(n-a)(3(n-a)-1)/2
を展開して整理して
① 2d=6an-3a^2-a
②a=(6n-1 (+ or -)sqrt((6n-1)^2-24d))/6
c=6n-1とすると
c^2-24d=b^2
c^2-b^2=24d
(c+b)(c-b)=24d
ここで
1からint(sqrt(24d))の範囲で24dを試し割して余りが出なければ
c-bをその割った数と確定ここから割り算のc+bが反対に出てくるので、ここでｃが確定する。
この確定したｃが条件①②を満たすか？
とやるとまあまあの計算量。
まだ計算量が多いがそれなりの時間で答えが出るのでまあいいかな。
多分条件を満たすc^2-b^2を連分数展開とかで求めるのではないかと思う。
以下効率の悪いコード。
手続き的で何かProlog的でないような気がします。
isPenta(N):-
N>0,
N1 is N*24+1,
sqrt(N1,N2),
NN is floor(N2),
NN =:=N2,
0 =:= (NN + 1) mod 6 .
penta(N,Result):-Result is (N*(3*N-1))//2.
calc(N,A,Result):-
T is 6*N-1+A,
T>0,
0=:=T mod 6,
Result is T //6.
calc(N,A,Result):-
T is 6*N-1-A,
T>0,
0=:=T mod 6,
Result is T //6.
calcN(D,A,Result):-
A>0,
T is 2*D+3*A*A+A,
0=:=T mod (6*A),
Result is T//(6*A).
main44:-
between(1,10000000,Seed),
penta(Seed,D),
D24 is D*24,
sqrt(D24,Up),
Up1 is floor(Up),
between(1,Up1,R),
0 =:= D24 mod R,
L is D24 // R,
0=:=(L+R) mod 2,
B is (L+R)//2,
0=:=(B+1) mod 6,
N is (B+1)//6,
A is B*B-D24,
sqrt(A,A1),
A2 is floor(A1),
A2=:=A1,
calc(N,A2,A3),
calcN(D,A3,N3),
N4 is N3-A3,
N4>0,
L1 is N3*(3*N3-1),
L2 is N4*(3*N4-1),
0=:=L1 mod 2,
0=:=L2 mod 2,
L11 is L1 // 2,
L22 is L2 // 2,
L33 is L11+L22,
isPenta(L33),
Ans is L11-L22,
write([L22,L11,Ans]).
*問い45
5角数と6角数が同じ値になるもののうち下から３つ目の数を答えよ。
penta(N,Result):-Result is (N*(3*N-1))//2.
hexa(N,Result):-Result is N*(2*N-1).
search(P,H):-
penta(P,PP),
hexa(H,HP),
PP=:=HP,!,
write(HP).
search(P,H):-
penta(P,PP),
hexa(H,HP),
PP<HP,!,
P1 is P+1,
search(P1,H).
search(P,H):-
penta(P,PP),
hexa(H,HP),
HP<PP,!,
H1 is H+1,
search(P,H1).
main45:-search(166,144).
*問い46
下から全部試すだけです
notPrime(N):-N<2,!.
notPrime(N):-
sqrt(N,N1),
N2 is floor(N1),
between(2,N2,N3),
0=:=N mod N3,
!.
isPrime(N):-not(notPrime(N)).
check(N):-
sqrt(N,N1),
N2 is floor(N1),
between(1,N2,N3),
N4 is N-N3*N3*2,
N4>0,
isPrime(N4),
!.
search(N):-
(isPrime(N);0=:=N mod 2),
!,
N1 is N+1,
search(N1).
search(N):-
(check(N)->
N1 is N+1,
!,
search(N1);
write(N)).
main46:-search(35).
*問い47
エレガントな解法はリンク先をご覧ください。
http://tsumuji.cocolog-nifty.com/tsumuji/2010/09/project-euler-p.html
ハクセル的発想というものを知らなければまず思いつけない重要な発想です。
こちらは私の書いた素朴そのもののコード。
下から愚直に数えてるだけ。
流石にこれでは恥ずかしいので４種類の素因数からなる合成数を小さい方から生成するルールでも考えてみます。
count_prime(1,_,List,Result):-!,sort(List,List1),length(List1,Result).
count_prime(N,R,List,Result):-
N<R*R,!,
count_prime(1,R,[N|List],Result).
count_prime(N,R,List,Result):-
0=:=N mod R,
!,
N1 is N//R,
count_prime(N1,R,[R|List],Result).
count_prime(N,R,List,Result):-
R1 is R+1,
!,
count_prime(N,R1,List,Result).
search(N,R):-
count_prime(N,2,[],4),
!,
N1 is N+1,
(R>2->
Ans is N-3,
write(Ans);
R1 is R+1,
search(N1,R1)).
search(N,_):-
N1 is N+1,
search(N1,0).
main47:-search(2,0).
*問い48
Σn^n(n=1..1000)の下１０ケタを答えよ
無制限の数整数が標準でサポートされている処理系ならこれは何の問題もない問題。
sum(1001,Sum,Result):-!,Result is Sum mod (10^10).
sum(N,Sum,Result):-Sum1 is Sum+N^N,N1 is N+1,sum(N1,Sum1,Result).
main48:-sum(1,0,Ans),write(Ans).
*問い49
３つの数列が４桁かつ等差数列かつ、どの数字も各桁を並べ替えたものとなる数列は２種類ある。
片方は問題文で提示されているのでもう片方を答えよ。
初項Nを各桁の数字のリストになおして並べ替えたものN1を1つ求める。
N<N1を確認してこれから３つ目の数字を２*N1-Nから求める。
この３つが等差数列になってて全て素数ならそれが答えです。
５桁で試してみてもこのコードはまあまあの速度で動きます。
notPrime(N):-N<2,!.
notPrime(N):-
sqrt(N,N1),
N2 is floor(N1),
between(2,N2,N3),
0=:=N mod N3,
!.
isPrime(N):-not(notPrime(N)).
toNum([],0):-!.
toNum(List,Result):-select(X,List,Rest),toNum(Rest,Num),Result is Num*10+X.
toList(0,[]):-!.
toList(Num,[X|Rest]):-X is Num mod 10,Num1 is Num //10,toList(Num1,Rest).
eqlist([],[]):-!.
eqlist([X|Rest],List):-select(X,List,Rest1),eqlist(Rest,Rest1).
search([N,N1,N2]):-between(1000,9997,N),
isPrime(N),
toList(N,List),
toNum(List,N1),
N<N1,
isPrime(N1),
N2 is 2*N1-N,
N2<10000,
toList(N2,List1),
eqlist(List,List1),
isPrime(N2).
main49:-findall(Set,search(Set),List),sort(List,List1),write(List1).
*問い50
連続した素数の和で表せる１００万以下の素数を考える。
この時、足す素数の個数が最長となる数を求めよ。
2～5万までの間にある素数リストを得る（問題で1000いかに２１項があると分かってるので答えが５万以上の素数から始まることはない）。
先頭を１つ消した素数リストとこれをたし（一つずらして足す）100万以下のもの求めてそこから素数を探してこれを記憶しsearch述語に再代入します。
後はこれを繰り返してリストが空になった時。
[2,3,5,7]なら
[2,3,5,7]
[3,5,7]
[5,8,12]
というリストを得ます、計算後 5,8,12に素数があればリストを先頭から調べて最初に見つかった素数5を暫定解として次へ行くわけです。
つまり全部検証してます。
配列があればもうすこし楽かつかなり高速に実装できるのですがピュアPrologには配列がなくリストしかないので処理はかなり遅いです。
配列があればもっと早いのですが。
notPrime(N):-N<2,!.
notPrime(N):-
sqrt(N,N1),
N2 is floor(N1),
between(2,N2,N3),
0=:=N mod N3,
!.
isPrime(N):-not(notPrime(N)).
prime_list(N):-
between(2,50000,N),
isPrime(N).
adds([],_,[]):-!.
adds(_,[],[]):-!.
adds([X|Rest],[Y|Rest1],[Z|Result]):-
Z is X+Y,Z<1000000,!,adds(Rest,Rest1,Result).
adds([_|Rest],[_|Rest1],Result):-!,adds(Rest,Rest1,Result).
prime_search([X|_],X):-isPrime(X),!.
prime_search([_|Rest],Result):-prime_search(Rest,Result).
search(_,[],Ans):-!,write(Ans).
search([],_,Ans):-!,write(Ans).
search(List,[_|List1],Temp):-
adds(List,List1,Next),
!,
(prime_search(Next,Result)->
search(Next,List1,Result);
search(Next,List1,Temp)).
