プロジェクトオイラーの問題を堀江伸一がPrologで解くページ。
問い51
http://projecteuler.net/problem=51
ピュアPrologには破壊的代入もなければ配列もstd::mapすらない。
リストで擬似的な11分木を構築して解いてみたけれど何か速度が出ない。
計算量の問題なのか、木を分解したり結合しなおしたりする作業が悪いのか？
こんな込み入ったことをするくらいなら単純な全部検証の方が早いとは思う。
まあリストで木構造を作ったらどの程度の速度が出るのか実験してみたかったというのが今回のコードです。
単純に全検証の方が早いよな絶対。
nums([a,0,1,2,3,4,5,6,7,8,9]).
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)).
add_a(X,[N|X]):-nums(A),
member(N,A).
add_list(List,[Y,X]):-
nums(A),
member(X,List),
member(Y,A).
create_tree(6,[0,0]):-!.
create_tree(N,Result):-
N1 is N+1,
findall(L,create_tree(N1,L),Re),
findall(L,add_list(Re,L),Result).
tree_add(_,M,[X,Y],[Re,Y1]):-
integer(X),
integer(Y),
!,
Y1 is Y+1,
((X=:=0;M<X)->
Re is M;
Re is X).
tree_add([],M,[[0,Rest]|Rest1],[[0,Result]|Rest1]):-
!,
tree_add([],M,Rest,Result).
tree_add([N|NumRest],M,[[N,Rest]|Rest1],[[N,Result]|Rest1]):-
!,
tree_add(NumRest,M,Rest,Result).
tree_add(NumList,M,[[N,Rest]|Rest1],[[N,Rest]|Result]):-
tree_add(NumList,M,Rest1,Result).
num_to_list(0,_,[]):-!.
num_to_list(N,R,[a|Result]):-
R =:= N mod 10,
N1 is N//10,
num_to_list(N1,R,Result).
num_to_list(N,R,[M|Result]):-
M is N mod 10,
N1 is N//10,
num_to_list(N1,R,Result).
check(List):-member(a,List),!.
perms(N,Result):-
between(0,9,R),
num_to_list(N,R,Result),
check(Result).
tree_adds(N,[],Tree):-!,search(N,Tree).
tree_adds(N,[L|List],Tree):-
N1 is N-1,
tree_add(L,N1,Tree,NextTree),
!,
tree_adds(N,List,NextTree).
tree_search([X,8]):-
integer(X),
!,
write([X,8]).
tree_search([[_,Rest]|Rest1]):-
member(X,[Rest,Rest1]),
tree_search(X).
search(1000000,Tree):-!,tree_search(Tree).
search(N,Tree):-
notPrime(N),
!,
N1 is N+1,
search(N1,Tree).
search(N,Tree):-
findall(L,perms(N,L),List),
N1 is N+1,
!,
tree_adds(N1,List,Tree).
main51:-create_tree(0,Tree),search(1,Tree).
問い52
この問題は1/7の割り算を行ったとき循環するまでに1/7,2/7,3/7,4/7,5/7,6/7という割り算が一度ずつ出てきてそれらは循環している。
だから1/7の循環節が一番小さな数字で答えとなる。
問題はこの数字が最小のものであるか？
これを証明しなくてはいけないが特に思いつかないので、結局いつも通り力技となりました
num_to_list(0,[]):-!.
num_to_list(N,[X|Result]):-
X is N mod 10,
N1 is N//10,
num_to_list(N1,Result).
eq_list([],[]):-!.
eq_list([X|Rest],List):-select(X,List,Rest1),eq_list(Rest,Rest1).
set(10,16).
set(100,166).
set(1000,1666).
set(10000,16666).
set(100000,166666).
check(N):-
num_to_list(N,List),
between(2,6,R),
N1 is N*R,
num_to_list(N1,List2),
( eq_list(List,List2)->
fail;!).
main52:-
set(A,B),
between(A,B,N),
not(check(N)),
!,
write(N).
問い53
そのまま計算して数えるだけ。
next_calc([1],[1]):-!.
next_calc([X,Y|Rest],[Z|Result]):-
Z1 is X+Y,
(1000000=<Z1 -> Z is 1000000;Z is Z1),
next_calc([Y|Rest],Result).
count100m([],0):-!.
count100m([X|Rest],Result):-
1000000=<X,!,count100m(Rest,Re),
Result is Re+1.
count100m([_|Rest],Result):-count100m(Rest,Result).
calc(100,_,Ans):-!,
write(Ans).
calc(N,List,Ans):-
N1 is N+1,
next_calc(List,NextList),
count100m(NextList,Re),
Ans1 is Ans + Re,
calc(N1,[1|NextList],Ans1).
main53:-calc(0,[1],0).
問い54
解法
自前の発想だと大きな役と照合していくだったけど、今回はリンク先の方が賢いのでリンク先を参考に記述（大体の問題で皆私より賢いのですがそこはおいといて）。
今回はカンニングしましたが、色々勉強になりました。
リンク先コードをPrologに翻訳してると物の考え方やPrologの欠点等結構勉強になった。
Prolog的に書こうと述語を小分けにするとコードサイズが膨らんでコードがスパゲッティ化するか、述語が増えて大きな役から試すのと変わらなくなったり記述量が増える。
またif文で書くとPrologは複雑なif文が記述がしにくいことなど色々勉強になった問題。
一気に書き上げて一か所の修正でアセプト。
to_list([],[]):-!.
to_list([[0,_]|Rest],Result):-
!,
to_list(Rest,Result).
to_list([[N,Num]|Rest],[Num|Result]):-
N1 is N-1,
to_list([[N1,Num]|Rest],Result).
calc_score(P,Score):-
sort(P,P1),
count_no(P1,[],NoCount),
type_count(P1,[],TypeCount),
length(NoCount,NoSize),
[[MaxCount,_]|_]=NoCount,
[[Min,_],_,_,_,[Max,_]]=P1,
Diff is Max-Min,
(NoSize=:=5 ->
(Diff>4 -> (TypeCount>1 -> T is 0 ; T is 5);
(TypeCount>1 -> T is 4;
(Min<10 -> T is 8 ;T is 9)));
(NoSize=:=4 -> T is 1;
(NoSize=:=3 -> (MaxCount=:=3->T is 3;T is 2)
;
(NoSize=:=2->(MaxCount=:=4 -> T is 7;T is 6)
;T is 0)))),
to_list(NoCount,Re1),
Score=[T|Re1].
read_one_set(5,[]):-!.
read_one_set(N,[[Num1,Type]|Result]):-
N1 is N+1,
get0(Num),
get0(Type),
get0(_),
!,
to_num(Num,Num1),
read_one_set(N1,Result).
p1Win([X|_],[Y|_]):-
Y<X,!.
p1Win([X|_],[Y|_]):-
X<Y,!,fail.
p1Win([_|RestP1],[_|RestP2]):-
p1Win(RestP1,RestP2).
reads(Result,Result):-peek_char(-1).
reads(Count,Result):-
read_one_set(0,P1),
read_one_set(0,P2),
calc_score(P1,P11),
calc_score(P2,P22),
(p1Win(P11,P22)->Count1 is Count+1;Count1 is Count),
reads(Count1,Result).
main54:-seen,see('e54.txt'),reads(0,Ans),write(Ans).
問い55
この計算必ず単調増加します。
配列があれば10000～1へ調べていきます。
ある数になった時そのあと何回以内で計算が終わるか、計算が終わらないか配列にメモっておけば５０回以内で計算が終わるかのチェックで少しだけ計算量を下げることが出来ます。
そこまでしなくてもすぐに答えが出るので気にしませんでした。
rev_num(0,Sum,Sum):-!.
rev_num(N,Sum,Result):-
Sum1 is Sum *10+ (N mod 10),
N1 is N //10,
rev_num(N1,Sum1,Result).
calc(50,_):-!.
calc(T,N):-
rev_num(N,0,N1),
T1 is T+1,
N2 is N+N1,
rev_num(N2,0,N3),
(N2=:=N3->
!,fail;
calc(T1,N2)).
search(N):-
between(1,10000,N),
calc(0,N).
main55:-findall(N,search(N),List),length(List,Len),write(Len).
問い56
そのまま計算するだけ。
keta_sum(0,0):-!.
keta_sum(N,Result):-
N1 is N//10,
M is N mod 10,
keta_sum(N1,Re),
Result is Re+M.
all_perm([Sum,A,B]):-
between(1,99,A),
between(1,99,B),
P is A^B,
keta_sum(P,Sum).
main56:-findall(L,all_perm(L),List),sort(List,List1),write(List1).
問い57
√２の連分数展開を一項ずつ計算していく時、１０００項目までで分子の桁数が分母の桁数を超える回数を答えよ。
計算速度を稼ぐために前の項の計算結果を再利用するくらいで特に工夫はしてません。
gcd(0, B, B).
gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G).
add(U1,D1,U2,D2,U4,D4):-
U3 is U1*D2+U2*D1,
D3 is D1*D2,
gcd(U3,D3,M),
U4 is U3//M,
D4 is D3//M.
isBig(0,0):-!,fail.
isBig(_,0):-!.
isBig(U,D):-U1 is U//10,D1 is D//10,isBig(U1,D1).
calc(1000,_,_,_,_):-!,fail.
calc(_,U,D,U1,D1):-
add(U,D,1,1,U1,D1),
isBig(U1,D1).
calc(N,U,D,U2,D2):-
add(U,D,2,1,U1,D1),
N1 is N+1,
calc(N1,D1,U1,U2,D2).
search(U):-
calc(0,1,2,U,_).
main57:-findall(U,search(U),List),length(List,Len),write(Len).
問い58
この問題はミラーラビン法というものを使うと高速に解けるらしい。
読んだけど一読しただけではよくわからなかったのでいつもの素朴な試し割で素数判定。
遅いけど一応答えは出てるからいいかな（向上心という言葉はないのかな俺？）。
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)).
search(_,ADD,_,PC,Count):-
PC>0,
PC1 is PC*10,
PC1<Count,
!,
Ans is ADD-1,
write([Ans,PC,Count]).
search(N,ADD,[],PC,Count1):-
!,
ADD1 is ADD+2,
search(N,ADD1,[ADD,ADD,ADD,ADD],PC,Count1).
search(N,ADD,[A1|Rest],PC,Count):-
!,
Count1 is Count+1,
N1 is N+A1,
( isPrime(N1)->
PC1 is PC+1;
PC1 is PC),
search(N1,ADD,Rest,PC1,Count1).
main58:-search(1,4,[2,2,2,2],0,1).
問い60
この問題ただいま効率的な実装がないか思考中。
今回どうしてもPrologで効率的に実装する方法がわからなくてC++で試しに実装。
試しがきなので関数の機能単位の分割とかは無視。
コードはちゃんと動いたけれどC++でコード実行時間４０秒は遅いかな。
素数判定の部分に大量の時間を取られてるのが問題だからここを改善すればいいとおもうのだけど、ちょっと遅すぎる？
C++ならコードを実行したら数秒もしない、理想的には一瞬のうちに答えが出るべきだしコード行数も１００行は多すぎる。
俺今最高にカッコ悪いコードを書いたぜって感じ。
まず組み合わせる素数の上限（今回は３万）を定めこの範囲内で答えがないか探します。
でてきた答えが３万以下ならこれより大きな素数を使った答えはないはずです。
コードの考え方は簡単。
5＝3+2です。
まず3万以下の全ての素数で、素数Piと結合して素数になる素数のリストLiを求めます。
Liの中の2つの要素a,bを選び、aとｂが結合可能か確認します。
L(Pi)を素数Piと結合可能な素数のリストを求めるものだとします。
L(Pi)∧L(a)∧L(b）となるリストを求め、このリストから2つ選びc,dとします。
ｃとdが結合可能ならPi,a,b,c,dは全部結合可能と判明します。
簡単な料理のレシピより短い発想です。
実際コーディングも簡単でした。
子供でも作れる簡単な料理のレシピを見ながら、これを食べれるレベルで作る方が私にとってはよほど難しいです。
試しがきなので書き終えてから多少無駄な処理を書いてたことに気付いたのは秘密です。
計算量大して変わらないから気にしません。
#include<stdio.h>
#include<map>
#include<set>
#include<vector>
#include<stdlib.h>
#include<iostream>
std::vector<__int64> primes;
const int up=100*1000;
const __int64 searchMax=30000;
void setPrime(){
bool isPrime[up];
memset(isPrime,true,sizeof(isPrime));
primes.push_back(2);
for(int i=3;i<up;i+=2){
if(isPrime[i]==false)continue;
primes.push_back(i);
int add=i*2;
for(int j=i*3;j<up;j+=add){
isPrime[j]=false;
}
}
}
bool isPrime(__int64 num){
for(int i=0;primes[i]*primes[i]<=num;i++){
if(num%primes[i]==0)return false;
}
return true;
}
__int64 numUnion(__int64 a,__int64 b,__int64 ten){
if(b % ten==b){
return a*ten+b;
}else{
return numUnion(a,b,ten*10);
}
}
std::set<__int64> numSet(std::set<__int64>& A,std::set<__int64>& B,std::set<__int64>& C){
std::set<__int64> result;
std::set<__int64>::iterator it;
__int64 a;
for(it=A.begin();it!=A.end();it++){
a=(*it);
if(B.find(a)!=B.end()&&C.find(a)!=C.end()){
result.insert(a);
}
}
return result;
}
int main(){
std::map<__int64,std::set<__int64> > cons;
std::map<__int64,std::set<__int64> >::iterator it;
std::set<__int64>::iterator itS,itS2,itS3,itS4;
setPrime();
for(int i=0;primes[i]<searchMax;i++){
__int64 a=primes[i];
for(int j=i+1;primes[j]<searchMax;j++){
__int64 b=primes[j];
if(isPrime(numUnion(a,b,1))&&isPrime(numUnion(b,a,1))){
cons[a].insert(b);
cons[b].insert(a);
}
}
std::cout<<a<<" "<<cons[a].size()<<" ";
}
printf("ok");
__int64 Ans=searchMax*5;
for(it=cons.begin();it!=cons.end();it++){
__int64 a,b,c,d;
for(itS=(*it).second.begin();itS!=(*it).second.end();itS++){
itS2=itS;
itS2++;
a=(*itS);
for(;itS2!=(*it).second.end();itS2++){
b=(*itS2);
if(cons[a].find(b)!=cons[a].end()&&cons[b].find(a)!=cons[b].end()){
std::set<__int64> A=numSet((*it).second,cons[a],cons[b]);
for(itS3=A.begin();itS3!=A.end();itS3++){
itS4=itS3;
itS4++;
for(;itS4!=A.end();itS4++){
c=(*itS3);
d=(*itS4);
if(cons[c].find(d)!=cons[c].end()&&cons[d].find(c)!=cons[d].end()){
if(Ans>a+b+c+d+(*it).first){
Ans=a+b+c+d+(*it).first;
}
}
}
}
}
}
}
}
std::cout<<"\nans="<<Ans<<"\n";
}
最終更新：2013年08月08日 17:14
