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

prolog勉強プロジェクトオイラー91~100 - (2013/09/12 (木) 11:53:32) の最新版との変更点

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

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

 プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。
 *問い91
-原点が特殊処理が必要なくらいで後は点OAのベクトルを90度回して最小公倍数でその長さ割ったベクトルの定数倍になる点を数えればOK。
+原点が特殊処理が必要なくらいで後は点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で高速に解くとなるとかなりの工夫が必要に見える。
-正直言うと解き方を全く思いつかない。
-原理的に言って無いのかもしれない。
+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).
+
+
+
+
+*問い96
+数独を解けという問題。
+なんかこの問題皆高速化に苦労してるようです。
+入る可能性のある数字が少ないところから選んでいっただけでそれなりの速度で解けました。
+数独データの末尾に改行を入れて少しコードを短縮しました。
+これがProlog的コードかといわれたらProlog的コードではなくC系統の発想です。
+また数独は驚くほど小さなコードで解けるようで、私の書いたような長ったらしくも素朴なコードはいらないようです。
+
+
+ %  Prolog
+ select_cell([],NowCell,_,NowCell):-!.
+ select_cell([[P,Nums]|Cells],NowCell,NowLen,ReCell):-
+ 	!,
+ 	length(Nums,Len),
+ 	(Len<NowLen ->NowCell1=[P,Nums],NowLen1 is Len;NowCell1=NowCell,NowLen1 is NowLen),
+ 	select_cell(Cells,NowCell1,NowLen1,ReCell).
+ 
+ 
+ 
+ dells_row(_,10,_,Cells,Cells):-!.
+ dells_row(X,Y,Num,Cells,Result):-
+ 	Y1 is Y+1,
+ 	select([[X,Y],Nums],Cells,Cells1),
+ 	([]=Nums->!,fail;true),
+ 	select(Num,Nums,Nums1),
+ 	([]=Nums1->!,fail;true),
+ 	!,
+ 	dells_row(X,Y1,Num,[[[X,Y],Nums1]|Cells1],Result).
+ dells_row(X,Y,Num,Cells,Result):-
+ 	Y1 is Y+1,
+ 	!,
+ 	dells_row(X,Y1,Num,Cells,Result).
+ 
+ 
+ 
+ dells_col(10,_,_,Cells,Cells):-!.
+ dells_col(X,Y,Num,Cells,Result):-
+ 	X1 is X+1,
+ 	select([[X,Y],Nums],Cells,Cells1),
+ 	([]=Nums->!,fail;true),
+ 	select(Num,Nums,Nums1),
+ 	([]=Nums1->!,fail;true),
+ 	!,
+ 	dells_col(X1,Y,Num,[[[X,Y],Nums1]|Cells1],Result).
+ dells_col(X,Y,Num,Cells,Result):-
+ 	X1 is X+1,
+ 	!,
+ 	dells_col(X1,Y,Num,Cells,Result).
+ 
+ 
+ 
+ dells_area(_,_,_,9,Cells,Cells):-!.
+ dells_area(X,Y,Num,DP,Cells,Result):-
+ 	DX is X+DP // 3,
+ 	DY is Y+DP mod 3,
+ 	select([[DX,DY],Nums],Cells,Cells1),
+ 	([]=Nums->!,fail;true),
+ 	select(Num,Nums,Nums1),
+ 	([]=Nums1->!,fail;true),
+ 	DP1 is DP+1,
+ 	!,
+ 	dells_area(X,Y,Num,DP1,[[[DX,DY],Nums1]|Cells1],Result).
+ dells_area(X,Y,Num,DP,Cells,Result):-
+ 	!,
+ 	DP1 is DP+1,
+ 	dells_area(X,Y,Num,DP1,Cells,Result).
+ 
+ write_ans(Ans):-
+ 	member([1,1,N1],Ans),
+ 	member([2,1,N2],Ans),
+ 	member([3,1,N3],Ans),
+ 	N4 is N1*100+N2*10+N3,
+ 	write(N4),nl.
+ 
+ search([],_,AnsN):-AnsN<81,!,fail.
+ search([],Ans,_):-!,sort(Ans,Ans1),write_ans(Ans1).
+ search(Cells,Ans,AnsN):-
+ 	select_cell(Cells,[],20,NowCell),
+ 	select(NowCell,Cells,Cells2),
+ 	[[X,Y],Nums]=NowCell,
+ 	member(Num,Nums),
+ 	Ans1=[[X,Y,Num]|Ans],
+ 	dells_row(X,1,Num,Cells2,Cells3),
+ 	dells_col(1,Y,Num,Cells3,Cells4),
+ 	X1 is ((X-1) // 3)*3+1,
+ 	Y1 is ((Y-1) // 3)*3+1,
+ 	dells_area(X1,Y1,Num,0,Cells4,Cells5),
+ 	AnsN1 is AnsN+1,
+ 	search(Cells5,Ans1,AnsN1).
+ create_cells([[X,Y],[1,2,3,4,5,6,7,8,9]]):-
+ 	between(1,9,X),
+ 	between(1,9,Y).
+ 
+ karayomi:-get0(N),!,(N=:=10->!;karayomi).
+ 
+ one_read(_,10,Cells,Ans,AnsN):-
+ 	%sort(Cells,Cells1),
+ 	%write(Cells1),nl,nl,write(Ans),nl,write(a),
+ 	!,search(Cells,Ans,AnsN).
+ 
+ one_read(_,Y,Cells,Ans,AnsN):-
+ 	peek_code(10),!,get0(_),
+ 	Y1 is Y+1,
+ 	one_read(1,Y1,Cells,Ans,AnsN).
+ one_read(X,Y,Cells,Ans,AnsN):-
+ 	peek_code(48),!,get0(_),
+ 	X1 is X+1,
+  	one_read(X1,Y,Cells,Ans,AnsN).
+ one_read(X,Y,Cells,Ans,AnsN):-
+ 	get0(N),
+ 	!,
+ 	Num is N-48,
+ 	select([[X,Y],_],Cells,Cells1),
+ 	dells_row(X,1,Num,Cells1,Cells2),
+ 	dells_col(1,Y,Num,Cells2,Cells3),
+ 	XX is ((X-1) // 3)*3+1,
+ 	YY is ((Y-1) // 3)*3+1,
+ 	dells_area(XX,YY,Num,0,Cells3,Cells4),
+ 	X1 is X+1,
+ 	AnsN1 is AnsN+1,
+ 	one_read(X1,Y,Cells4,[[X,Y,Num]|Ans],AnsN1).
+ 
+ 
+ 
+ main96:-findall(Cell,create_cells(Cell),Cells),
+ 	seen,
+ 	see('e96.txt'),
+ 	between(1,51,_),
+ 	karayomi,
+ 	one_read(1,1,Cells,[],0),fail.
+
+
+
+*問い97
+28433×2^27830457+1の末尾10ケタを答えよという問題。
+簡単なのですぐ答えが出る。
+高校生でもBigO(100)以下で答えが出る簡単な問題。
+
+ mode(10000000000). 
+ calc(Ans,0,_,Ans,_):-!.
+ calc(Ans,N,N2,Result,M):-
+ 	N22 is (N2*N2) mod M,
+ 	N_Herf is N//2,
+ 	(0=:=N mod 2->Ans1 is Ans;Ans1 is (Ans*N2) mod M),
+ 	calc(Ans1,N_Herf,N22,Result,M).
+ 
+ main97:-mode(M),
+ 	calc(28433,7830457,2,Ans,M),Ans1 is (Ans+1) mod M,write(Ans1).
+ 
+
+
+*問い98 Problem 98 「アナグラム平方数」 †
+http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2098
+うーん?
+高速化を目指しても良かったけれど少しめんどくさい実装になりそうだったので手を抜いた。
+高速化を目指すと結構やりがいはありそうだったんだけど。
+そこまでしなくても答えが出たのでいいかな。
+
+ is_square(N):-
+ 	N=:=floor(sqrt(N))*floor(sqrt(N)).
+ 
+ one_word(C):-
+ 	repeat,
+ 	get0(C),
+ 	(C=:=34->!,fail;true).
+ word_list([Len,Word,Word1]):-
+ 	repeat,
+ 	get0(C),
+ 	(C=:= -1->!,fail;true),
+  	C=:=34,
+ 	findall(C1,one_word(C1),Word),
+ 	msort(Word,Word1),
+ 	length(Word,Len).
+ 
+ to_num(_,[],AnsN,AnsN):-!,is_square(AnsN).
+ to_num(ToN,[C|Word],AnsN,Result):-
+ 	!,
+ 	member([C,N],ToN),
+ 	AnsN1 is AnsN*10+N,
+ 	0<AnsN1,
+ 	to_num(ToN,Word,AnsN1,Result).
+ 
+ search([],Word1,ToN,_,AnsN,Result):-!,
+ 	is_square(AnsN),
+ 	to_num(ToN,Word1,0,AnsN1),
+ 	(AnsN<AnsN1->Result is AnsN1;Result is AnsN).
+ search([C|Word],Word1,ToN,Nums,AnsN,Result):-
+ 	!,
+ 	(member([C,N],ToN) ->
+ 	ToN1 = ToN,Nums1 = Nums;
+         select(N,Nums,Nums1),ToN1=[[C,N]|ToN]
+ 	),
+ 	AnsN1 is AnsN*10+N,
+ 	0<AnsN1,
+ 	search(Word,Word1,ToN1,Nums1,AnsN1,Result).
+ 
+ 
+ roop2(_,[],AnsN,AnsN):-!.
+ roop2([Len,_,_],[[Len1,_,_]|_],AnsN,AnsN):-Len<Len1,!.
+ roop2([Len,Word,SortWord],[[_,Word1,SortWord]|Words],AnsN,Result):-
+ 	!,
+  	search(Word,Word1,[],[0,1,2,3,4,5,6,7,8,9],0,AnsN1),
+ 	(AnsN<AnsN1 -> AnsN2 is AnsN1;AnsN2 is AnsN),
+ 	roop2([Len,Word,SortWord],Words,AnsN2,Result).
+ roop2(Word,[_|Words],AnsN,Result):-!,
+ 	roop2(Word,Words,AnsN,Result).
+ 
+ 
+ roop1([],AnsN,AnsN):-!.
+ roop1([WordSet|Words],AnsN,Result):-
+ 	roop2(WordSet,Words,AnsN,AnsN1),
+ 	!,
+ 	roop1(Words,AnsN1,Result).
+ roop1([_|Words],AnsN,Result):-
+ 	!,
+ 	roop1(Words,AnsN,Result).
+ 
+ main98:-seen,
+ 	see('e98.txt'),
+ 	findall(Word,word_list(Word),Words),
+ 	sort(Words,Words1),
+ 	roop1(Words1,0,Ans),write(Ans).
+
+
+
+*問い99 Problem 99 「最大のべき乗」 
+http://projecteuler.net/problem=99
+数万、数万と2つの数字が一行ずつ書かれたテキストがある。
+各行を数万の数万乗と解釈して最も大きな値となる行を答えよ。
+非常に微妙な差があると困りますがlog計算で答えが出ました。
+読み込み処理を簡素化するためにデータテキストの先頭と末尾に[ ]を加え改行の手前に , を追加しております。
+
+ search([],_,AnsRow,_):-write(AnsRow).
+ search([A,B|Nums],AnsN,AnsRow,Row):-
+	C is log(A)*B,
+  	(C<AnsN->AnsN1 is AnsN,AnsRow1 is AnsRow;
+ 	AnsN1 is C,AnsRow1 is Row),
+ 	Row1 is Row+1,
+ 	search(Nums,AnsN1,AnsRow1,Row1).
+ 
+ main99:-
+ 	seen,
+ 	see('e99.txt'),
+ 	read(Nums),
+ 	X=Nums,
+ 	search(X,0,0,1).
+
+
+
+*問い100 Arranged probability
+http://projecteuler.net/index.php?section=problems&id=100
+整数論も行列も人生で一度も授業を受けたことがないので力技で解きました。
+問題を式変形していくと青の数をb、全体の数をaとします。
+2(2b-1)^2=(2a-1)^2+1が出てきます。
+b1=2b-1,a1=2a-1と置いてExcelでこの式を小さい方から試すと
+(1,1),(7,5),(41,29)という答えのセットが得られます。
+(a b) (1) (7)
+(c d)*(1)=(5)
+
+(a b) (7) (41)
+(c d)*(5)=(29)
+という行列式を立てることが出来ます。
+これを解いて後は(a1+1)/2が整数かつ1兆を超えたところが答えです。
+
+
+ calc(A,B):-
+ 	0=:=(A+1) mod 2,
+ 	0=:=(B+1) mod 2,
+ 	A1 is (A+1)//2,
+  	B1 is (B+1)//2,
+ 	Up is 10^12,
+ 	A1>Up,
+ 	!,
+ 	write([A1,B1]).
+ calc(A,B):-
+ 	A1 is A*3+B*4,
+ 	B1 is A*2+B*3,
+ 	!,
+ 	calc(A1,B1).
+ main100:-calc(1,1).