Project Euler Problem 31
問題31: イギリス硬貨 1ペンス〜2ポンドを使い 2ポンドをつくる方法は何通りあるか。[=>問題文]
問題31の解答
;; Emacs lisp (require 'cl) (defun problem031a (n) (labels ((repeat (coins acc) (if (null coins) 0 (do ((add (car coins)) (next-coins (cdr coins)) (count 0)) ((>= acc n) (if (= acc n) (1+ count) count)) (incf count (repeat next-coins acc)) (incf acc add))))) (repeat '(200 100 50 20 10 5 2 1) 0))) (problem031a 200)
組み合わせパターンの生成は簡単でしたが、count の置き方、受け渡し方に迷いました。do の終了処理周りがいまいちスッキリしませんが･･･とりあえず答えは出たのでこれで・・・。
ところでいろいろコードをいじってて思ったんですが、なんか labels って遅くないですか？ というわけで実験。
ロジックは全部同じです
(require 'cl) ;; flet 版 (defun problem031b (n) (flet ((repeat (coins acc) (if (null coins) 0 (do ((add (car coins)) (next-coins (cdr coins)) (count 0)) ((>= acc n) (if (= acc n) (1+ count) count)) (incf count (repeat next-coins acc)) (incf acc add))))) (repeat '(200 100 50 20 10 5 2 1) 0))) ;; let + lambda 版 (defun problem031c (n) (let ((repeat (lambda (coins acc) (if (null coins) 0 (do ((add (car coins)) (next-coins (cdr coins)) (count 0)) ((>= acc n) (if (= acc n) (1+ count) count)) (incf count (funcall repeat next-coins acc)) (incf acc add)))))) (funcall repeat '(200 100 50 20 10 5 2 1) 0)))
例によって Core2 Quad 6600 にて計測。
Function Avg (sec) =========== ========== problem031a 3.822200 labels problem031b 0.959400 flet problem031c 1.037200 let + labmda
やっぱり labels は遅かった。labels と flet はどちらも cl パッケージのマクロです。レキシカルスコープの実装がオーバーヘッドになっているんでしょうか？
コンパイル時の警告
let とflet はコンパイル時(compile-defun)に警告が出ます。
;; flet Warning: repeat called with 2 arguments, but accepts only 1 Warning: repeat called with 2 arguments, but accepts only 1 ;; let + lambda Warning: reference to free variable `repeat'
後者は、let + lambda で再帰を書くといつもでる警告ですが、前者はどういうことなんでしょう？ 引数が 1 つしか accept されない？ よくわかりません。
追記： repeat は Emacs Lisp で定義済みの関数でした。直前のコマンドを繰り返す機能です。上記 3 つの関数では flet だけが、グローバルな関数を参照するため、既存の repeat と引数が合わないという警告を出しているようです。
Project Euler Problem 30
問題30: 2以上の自然数について、各桁を5乗した数の和が元の数と一致する数をすべて求め、その和を答えよ。[=>問題文]
上限が与えられない問題です。左辺の値いくつまで調べればよいのか。注目するのは「桁数」です。
k 桁の数を考えます。各桁の最大値は 9 なので右辺 a のとりうる最大値 a-max は次のようになります。
[右辺 a 最大値] a-max(1): 9^5 a-max(2): 9^5 + 9^5 a-max(3): 9^5 + 9^5 + 9^5 : a-max(k-1): 9^5 + 9^5 + 9^5 ... 9^5 a-max(k) : 9^5 + 9^5 + 9^5 ... 9^5 + 9^5
つまり、a の最大値は桁が増える毎に 9^5 ずつ増えていきます。次は左辺 b の最小値 b-min について考えます。
[左辺 b 最小値] b-min(1): 1 (0) b-min(2): 10 b-min(3): 100 : b-min(k-1): 1000...00 (k-1桁) b-min(k) : 1000...000 (k桁)
k=2 の時は2桁の最小値なので10、k=3の bの最小値は 100 です。b の最小値は桁が増えるごとに 10 倍になります。1桁の最小値は 0 ですが便宜的にここでは 1 とします。
もう一度 b-min と a-max を並べて書いてみます。
k b-min a-max 1 1 59049 (9^5) 2 10 118098 (9^5 + 9^5) 3 100 177147 (9^5 + 9^5 + 9^5) 4 1000 236196 (9^5 + 9^5 + 9^5 + 9^5) :
このように、はじめのうちは b-min(左辺最小）< a-max (右辺最大)ですが、桁が増える毎の値の増加量は b-min が指数的増加量（というか指数関数そのものですが）であるのに対して、a-max は単調増加です。だから、k を増やしていけばいずれどこかで b-min > a-max となる瞬間があるはずです。そして大小が逆転した後はその差は開いていく一方です。この問題は b = a となる数を探すわけですから、b-min > a-max となった後は調べる必要がありません。
b-min > a-max となった瞬間の a-max が、調べる値の上限となります。
問題30の解答その1
;;Emacs Lisp (require 'cl) (defun problem030a (n) (let ;; 0^n 〜 9^n を格納した配列。 ;; n乗計算を省略するためにメモ化しておく。 ((expt-n-vect (do ((v (make-vector 10 0)) (i 1 (1+ i))) ((= i 10) v) (aset v i (expt i n)))) ;; 計算対象の上限。 ;; k 桁での左辺最小値を bk。 右辺最大値を ak とし、 ;; bk が ak を超えたとき、その時点の ak の値を ;; 計算すべき上限値とする。 (limit (do ((k 1 (1+ k)) ; 実は不要だけど、一応 (a-max (expt 9 n)) ; a(k)の最大値 (b-min 1) ; b(k)の最小値 (a-inc (expt 9 n))) ; a増加量 ((<= a-max b-min) a-max) (incf a-max a-inc) (incf b-min (* b-min 10))))) (let* ;; 桁分割リスト作成関数。 123 => (1 2 3) ((split-digit (lambda (x &optional lis) (if (zerop x) lis (funcall split-digit (/ x 10) (cons (mod x 10) lis))))) ;; 左辺 b から 右辺 a を計算する関数 (calc-a (lambda (b) (reduce #'(lambda (acc y) (+ acc (aref expt-n-vect y))) (funcall split-digit b) :initial-value 0)))) ;; メインルーチン ;; 左辺 b (limit 〜 2) に対する 右辺 a を計算し ;; b = a となる b をリストに保存し、合計する。 (do ((lis) (b limit (1- b))) ((= b 1) (apply #'+ lis)) (when (= b (funcall calc-a b)) (push b lis)))))) (problem030a 5)
やってることは、最初の説明のとおりです。最後のメインルーチンでは、左辺 b の値を上限 limit から 2 までデクリメントし、すべての b に対する a を計算して b = a となるものを得ています。
このやり方だとn=5の場合、約47万回の右辺計算をしています。しかもその中には b = 1234,1342,2341... のように右辺 a の値が同じになるものが多数含まれています。この無駄をなんとかしてみましょう。
問題30の解答その2
(require 'cl) (defun problem030b (n) ;; --- 前半は、problem030a と同じ ---------- (let ((expt-n-vect ; 0^n 〜 9^n を格納した配列 (do ((v (make-vector 10 0)) (i 1 (1+ i))) ((= i 10) v) (aset v i (expt i n)))) ;; 計算対象の上限。 (limit (do ((k 1 (1+ k)) ; 実は不要だけど、一応 (a-max (expt 9 n)) ; a(k)の最大値 (b-min 1) ; b(k)の最小値 (a-inc (expt 9 n))) ; a増加量 ((<= a-max b-min) a-max) (incf a-max a-inc) (incf b-min (* b-min 10))))) (let ;; 桁分割リスト作成関数。 123 => (1 2 3) ((split-digit (lambda (x &optional lis) (if (zerop x) lis (funcall split-digit (/ x 10) (cons (mod x 10) lis))))) ;; --- ここまで、problem030a と同じ ---------- ;; リスト要素のn乗和計算関数。 (1 2 3) => 1^n+2^n+3^n (sum-of-expt-n (lambda (lis) (reduce #'(lambda (acc y) (+ acc (aref expt-n-vect y))) lis :initial-value 0)))) (let* ;; limitの桁数 ((limit-digits (1+ (truncate (log limit 10)))) ;; 2〜limit までの数を桁分割し sortしたリストを key として ;; 保存した hash tableを作成。hash table の key なので重複は ;; 除去される。 (ht (do ((ht (make-hash-table :test #'equal)) (x limit (1- x)) (border (expt 10 (1- limit-digits))) (zero-list)) ((= x 1) ht) ;; すべてのリストが同じ長さになるよう zero-list を ;; append して調整する。 (let ((lis (append zero-list (funcall split-digit x)))) (puthash (sort lis #'<) t ht)) ;; x の桁が減ったら、zero-list に 0 を追加 (when (= x border) (push 0 zero-list) (setq border (/ border 10))))) ;; 題意に該当する数を保存する為のリスト (correspond-numbers)) ;; 作った hash table (ht) のすべてのキーに対して題意のテスト ;; を行い合致するものを correspond-numbers に保存する (maphash #'(lambda (a-list v) ;; a-list は右辺の基数となる数値のリスト。 ;; a-list の要素の n乗和を 左辺 bとし、 ;; bを桁にばらして sortしたリストを b-list とする. (let* ((b (funcall sum-of-expt-n a-list)) (b-list (sort (funcall split-digit b) #'<)) (b-list-len (length b-list))) ;; a-list の桁数はすべて limit-digits なので ;; b-list の桁数も 0 埋めしてあわせる (if (< b-list-len limit-digits) (setq b-list (append (make-list (- limit-digits b-list-len) 0) b-list))) ;; a-list と b-list が一致したら、それが求める数。 (if (equal a-list b-list) (push b correspond-numbers)))) ht) ;; correspond-numbers には 1 も含まれるので、 ;; sort して car を除いてから合算。 (apply #'+ (cdr (sort correspond-numbers #'<)))))))) (problem030b 5)
1234,1342,2341...のような重複要素を除去してからメインの計算ループに入ろうという発想です。上限 limit までの数をすべて (1 2 3 4) (1 3 4 2) (2 3 4 1) のようなリストに変換し、sort した上で hash table に key として保存しています。ソートすることにより、上記 3 つのような数はすべて同じパターンとなり hash table によって重複要素は除去されます。
n=5の場合、すべてのリストを ht に登録し重複除去された時点でパターン数は 4542個まで絞り込まれました。右辺 a の計算回数は、problem030a の 100分の1にまで減っています。その代わり前準備がかなり複雑になっています。
プロファイル
Function Avg (sec) =========== ========= problem030a 4.156500 problem030b 2.421500
hash table 登録や リストのソートなど、オーバーヘッドになりそうな処理をしているにもかかわらず、後者の方式のほうが速いようです。ループ回数 100倍の差が効いているのでしょうか？
Project Euler Problem 29
問題29: 2<=a<=100, 2<=b<=100 において、a^b で作成される数は重複を除いて何個あるか[=>問題文]
問題29の解答 (list + delete-dups)
;; Emacs Lisp (require 'cl) (require 'bigint) ;integer-to-bigint, bigint* (defun problem029a (limit) (length (delete-dups (do ((a 2 (1+ a)) (lis)) ((> a limit) lis) (let ((n (integer-to-bigint a)) (acc (integer-to-bigint a))) (do ((b 2 (1+ b))) ((> b limit)) (setq acc (bigint* acc n)) (push acc lis))))))) (problem029a 100)
特別なことはしていません。作成される数をリストにすべて保存し、delete-dups で重複を除去しています。100の100乗は明らかに多倍長数なので、bigintライブラリを使っています。
bigintライブラリに新たに integer-to-bigint を追加しました。整数からbigintへ文字列を経由することなくダイレクトに変換できるようになりました。また掛け算関数 bigint* にバグがあったので修正しました。
ところで、delete-dups は遅いらしいです。対象データが多いときは hash-table を使ったほうが良いらしい。というわけで。
問題29の解答その２ (hash-table + hash-table-count)
(defun problem029b (limit) (hash-table-count (do ((a 2 (1+ a)) (ht (make-hash-table :test #'equal))) ((> a limit) ht) (let ((n (integer-to-bigint a)) (acc (integer-to-bigint a))) (do ((b 2 (1+ b))) ((> b limit)) (setq acc (bigint* acc n)) (puthash acc t ht)))))) (problem029b 100)
hash-table はデフォルトではキーの比較に eq を使いますが、今回はbigint、つまりリストをキーにするので equal で比較するよう、make-hash-table に指定しています。ちなみに、delete-dups は equal 比較で固定されています。
どれくらい違うか
Function Avg (sec) =========== ========= problem029a 1.615900 delete-dups 版 problem029b 0.196600 hash-table 版
これはまた、想像以上に差がつきました。
問題29の解答その３
実際に指数計算しなくても、素因数分解結果を保持していれば生成される数のパターンは数えられることに気づきました。これならば多倍長演算不要です。
(require 'cl) (require 'math-lib) ;factorize (defun problem029c (limit) (let* ((primes (primes-sieve (truncate (sqrt limit)))) (parse (lambda (a b) (mapcar #'(lambda (pair) (cons (car pair) (* b (cdr pair)))) (factorize a primes))))) (hash-table-count (do ((a 2 (1+ a)) (ht (make-hash-table :test #'equal))) ((> a limit) ht) (do ((b 2 (1+ b))) ((> b limit)) (puthash (funcall parse a b) t ht)))))) (problem029c 100)
parse にて (a b) の組から a^b の計算結果を素因数分解したリストを作成します。factorize は素因数分解の結果を 素因数の小さいほうから並べたリストで返しますので sort 不要です。（順序が保証されなければここで sort が必要になります。
Project Euler Problem 28
問題28:渦巻状に自然数を並べた1001x1001数方陣の対角線上にある数の合計を求める。方陣の詳細はリンク先参照[=>問題文]
問題28の解答
;; Emacs Lisp (require 'cl) (defun problem028 (N) (do ((n 3.0 (+ n 2)) (result 1)) ((> n N) result) (let ((d (1- n)) (x (* n n))) (incf result x) (incf result (decf x d)) (incf result (decf x d)) (incf result (decf x d))))) (problem028 1001)
これは簡単ですね。外周の一辺のマス目の数を n とすると右上隅は n の二乗。その他の隅は n-1 ずつ減じた数になる。これを合わせれば外周上の四隅の合計が求まり、外周を増やしながらそれらを加算していけば、最終的には対角線上の総計が求まります。
Project Euler Problem 27
問題27: 二次式 n^2 + an + n = p (|a| < 1000, |b| < 1000, p = 素数) において、n に 0から連続する整数を代入したときもっとも大きな n まで式が成り立つ a b の組を求め、その積を答えよ。[=>問題文]
効率的な素数判定方法を考えていたら2日経ってしまいました。最終的に、処理時間は 1.89秒になりました。これは速いなのか遅いほうなのか、これから他の人の解答を見るのが楽しみです。
問題27の解答
;; Emacs Lisp (require 'cl) (require 'math-lib) (defun problem027 (limit) (let* ((b-list (primes-sieve limit)) (prime-p (make-prime-number-p b-list)) (pairs (let ((lis) (b)) (dolist (b (reverse b-list) lis) (do ((a (1- limit) (1- a)) (end (* -1 limit))) ((= a end)) (push (list a b) lis)))))) (do ((n 1 (1+ n))) ((>= 1 (length pairs)) (apply #'* (car pairs))) (let ((lis) (n2 (* n n)) (continue t)) (while (and (consp pairs) continue) (let ((pair (pop pairs))) (destructuring-bind (a b) pair (if (>= n b) (setq continue nil) (let ((x (+ n2 (* a n) b))) (if (and (> x 1) (funcall prime-p x)) (push pair lis))))))) (setq pairs lis))))) (problem027 1000)
n^2+an+b=p に n=0 を代入すると b=p となり、bは素数で無ければならないことがわかります。また、n=b を代入すると左辺は b^2+ab+b となりこれは明らかに b で割り切れます。したがって 連続する nにおいて、式が成り立つ n の最大値は高々 b - 1 であることがわかります。
problem027 では、はじめに a, b の組み合わせをすべて作成 (pairs) しています。但しbは1000未満の素数だけです。メインのループでは n の値を増やしながら、成り立たないものを除去ししてゆき リスト pairs の要素が 1 つになったところで、処理を完了としています。
なお素数判定関数 prime-p は「素数判定関数作成関数」によって取得しています。
素数判定関数作成関数
(defun make-prime-number-p (primes) (lexical-let* ((primes (if (> 2 (length primes)) [2 3] (vconcat primes))) (capacity (length primes)) (index (1- capacity)) (expand-capacity (lambda () (setq primes (vconcat primes (make-vector capacity 0))) (setq capacity (* 2 capacity)))) (iterate) (iterate-2 (lambda () (when (= index (1- capacity)) (funcall expand-capacity)) (let ((n (aref primes index)) (r 0)) (while (zerop r) (incf n 2) (setq r 1) (do ((i 0 (1+ i)) (p 0) (end (truncate (sqrt n)))) ((or (>= p end) (zerop r))) (setq p (aref primes i)) (setq r (mod n p)))) (aset primes (incf index) n)))) (iterate-1 (lambda () (if (< index (1- capacity)) (aref primes (incf index)) (progn (setq iterate iterate-2) (funcall iterate-2))))) (binary-search (lambda (x) (let ((i 0) (j index) (exists)) (while (and (<= i j) (not exists)) (let* ((k (/ (+ i j) 2)) (m (aref primes k))) (cond ((< m x) (setq i (1+ k))) ((> m x) (setq j (1- k))) (t (setq exists t))))) exists))) (prime-number-p (lambda (x) (let ((last-prime (aref primes index))) (cond ((< x last-prime) (funcall binary-search x)) ((> x last-prime) (do ((p (funcall iterate) (funcall iterate))) ((>= p x) (= p x)))) (t t)))))) (setq iterate iterate-1) (lambda (x) (if (evenp x) nil (funcall prime-number-p x)))))
最大値のわからない素数判定に使える関数が欲しくて考えた末できたのがこれです。はじめに 素数を格納した vector primes を作成し、これに収まる整数 x についてはバイナリ探索で該当する素数があるかを判定します。primes よりもおおきな素数が渡されたときは、primes を拡張して対応できるようにしています。一度拡張するとそれよりさらに大きい数がくるまではバイナリ探索ですむので、大量のランダムな数を素数判定するには有利なんじゃないかと。(でも条件によるか・・・小さな素因数を持つ合成数の判定はむしろ遅くなるな、すごくダメな気がしてきた。普通に試し割りの方式に修正するかも）。
無駄に複雑に見えるかもしれませんが、実はもう少し発展させて便利なものをつくろうと画策しています。うまくいくかな。
Project Euler Problem 26
問題26: 1000未満のdにおいて 1/d を小数にした時、循環部の長さがもっとも長くなる d を求めよ。[=>問題文]
問題26の解答
気がつけば do ばっかり使ってる自分がいる。
;; Emacs Lisp (defun problem026 (limit) (let ((recurring-length ; 1/d の循環部の長さを求める関数 (lambda (d) (do ((table (make-hash-table)) (r 1) (count 0 (1+ count)) (result)) ((not (null result)) result) (cond ((zerop r) (setq result 0)) ((gethash r table) (setq result (- count (gethash r table)))) (t (puthash r count table) (setq r (mod (* r 10) d)))))))) ;; メインループ (do ((max-length 0) (result) (i (1- limit) (1- i))) ((<= i max-length) result) (let ((len (funcall recurring-length i))) (when (> len max-length) (setq max-length len) (setq result i)))))) (problem026 1000)
考え方のベースは割り算の筆算です。ここで重要なのは「商」として並ぶ数字列ではなく「余り」の部分。無限小数になるということはいつまでも割り切れない･･･つまり「余り」が０にならないことが条件になります。
筆算をしていった時、以前に出現した「余り」が再び現れると、それ以降は以前と同じ計算が繰り返されることになり、結果として循環した小数になります。同じ「余り」が出現するまでのループ回数を数えれば、それが循環部の長さになります。
recurring-length では、hash table に、「キー:余り、値:ループカウンタ」を保存してゆき、計算した「余り」がすでに登録されていた場合、保存されているカウンタ値と、現在のカウンタ値の差をとることで、求める循環部の長を得ています。
割り算の「余り」はその性質上「割る数」より小さい値にしかなりえません。したがって循環するための「余り」のバリエーションは「割る数」以上にはならないため、循環小数の循環部も「割る数」より長くはなりません。だから 1/d の d は大きい方からチェックを行い、既に計算済みの中での最大長(max-length)より d が小さくなったらそれ以下の d は調べる必要がありません。
Project Euler Problem 25
問題25: フィボナッチ数列の項 F(n) が千桁になる最初の n を求めよ。ただし F(1)=F(2)=1とする。[=>問題文]
問題25の解答
bigint ライブラリの出番です。
;; Emacs Lisp (require 'bigint) (defun problem025 (digit) (do ((fn (string-to-bigint "0")) (fn-1 (string-to-bigint "1")) (fn-2 (string-to-bigint "1")) (n 2 (1+ n))) ((= digit (bigint-length fn)) n) (setq fn (bigint+ fn-1 fn-2)) (setq fn-1 fn-2 fn-2 fn))) (problem025 1000)
bigint ライブラリに新たに桁数を取得する関数 bigint-length を追加しました。bigint-to-stringで文字列にして length を使えば長さは求まりますが処理コストに差が出ます。
多倍長整数(bigint)の桁数を取得する関数
(defun bigint-length (n) "桁数を取得" (let ((m (reverse n))) (+ (length (number-to-string (car m))) (* 4 (length (cdr m))))))
problem025 は、bigint-to-string & length で作ると 2.5秒、bigint-length で作ると 0.6 秒となりました。
