「prolog勉強プロジェクトオイラー141～150」の編集履歴（バックアップ）一覧に戻る
プロジェクトオイラーの問題を堀江伸一さんがProlog言語でとくページ。
整数論の知識はないので早いコードはあまりかけない。
*Problem 142 「完全平方数コレクション」 †
x + y, x - y, x + z, x - z, y + z, y - zが全て平方数となる整数の組 x > y > z > 0で, 最小の x + y + z を求めよ.
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20142
*解法
最新のパソコンで1秒くらいで答えが出る方法。
x+y=s1
x-y=s2
y+z=s3
y-z=s4
x+z=s5
x-z=s6
として
x-s1=-y
x-s2=y
y-s3=-z
y-s4=z
として
s1+s2の半分をx
x-s2=y
としｓ1の上限を適当に定めてまずこのｘ、ｙの対応リストを求めこれを使って全組み合わせを検討します。
(x,ys)を値ｘが決まった時の対応するｙのリストとします。
ysの要素ｙを一つずつ検証します。
(x,ys)のｘをｙに変えて代入すると(y,zs)とｙに対応するｚのリストが得られます。
ｚｓから一つずつｚを取り出すと
この時点で
x+y,x-y,y+z,y-zは平方数となることがわかっているので
残りはｘ+ｚ、ｘ-ｚが平方数となれば、ｘ、ｙ、ｚがすべての条件を見たします。
あとは条件を満たしたｘ、ｙ、ｚの組でｘ+ｙ+ｚが最小となるものを探せば答えとなります。
答えが見つからなければｓ1の値を大きくして再度コードを実行します。
たぶん数式で解くのが正道でしょうが私には無理でした。
#include<stdio.h>
#include<set>
#include<map>
#include<vector>
int main(){
std::set<int> is_square;
std::map<int,std::vector<int> > memo;
std::map<int,std::vector<int> >::iterator it;
std::vector<int>::iterator vIt,vIt2;
int ans=10000*100000;
for(int i=1;i<2000;i++)is_square.insert(i*i);
for(int i=2;i<1500;i++){
for(int j=1;j<i;j++){
int s1=i*i;
int s2=j*j;
int s3=s1+s2;
if(s3%2==1)continue;
int x=s3/2;
int y=x-s2;
memo[x].push_back(y);
}
}
for(it=memo.begin();it!=memo.end();it++){
int x=(*it).first;
for(vIt=(*it).second.begin();vIt!=(*it).second.end();vIt++){
int y=(*vIt);
if(memo.find(y)==memo.end())continue;
for(vIt2=memo[y].begin();vIt2!=memo[y].end();vIt2++){
int z=(*vIt2);
if(is_square.find(x+z)==is_square.end())continue;
if(is_square.find(x-z)==is_square.end())continue;
if(x+y+z<ans){
ans=x+y+z;
}
}
}
}
printf("ans=%d",ans);
}
*Problem 144 「レーザービームの多重反射について調べ上げよ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20144
レーザー物理学では, "white cell"とはレーザー光線を遅延させるための鏡の装置のことである. 光線はcellに入り, 鏡で反射して, 最終的には飛び出す.
楕円のwhite cellを題材にした問題。
解法
光線の向きを線対称変換する行列と楕円とレーザーの交点を求める関数を作るだけの実装ゲーム。
特に工夫するところはないです。
ｙ＝0近辺で光線が反射する場合を考えベクトルで楕円の接線の傾きを現して線対称変換しなくてはいけないが、
楕円の性質と最初の入射角から、ｙ＝0近辺にレーザーが来ることはないので問題で指定された傾きの式をそのまま使う。
今日風邪で頭が痛いが、プロのPGなら風邪ひき3日徹夜でも納期が迫れば働かざるおえない事もあることを考えるとこれくらいなんてことはないか。
max(A,B,A):-A>B,!.
max(_,B,B):-!.
is_esc(X,Y):-
-0.01=<X,
X=<0.01,
Y>0.
next_point(X1,Y1,Dx,Dy,NextX,NextY):-
S1 is 4*Dx^2+Dy^2,
S2 is 8*X1*Dx+2*Y1*Dy,
S3 is 4*X1^2+Y1^2-100,
M is sqrt(S2^2-4*S1*S3),
Re1 is (-S2+M)/(2*S1),
Re2 is (-S2-M)/(2*S1),
max(Re1,Re2,T),
NextX is X1+Dx*T,
NextY is Y1+Dy*T.
next_vector(X1,Y1,Dx,Dy,NextDx,NextDy):-
M is -4*X1/Y1,
NextDx is Dx*((1-M^2)/(M^2+1))+Dy*((2*M)/(M^2+1)),
NextDy is Dx*(2*M/(M^2+1))+Dy*((M^2-1)/(M^2+1)).
search(_,_,_,_,1000):-write(ok),!.
search(X,Y,Dx,Dy,Turn):-
next_point(X,Y,Dx,Dy,NextX,NextY),
next_vector(NextX,NextY,Dx,Dy,NextDx,NextDy),
write([NextX,NextY,NextDx,NextDy]),nl,
Turn1 is Turn+1,
(is_esc(NextX,NextY)->write(Turn);
search(NextX,NextY,NextDx,NextDy,Turn1)).
*Problem 145 「10億未満に存在するreversibleな数はいくつか?」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20145
ある正の整数nについて, [n + reverse(n)]が奇数のみで表されるようなnが存在する. 例えば, 36 + 63 = 99, 409 + 904 = 1313 のように. この性質を持つ数を, reversibleと呼ぶことにする. つまり, 36, 63, 409, 904はrevesibleである. 先頭の0はnでもreverse(n)でも許されない.
1000未満には120個のreversibleな数が存在する.
10億(10^9)未満では, いくつのreversibleな数が存在するか.
小学生でもわかる考え方で自分なりにコードを考えてみました。
メモ化計算的に末尾から数字を確定していきます。
-何桁の数字を計算しようとするか、
-何桁目の処理か、
-前の桁が奇数だったか偶数だったか、
-繰り上がりはあったかなかったか。
この4種類のデータを変数として漸化式を適用すれば10＾100くらいは一瞬でもとまります。
実際一瞬で答えが出てそれぞれ下記の通りとなりました。
10＾9で正しい答えが出ているので下の答えもたぶん正しいとは思いますが確認してくださる方がいらっしゃれば幸いですね
解法
ｎ桁の数字を各桁に分解しA1A2,,,Anとします。
A1A2,,,,An
+AnAn-1,,A1
を考えます。
まずAn+A1はどう考えても奇数です。
よって合計が奇数になり繰り上がる集合B1と、合計が奇数になり繰り上がらない組の集合B2を考えます。
B1だと繰り上がるので下２桁目はA2+An-1の合計は偶数となります。
もしA2+An-1が繰り上がる偶数だと左側の足し算でA2+An-1が繰り上がりA1+Anが偶数になってしまいます。
よってB1だと2桁目は偶数で繰り上がらないA2,An-1の組み合わせをB1の個数に掛ければ二ケタ目の組み合わせ数が出ます。
またB2だと2桁目は合計が奇数で上と同じ理由で繰り上がらない数となります。
こうして2桁目が完全に決まり、3桁目も同様の発想で組み合わせ数が決まります。
最後は足し算の組み合わせ数の検討が真ん中まで到達したらそれより左の足し算は確定しています。
ｎ桁が奇数桁と偶数桁の場合に注意して組み合わせを計算すれば答えとなります。
この発想を漸化式に落とし込めばｎが100桁でも200桁でも答えが一瞬で出ます。
学術系のプログラム言語なら数万桁程度を普通に扱えるのが基本機能だったりするのでここで答えの桁数は問題になりません。
高校生でもわかる簡単な問題でした。
10^9をほかの値に変えた答え。
10^20 413497007808720
10^21 413497007808720
10^22 12223297007808720
10^23 15348297007808720
10^24 369642297007808720
10^25 369642297007808720
10^26 10998462297007808720
10^30 9897363562297007808720
10^50 5843418646273505631953562297007808720
10^100 49510212020123267770652667541980512749637205445373505631953562297007808720
10^2000 9117729789522804392347967308635478385968427812083780466654623229150542050666930341660540095521434910095929501357741947282420376489501300035868050706642098037447697561586566995577403565558615496562693222576482400251763331720240551407113439652548775147282842280434336921765061097857139782702862670626408509379391995047291050086489079555723810628666596103200756251762790211620833179666617305416643649319878976157041007994301058589208452099736862822577396048921353188422953957469022513993504249879068481791168543984520765669269573630018658005666505424642388224725312694354225692764840024877340888673899523184299633750259138967590353120033169787851565199364245732845000345518623453804160044226383802086932485660977126667127358164605072213392301845069449243314214636168889503144219473429617856402460092598991085619514891852670858959297906157141869946790131988114159353189136894478612397208209522493262386842650818879137585515859304816529610946029991016515790201091838850114021145739755372814594706654688687720268122451800152028194319673830419459608872918250293690829935733536037592426231773892612811830557667058254439914311381383456568309031856817082440743556077672586552415175177942091078709142422776587658074770230115403220233570589454771612189897035450210766360306820537626978094119273028816253196047266947688480409094050169304125492364038421670928063022596917973878792066892405500656485384562227904084030129223965171722755856540667541980512749637205445373505631953562297007808720
sum([],0):-!.
sum([[_,_,C]|Rest],Result):-sum(Rest,Re),Result is Re+C.
union_sum([],Set,[Set]):-!.
union_sum([[P,U,C]|Rest],[P,U,C1],Result):-
!,
C2 is C+C1,
union_sum(Rest,[P,U,C2],Result).
union_sum([A|Rest],Set,[Set|Result]):-
union_sum(Rest,A,Result).
% KiAdd0 30
% ki Add1 20
% gu Add0 25
% gu Add1 25
next_calc([0,1,C],[0,1,C1]):-!,C1 is C*25.
next_calc([0,0,C],[1,1,C1]):-!,C1 is C*20.
next_calc([1,1,C],[0,0,C1]):-!,C1 is C*25.
next_calc([1,0,C],[1,0,C1]):-!,C1 is C*30.
next_calc_w(Memo,Result):-
member(M,Memo),
next_calc(M,Result).
last1([1,1,C],[0,0,C1]):-!,C1 is C*5.
last1([0,1,C],[0,1,C1]):-!,C1 is C*5.
last0([0,1,C],[0,1,C1]):-!,C1 is C*25.
last0([1,0,C],[0,0,C1]):-!,C1 is C*30.
seed([1,U,1]):-
between(1,9,A),
between(1,9,B),
1 =:= (A+B) mod 2,
(A+B<10->U is 0;U is 1).
calc(_,2,_,[0,0,20]).
calc(Limit,KetaP,_,_):-
Limit<KetaP*2+1,
!,fail.
calc(_,_,Memo,Result):-
member(M,Memo),
last1(M,Result).
calc(Limit,KetaP,_,_):-
Limit<KetaP*2+2,
!,fail.
calc(_,_,Memo,Result):-
member(M,Memo),
last0(M,Result).
calc(Len,KetaP,Memo,Result):-
!,
findall(M,next_calc_w(Memo,M),Memo1),
msort(Memo1,Memo2),
[Top|Memo3]=Memo2,
union_sum(Memo3,Top,Memo4),
KetaP1 is KetaP+1,
calc(Len,KetaP1,Memo4,Result).
seed_create(Result):-
findall(Seed,seed(Seed),Seeds),
msort(Seeds,Seeds1),
[Top|Seeds2]=Seeds1,
union_sum(Seeds2,Top,Result).
calc_start(Seed,Sum):-
read(Limit),
calc(Limit,1,Seed,Sum).
main145:-seed_create(Seed),
findall(Sum,calc_start(Seed,Sum),Sums),
sum(Sums,Ans),write([ans,Ans]).
*Problem 147 「斜交平行格子内の長方形」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20147
長方形の中に線が描かれており何個の長方形が隠れているか答える問題。
とりあえずｎ*ｍサイズでｈ=<wのときのみ正しい答えをだす関数を動的計画法で作りまあまあの速度計算時間７秒で正答。
最初全探索する遅いコードで正答。
他人のコードをチラ見するもなんだか漸化式を使った方法で見当がつかなかったので自前の動的計画法となりました。
カンニングするのもなんだし動的計画法でもう少し攻めて、それが十分早くなってから他人の漸化式を理解するところに行きたいと思う。
たぶん動的計画法を極めた先に漸化式があると思われます。
現状あるサイズを計算するときはまあまあ効率的になってるけど、別のサイズを計算するときと計算がかぶってて同じ計算をしている。
ここを改善すれば速度は出るかどうか微妙なところです。
sum([],0):-!.
sum([P|Rest],Result):-sum(Rest,Re),Result is Re+P.
sum_up1([],P1,Sum,Result):-!,Result is Sum+P1.
sum_up1([P|Rest],P1,Sum,Result):-
(P<P1->Sum1 is Sum+P;Sum1 is Sum+P1),
sum_up1(Rest,P1,Sum1,Result).
sum_up(_,[],Sum,Sum):-!.
sum_up(Left,[[_,_,P]|Right],Sum,Result):-
sum_up1(Left,P,Sum,Sum1),
sum_up([P|Left],Right,Sum1,Result).
in_area(X,Y,H,W):-
(1-(Y mod 2))=<X,
X=<(2*(W-1)-((Y+1) mod 2)),
Y=<(H-1)*2,
0=<Y.
seed([0,1,1],H,W):-in_area(0,1,H,W).
seed([1,0,1],H,W):-in_area(1,0,H,W).
next_cell([[X,Y,P]|_],[X1,Y1,P1],H,W):-
X1 is X+1,
Y1 is Y+1,
P1 is P+1,
in_area(X1,Y1,H,W).
next_cell([[X,Y,_]],[X1,Y1,1],H,W):-
X1 is X+2,
Y1 is Y,
in_area(X1,Y1,H,W).
next_cell([_|Cells],Result,H,W):-
next_cell(Cells,Result,H,W).
top_add([],[],_,_):-!.
top_add([[X,Y,1]],[[X1,Y1,1]],H,W):-
X1 is X,Y1 is Y+2,
in_area(X1,Y1,H,W),!.
top_add([[X,Y,P]|Rest],[[X1,Y1,1],[X,Y,P]|Rest],H,W):-
X1 is X-1,Y1 is Y+1,
in_area(X1,Y1,H,W),!.
top_add(A,A,_,_):-!.
calc([],Sum,Sum1,H,W):-!,Sum1 is Sum+(H*(H+1)*W*(W+1))//4.
calc(Cells,Sum,Result,H,W):-
sum_up([],Cells,Sum,Sum1),
findall(Cell,next_cell(Cells,Cell,H,W),Cells1),
top_add(Cells1,Cells2,H,W),
calc(Cells2,Sum1,Result,H,W).
test(Result1):-
between(1,43,H),
between(H,47,W),
findall(Cell,seed(Cell,H,W),Cells),
calc(Cells,0,Result,H,W),
write([H,W,Result]),nl,
((H=:=W;43<W)->Result1 is Result;Result1 is Result*2).
main147:-findall(X,test(X),Xs),
sum(Xs,Ans),write(Ans).
*Problem 148 「パスカルの三角形を探し当てよ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20148
パスカルの三角形の最初の7列には7で割り切れる要素は一つもないことが簡単に分かる.
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
しかし最初の100列を調べると, 5050個の要素の内, 7で割り切れないものは2361個しかない.
パスカルの三角形の最初の10億列 (10^9列) の要素で7で割り切れないものの数を答えよ.
解法
大きさの違う三角形の個数から次の段の三角形の個数が決まるオートマトンとみて解いてみたのですが。
数分で解けるものの余り早くありません。
正答率が高いということは何か整数論の知識があれば簡単に解ける何かがあるのでしょうか？
#include<iostream>
#include<math.h>
const __int64 R=7;
__int64 calc_wide(__int64 row,__int64 maxWide){
return maxWide-((row+maxWide-1) % maxWide)-1;
}
int main(){
__int64 memo[10]={0},add,size[10],w,split,ans=0;
for(int i=0;i<=9;i++)size[i]=pow(R,10-i);
for(__int64 row=1;row<=pow(10,9);row++){
if(row<=R)continue;
if(row<=R*R+1){
memo[9]+=((row%R)==1);
}else{
add=0;
for(int i=0;i<9;i++){
w=calc_wide(row,size[i]);
if((size[i]-w)%R==1){
add=1;
break;
}
}
memo[9]+=add;
}
add=0;
for(int i=9;i>=0;i--){
memo[i]+=add;
if(memo[i]>=R){
memo[i]=0;
add=1;
}else{
add=0;
}
}
split=0;
for(int i=0;i<=9;i++){
if(split>0&&memo[i]>0){
ans+=calc_wide(row,size[i])*split*memo[i];
split*=(memo[i]+1);
}
if(split==0&&memo[i]>0){
ans+=calc_wide(row,size[i])*memo[i];
split=memo[i]+1;
}
}
if(row%100000==0)std::cout<<row<<" "<<ans<<"\n";
}
__int64 L1=pow(10,9),L2=pow(10,9)+1;
std::cout<<"\n"<<(L1*L2)/2-ans<<"\n";
}
*Problem 149 「和が最大となる部分列を探し出せ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20149
2000*2000マスの中の縦横斜めの和の最大値を探す問題。
解法
動的計画法とPrologがどんなに相性が悪いか見本みたいな問題。
私の頭の悪さを証明するようなコード。
C++なら素直にループで書けばいいところを、再帰で書くとこうなる。
もしかしてエレガントに書く方法もあるかな？
Dammyデータを入れれば述語の数を減らせる可能性はあるけれど？
lfg(K,Seeds,Result,S):-
K=<55,
!,
S is ((100003 - 200003*K + 300007*(K^3)) mod 1000000) - 500000,
append(Seeds,[S],Result).
lfg(_,[S0|Seeds],Result,S):-
!,
nth1(31,Seeds,S31),
S is ((S31+S0+1000000) mod 1000000)-500000,
append(Seeds,[S],Result).
max([],Max,Max):-!.
max([T|Ts],Max,Result):-Max>T,!,max(Ts,Max,Result).
max([T|Ts],_,Result):-max(Ts,T,Result).
search(4000001,_,_,_,
_,_,_,
Ans):-!,write(Ans).
search(N,Seeds,Row,2001,
UL,Tate,UR,Ans):-
!,
write([d,Row]),
Row1 is Row+1,
write(e),
reverse(UL,UL1),
write([f]),
reverse(Tate,Tate1),
write(g),
reverse(UR,UR1),
write([c,N,Ans]),nl,
search(N,Seeds,Row1,1,
UL1,Tate1,UR1,
Ans).
search(N,Seeds,Row,1,
ULs,[T|Tate],[UR|URs],
Ans):-!,
write(a),
N1 is N+1,
lfg(N,Seeds,Seeds1,S),
T1 is T+S,
UR1 is UR+S,
max([UR1,T1,S],Ans,Ans1),
max([T1],S,T2),
search(N1,Seeds1,Row,2,
ULs,Tate,URs,
[S], [T2], [],S,Ans1).
search(N,Seeds,Row,2000,
[UL] ,[T], [],
NUL ,NTate,NUR,Yoko,
Ans):-!,
N1 is N+1,
lfg(N,Seeds,Seeds1,S),
UL1 is UL+S,
T1 is T+S,
Yoko1 is Yoko+S,
max([UL1,T1,Yoko1,S],Ans,Ans1),
max([T1],S,T2),
search(N1,Seeds1,Row,2001,
NUL,[T2|NTate],[S|NUR],
Ans1).
search(N,Seeds,Row,Col,
[UL|ULs],[T|Tates],[UR|URs],
NULs,NTate,NURs,Yoko,
Ans):-
lfg(N,Seeds,Seeds1,S),
UL1 is UL+S,
Tate1 is T+S,
UR1 is UR+S,
Yoko1 is Yoko+S,
N1 is N+1,
(Yoko1<S -> Yoko2 is S;Yoko2 is Yoko1),
Col1 is Col+1,
max([UL1,Tate1,UR1,Yoko2,S],Ans,Ans1),
max([UL1],S,UL2),
max([Tate1],S,Tate2),
max([UR1],S,UR2),
search(N1,Seeds1,Row,Col1,
ULs,Tates,URs,
[UL2|NULs],[Tate2|NTate],[UR2|NURs],Yoko2,
Ans1).
dammy(N,0):-
between(1,N,_).
main149:-
findall(S,dammy(2000,S),Tate),
findall(S1,dammy(1999,S1),ULs),
findall(S2,dammy(1999,S2),URs),
search(1,[],1,1,
ULs,Tate,URs,
0).
*Problem 150 「三角形配列から和が最小となる部分三角形を探し出せ」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20150
問題要約
１０００段のパスカルの三角形風のものに線形合同法で作った数を入れる。
この三角形のある範囲を三角形で囲むとき、囲まれた範囲の合計値が最小になる場所がある。
その値をこたえよ。
解法
一段ずつの動的計画法でアセプト。
遅いのはＰｒｏｌｏｇのせいか私の発想のせいか？
reverseが遅いということはないと思うけど。
うれしいのは問149よりは再帰的な処理がきれいに決まった点くらい。
rand(T,S,T1):-
T1 is (T*615949+797807) mod 2^20,
S is T1-2^19.
min([],Min,Min):-!.
min([T|Ts],Min,Result):-Min<T,!,min(Ts,Min,Result).
min([T|Ts],_,Result):-min(Ts,T,Result).
add([],_,[]):-!.
add([X|Rest],Y,[Z|Result]):-Z is X+Y,add(Rest,Y,Result).
one_calc([],[],[]):-!.
one_calc([X|Rest],[Y|Rest1],[Z|Result]):-Z is X+Y,one_calc(Rest,Rest1,Result).
next_row_calc(_,_,1001,1,_,Ans):-
!,
write(Ans).
next_row_calc(N,T,Row,1,Memo,Ans):-
!,
N1 is N+1,
rand(T,S,T1),
min([S],Ans,Ans1),
next_row_calc(N1,T1,Row,2,Memo,[[S]],[S],Ans1).
next_row_calc(N,T,Row,Col,_,NextMemo,_,Ans):-
Row<Col,
!,
write([Row,Ans]),nl,
Row1 is Row+1,
reverse(NextMemo,NextMemo1),
next_row_calc(N,T,Row1,1,NextMemo1,Ans).
next_row_calc(N,T,Row,Col,[Sums|Rest],NextMemo,Yoko,Ans):-
!,
rand(T,S,T1),
add(Yoko,S,Yoko1),
one_calc(Sums,Yoko1,Sums1),
Col1 is Col+1,
Sums2 =[S|Sums1],
min(Sums2,Ans,Ans1),
next_row_calc(N,T1,Row,Col1,Rest,[Sums2|NextMemo],[S|Yoko1],Ans1).
