「prolog勉強プロジェクトオイラー121～130」の編集履歴（バックアップ）一覧に戻る
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ
どうでもいい雑記。
今日宮沢賢治のやまなし という作品を読んだ。
自然描写力が圧倒的、作画が一流のアニメを見た後のような読後感があった。
*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を直方体の高さ、ｈを直方体の横幅とし、直方体の縦Zを求める関数ｆ(i,w,h)とすると
all=使う直方体の数
Z=(all-2*(w*h+2wi+2hi+2i(i-1)))/(2w+2h+4i)のあまりが0かつｗｈ＜allかつW=<H=<Zであるならひとつ条件を満たすと考えればうまくいきそうに思える。
配列のある言語なら適当に上限を定めてN*(2w+2h+4i)+2*(w*h+2wi+2hi+2i(i-1)と安易に計算すればよいだけの話です。
Prologだと配列がないので愚直に実装するとためし割が遅いです。
何か工夫が必要ですが、うーん？
とりあえず答えが見たくてためし割ったコード。
とても遅い。
答えが思ったより小さい範囲で出たので、findallとsortを使った富豪プログラムでも解決はできそうですが？
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).
