「prolog勉強プロジェクトオイラー91～100」の編集履歴（バックアップ）一覧はこちら
追加された行は緑色になります。
削除された行は赤色になります。
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。
*問い91
原点が特殊処理が必要なくらいで後は点OAのベクトルを90度回してOA'(a,b)としgcd(a,b)でその長さ割ったベクトルを求め、そのベクトルの定数倍地点になる点の数を数えればOK。
gcd(0, B, B).
gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G).
min(X,Y,X):-X<Y,!.
min(_,Y,Y):-!.
point_count(_,0,1000):-!.
point_count(L,D,Result):-
Result is L//D.
r90(0,_,1,0):-!.
r90(_,0,0,1):-!.
r90(DY,DX,DY1,DX1):-
gcd(DY,DX,G),
DY1 is DX//G,
DX1 is DY//G.
calc_perm(Y,X,Ans,Ans1):-
r90(Y,X,DY,DX),
X1 is 50-X,
Y1 is 50-Y,
point_count(X ,DX,T1),
point_count(Y1,DY,T2),
point_count(X1,DX,T3),
point_count(Y,DY,T4),
min(T1,T2,Re1),
min(T3,T4,Re2),
Ans1 is Ans+Re1+Re2.
search(51,0,Ans):-!,write(Ans).
search(Y,51,Ans):-
!,
Y1 is Y+1,
search(Y1,0,Ans).
search(Y,X,Ans):-
calc_perm(Y,X,Ans,Ans1),
X1 is X+1,
search(Y,X1,Ans1).
main91:-search(0,1,2500).
*問い92 Square digit chains
桁の2乗による連鎖
http://projecteuler.net/problem=92
この問題はC++なら配列を使えば結構な速度で答えが出る。
Prologでは遅い。
少しだけ工夫が必要だ問題の関数をｆと定義すると。
f(102345)=f(201345)
のように数字を並べ替えたものは同じ値になることとどうやら
f(1～9999999)<=567
となることを利用して計算量を下げれそうに見える。
これで試してみることにしよう。
計算を簡単にするために0012407のような0づめを許して組み合わせ数を求めることにする。
動作が遅いPrologだがほとんど瞬きする間に答えが出た。
こんなものだろう。
C++に翻訳して試しに一千万を1兆に変えてみたところやはりまたたきする間に答えが出た。
to_perm(N,Result):-nth0(N,[1,1,2,6,24,120,720,5040,40320],Result).
search(Num,7,Div,NowDiv,_,[Num,Result]):-
!,
Num>0,
to_perm(NowDiv,NowDiv1),
Result is 5040//(Div*NowDiv1).
search(Num,Deep,Div,NowDiv,[N|Nums],Result):-
NowDiv1 is NowDiv+1,
Num1 is Num+N*N,
Deep1 is Deep+1,
search(Num1,Deep1,Div,NowDiv1,[N|Nums],Result).
search(Num,Deep,Div,NowDiv,[_|Nums],Result):-
to_perm(NowDiv,NowDiv1),
Div1 is Div*NowDiv1,
search(Num,Deep,Div1,0,Nums,Result).
sum([],[N,P],[[N,P]]):-!.
sum([[N,P]|Rest],[N,P1],Result):-
!,
P2 is P+P1,
sum(Rest,[N,P2],Result).
sum([[N,P]|Rest],[N1,P1],[[N1,P1]|Result]):-
sum(Rest,[N,P],Result).
calc(0,Sum,Sum):-!.
calc(N,Sum,Result):-
M is N mod 10,
Sum1 is Sum+M*M,
N1 is N//10,
!,
calc(N1,Sum1,Result).
search(1):-!,fail.
search(89):-!.
search(N):-
calc(N,0,N1),
search(N1).
seed(N):-
between(1,567,N),
search(N).
sum2([],_,Sum,Sum):-!.
sum2(_,[],Sum,Sum):-!.
sum2([[N,P]|Perm],[N|OKs],Sum,Result):-
!,
Sum1 is Sum+P,
sum2(Perm,OKs,Sum1,Result).
sum2([[N,_]|Perm],[N1|OKs],Sum,Result):-
N<N1,
!,
sum2(Perm,[N1|OKs],Sum,Result).
sum2(Perm,[_|OKs],Sum,Result):-
!,
sum2(Perm,OKs,Sum,Result).
main92:-Nums=[0,1,2,3,4,5,6,7,8,9],
findall(P,search(0,0,1,0,Nums,P),Perm),
msort(Perm,Perm1),
[[N,_]|_]=Perm1,
sum(Perm1,[N,0],Perm2),
findall(N1,seed(N1),OKs),
sum2(Perm2,OKs,0,Ans),write(Ans).
*問い92 Bcc5.5 C++ Code
1000万を一兆にしてみて計算、またたきする間に答えが出る。
答えは 850878696414通り
#include<stdio.h>
#include<iostream>
#include<vector>
std::vector<__int64> to_perm;
int oks[1001]={0};
const int up=1000;
const int keta=12;//7 1000万 8 1億 9 10億 10 100億 11 10000億、12 一兆
int calc(int n){
if(n==0){
return 0;
}else{
int m=n % 10;
return calc(n/10)+m*m;
}
}
int calc_chain(int n){
if(oks[n]==0){
oks[n]=calc_chain(calc(n));
}
return oks[n];
}
__int64 search(int num,int deep,__int64 div,int nowDiv,int p){
__int64 ans=0;
if(deep==keta){
return to_perm[deep]/(div*to_perm[nowDiv])*oks[num];
}else{
ans+=search(num+p*p,deep+1,div,nowDiv+1,p);
if(p<9){
ans+=search(num,deep,div*to_perm[nowDiv],0,p+1);
}
}
return ans;
}
int main(){
oks[89]=1;
oks[1]=-1;
for(int i=1;i<=up;i++){
calc_chain(i);
}
for(int i=1;i<=up;i++){
if(oks[i]!=1)oks[i]=0;
}
__int64 p=1;
to_perm.push_back(p);
for(int i=1;i<=12;i++){
p*=i;
to_perm.push_back(p);
}
__int64 ans=search(0,0,1,0,0);
std::cout<<ans;
}
*Problem 94 「擬正三角形」 †
http://odz.sakura.ne.jp/projecteuler//index.php?cmd=read&page=Problem%2094
色々ショートカットテクニックを考えたのですが結局原始ピタゴラス数でNを１ずつ増加させて、Nから斜辺と底辺*２が１差になるMを2次方程式の解法から求めて解きました。
計算の終了条件は少し手抜きです。
calc(N,Len):-述語の適用式のmember化は少し楽しいテクニックを思いつけたと思っています。
多分Prologをわかってる人からみたら初歩だと思いますが自力でこういうのを思いつけたとこの楽しさは格別ですね。
以下コード。
over(N):-1000000000<N.
is_square(N,Result):-sqrt(N,N1),
Result is floor(N1),
N=:=Result*Result.
gcd(0, B, B).
gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G).
add([],Ans,Ans):-!.
add([X|Lens],Ans,Result):-
!,
Ans1 is Ans+X,
add(Lens,Ans1,Result).
calc(N,Len):-
member([S1,S2,S3]
,[[3*N*N-1,2*N+T1,2*M*M+2*N*N+4*M*N],
[3*N*N+1,2*N+T1,2*M*M+2*N*N+4*M*N],
[3*N*N+1,T1,4*M*M],
[3*N*N-1,T1,4*M*M]]),
T is S1,
is_square(T,T1),
M is S2,
M>N,
(M-N) mod 2=:=1,
gcd(M,N,G),
G=:=1,
Len is S3,
not(over(Len)).
search(N,Ans,Ans):-
T is 4*N*N,over(T),!.
search(N,Ans,Result):-
N1 is N+1,
findall(Len,calc(N,Len),Lens),
add(Lens,Ans,Ans1),
search(N1,Ans1,Result).
main94:-search(1,0,Ans),write(Ans).
*問い95
この問題はC++なら中学生でも解ける問題である。
Prologで高速に解くとなるとすこしコードがめんどくさい。
まずはメモ化を使ったC++コード、素直に書くと中学生でも解けるので使用するメモリ量抑えてみた。
未探索の数字は次のチェイン数字、探索済みは0、今探索中はマイナスを割り当て、戻り値が-1を返すループになったと判明したｎを返すかで答えの更新を行うかを分岐している。
メモリの安くなった今時こんなコードを書く理由もないので多分に趣味的なコードになっております。
Prologで書くなら木構造もどきを使うことになり少しコードがめんどくさいので後日作ります。
#include<stdio.h>
#include<algorithm>
const int up=1000000;
int chain[up+1]={0};
int maxLen=0;
int ans=up;
int search(int n,int deep){
if(up<n||chain[n]==0){
return -1;
}
if(chain[n]<0){
int sa=-deep+chain[n];
if(maxLen<sa){
maxLen=sa;
ans=n;
return n;
}else if(maxLen==sa){
ans=std::min(ans,sa);
return n;
}else{
return -1;
}
}
int next=chain[n];
chain[n]=deep;
int re=search(next,deep-1);
chain[n]=0;
if(re==-1){
return re;
}else if(re==n){
return -1;
}else{
ans=std::min(ans,n);
return re;
}
}
int main(){
for(int i=1;i<=up/2;i++){
for(int j=i+i;j<=up;j+=i){
chain[j]+=i;
}
}
for(int i=2;i<=up;i++){
if(chain[i]>0){
search(i,-1);
chain[i]=0;
}
}
printf("%d",ans);
}
プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。
*問い91
原点が特殊処理が必要なくらいで後は点OAのベクトルを90度回してOA'(a,b)としgcd(a,b)でOA'割ったベクトルを求め、そのベクトルの定数倍地点になる点の数を数えればOK。
gcd(0, B, B).
gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G).
min(X,Y,X):-X<Y,!.
min(_,Y,Y):-!.
point_count(_,0,1000):-!.
point_count(L,D,Result):-
Result is L//D.
r90(0,_,1,0):-!.
r90(_,0,0,1):-!.
r90(DY,DX,DY1,DX1):-
gcd(DY,DX,G),
DY1 is DX//G,
DX1 is DY//G.
calc_perm(Y,X,Ans,Ans1):-
r90(Y,X,DY,DX),
X1 is 50-X,
Y1 is 50-Y,
point_count(X ,DX,T1),
point_count(Y1,DY,T2),
point_count(X1,DX,T3),
point_count(Y,DY,T4),
min(T1,T2,Re1),
min(T3,T4,Re2),
Ans1 is Ans+Re1+Re2.
search(51,0,Ans):-!,write(Ans).
search(Y,51,Ans):-
!,
Y1 is Y+1,
search(Y1,0,Ans).
search(Y,X,Ans):-
calc_perm(Y,X,Ans,Ans1),
X1 is X+1,
search(Y,X1,Ans1).
main91:-search(0,1,2500).
*問い92 Square digit chains
桁の2乗による連鎖
http://projecteuler.net/problem=92
この問題はC++なら配列を使えば結構な速度で答えが出る。
Prologでは遅い。
少しだけ工夫が必要だ問題の関数をｆと定義すると。
f(102345)=f(201345)
のように数字を並べ替えたものは同じ値になることとどうやら
f(1～9999999)<=567
となることを利用して計算量を下げれそうに見える。
これで試してみることにしよう。
計算を簡単にするために0012407のような0づめを許して組み合わせ数を求めることにする。
動作が遅いPrologだがほとんど瞬きする間に答えが出た。
こんなものだろう。
C++に翻訳して試しに一千万を1兆に変えてみたところやはりまたたきする間に答えが出た。
to_perm(N,Result):-nth0(N,[1,1,2,6,24,120,720,5040,40320],Result).
search(Num,7,Div,NowDiv,_,[Num,Result]):-
!,
Num>0,
to_perm(NowDiv,NowDiv1),
Result is 5040//(Div*NowDiv1).
search(Num,Deep,Div,NowDiv,[N|Nums],Result):-
NowDiv1 is NowDiv+1,
Num1 is Num+N*N,
Deep1 is Deep+1,
search(Num1,Deep1,Div,NowDiv1,[N|Nums],Result).
search(Num,Deep,Div,NowDiv,[_|Nums],Result):-
to_perm(NowDiv,NowDiv1),
Div1 is Div*NowDiv1,
search(Num,Deep,Div1,0,Nums,Result).
sum([],[N,P],[[N,P]]):-!.
sum([[N,P]|Rest],[N,P1],Result):-
!,
P2 is P+P1,
sum(Rest,[N,P2],Result).
sum([[N,P]|Rest],[N1,P1],[[N1,P1]|Result]):-
sum(Rest,[N,P],Result).
calc(0,Sum,Sum):-!.
calc(N,Sum,Result):-
M is N mod 10,
Sum1 is Sum+M*M,
N1 is N//10,
!,
calc(N1,Sum1,Result).
search(1):-!,fail.
search(89):-!.
search(N):-
calc(N,0,N1),
search(N1).
seed(N):-
between(1,567,N),
search(N).
sum2([],_,Sum,Sum):-!.
sum2(_,[],Sum,Sum):-!.
sum2([[N,P]|Perm],[N|OKs],Sum,Result):-
!,
Sum1 is Sum+P,
sum2(Perm,OKs,Sum1,Result).
sum2([[N,_]|Perm],[N1|OKs],Sum,Result):-
N<N1,
!,
sum2(Perm,[N1|OKs],Sum,Result).
sum2(Perm,[_|OKs],Sum,Result):-
!,
sum2(Perm,OKs,Sum,Result).
main92:-Nums=[0,1,2,3,4,5,6,7,8,9],
findall(P,search(0,0,1,0,Nums,P),Perm),
msort(Perm,Perm1),
[[N,_]|_]=Perm1,
sum(Perm1,[N,0],Perm2),
findall(N1,seed(N1),OKs),
sum2(Perm2,OKs,0,Ans),write(Ans).
*問い92 Bcc5.5 C++ Code
1000万を一兆にしてみて計算、またたきする間に答えが出る。
答えは 850878696414通り
#include<stdio.h>
#include<iostream>
#include<vector>
std::vector<__int64> to_perm;
int oks[1001]={0};
const int up=1000;
const int keta=12;//7 1000万 8 1億 9 10億 10 100億 11 10000億、12 一兆
int calc(int n){
if(n==0){
return 0;
}else{
int m=n % 10;
return calc(n/10)+m*m;
}
}
int calc_chain(int n){
if(oks[n]==0){
oks[n]=calc_chain(calc(n));
}
return oks[n];
}
__int64 search(int num,int deep,__int64 div,int nowDiv,int p){
__int64 ans=0;
if(deep==keta){
return to_perm[deep]/(div*to_perm[nowDiv])*oks[num];
}else{
ans+=search(num+p*p,deep+1,div,nowDiv+1,p);
if(p<9){
ans+=search(num,deep,div*to_perm[nowDiv],0,p+1);
}
}
return ans;
}
int main(){
oks[89]=1;
oks[1]=-1;
for(int i=1;i<=up;i++){
calc_chain(i);
}
for(int i=1;i<=up;i++){
if(oks[i]!=1)oks[i]=0;
}
__int64 p=1;
to_perm.push_back(p);
for(int i=1;i<=12;i++){
p*=i;
to_perm.push_back(p);
}
__int64 ans=search(0,0,1,0,0);
std::cout<<ans;
}
*Problem 94 「擬正三角形」 †
http://odz.sakura.ne.jp/projecteuler//index.php?cmd=read&page=Problem%2094
色々ショートカットテクニックを考えたのですが結局原始ピタゴラス数でNを１ずつ増加させて、Nから斜辺と底辺*２が１差になるMを2次方程式の解法から求めて解きました。
計算の終了条件は少し手抜きです。
calc(N,Len):-述語の適用式のmember化は少し楽しいテクニックを思いつけたと思っています。
多分Prologをわかってる人からみたら初歩だと思いますが自力でこういうのを思いつけたとこの楽しさは格別ですね。
以下コード。
over(N):-1000000000<N.
is_square(N,Result):-sqrt(N,N1),
Result is floor(N1),
N=:=Result*Result.
gcd(0, B, B).
gcd(A, B, G) :- A > 0, R is B mod A, gcd(R, A, G).
add([],Ans,Ans):-!.
add([X|Lens],Ans,Result):-
!,
Ans1 is Ans+X,
add(Lens,Ans1,Result).
calc(N,Len):-
member([S1,S2,S3]
,[[3*N*N-1,2*N+T1,2*M*M+2*N*N+4*M*N],
[3*N*N+1,2*N+T1,2*M*M+2*N*N+4*M*N],
[3*N*N+1,T1,4*M*M],
[3*N*N-1,T1,4*M*M]]),
T is S1,
is_square(T,T1),
M is S2,
M>N,
(M-N) mod 2=:=1,
gcd(M,N,G),
G=:=1,
Len is S3,
not(over(Len)).
search(N,Ans,Ans):-
T is 4*N*N,over(T),!.
search(N,Ans,Result):-
N1 is N+1,
findall(Len,calc(N,Len),Lens),
add(Lens,Ans,Ans1),
search(N1,Ans1,Result).
main94:-search(1,0,Ans),write(Ans).
*問い95
この問題はC++なら中学生でも解ける問題である。
Prologで高速に解くとなるとすこしコードがめんどくさい。
まずはメモ化を使ったC++コード、素直に書くと中学生でも解けるので使用するメモリ量抑えてみた。
未探索の数字は次のチェイン数字、探索済みは0、今探索中はマイナスを割り当て、戻り値が-1を返すループになったと判明したｎを返すかで答えの更新を行うかを分岐している。
メモリの安くなった今時こんなコードを書く理由もないので多分に趣味的なコードになっております。
Prologで書くなら木構造もどきを使うことになり少しコードがめんどくさいので後日作ります。
#include<stdio.h>
#include<algorithm>
const int up=1000000;
int chain[up+1]={0};
int maxLen=0;
int ans=up;
int search(int n,int deep){
if(up<n||chain[n]==0){
return -1;
}
if(chain[n]<0){
int sa=-deep+chain[n];
if(maxLen<sa){
maxLen=sa;
ans=n;
return n;
}else if(maxLen==sa){
ans=std::min(ans,sa);
return n;
}else{
return -1;
}
}
int next=chain[n];
chain[n]=deep;
int re=search(next,deep-1);
chain[n]=0;
if(re==-1){
return re;
}else if(re==n){
return -1;
}else{
ans=std::min(ans,n);
return re;
}
}
int main(){
for(int i=1;i<=up/2;i++){
for(int j=i+i;j<=up;j+=i){
chain[j]+=i;
}
}
for(int i=2;i<=up;i++){
if(chain[i]>0){
search(i,-1);
chain[i]=0;
}
}
printf("%d",ans);
}
