プロジェクトオイラーの問題を堀江伸一さんがPrologで解くページ。 *Problem 164 「どの連続した3桁の和も与えられた数を超えない数」 † http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20164 教科書的な動的計画法で何も工夫するところがありません。 破壊的代入がないので段階を置いて集計する必要があるくらいです。 配列や破壊的代入がない言語なのでこの一手間の面倒くさいこと。 seed([N1,N2,N3,1]):- between(1,9,N1), between(0,9,N2), between(0,9,N3), N1+N2+N3=<9. sum([],0):-!. sum([[_,_,_,P]|Ps],Result):-sum(Ps,Re),Result is Re+P. union_sum([],Set,[Set]):-!. union_sum([[N1,N2,N3,P]|Memo],[N1,N2,N3,P1],Result):- !, P2 is P1+P, union_sum(Memo,[N1,N2,N3,P2],Result). union_sum([Set|Memo],Set1,[Set1|Result]):- union_sum(Memo,Set,Result). calc_next_B(Memo,[N2,N3,N,P]):- member([_,N2,N3,P],Memo), between(0,9,N), N2+N3+N=<9. calc_next_A(Memo,Memo4):- findall(Set,calc_next_B(Memo,Set),Memo1), msort(Memo1,Memo2), [Top|Memo3]=Memo2, union_sum(Memo3,Top,Memo4). calc(20,Memo):-!,sum(Memo,Ans),write(Ans). calc(N,Memo):- calc_next_A(Memo,Memo4), N1 is N+1, calc(N1,Memo4). main164:- findall(Set,seed(Set),Memo), calc(3,Memo). *Problem 165 「交点」 † http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20165 全部の線分の交差判定をしました、遅いです。 交点は分数座標で保持することで同じ点かを判断します。 グローバルスタックが足らないのでコードの実行時に拡張してアセプト。 std::setがないのでメモリや速度が無駄になる。 線分の生成規則から何か交点にの発生に性質があるのだろうか？ gcd(0, B, G) :- G is abs(B). gcd(A, B, G) :- A =\= 0, R is B mod A, gcd(R, A, G). seed(20000,_,Line,[Line]):-!. seed(N,S,Line,[Line1|Result]):- 0=:= N mod 4,N>0, !, S1 is (S*S) mod 50515093, N1 is N+1, T is S1 mod 500, reverse(Line,Line1), seed(N1,S1,[T],Result). seed(N,S,Line,Result):- S1 is (S*S) mod 50515093, N1 is N+1, T is S1 mod 500, seed(N1,S1,[T|Line],Result). cross(X1,Y1,X2,Y2,Result):- Result is X1*Y2-X2*Y1. check([X1,Y1,X2,Y2],[X3,Y3,X4,Y4]):- SaX is X2-X1, SaY is Y2-Y1, SaX3 is X3-X1, SaY3 is Y3-Y1, SaX4 is X4-X1, SaY4 is Y4-Y1, 0>(SaX*SaY3-SaY*SaX3)*(SaX*SaY4-SaY*SaX4). calc_point([X1,Y1,X2,Y2],[X3,Y3,X4,Y4],[AnsX2,AnsY2,GcdX1,GcdY1]):- SaX2 is X2-X1, SaY2 is Y2-Y1, SaX4 is X4-X3, SaY4 is Y4-Y3, cross(SaX4,SaY4,SaX2,SaY2,Down), SaX31 is X3-X1, SaY31 is Y3-Y1, cross(SaX4,SaY4,SaX31,SaY31,Up), AnsX1 is X1*Down+SaX2*Up, AnsY1 is Y1*Down+SaY2*Up, gcd(AnsX1,Down,GcdX), gcd(AnsY1,Down,GcdY), !, AnsX2 is abs(AnsX1//GcdX), AnsY2 is abs(AnsY1//GcdY), GcdX1 is abs(Down//GcdX), GcdY1 is abs(Down//GcdY). roopB(L1,[L2|_],Ans):- check(L1,L2), check(L2,L1), calc_point(L1,L2,Ans). roopB(L1,[_|Lines],Result):- !, roopB(L1,Lines,Result). roopA([L1|Lines],Result):- roopB(L1,Lines,Result). roopA([_|Lines],Result):- !, roopA(Lines,Result). main165:-X is 512*1024*1024,set_prolog_stack(global,limit(X)), seed(0,290797,[],Lines),nl, findall(Ans,roopA(Lines,Ans),Answers), sort(Answers,Answers1), length(Answers1,Len),write(Len). *Problem 167 「Ulam数列について調べ上げよ」 † http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20167%E3%81%88 まずは問題を理解するところから始めてみよう。 何か問題に取り掛かるとき丁寧に列挙すること、丁寧に整理整頓するところから始めるべきだ。 数学の問題だって同じだ。 まずは問題が何を言っているか考えてみる。 ∑U(2,2n+1)1000億,2=<n=<10ということは ΣU(2,5)1000億 ΣU(2,7)1000億 ΣU(2,9)1000億 ΣU(2,11)1000億 ΣU(2,13)1000億 ΣU(2,15)1000億 ΣU(2,17)1000億 ΣU(2,19)1000億 ΣU(2,21)1000億 こういう問題はおそらく小さな値を見て規則性を発見するのが一番素朴な方法だと思う。 まずはそれを試してみよう。 ＊Problem 168 「数の循環」 † 自然数 142857 を考える. 最後の桁 (7) を一番前に持っていく, すなわち, 右に循環させると 714285 を得る. 714285 = 5 × 142857 が確認できる. これは142857の珍しい性質を示している. つまり, 右に循環させた数の約数になっている. 10 < n < 10^100 の n について, この性質をもつ数の総和の下 5 桁を答えよ. http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20168 まあ数秒で答えが出るところまで進んだ。 下一桁を決めたらあとは全自動で残りの桁すべてが決まるとわかれば結構簡単でした。 試行錯誤して解いたので計算に要らない要素がコードに残ってたり、もう少し計算量を抑える方法はあるけれどめんどくさいのでここで終了。 sum([],0):-!. sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X. calc_next(Memo,KetaP,Mult,Result):- member([FNum,LNum,Num,Up,Count,Ans],Memo), Temp is Num*Mult+Up, NextUp is Temp//10, NextNum is Temp mod 10, (KetaP=<10 -> Ans1 is Ans+Count*NextNum*10^(KetaP-1); Ans1 is Ans), Result=[FNum,LNum,NextNum,NextUp,Count,Ans1]. last_calc(Memo,Mult,_,Ans):- member([FNum,_,Num,Up,_,Ans],Memo), Num>0, FNum=:=Num*Mult+Up. search(Limit,Limit,Mult,Memo,Result):- !, findall(E,last_calc(Memo,Mult,Limit,E),Es), sum(Es,Result). search(Limit,KetaP,Mult,Memo,Result):- !, findall(M,calc_next(Memo,KetaP,Mult,M),Memo1), KetaP1 is KetaP+1, search(Limit,KetaP1,Mult,Memo1,Result). seed(Seed,Mult):- between(1,9,FNum), LNum is (FNum*Mult) mod 10, Seed =[FNum,LNum,FNum,0,1,FNum]. calc(Result):- between(3,101,Limit), between(1,9,Mult), findall(Seed,seed(Seed,Mult),Seeds), search(Limit,2,Mult,Seeds,Result). main168:-findall(Sum,calc(Sum),Sums), sum(Sums,Ans), Ans1 is Ans mod 100000, write(Ans1). *Problem 169 「ある数を2のべき乗の和で表せる方法の数を探し当てよ」 † http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20169 何も考えない動的計画法で解けます。 % 答えは１２桁 to_bit2(0,[]):-!. to_bit2(N,[B|Result]):- N1 is N//2, B is N mod 2, to_bit2(N1,Result). check(0,0,0). check(0,0,1). check(0,1,1). check(1,0,0). check(1,1,0). check(1,1,1). union_sum([],Set,[Set]):-!. union_sum([[Add,P]|Memo],[Add,P1],Result):- !, P2 is P1+P, union_sum(Memo,[Add,P2],Result). union_sum([Set|Memo],Set1,[Set1|Result]):- union_sum(Memo,Set,Result). calc_nextB(Memo,Bit,[NextAdd,P]):- member([Add,P],Memo), check(Bit,Add,NextAdd). calc_nextA(Memo,Bit,Memo4):- findall(Set,calc_nextB(Memo,Bit,Set),Memo1), msort(Memo1,Memo2), [Top|Memo3]=Memo2, union_sum(Memo3,Top,Memo4). calc(Memo,[1]):-!,member([0,P0],Memo),member([1,P1],Memo), Ans is P0+P1, write(Ans). calc(Memo,[Bit|Bits]):- !, calc_nextA(Memo,Bit,Memo4), write([Bit,Memo4]),nl, calc(Memo4,Bits). main169:- N is 10^25, to_bit2(N,Bits), calc([[0,1]],Bits).
