※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

「prolog勉強プロジェクトオイラー91~100」の編集履歴(バックアップ)一覧に戻る

prolog勉強プロジェクトオイラー91~100 - (2013/09/15 (日) 05:46:30) の1つ前との変更点

追加された行は青色になります

削除された行は赤色になります。

 プロジェクトオイラーの問題を堀江伸一さんが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と定義すると。
 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 93
 この問題はまったく自信が持てません。
 a~dが1~9の間にあるんじゃないかなと考え最も長い数列を生成した組み合わせを答えとしました。
 不正回だったら探す数を増加させる予定で、たまたま4桁の数で運良く答えでしたが、なぜこの数が最小になるのか?
 原理を全く理解できてない状態です。
 難しい。
 
 
 
  add([U1,D1],[U2,D2],[U3,D3]):-
  	D3 is D1*D2, 
  	U3 is U1*D2+U2*D1.
  dell([U1,D1],[U2,D2],[U3,D3]):-
  	D3 is D1*D2,
  	U3 is U1*D2-U2*D1.
  mult([U1,D1],[U2,D2],[U3,D3]):-
  	D3 is D1*D2,
  	U3 is U1*U2.
  div([U1,D1],[U2,D2],[U3,D3]):-
  	U2 >0,
  	D3 is D1*U2,
  	U3 is U1*D2.
  
  calc([X|Rest],[X|Result]):-
  	calc(Rest,Result).
  calc([X1,X2|Rest],[X3|Rest]):-
  	!,
  	member(Siki,[add(X1,X2,X3),
  		     dell(X1,X2,X3),
  		     mult(X1,X2,X3),
  		     div(X1,X2,X3)]),
 	Siki. 
  
  selects([],_):-!.
  selects([X|Rest],List):-select(X,List,List1),
  	selects(Rest,List1). 
  
  search([[U,D]],Result):-!,D>0,0=:=U mod D,Result is U//D,Result>0.
  search(List,Result):-
  	calc(List,List1),
  	search(List1,Result).
  search_w(List,Result):-
  	selects([X1,X2,X3,X4],List),
  	search([X1,X2,X3,X4],Result).
  
  perm(_,List,List):-length(List,4),!.
  perm([X|Rest],List,Result):-
  	perm(Rest,[[X,1]|List],Result).
  perm([_|Rest],List,Result):-
  	perm(Rest,List,Result).
  count(N,[N|List],Result):-!,
  	N1 is N+1,
  	count(N1,List,Result).
  count(N,_,N):-!.
  
  test([Len,Set]):-findall(Set,perm([1,2,3,4,5,6,7,8,9],[],Set),Sets),
  	member(Set,Sets),
  	findall(N,search_w(Set,N),Ans),
  	sort(Ans,Ans1),
          count(1,Ans1,Len).
  main93:-findall(Len,test(Len),Ans),
  	sort(Ans,Ans1),
  	write(Ans1).
 
 
 
 
 *Problem 94 「擬正三角形」 †
 http://odz.sakura.ne.jp/projecteuler//index.php?cmd=read&page=Problem%2094
 色々ショートカットテクニックを考えたのですが結局原始ピタゴラス数でNを1ずつ増加させて、Nから斜辺と底辺*2が1差になる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++コード、素直に書くと中学生でも解けるので使用するメモリ量抑えてみた。
 search関数はchain[i]がiが未探索の数字なら次のチェイン数字を取り出し、探索済みは0、今探索中はマイナスへと更新しこの値はループの終了条件判定で使っております。
 探索の戻り値はループになったと判明したらnを返し、でないなら-1を返す。
 search関数の後半では戻り値によって答えの更新を行うかを分岐しています。
 n=戻り値ならループ部分は終了したと判定し以後戻り値が-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;
  		}
  	}
  	int count_Abundant=0;
  	for(int i=2;i<=up;i++){
  		if(chain[i]>i){
  			//ループは過剰数を含むのでループを期待して過剰数から始める
   			count_Abundant++;
  			search(i,-1);
  			chain[i]=0;
   		}
  	}
  	printf("%d %d %d",ans,maxLen,count_Abundant);
  }
