「prolog勉強プロジェクトオイラー121～130」の編集履歴（バックアップ）一覧に戻る
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ
どうでもいい雑記。
今日宮沢賢治のやまなし という作品を読んだ。
自然描写力が圧倒的、作画が一流のアニメを見た後のような読後感があった。
コードはページの下のほうに掲載。
コードを見てもらう前に競技プログラムの世界の話を少し書いておきます。
競技プログラムとは数学の問題のようにプログラムの問題が出題されそれを解くことを競う競技で国際的な大会も常時開催されていたりします。
ふつうはN月M日の１６時に出題され１８時までに回答を出題したかたのみ採点のように制限時間付となりますが、過去問題や問題集だけを扱ったサイトもあります。
プロジェクトオイラーは問題集サイトになります。
ここで問題を解くプログラムの話をします。
プログラムの世界では書いたプログラムをコードと呼びます。
競技プログラムの世界ではプログラムは。
複雑でメモリ使用量多く速度が遅いコード＜シンプルでメモリ使用量が多く速度が遅いコード＜
複雑でメモリ使用量が少なく速度が遅いコード＜シンプルでメモリ使用量が少なく速度が遅いコード＜
複雑でメモリ使用量が少なく速度が速いコード＜シンプルでメモリ使用量が少なく速度が速いコード＜エレガントなコード。
右側ほど素晴らしいコードとされています。
素人から見ると
複雑なコード（＝かっこいい）＞＞＞＞越えられない壁＞＞＞シンプルなコード（小学生レベルかよ笑いもの）
のように見えるようです。
私のサイトなどは、基本的に右側を最初に思いつけたら最良、だめなら左側から初めて、だんだん右側へ進化させていくことになります。
ここで書いておきますがシンプルなコードは問題の本質だけをつかんだ立派なコードです。
不必要なものがないから処理がシンプルになります。
物理でいうところのE=MC^2的な本質だけをつかんだコードであり、熟練の職人が最小の動作で素晴らしい仕事をこなすのに似ています。
無駄がなく本質をつかんでるからコードがシンプルになるのです。
*Problem 121 「円盤ゲームの賞金額」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20121
問題はリンク先。
無意味に動的計画法で計算量を百回程度に抑えたりしてみましたが、全探索してもたった数万回の探索量。
この程度では速度差は出ません。
プロジェクトオイラーのことですから何か綺麗な数式的答えがあるのかもしれませんね。
fact(0,1):-!.
fact(N,Result):-N1 is N-1,fact(N1,Re),Result is Re*N.
calc_next(Sum,[],[],Sum):-!.
calc_next(Sum,[[M,Count]|Counts],[[M,Sum]|Counts1],Result):-
Sum1 is Sum+M*Count,
!,
calc_next(Sum1,Counts,Counts1,Result).
calc(8,_,Ans):-fact(16,Up),Ans1 is Up//Ans,write(Ans1).
calc(N,Counts,Ans):-
[[M,Count]|Counts1]=Counts,
Sum is M*Count,
calc_next(Sum,Counts1,Counts2,Add),
Ans1 is Ans+Add,
N1 is N-1,
!,
calc(N1,Counts2,Ans1).
seed([N,1]):-
between(1,15,N).
main121:-findall(CSet,seed(CSet),Counts),
calc(15,Counts,1).
*Problem 122 「効率的なべき乗計算」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20122
この問題厳密解を高速に求める方法は未解決問題らしいです。
値が小さい範囲を求めるだけならなんとかなるということだと思います。
でもどうやって？
考えてもよくわからないので色々カンニングして解くことにします。
只今調査中。
この問題は考えれば考えるほど難しい、未解決なわけだ。
２数の積は求められた数字列のうち一番大きな数を片方に取ればよい。
これはたとえばA=6+3とB=6+2の場合A'=6+6とB'=３+２とみてA+B=A'+B'を作ってもおなじ。
これで探索の量を少し減らしている。
発想ひとつで一位が狙えるAOJと違ってプロジェクトオイラーは整数論の知識が要求される。
整数論の知識がない私にはこの問題は今のところお手上げだったのだが。
海外のコードも結構いじましい努力してるな。
深さを制限して1～200まで探索して制限内である数Nを作れたら確定で探索をまとめている。
もっと賢い方法を期待してたけど、こういう小手先のテクニックでも結構速くなるらしい。
プロジェクトオイラーのほうはどうせ整数論の知識もない私のことだし、一位を目指すとか上位を目指すとかそういうことはしない、まあできないけど。
深さ制限探索してある数Nが作れたときそのNをそれまでにもっと低い回数で作れていたら探索しない。
ということでまあ何とか見れる速度になった。
海外のコードを参考にしたし、なぜそれでうまくいくかの原理もよくわかってない。
参考サイトでも探そう。
search(_,_,Count,Up,_):-Count>Up,!,fail.
search([N|_],_,_,_,_):-N>200,!,fail.
search([Commit|_],Commits,Count,_,Commit):-
member([Commit,Up1],Commits),
(Up1<Count -> (!,fail);true),
Count<Up1.
search([N|Nums],Commits,Count,Up,Result):-
member(A,[N|Nums]),
N1 is N+A,
Count1 is Count+1,
search([N1,N|Nums],Commits,Count1,Up,Result).
commit_sub(Commit,Adds,UpA,UpB,UpRe):-
member(Commit,Adds),!,( UpA<UpB -> UpRe is UpA;UpRe is UpB).
commit_sub(_,_,UpA,_,UpA):-!.
commit(Commits,Adds,UpB,[Commit,UpRe]):-
member([Commit,UpA],Commits),
commit_sub(Commit,Adds,UpA,UpB,UpRe).
calc(12,_,Ans):-
!,
write(Ans).
calc(Up,Commits,Ans):-
findall(Commit,search([1],Commits,0,Up,Commit),CommitAdds),
sort(CommitAdds,CommitAdds1),
length(CommitAdds1,Len),
Ans1 is Ans+Len*Up,
Up1 is Up+1,
findall(NextCommit,commit(Commits,CommitAdds1,Up,NextCommit),NextCommits),
calc(Up1,NextCommits,Ans1).
seed([N,100]):-
between(1,200,N).
main122:-findall(S,seed(S),Seeds),calc(0,Seeds,0).
*Problem 123 「素数の自乗で割った余り」 †
Pn-1とPn+1のｎ乗を展開した末尾に注目すればPn^2で割ると余りは末尾の２項以外のこりません。
それもｎが奇数以外計算する必要がありませんしｎが奇数のときはPn^1の項以外のこりません。
あとは地道に小さいほうから計算するだけです。
not_prime(P,[P1|_]):-
P2 is P1*P1,
P<P2,
!,fail.
not_prime(P,[P1|_]):-
0=:=P mod P1,!.
not_prime(P,[_|Ps]):-
not_prime(P,Ps).
is_prime(P,Ps):-not(not_prime(P,Ps)).
calc(P,_,N1):-N2 is (P*2*N1) mod (P*P),
N2>10^10,
!,
write(N1).
calc(P,Ps,N1):-
P1 is P+2,
next_prime(P1,Ps,N1).
next_prime(P,Ps,N):-
is_prime(P,Ps),
!,
(P=<10^5+3->append(Ps,[P],Ps1);Ps1=Ps),
N1 is N+1,
P1 is P+2,
(N1 mod 2=:=1 ->
calc(P,Ps1,N1);
next_prime(P1,Ps1,N1)).
next_prime(P,Ps,N):-
P1 is P+2,
next_prime(P1,Ps,N).
main123:-calc(5,[2,3,5],3).
*Problem 124 Ordered radicals
http://projecteuler.net/problem=124
1～10万での値をrad(n)という値で変換し小さい方から10000個目になる数を答える問題。
今回はPrologのように特殊な言語ばかりやってたら勘が狂うのでちょっと先にC++で書く。
C++なら配列で楽勝、計算量を無意味に下げる楽しみも味わえました。
さてこの処理をPrologで同じことをしようとすると？
リスト処理が間に入ってくるので高速化がちょっと難しいですね。
リストの利点は何でも入るですがそれが逆に配列的利点が弱くなる。
javascriptはそういう意味ではリストと配列のいいとこどり言語なのだなと実感します。
Prologコードのほうは真面目にためし割になりました。
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<vector>
struct Rad{
int rad,n;
bool operator<(const Rad& r)const{
if(rad!=r.rad)return rad<r.rad;
return n<r.n;
}
};
const int up=100001;
bool ps[up];
std::vector<Rad> rads;
void calc_rad(){
memset(ps,true,sizeof(ps));
Rad rad;
for(int i=0;i<up;i++){
rad.n=i;
rad.rad=1;
rads.push_back(rad);
}
rads[1].n=rads[1].rad=1;
for(int i=2;i<up;i+=1){
if(ps[i]==false)continue;
for(int j=i;j<up;j+=i){
ps[j]=false;
rads[j].rad*=i;
}
}
}
int main(){
calc_rad();
std::sort(rads.begin(),rads.end());
printf("(%d %d)",rads[10000].n,rads[10000].rad);
}
124のProlog版。
not_prime(P,[P1|_]):-
P2 is P1*P1,
P<P2,
!,fail.
not_prime(P,[P1|_]):-
0=:=P mod P1,!.
not_prime(P,[_|Ps]):-
not_prime(P,Ps).
is_prime(P,Ps):-not(not_prime(P,Ps)).
prime_list(351,Ps,Ps):-!.
prime_list(P,Ps,Result):-
is_prime(P,Ps),
!,
P1 is P+2,
append(Ps,[P],Ps1),
prime_list(P1,Ps1,Result).
prime_list(P,Ps,Result):-
P1 is P+2,
prime_list(P1,Ps,Result).
calc_divs(N,P,N):-T is N mod P,0<T,!.
calc_divs(N,P,Result):-N1 is N//P,
calc_divs(N1,P,Result).
calc_rad(1,_,1):-!.
calc_rad(N,[],N):-!.
calc_rad(N,[P|Ps],Result):-
0=:=N mod P,
!,
calc_divs(N,P,N1),
calc_rad(N1,Ps,Re),
Result is Re*P.
calc_rad(N,[_|Ps],Result):-
calc_rad(N,Ps,Result).
rads(Ps,[Rad,N]):-
between(1,100000,N),
calc_rad(N,Ps,Rad).
search(10000,[[R,Num]|_]):-
write([R,Num]),!.
search(No,[_|Sets]):-
No1 is No+1,
search(No1,Sets).
main124:-prime_list(5,[2,3],Ps),
findall(Set,rads(Ps,Set),Sets),
sort(Sets,Sets1),
search(1,Sets1).
*Problem 125 「回文数の和」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20125
この問題、多少奥深い話があるかと思って他の人のコードを検索。
日本語では単純な全検証くらいしか引っかからなかった。
自分も同じコードで解く。
ループと再起の置き換えはとても機械的なので特に工夫するところなし。
sum([],0):-!.
sum([X|Rest],Result):-sum(Rest,Re),Result is Re+X.
to_list(0,[]):-!.
to_list(N,[X|Result]):-
X is N mod 10,
N1 is N // 10,
to_list(N1,Result).
is_palin(N):-to_list(N,List),
reverse(List,List).
sum_d(_,N,[]):-10^8=<N*N+(N-1)*(N-1),!.
sum_d(Sum,N,[Sum|Result]):-
Sum1 is Sum+N*N,
N1 is N+1,
sum_d(Sum1,N1,Result).
roop2(_,[],_):-!,fail.
roop2(X,[X1|_],_):-Sa is X1-X,10^8=< Sa,!,fail.
roop2(X,[X1|_],Sa):-
Sa is X1 - X,
is_palin(Sa).
roop2(X,[_|SumsD],Result):-
!,
roop2(X,SumsD,Result).
roop1([X,_|SumsD],Result):-
roop2(X,SumsD,Result).
roop1([_|SumsD],Result):-
roop1(SumsD,Result).
test:-sum_d(0,1,SumsD),
findall(X,roop1(SumsD,X),Xs),sort(Xs,Xs1),
sum(Xs1,Ans),write(Ans).
*Problem 126 「直方体層」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20126
この問題直方体のｎ層コーティングを考えるとｎ層目を計算の簡便のためにI=n-1,wを直方体の高さ、ｈを直方体の横幅とし、ｎ層目のコーティングに称される直方体個数をallとし
もとの直方体の縦Zを求める関数ｆ(all,i,w,h)とすると
この問題はｎ層コーティングのｎ層目を平行な平面で卵焼きよろしく輪切りにしたときの輪切りにした部分の直方体の個数を数えると簡単に出てくる式があります。
特に輪の両端と胴体部分の２つに分けると一瞬で考えるまでもなく下記の式が出てきます。
Z=(all-2*(w*h+2wi+2hi+2i(i-1)))/(2w+2h+4i)のあまりが0かつｗｈ＜allかつW=<H=<Zであるならひとつ条件を満たすと考えればうまくいきそうに思えます。
試すと実際うまくいきました。
配列のある言語なら適当に上限を定めてZ*(2w+2h+4i)+2*(w*h+2wi+2hi+2i(i-1)),ZはH以上の整数と安易に計算すればよいだけの話だと判明します。
C++なら一瞬で答えが出ます。
Prologだと配列がないので愚直に実装するとためし割が遅いです。
何か工夫が必要ですが、うーん？
#include<stdio.h>
#include <string.h>
const int up =50000;
int main(void){
int memo[up+2];
memset(memo,0,sizeof(memo));
for(int w=1;w*w<up;w++){
for(int h=w;h*h<up;h++){
int sum1=2*w+2*h+2*w*h;
if(sum1>up)break;
for(int i=0;i<up;i++){
int sumL=2*w+2*h+4*i;
int sumR=2*(w*h+2*w*i+2*h*i+2*i*(i-1));
if(sumL+sumR>up)break;
for(int A=h;A<up;A++){
int Ans = sumL*A+sumR;
if(Ans>up)break;
memo[Ans]++;
}
}
}
}
for(int i=6;i<up;i++){
if(memo[i]==1000){
printf("ans=%d",i);
break;
}
}
}
Prologでためし割したコード、配列が使えないためにとても遅い。
とりあえず２分木が手っ取り早い代価手段ではあるがワンパターンすぎてそれも楽しくないよなと思う。
calcI(All,W,H,I,Z):-
between(0,All,I),
UpN is All-2*(W*H+2*W*I+2*H*I+2*I*(I-1)),
(UpN<0 -> !,fail;true),
Z is UpN // (2*W+2*H+4*I),
(Z<H->!,fail;true),
0=:=UpN mod (2*W+2*H+4*I),
W=<H,
H=<Z.
calc(All,[W,H,Z,I]):-
between(1,All,W),
HUp is (All-1)//W,
(HUp<W -> !,fail;true),
between(W,HUp,H),
calcI(All,W,H,I,Z).
search(All):-
findall(Ans,calc(All,Ans),Answers),
length(Answers,Len),
write([All,Len]),nl,
Len=:=1000,!,nl,
write(All).
search(All):-All1 is All +2,
search(All1).
*Problem 127 「abc-hit」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20127
この問題の難しいところは、a,bの組み合わせを検討した時点で、ｃを検討したくなるってこと。
12万＊１２万/4*gcdとradの計算量となる。
のだけどもっと賢いほうがあるのではないかと？
なにしろプロジェクトオイラーのことですからなにかありそうですが、どうでしょう。
とりあえず何も考えずにradのテーブルを用意して、全探索してみたコード。
最新のパソコンで３０秒も時間がかかるなんて話にならない遅さ。
論外である。
多分行列演算とか特性方程式とか整数論の奥深い知識を使うのだと思うけど？
ギブアップしてほかの人のコードを眺めてみようかな？
http://www.mathblog.dk/project-euler-127-abc-hits/
英語なので読めないのだが、コードを斜め読みした感じCを基準にループしましょう程度でこの工夫は私もしてみて確かに速度は少し上がった。
うーん、リンク先の人ですらあまり効率的にならないということはこの問題、効率的な方法はかなり難しい？
#include<stdio.h>
#include <iostream>
int gcd(int x, int y)
{
int t;
while(y != 0)
{
t = x % y;
x = y;
y = t;
}
return x;
}
const int up =120000;
__int64 rads[up+2];
int main(void){
__int64 ans=0,r;
for(int i=1;i<=up;i++)rads[i]=1;
for(int i=2;i<=up;i++){
if(rads[i]>1)continue;
for(int j=i;j<=up;j+=i){
rads[j]*=i;
}
}
for(int a=1;a<=up;a++){
for(int b=a+1;a+b<=up;b++){
int c =a+b;
r=rads[a]*rads[b]*rads[c];
if(c<=r)continue;
if(gcd(a,b)*gcd(a,c)*gcd(b,c)==1)ans+=c;
}
}
std::cout<<ans;
}
Problem 128 「六角形タイルの差分」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20128
へクス格子の隣接するマスの値との差分が素数になるものの２０００番目をこたえよという問題。
解法
図を書けばすぐにわかりますが、N層目の六角形を描いたとき、六角形の頂点に来るマスと層のラストマス以外検討する必要がないことがわかります。
２層目と中心１は特別に手計算で求めてあとはは頂点と、最後のマスの７か所で差分が素数になる隣接マスとで差分を出す式を算出すれば答えとなりますが。
少し計算時間がかかるのでもう少し賢い方法がありそうです。
not_prime(N):-
Limit is floor(sqrt(N)),
between(2,Limit,Div),
0=:=N mod Div,
!.
is_prime(N):-not(not_prime(N)).
%4はdammyData
calc_d(N,0,D1,D2,D3,D4):-
!,
D1 is 6*N-1,
D2 is 6*N+1,
D3 is 12*N+5,
D4 is 4.
calc_d(N,6,D1,D2,D3,D4):-
!,
D1 is 6*N-1,
D2 is 12*N-7,
D3 is 6*N+5,
D4 is 4.
calc_d(N,M,D1,D2,D3,D4):-
!,
D1 is 6*N+(M-1),
D2 is 6*N+M,
D3 is 6*N+M+1,
D4 is 6*(N-1)+M.
prime_count([],0):-!.
prime_count([D|Ds],Result):-
is_prime(D),
!,
prime_count(Ds,Re),
Result is Re+1.
prime_count([_|Ds],Result):-
prime_count(Ds,Result).
ans_calc(N,M,Result):-
Z is 6*(N*(N-1))//2+2+N*M,
(M=:=6->Result is Z-1;Result is Z).
roopNM(N,M,2000):-
!,
M1 is M-1,
ans_calc(N,M1,Result),
write([ans,Result]).
roopNM(N,7,Ans):-
!,
N1 is N+1,
roopNM(N1,0,Ans).
roopNM(N,M,Ans):-
calc_d(N,M,D1,D2,D3,D4),
prime_count([D1,D2,D3,D4],PC),
PC=:=3,
!,
Ans1 is Ans+1,
M1 is M+1,
roopNM(N,M1,Ans1).
roopNM(N,M,Ans):-
M1 is M+1,
roopNM(N,M1,Ans).
main128:-roopNM(2,0,2).
*問129 Repunit divisibility
http://projecteuler.net/problem=129
レピュニット数を題材にした問題。
解法
単純に割り算で解きました。
鳩ノ巣原理より100万以下に答えがないはずなので、100万から探索を開始。
すると答えがすぐに見つかりました。
これは邪道でシンプルなコードだけどずるいコードです。
もっと本質的なコードを書きたいと感じた次第。
シンプルならいいってものでもないですね。
calc_len(0,_,_,Len,Len):-!.
calc_len(R,N,Add,Len,Result):-
R<N,
!,
R1 is R*10+1,
Add1 is Add+1,
calc_len(R1,N,Add1,Len,Result).
calc_len(R,N,Add,Len,Result):-
Len1 is Len+Add,
R1 is R mod N,
calc_len(R1,N,0,Len1,Result).
main129:-between(1000000,2000000,N),N mod 2>0,N mod 5>0,
calc_len(1,N,0,1,Re),write([N,Re]),nl,1000000=<Re.
*Problem 130 「素数桁レピュニットと同じ性質を持つ合成数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20130
レピュニット数を題材にした問題その２．
解法
とりあえず問129のコードをそのまま流用してもそれなりの時間で答えは出た。
数字とレピュニット数の桁数の関係をみているのだがいくら見ても規則性が見えてこない。
難しい。
not_prime(N):-
Limit is floor(sqrt(N)),
between(2,Limit,R),
N mod R=:=0,
!.
calc_len(0,_,_,Len,Len):-!.
calc_len(R,N,Add,Len,Result):-
R<N,
!,
R1 is R*10+1,
Add1 is Add+1,
calc_len(R1,N,Add1,Len,Result).
calc_len(R,N,Add,Len,Result):-
Len1 is Len+Add,
R1 is R mod N,
calc_len(R1,N,0,Len1,Result).
search(25,_,Sum):-!,write(Sum).
search(Hit,N,Sum):-
N mod 2>0,N mod 5>0,
not_prime(N),
calc_len(1,N,0,1,Re),(N-1) mod Re =:=0 ,
!,
N1 is N+1,
Hit1 is Hit+1,
Sum1 is Sum+N,
write([N,Re]),nl,
search(Hit1,N1,Sum1).
search(Hit,N,Sum):-
N1 is N+1,
search(Hit,N1,Sum).
main130:-search(0,4,0).
