プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ
Problem 61 「巡回図形数」 †
三角数 P3,n=n(n+1)/2 1, 3, 6, 10, 15, ...
四角数 P4,n=n^2 1, 4, 9, 16, 25, ...
五角数 P5,n=n(3n-1)/2 1, 5, 12, 22, 35, ...
六角数 P6,n=n(2n-1) 1, 6, 15, 28, 45, ...
七角数 P7,n=n(5n-3)/2 1, 7, 18, 34, 55, ...
八角数 P8,n=n(3n-2) 1, 8, 21, 40, 65, ...
3つの4桁の数の順番付きの集合 (8128, 2882, 8281) は以下の面白い性質を持つ.
- この集合は巡回的である. 最後の数も含めて, 各数の後半2桁は次の数の前半2桁と一致する
- それぞれ多角数である: 三角数 (P3,127=8128), 四角数 (P4,91=8281), 五角数 (P5,44=2882) がそれぞれ別の数字で集合に含まれている
4桁の数の組で上の2つの性質をもつはこの組だけである.
三角数, 四角数, 五角数, 六角数, 七角数, 八角数が全て表れる6つの巡回する4桁の数からなる唯一の順序集合の和を求めよ.
この問題は探索でといても十分間に合いますが。
スタック操作だけで動的計画法で解くというチャレンジをしてみました。
どの部分の処理も線形時間か n log(n)で終わるようになっています。
appendを使ってるところだけが効率が悪いくらいです。
calcP3(N1):-between(1,200,N),N1 is (N*(N+1))//2,999<N1,N1<10000.
calcP4(N1):-between(1,100,N),N1 is N*N,999<N1,N1<10000.
calcP5(N1):-between(1,100,N),N1 is (N*(3*N-1))//2,999<N1,N1<10000.
calcP6(N1):-between(1,100,N),N1 is (N*(2*N-1)),999<N1,N1<10000.
calcP7(N1):-between(1,100,N),N1 is (N*(5*N-3))//2,999<N1,N1<10000.
calcP8(N1):-between(1,100,N),N1 is (N*(3*N-2)),999<N1,N1<10000.
mysweep([],[]):-!.
mysweep([[NNum,FNum,Sum]],[]):-integer(NNum),integer(FNum),integer(Sum),!.
mysweep([[NNum,FNum,_],[NNum,FNum,Sum]|Rest],Result):-!,
mysweep([[NNum,FNum,Sum]|Rest],Result).
mysweep([_,Set|Rest],[Set|Result]):-
!,
mysweep([Set|Rest],Result).
calc2([NN,FN,Sum],[[NN,NB]|_],[NB,FN,Sum1]):-
Sum1 is Sum+NN*100+NB.
calc2(Set,[_|Rest],Result):-
calc2(Set,Rest,Result).
calc([],_,Ans,Perm):-!,search(Ans,Perm).
calc(_,[],Ans,Perm):-!,search(Ans,Perm).
calc([[NN,FN,Sum]|Rest],[[NN,NB]|Rest1],Ans,Perm):-
!,
findall(Re1,calc2([NN,FN,Sum],[[NN,NB]|Rest1],Re1),List),
append(Ans,List,Ans1),
calc(Rest,[[NN,NB]|Rest1],Ans1,Perm).
calc([[NN,_,_]|Rest],[[NF,NB]|Rest1],Ans,Perm):-
NN<NF,!,
calc(Rest,[[NF,NB]|Rest1],Ans,Perm).
calc(List,[_|Rest],Ans,Perm):-
!,
calc(List,Rest,Ans,Perm).
check([[A,A,Ans]|_],Ans):-!.
check([_|Rest],Ans):-
check(Rest,Ans).
search([],[]):-!,fail.
search(Ans,[]):-check(Ans,Ans1),write([ansis,Ans1]).
search(Ans,Perm):-
sort(Ans,Ans1),
[Top|_]=Ans1,
mysweep(Ans1,Re1),
Ans2=[Top|Re1],
select(L,Perm,Rest1),
calc(Ans2,L,[],Rest1).
to_date3([],[]):-!.
to_date3([Y|Rest],[[N2,N1,Y]|Result]):-
N1 is Y // 100,
N2 is Y mod 100,
to_date3(Rest,Result).
to_date2([],[]):-!.
to_date2([Y|Rest],[[N2,N1]|Result]):-
N2 is Y // 100,
N1 is Y mod 100,
to_date2(Rest,Result).
to_date([],[]):-!.
to_date(List,Re2):-
member(L,List),
to_date2(L,Re1),
sort(Re1,Re2).
main61:-
findall(Y,calcP3(Y),L3),
findall(Y,calcP4(Y),L4),
findall(Y,calcP5(Y),L5),
findall(Y,calcP6(Y),L6),
findall(Y,calcP7(Y),L7),
findall(Y,calcP8(Y),L8),
findall(Re,to_date([L3,L4,L5,L6,L7],Re),Perm),
to_date3(L8,L88),
search(L88,Perm).
Problem 62 「立方数置換」 †
Nを小さい数字から始めN^3の各桁の数字の出現数でカウントしていき大きな数字へ続けます。
条件を満たせばそれを暫定解として探索していき暫定解*10＜N＾３となればそれより小さい答えはありません。
ｃ++のstd::mapにあたる機能がピュアPrologにはないので線形探索をしておりコード実行速度は遅いです。
std::mapに当たる機能があったら非常に高速に解ける問題です。
mycount(0,Temp,Result):-sort(Temp,Result),!.
mycount(N,Temp,Result):-
N1 is N//10,
M is N mod 10,
(select([M,C],Temp,Rest) ->
C1 is C+1,
mycount(N1,[[M,C1]|Rest],Result);
mycount(N1,[[M,1]|Temp],Result)).
search(N,Ans,_):-
0<Ans,
T is Ans*10,
N3 is N*N*N,
T<N3,
!,
write(Ans).
search(N,Ans,Counts):-
N3 is N*N*N,
N1 is N+1,
mycount(N3,[],Re1),
select([Re1,C,Min],Counts,Rest),!,
(4 =< C ->
((Ans>Min;Ans=:=0) -> search(N1,Min,Counts);
search(N1,Ans,Counts));
C1 is C+1,
search(N1,Ans,[[Re1,C1,Min]|Rest])
).
search(N,Ans,Counts):-
N3 is N*N*N,
N1 is N+1,
mycount(N3,[],Re1),
search(N1,Ans,[[Re1,1,N3]|Counts]).
Problem 63 「べき乗の桁の個数」 †
keta_size(0,0):-!.
keta_size(N,Result):-
N1 is N//10,
keta_size(N1,Re),
Result is Re+1.
search2(A,B,N):-N is A^B,keta_size(N,L),B>L,!,fail.
search2(A,B,N):-N is A^B,keta_size(N,L),B=:=L.
search2(A,B,N):-B1 is B+1,search2(A,B1,N).
search(N):-
between(1,9,A),
search2(A,1,N).
main63:-findall(N,search(N),List),sort(List,List1),length(List1,Len),write(Len).
問い64
gcd(0, B, B) :-!.
gcd(A, B, G) :- A \== 0, R is B mod A, gcd(R, A, G).
calc(N,NSqrt,U,DR,Deep,List,Result):-
D is N-DR*DR,
gcd(D,U,G),
D1 is D//G,
U1 is U//G,
L1 is floor((NSqrt+DR)*U1/D1),
L2 is L1*D1,
UR2 is L2-DR,
Deep1 is Deep+1,
(member([L1,UR2,Deep2,D1],List) ->
!,Result is (Deep1-Deep2) mod 2
;calc(N,NSqrt,D1,UR2,Deep1,[[L1,UR2,Deep1,D1]|List],Result)).
search(10001,Ans):-!,write(Ans).
search(N,Ans):-
sqrt(N,N1),
N=:= floor(N1)*floor(N1),
!,
N2 is N+1,
search(N2,Ans).
search(N,Ans):-
sqrt(N,NSqrt),
FloorSQ is floor(NSqrt),
calc(N,NSqrt,1,FloorSQ,0,[],Result),
Ans1 is Ans+Result,
N1 is N+1,
search(N1,Ans1).
Problem 65 「e の近似分数」 †
解法
問題の指定通り求めてみました。
gcd(0, B, B) :-!.
gcd(A, B, G) :- A \== 0, R is B mod A, gcd(R, A, G).
add(U,D,U1,D1,ReU,ReD):-
U2 is U*D1+U1*D,
D2 is D*D1,
gcd(U2,D2,G),
ReU is U2//G,
ReD is D2//G.
calc([X],U,D):-
integer(X),
U is 1,
D is X.
calc([X|Rest],U1,D1):-
calc(Rest,U,D),
add(U,D,X,1,D1,U1).
calc_list(M):-
between(0,98,N),
(1=:=N mod 3->
M is N //3*2+2;
M is 1).
keta_sum(0,0):-!.
keta_sum(N,Result):-
N1 is N//10,
keta_sum(N1,Re),
Result is Re+N mod 10.
main65:-
findall(M,calc_list(M),List),calc(List,U,D),add(U,D,2,1,U1,_),
keta_sum(U1,Ans),write(Ans).
Problem 66 「ディオファントス方程式」 †
解法
この問題は色々式変形したり漸化式が立たないか考えたけど自力では解けなかった問題。
ネットに転がっていたjavascriptコードをＰｒｏｌｏｇコードに翻訳して100％カンニングして解いた。
驚くほど高速に答えが出てすこし驚いた。
この手法は連分数展開というものらしい。
大雑把な考え方までは理解したけれど、細かい計算は理解してない状態。
D=(ｘ＾２-1)/y^2はｙが大きくなると
x^2/y^2と見分けがつかなくなってくる。
√Ｄ＝x/yを満たす分数を求める問題で近似するという発想らしい。
聞けばなるほどの着想だが厳密解といいう考え方にこだわってると中々出てこない発想だな。
一つ勉強になった。
calc(_,_,_,_,E,F,K,N,_,P2,_,_,X,Result):-
B1 is (X-E*E)/F,
C1 is floor((K+E)/B1),
A1 is C1,
A1=:=2*K,0=:=N mod 2,!,Result is P2.
calc(_,_,_,_,E,F,K,N,P1,P2,Q1,Q2,X,Result):-
B1 is (X-E*E)/F,
C1 is floor((K+E)/B1),
D1 is C1*B1-E,
A1 is C1,
E1 is D1,
F1 is B1,
P3 is A1*P2+P1,
Q3 is A1*Q2+Q1,
N1 is N+1,
calc(A1,B1,C1,D1,E1,F1,K,N1,P2,P3,Q2,Q3,X,Result).
calc_first([Result,D]):-
between(2,1000,D),
sqrt(D,T),
K is floor(T),
K2 is K*K,
D\==K2,
A is K,
E is K,
F is 1,
N is 1,
P1 is 1,
P2 is K,
Q1 is 0,
Q2 is 1,
calc(A,_,_,_,E,F,K,N,P1,P2,Q1,Q2,D,Result).
main66_2:-
findall([D,X],calc_first([D,X]),List),
sort(List,List1),write(List1).
Problem 67 「最大経路の和 その2」 †
問題自体は簡単。
ただPrologには破壊的代入がないせいで少しめんどくさいので只今効率的なコードとは何か熟慮中。
とりあえずProlog風味なファイル読み込みとは何かと考えてコーディングをしてみました。
データにゴミや入力ミスが入ってないから練習問題って楽でいいですね。
次の行の計算が少し混乱してるがまあ正しい答えが出てるからいいかな？
one_read(Result,State):-get0(C),((C =:= -1 ; C=:=10;C=:=32)->Result=[],State is C;
one_read(Re,State),Result=[C|Re]).
row_read([Re1|Result],ReState):-
one_read(Re,State),(Re==[]->Re1=[];to_num(Re,0,Re1)),((State=:= -1;State=:=10)->ReState is State,Result=[];
row_read(Result,ReState)).
all_read([Result|Re]):-row_read(Result,State),(State=:= -1 ->Re=[];
all_read(Re)).
to_num([],Sum,Sum):-!.
to_num([X|Rest],Sum,Result):-
Sum1 is Sum*10+(X-48),to_num(Rest,Sum1,Result).
calc_row(_,[],Re,Result):-!,reverse(Re,Result).
calc_row([L,R|Rest],[Z|Sums],[RS|Rest1],Result):-
!,
L1 is L+Z,
R1 is R+Z,
(L1<RS -> R2 is RS;R2 is L1),
calc_row([R|Rest],Sums,[R1,R2|Rest1],Result).
calc_row_first([L,R|Rest],[Z|Sums],Result):-
L1 is L+Z,
R1 is R+Z,
calc_row([R|Rest],Sums,[R1,L1],Result).
calc([[[]]],Sums):-!,write(a),sort(Sums,Sums1),write(Sums1).
calc([Row|Rest],Sums):-
calc_row_first(Row,Sums,Sums1),
calc(Rest,Sums1).
main67:-seen,see('e67.txt'),all_read(Re),
[Sums|Rest]=Re,
calc(Rest,Sums).
Problem 68 「Magic 5-gon ring」 †
解法
こういう探索系の問題はPrologは強いですね。
いかにもPrologらしいという感じで処理をかけた気がします。
search([],[],[],Sum,Sum,Ans):-!,write(Ans).
search([],[L|Perm],Nums,Sum,Sum,Ans):-!,
[Xi|L1]=L,
[X1|_]=Ans,
select(Xi,Nums,NumsRest),
(X1=:=10;X1<Xi),
search(L1,Perm,NumsRest,Xi,Sum,Ans).
search([Xi|Rest],Perm,Nums,NowSum,Sum,Ans):-
integer(Xi),
!,
Xi<10,
NowSum1 is NowSum+Xi,
search(Rest,Perm,Nums,NowSum1,Sum,Ans).
search([Xi|Rest],Perm,Nums,NowSum,Sum,Ans):-
select(Xi,Nums,RestNum),
NowSum1 is NowSum+Xi,
search(Rest,Perm,RestNum,NowSum1,Sum,Ans).
main68:-
Ans=[X1,X2,X3, X4,X3,X5, X6,X5,X7, X8,X7,X9, X10,X9,X2],
Perm=[[X4,X3,X5],[X6,X5,X7],[X8,X7,X9],[X10,X9,X2]],
select(X1,[9,8,7,6,5,4,3,2,1,10],Rest),
select(X2,Rest,Rest1),
select(X3,Rest1,Rest2),
X2<10,
X3<10,
Sum is X1+X2+X3,
search([],Perm,Rest2,Sum,Sum,Ans).
problem 69 「トーティエント関数の最大値」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2069
この問題はわざわざプログラムを書くまでもありません。
φ(n)=n*(1-1/p1)*(1-1/p2),,,*(1-1/pn)なのですから
A=n/φ(n)=1/(1-1/p1)*(1-1/p2),,,*(1-1/pn)です。
Aが最大化されるのは右辺の分母が最小になるとき。
素数を小さい方からかけていき、2,3,5,7,11,13、、、とかけていき100万を超える手前が答えとなります。
Problem 70 「トーティエント関数の置換」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2070
オイラーのトーティエント関数 φ(n) (ファイ関数とも呼ばれる) とは, n 未満の正の整数で n と互いに素なものの個数を表す. 例えば, 1, 2, 4, 5, 7, 8 は9未満で9と互いに素であるので, φ(9) = 6 となる.
1 は全ての正の整数と互いに素であるとみなされる. よって φ(1) = 1 である.
面白いことに, φ(87109)=79180 であり, 87109は79180を置換したものとなっている.
1 < n < 10^7 で φ(n) が n を置換したものになっているもののうち, n/φ(n) が最小となる n を求めよ.
解法
１素因数は一ケタ目が必ず違うのであり得ない。
３素因数のこれ以上下は絶対ないと言える最小値は1.014
試しに１０万までの答えを探ると
75841で1.0087、、、という答えが出てくる。
４素因数以上は３より大きくなることは確実。
よって２素因数だけを調べればよい。
とりあえず試しに計算したらこうなっただけで何で２素因数だけを調べればよいのかさっぱり理屈がわからない。
直感的には√1000万に近い値同士の掛け算のあたりに答えがあるのでそのあたりだけを調べればよいとわかるのですが。
厳密に説明しろと言われたら困る状態です。
is_prime(N):-
between(1,N,N1),
N2 is 2*N1+1,
N3 is N2*N2,
(N<N3->!,true;
0=:=N mod N2,!,fail).
prime_list(N1):-between(5,4999,N),N1 is N*2+1,is_prime(N1).
countN(0,Result,Result):-!.
countN(N,Counts,Result):-
M is N mod 10,
N1 is N//10,
Counts1 is Counts+8^M,
countN(N1,Counts1,Result).
check(N,M):-
countN(N,0,Re1),
countN(M,0,Re2),
Re1=:=Re2.
calc2(_,[],Ps,Ans,AnsN):-
!,
calc1(Ps,Ans,AnsN).
calc2(P,[P1|_],Ps,Ans,AnsN):-
T is P*P1,
10000000=<T,
!,
calc1(Ps,Ans,AnsN).
calc2(P,[P1|Rest],Ps,Ans,AnsN):-
T is (P*P1)/((P-1)*(P1-1)),
Ans=<T,
!,
calc2(P,Rest,Ps,Ans,AnsN).
calc2(P,[P1|Rest],Ps,_,_):-
N is P*P1,
M is (P-1)*(P1-1),
check(N,M),
!,
Ans is (P*P1)/((P-1)*(P1-1)),
AnsN is P*P1,
calc2(P,Rest,Ps,Ans,AnsN).
calc2(P,[_|Rest],Ps,Ans,AnsN):-
calc2(P,Rest,Ps,Ans,AnsN).
calc1([],Ans,AnsN):-!,write([Ans,AnsN]).
calc1(Ps,Ans,AnsN):-
[P|Ps1]=Ps,
T is (P+1)/P,
Ans>T,
!,
calc2(P,Ps,Ps1,Ans,AnsN).
calc1([_|Ps],Ans,AnsN):-
calc1(Ps,Ans,AnsN).
main70:-findall(N,prime_list(N),Ps),calc1(Ps,2,2).
最終更新：2013年12月26日 17:00