+
+Prologで解いてみた結果。
+非常に低速なコードとなった。
+どうやって高速化したらよいかよくわからない状態。
+
+ is_prime(N):-
+ 	between(2,N,N1),
+ 	T is N1*N1,
+ 	(N<T->!;0=:= N mod N1,!,fail).
+ 
+ yakusu_sum(0,_,_,_,0):-!.
+ yakusu_sum(N,[],Sum,_,Result):-
+ 	!,
+ 	(N>1->Result is Sum*(N+1);Result is Sum).
+ yakusu_sum(N,[P|Ps],Sum,A,Result):-
+ 	0=:=N mod P,
+ 	!,
+ 	A1 is A+1,
+ 	N1 is N // P,
+ 	yakusu_sum(N1,[P|Ps],Sum,A1,Result).
+ yakusu_sum(N,[P|Ps],Sum,A,Result):-
+ 	(A>0->Sum1 is Sum*((P^(A+1)-1)//(P-1));Sum1 is Sum),
+ 	yakusu_sum(N,Ps,Sum1,0,Result).
+ yakusu_sum_w(N,Ps,Result):-
+ 	yakusu_sum(N,Ps,1,0,Re),
+ 	Result is Re-N.
+ tree_create(Down,Up,_,[]):-
+  	Sa is Up-Down,
+ 	Sa<2,
+ 	!.
+ tree_create(Down,Up,Ps,[Ls,N,Rs]):-
+ 	!,
+ 	M is (Up+Down)//2,
+ 	yakusu_sum_w(M,Ps,N),
+ 	tree_create(Down,M,Ps,Ls),
+ 	tree_create(M,Up,Ps,Rs).
+ 
+ tree_update(P,Down,Up,[Ls,M,Rs],[Ls,0,Rs],M):-
+ 	P =:= (Down+Up)//2,
+ 	!.
+ tree_update(P,Down,Up,[Ls,M,Rs],[Ls1,M,Rs],Result):-
+ 	P1 is (Down+Up)//2,
+ 	P<P1,
+ 	!,
+ 	tree_update(P,Down,P1,Ls,Ls1,Result).
+ tree_update(P,Down,Up,[Ls,M,Rs],[Ls,M,Rs1],Result):-
+ 	!,
+ 	P1 is (Down+Up)//2,
+ 	tree_update(P,P1,Up,Rs,Rs1,Result).
+ 
+ min_get([[N,_]|_],N,Min,Min):-!.
+ min_get([[N,_]|Log],N1,Min,Result):-
+ 	N<Min,
+ 	!,
+ 	min_get(Log,N1,N,Result).
+ min_get([_|Log],N,Min,Result):-
+ 	min_get(Log,N,Min,Result).
+ 
+ chain(N1,N,_,_,Tree,Up,MaxLen,Ans):-
+ 	tree_update(N,-1,Up,Tree,NextTree,NN),
+ 	(NN=:=0;Up=<NN),
+ 	!,
+ 	search(N1,NextTree,MaxLen,Ans,Up).
+ chain(N1,N,Log,Deep,Tree,Up,MaxLen,Ans):-
+ 	tree_update(N,-1,Up,Tree,NextTree,NN),
+ 	member([NN,D1],Log),
+ 	!,
+ 	Sa is Deep-D1+1,
+ 	min_get(Log,NN,NN,ReMin),
+ 	(MaxLen<Sa->search(N1,NextTree,Sa,ReMin,Up);
+ 	(MaxLen>Sa->search(N1,NextTree,MaxLen,Ans,Up);
+ 	(Ans=<ReMin->search(N1,NextTree,MaxLen,Ans,Up);
+ 	 search(N1,NextTree,MaxLen,ReMin,Up)))).
+ chain(N1,N,Log,Deep,Tree,Up,MaxLen,Ans):-
+ 	!,
+ 	tree_update(N,-1,Up,Tree,NextTree,NN),
+ 	Deep1 is Deep+1,
+ 	chain(N1,NN,[[N,Deep]|Log],Deep1,NextTree,Up,MaxLen,Ans).
+ 
+ search(Up,_,MaxLen,Ans,Up):-!,write([MaxLen,Ans]).
+ search(N,Tree,MaxLen,Ans,Up):-
+ 	!,
+ 	tree_update(N,-1,Up,Tree,NextTree,NN),
+ 	N1 is N+1,
+ 	((NN=:=0;Up=<NN)->!,search(N1,NextTree,MaxLen,Ans,Up);
+ 	!,chain(N1,NN,[[N,0]],1,NextTree,Up,MaxLen,Ans)).
+ 
+ primes(N):-
+ 	between(2,1000,N),
+ 	is_prime(N).
+ main95:-findall(P,primes(P),Ps),
+ 	Up is 1000001,
+ 	tree_create(-1,Up,Ps,Tree),
+  	search(1,Tree,0,Up,Up).