「prolog勉強プロジェクトオイラー61～70」の編集履歴（バックアップ）一覧に戻る
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ
*Problem 61 「巡回図形数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2061
三角数, 四角数, 五角数, 六角数, 七角数, 八角数は多角数であり, それぞれ以下の式で生成される.
三角数 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 「立方数置換」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2062
立方数 41063625 (345^3) は, 桁の順番を入れ替えると2つの立方数になる: 56623104 (384^3) と 66430125 (405^3) である. 41063625は, 立方数になるような桁の置換をちょうど3つもつ最小の立方数である.
立方数になるような桁の置換をちょうど5つもつ最小の立方数を求めよ.
解法
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 「べき乗の桁の個数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2063
5桁の数 16807 = 7＾5は自然数を5乗した数である. 同様に9桁の数 134217728 = 8＾9も自然数を9乗した数である.
自然数を n 乗して得られる n 桁の正整数は何個あるか?
解法
この問題は10より大きな数を考えると10＾２＞2桁なのでそれより大きな数の乗数は考慮する必要はありません。
1～9まで検証すれば十分です。
Logを使った方がすっきりしますが桁数を簡単に数えてときました。
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
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2064
10000以下の平方数でない数を平方根を連分数展開で計算する。
これを連分数展開した時、連分数展開が奇数周期になる物の個数を答えよ。
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 の近似分数」 †
e についての連分数である近似分数の100項目の分子の桁の合計を求めよ.
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2065
詳細はリンク先を参照のこと。
解法
問題の指定通り求めてみました。
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 「ディオファントス方程式」 †
http://projecteuler.net/problem=66
2次のディオファントス方程式に関する問題
詳細はリンク先を参照のこと。
解法
この問題は色々式変形したり漸化式が立たないか考えたけど自力では解けなかった問題。
ネットに転がっていた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」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2067
パスカルの三角形風味に数字が並んでいるので上から下までの経路の数字の和が最大になる経路の合計値をこたえよという問題。
問題自体は簡単。
ただ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」 †
5角形の特殊な魔法陣を解く問題。
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2068
解法
こういう探索系の問題は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).
