プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。 多分このページをIT系の仕事している私の弟に見られたらPrologなんて仕事場でほとんど使われていない言語なんか勉強してどうするの？ とかプロジェクトオイラーのサイトの問題が解けても実務にはつながらないだろ、もっと役立つことを勉強しろとか言われる気はする。 言われたら反論が出来ないのがこまる。 問い71 http://projecteuler.net/problem=71 ファレイ数列の定義そのままです。 ファレイ数列が隙間を埋めていくというのは分かりますが、出てくる数字が既約というのは何だか不思議。 整数論は初歩の話でも面白いことがおおいですね、私の性根はこういう初等的な数学の話を知ることが好きなおとなしい人間だったりします。 farey(U0,D0,_,D1):- T is D0+D1, T>=1000000, !, write([U0,D0]). farey(U0,D0,U1,D1):- U2 is U0+U1, D2 is D0+D1, farey(U2,D2,U1,D1). main71:- farey(2,5,3,7). *問い72 Counting fractions この問題は一番素朴な方法は試し割でPhi関数を実装、これだと遅いですがメモリ使用量は抑えられますしまあいいかなと。 試し割は無駄な処理が入りますが試し割の無駄を省こうとすると素数リストが必要となり素数リスト関係の計算コストやメモリ使用量が上がります。 整数論に通じてない素朴な身としては試し割が一番シンプルでいいという結論になってしまいました。 なにやら高速に計算する関数があるようですが原理がよくわからないので飛ばしています。 is_prime(N):- between(2,N,N1), N2 is N1*N1, (N<N2->!,true; 0=:=N mod N1,!,fail). prime_list(N):-between(2,1100,N),is_prime(N). calc(1000000,Ans,_):-!,write(Ans). calc(N,Ans,Ps):-N1 is N+1, calc_phi(N,Phi,Ps), Ans1 is Ans + Phi, calc(N1,Ans1,Ps). main72:-findall(N,prime_list(N),Ps), calc(1,0,Ps). div(N,P,Result):-0=:=N mod P,!, N1 is N//P,div(N1,P,Result). div(N,_,N):-!. calc_phi(N,Result,[P|_]):- P2 is P*P, N<P2, !, (N>1 -> Result is N-1; Result is 1). calc_phi(N,Result,[P|Ps]):-!, div(N,P,N1), calc_phi(N1,Re,Ps), (N>N1-> Result is Re*(N//(N1*P))*(P-1); Result is Re). *問い73 遅いけどファレイ数列でとりあえず解。 もちょっとましな方法はありそう。 calc(D0,D1,_):- D2 is D0+D1, D2>12000,!,fail. calc(D0,D1,D2):- D2 is D0+D1. calc(D0,D1,D3):- D2 is D0+D1, calc(D0,D2,D3). calc(D0,D1,D3):- D2 is D0+D1, calc(D2,D1,D3). test:-findall(N,calc(3,2,N),List),length(List,Len),write(Len). *問い74 http://projecteuler.net/problem=74 この問題はリスト処理だけで高速に解くのはひと工夫必要でした。 一回目の変換で123456も654321も152346も数字を並び変えたものは全部同じ値になるということで纏めて計算します。 これを利用して計算量を下げました。 纏めをもっと上手にやればもっと計算量が下がりますが、そのためにはコードが膨らむのでこのへんで手打ちとしました。 kaizyo(0,1):-!. kaizyo(1,1):-!. kaizyo(2,2):-!. kaizyo(3,6):-!. kaizyo(4,24):-!. kaizyo(5,120):-!. kaizyo(6,720):-!. kaizyo(7,5040):-!. kaizyo(8,40320):-!. kaizyo(9,362880):-!. calc(0,0):-!. calc(N,Result):-M is N mod 10,N1 is N//10,calc(N1,Re), kaizyo(M,Re1),Result is Re+Re1. calc2([],1):-!. calc2([P|Rest],Result):- calc2(Rest,Re), kaizyo(P,P1), Result is Re*P1. sum([],Result,Result):-!. sum([X|Rest],Sum,Result):-Sum1 is Sum+X,sum(Rest,Sum1,Result). perm_calc(Deep,Perm,Result):- Deep1 is Deep-1, kaizyo(Deep1,P1), reverse(Perm,RevPerm), [_|Perm1]=RevPerm, select(P2,Perm1,Perm2), P3 is P2-1, calc2([P3|Perm2],Div1), Result is P1//Div1. perm_calc_w(Deep,Perm,Result):- findall(Re,perm_calc(Deep,Perm,Re),List), sum(List,0,Result). perm(_,_,_,_,7,_):-!,fail. perm(N,Down,K,Perm,Deep,[N,Perm1]):- Down>0, perm_calc_w(Deep,[K|Perm],Perm1). perm(N,Down,K,Perm,Deep,Result):- K1 is K+1, Deep1 is Deep+1, kaizyo(Down,P1), N1 is N+P1, perm(N1,Down,K1,Perm,Deep1,Result). perm(N,Down,K,Perm,Deep,Result):- Deep1 is Deep+1, Down1 is Down+1, between(Down1,9,Down2), kaizyo(Down2,P1), N1 is N+P1, perm(N1,Down2,1,[K|Perm],Deep1,Result). test(Result):-perm(0,0,0,[],0,Result). calc3(_,_,58):-!. calc3(N,List,Deep):- Deep1 is Deep+1, calc(N,N1), not(member(N1,List)), calc3(N1,[N1|List],Deep1). search(List,P):- member([N,P],List), calc3(N,[N],0). main74:-findall(Re,test(Re),List), findall(P,search(List,P),Ans), sum(Ans,0,AnsN),write(AnsN). 独り言 http://share-wis.com/ この教育系サイトに誰かスポンサーやパトロンや協力者が増えたらいいなと思う、社会的責任投資の一種としてそんなに悪くないと思う。 勉強サイトだけど評価がゲームに近く、一動画が短くサクサク見れて見た動画の数だけプチ勉強した気分になれるというのがいい。 独り言２ ここら辺までのチュートリアル問題と違い問75あたりから難易度がすこしずつ上がり始める。 流石に基本リスト処理しかないPrologでは計算量的に厳しいものが増えてきた。 配列とstd::mapと優先順位付きキューがPrologにもあればまだ楽なんだけどな。 *問い75 Singular integer right triangles http://projecteuler.net/problem=75 この問題は配列があれば高校生でも簡単に解ける。 リストしかないピュアPrologで解くには少しばかり骨を折る必要がある。 スタック操作だけでこの問題を効率よく解くのはひと手間でした。 今考えているのは以下のような方法。 二つ以上の折り曲げ方がある直角三角形は複数の原始ピタゴラス数の周長を約数に持つはずであり約数の個数は少ないということを利用して計算量とメモリ使用量を下げる。 周長を50000単位でなど適当な長さで区切って計算し結果を集計します。 1～50000までの数をLとし、[Bi,L]をBiをLの約数としこれを全て集めたリストBを作ります。 リストBをBiを基準にBiが小さい順にソートします。 次に原始ピタゴラス数から求まる周長を小さい方から並べて重複する長さがあるものを排除したリストCを作ります。 BとCの先頭要素CiでBiと周長を比べます。 Biの方が小さいならBの先頭を取り除きます。 BiとCiが同じならBの先頭を取り除きLをリストDに追加します。 Ciのほうが小さいならCの要素を取り除きます。 CかBが空になったら計算終了でリストDを返します。 リストDを小さい順に並べLが２個以上あるものを取り除きます。 残ったLの数が１～50000までの周長で一つしか折り曲げ方がないものの答えとなります。 後はこれを50001～100000、、、700001～750000までと集計していけばそこそこの時間と現実的なメモリ使用量で答えが出ます。 という発想で下記のコードを書いてみましたが答えが出るまで５分かかり思ったより速度が出ませんでした。 言語的に言ってスタック操作以外禁止という条件ですからこれくらいの速度が妥当なのかもしれません。 リストで木を表現してstd::mapもどきを作ればもう少し何とかなったかな。 gcd(0, B, B). gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G). check(L,R,R):-0=:=L mod R,R>5. check(L,R,V):-0=:=L mod R, V is L//R,V>R,V>5. calc2(L,B1):- between(1,L,B), B2 is B*B, (L<B2->!,fail;true), check(L,B,B1). calc1(Down,Up,[B,L]):- between(Down,Up,L), calc2(L,B). calcL(M,L):- between(1,M,N), N<M, L is M*(M+N), (750000<L->!,fail;true), 1=:=(M-N) mod 2, gcd(M,N,G), G=:=1. calcM(L):- between(2,866,M), T is M*M+M, (750000<T->!,fail;true), calcL(M,L). single([],L,1,[L]):-!. single([],_,_,[]):-!. single([L|Rest],L,K,Result):- !, K1 is K+1, single(Rest,L,K1,Result). single([L|Rest],L1,1,[L1|Result]):- L\==L1, !, single(Rest,L,1,Result). single([L|Rest],L1,_,Result):- L\==L1, !, single(Rest,L,1,Result). merge_count([],_,[]):-!. merge_count(_,[],[]):-!. merge_count([L|Rest],[[L,L1]|Rest1],[L1|Result]):- !, merge_count([L|Rest],Rest1,Result). merge_count([L|Rest],[[L1,L2]|Rest1],Result):- L<L1,!, merge_count(Rest,[[L1,L2]|Rest1],Result). merge_count([L|Rest],[[L1,_]|Rest1],Result):- L1<L,!, merge_count([L|Rest],Rest1,Result). calc3(List2,Len):- between(0,14,N), Down is N*50000+1, Up is (N+1)*50000, findall(S,calc1(Down,Up,S),List3), sort(List3,List4), merge_count(List2,List4,Ans), msort(Ans,Ans1), single(Ans1,0,0,Ans2), length(Ans2,Len),write([ans,Len]),nl. sum([],Sum,Sum):-!. sum([X|Rest],Sum,Result):-Sum1 is Sum+X, sum(Rest,Sum1,Result). main75:-findall(L,calcM(L),List), msort(List,List1), single(List1,0,0,List2), findall(X,calc3(List2,X),Sums),nl, write(Sums), sum(Sums,0,Ans),write(Ans). *問い76 100の分割数-1の数を求めよという問題。 Wikiの定義通り実装。 calc([],_,_,_,0):-!. calc([X|Rest],K,Deep,Deep,Result):- !, (K<0->K1 is -K+1;K1 is -K), M is (K1*(3*K1-1))//2, Deep1 is Deep+1, calc(Rest,K1,Deep1,M,Re), (1=:= K mod 2->Add is X;Add is -X), Result is Re+Add. calc([_|Rest],K,Deep,M,Result):- !, Deep1 is Deep+1, calc(Rest,K,Deep1,M,Result). calc1(List):- calc(List,1,1,1,Re), length(List,Len), write([Len,Re]),nl, Ans is Re-1, (Len=:=100->!,write([ans,Ans]);calc1([Re|List])). *問い77 http://projecteuler.net/problem=77 ちょっとよんで即座に速度の出る方法を思いつかなかったので他の人の考えを検索。 素数和で素数を小さい順から並べた時末尾の素数と同じより大きな数字を足す動的計画法でアセプト。 という一行を読んで発想法を3秒で理解。 なるほどなるほどとPrologで実装したらすこしめんどくさかった。 std::mapがないから、ソートして集計し直す処理でそれを代用するあたりがめんどくさい。 is_prime(N):- between(2,N,N1), M is N1*N1, (N<M ->!,true; 0=:=N mod N1,!,fail). sums([],[N,LP,C],[[N,LP,C]]):-!. sums([[N,LP,C]|Rest],[N,LP,C1],Result):- !, C2 is C+C1, sums(Rest,[N,LP,C2],Result). sums([X|Rest],[N1,LP,C1],[[N1,LP,C1]|Result]):- !, sums(Rest,X,Result). n_sum(_,[],Result,Result):-!. n_sum(N,[[N,_,C]|Rest],Sum,Result):- !, Sum1 is Sum+C, n_sum(N,Rest,Sum1,Result). n_sum(N,[_|Rest],Sum,Result):- n_sum(N,Rest,Sum,Result). calc3(_,[],_,[]):-!. calc3(_,_,[],[]):-!. calc3(N,[[N1,LP,C]|List],[P|Ps],[[N,P,C]|Result]):- LP=<P, N=:=N1+P, !, calc3(N,List,[P|Ps],Result). calc3(N,[[N1,_,_]|List],[P|Ps],Result):- T is N1+P, N>T, !, calc3(N,List,[P|Ps],Result). calc3(N,List,[_|Ps],Result):- !, calc3(N,List,Ps,Result). calc2(N,Ps,List):- calc3(N,List,Ps,List1), append(List,List1,List2), msort(List2,List3), [Top|List4]=List3, sums(List4,Top,List5), n_sum(N,List5,0,Perm), write([N,Perm]),nl, N<72, N1 is N+1, calc(N1,Ps,List5). calc(N,Ps,List):- is_prime(N), !, calc2(N,[N|Ps],List). calc(N,Ps,List):- !, calc2(N,Ps,List). main77:-calc(2,[],[[0,0,1]]).
