« 2009年12月 | トップページ | 2010年2月 »

2010年1月

2010年1月21日 (木)

Project Euler : Problem 12 ~ 三角数の約数の数 (2)

 ちょっと間が空きましたが、Problem 12 の続きです。

 問題はこちらをご覧ください。

 

 実際の解法を示す前に、前提となる条件をいくつか挙げておきます。

  1. n 番目の三角数は Tn = n * (n + 1) / 2 で表せる。
  2. 隣り合う正の整数 n と n + 1 は「互いに素」である。ここで、n が偶数の場合は n / 2 と n + 1 が、n が奇数の場合は n と (n + 1) / 2 がそれぞれ「互いに素」となる。
  3. 正の整数 x の約数の数を f(x) とすると、a と b が「互いに素」であれば、f(a * b) = f(a) * f(b) となる。
  4. 以上のことより、n が偶数の場合は f(Tn) = f(n / 2) * f(n + 1) となり、n が奇数の場合は f(Tn) = f(n) * f((n + 1) / 2) となる。
 ここで、f(Tn) を小さい方から見ていくと、次のようになります。
  f(T1) = f(1) * f(2/2)
  f(T2) = f(2/2) * f(3)
  f(T3) = f(3) * f(4/2)
  ...
 よく見ると、前の式の一部が必ず次の式に出てきています。ということは、小さい方から三角数の個数を調べていく場合、計算結果の一部を次に引き継ぐことが出きるので、計算量の節約になります。

 以上のことを踏まえて、Haskell のコードを書いてみました。

import Data.List (group) -- 素因数分解 factorize :: Integral a => a -> [(a, Int)] factorize 1 = [(1, 0)] factorize n = format $ factorize' n divs where divs = 2 : 3 : [x + y | x <- [6, 12 ..], y <- [-1, 1]] factorize' n ps@(p : ps') | n < p * p = [n] | rem n p == 0 = p : factorize' (div n p) ps | otherwise = factorize' n ps' format ps = [(head xs, length xs) | xs <- group ps] -- 約数の個数 countOfDivs :: Integral a => a -> Int countOfDivs n = product [b + 1 | (_, b) <- factorize n] -- 多角数の一般項 from 『Wikipedia』, modified by Tsumuji polyNum :: Integral a => a -> a -> a polyNum p n = div (p * n * (n - 1)) 2 - n * (n - 2) problem012 :: Int -> Int problem012 num = loop (countOfDivs 1) 2 where loop d1 n | d1 * d2 >= num = polyNum 3 (n - 1) | otherwise = loop d2 (n + 1) where d2 = countOfDivs (if even n then (div n 2) else n) main = print $ problem012 501
 最適化オプションを付けてコンパイルしたものを、 Pentium M, 1500MHz のノートパソコンで実行すると、次のようになりました。
$ time ./problem012 76576500 real 0m0.036s user 0m0.020s sys 0m0.000s
 前回のコードの約 5 倍の速度です。自分でもちょっとびっくりでした。

2010年1月 7日 (木)

Project Euler : Problem 12 ~ 三角数の約数の数 (1)

 問題はこちらをご覧ください。

 

 今回は三角数の数列の問題です。
 問題では「三角数は自然数の和の形で表わせる」と説明されています。これを基に三角数を定義すると

triNum :: Integral a => a -> a triNum n = div ((1 + n) * n) 2
ということになります。
 しかし今回は、私が作った多角数の無限リストを返す関数を使います。
-- 多角数のリスト polyNumList :: Integral a => a -> [a] polyNumList n = poly 1 (n - 1) (n - 2) where poly v add d = v : poly (v + add) (add + d) d
 例えば n を 3 にすると三角数の無限リストを返します。
*Main> take 10 $ polyNumList 3 [1,3,6,10,15,21,28,36,45,55]

 約数の個数は、実際の約数を調べなくても、素因数分解の結果を利用すると簡単に出せます。

-- 約数の個数 numOfDivs :: Integral a => a -> Int numOfDivs n = product [b + 1 | (_, b) <- factorize n]

 これらをまとめると、次のようになりました。

import Data.List (group) -- 素因数分解 (要 : priems) factorize :: Integer -> [(Integer, Int)] factorize 1 = [(1, 0)] factorize n = format $ factorize' n divs where factorize' n xs'@(x : xs) | n < x * x = [n] | rem n x == 0 = x : factorize' (div n x) xs' | otherwise = factorize' n xs format xss = [(head xs, length xs) | xs <- group xss] divs :: [Integer] divs = 2 : 3 : [x + y | x <- [6, 12 ..], y <- [-1, 1]] -- 約数の個数 numOfDivs :: Integer -> Int numOfDivs n = product [b + 1 | (_, b) <- factorize n] -- 多角数のリスト polyNumList :: Integer -> [Integer] polyNumList n = poly 1 (n - 1) (n - 2) where poly v add d = v : poly (v + add) (add + d) d problem012 :: Integer problem012 = head $ dropWhile ((<= 500) . numOfDivs) triNums where triNums = polyNumList 3 main = print problem012

 最適化オプションを付けてコンパイルしたものを、 Pentium M, 1500MHz のノートパソコンで実行すると、約 0.15 秒で答えが出ました。

 今回の方法でもそこそこ速いのですが、実はもっと速く答えが出せるやり方があります(同じノートパソコンで約 0.05 秒)。
 その具体的な方法は、次回ということで……。

 

追記:
 rst76 さんのところで、素因数分解のコードの「Haskell っぽいコードってどんなだろう」という話題が出ていました。
 rst76 さんは「高階関数」を挙げていらっしゃいますが、私自身は以前、Scheme をいじっていたので、「高階関数」は普通に使っていたしなあ……。
 ということで、Haskell 初心者の私が Ruby や Scheme のコードを書く場合と比べて考えるに、「if の使用頻度」と「リスト内包表記」かなと思います。
 Haskell のコードを書くときには、「パターンマッチ」や「ガード」が便利過ぎて、関数定義をする際に if を使う機会がめっきり減りました。
 しかも、「パターンマッチ」、「ガード」、「リスト内包表記」を使うとすっきりして見やすくなるような気がするので、つい多用してしまいます。

2010年1月 2日 (土)

Project Euler : Problem 11

 問題はこちらをご覧ください。

 

 このブログを見た人の中には、「なぜ、すでに Scheme や Ruby で解いた問題を、また Haskell で解き直しているのだろう?」と思われる方もいらっしゃるかもしれません。
 その理由は、「Haskell らしい考え方や、Haskell らしいコードの書き方」を探すことにあります。
 これまで Haskell で解いた問題では、解法のアルゴリズムは以前のものと変わっていないものがほとんどです。それでも、特定のアルゴリズムを Haskell でコーディングするには、当然 Scheme や Ruby とは違ったアプローチが必要になる事も多いです。
 また、「Haskell らしい書式」に関して、試行錯誤している段階でもあります。

 

 今回の問題ですが、まずは手っ取り早くデータを二次元リストに放り込んで解いてみました。

nums :: [[Int]] nums = [[08,02,22,97,38,15,00,40,00,75,04,05,07,78,52,12,50,77,91,08], [49,49,99,40,17,81,18,57,60,87,17,40,98,43,69,48,04,56,62,00], [81,49,31,73,55,79,14,29,93,71,40,67,53,88,30,03,49,13,36,65], [52,70,95,23,04,60,11,42,69,24,68,56,01,32,56,71,37,02,36,91], [22,31,16,71,51,67,63,89,41,92,36,54,22,40,40,28,66,33,13,80], [24,47,32,60,99,03,45,02,44,75,33,53,78,36,84,20,35,17,12,50], [32,98,81,28,64,23,67,10,26,38,40,67,59,54,70,66,18,38,64,70], [67,26,20,68,02,62,12,20,95,63,94,39,63,08,40,91,66,49,94,21], [24,55,58,05,66,73,99,26,97,17,78,78,96,83,14,88,34,89,63,72], [21,36,23,09,75,00,76,44,20,45,35,14,00,61,33,97,34,31,33,95], [78,17,53,28,22,75,31,67,15,94,03,80,04,62,16,14,09,53,56,92], [16,39,05,42,96,35,31,47,55,58,88,24,00,17,54,24,36,29,85,57], [86,56,00,48,35,71,89,07,05,44,44,37,44,60,21,58,51,54,17,58], [19,80,81,68,05,94,47,69,28,73,92,13,86,52,17,77,04,89,55,40], [04,52,08,83,97,35,99,16,07,97,57,32,16,26,26,79,33,27,98,66], [88,36,68,87,57,62,20,72,03,46,33,67,46,55,12,32,63,93,53,69], [04,42,16,73,38,25,39,11,24,94,72,18,08,46,29,32,40,62,76,36], [20,69,36,41,72,30,23,88,34,62,99,69,82,67,59,85,74,04,36,16], [20,73,35,29,78,31,90,01,74,31,49,71,48,86,81,16,23,57,05,54], [01,70,54,71,83,51,54,69,16,92,33,48,61,43,52,01,89,19,67,48]] -- 数列の生成 same :: Int -> [Int] same i = [i, i, i, i] plus :: Int -> [Int] plus i = [i + n | n <- [0 .. 3]] minus :: Int -> [Int] minus i = [i - n | n <- [0 .. 3]] -- 4 個の数の積 calc :: [Int] -> [Int] -> Int calc xs ys = product $ map ref $ zip xs ys where ref (x, y) = nums !! y !! x -- 横の積の集合 line1 :: [Int] line1 = [calc (plus x) (same y) | x <- [0 .. 16], y <- [0 .. 19]] -- 縦の積の集合 line2 :: [Int] line2 = [calc (same x) (plus y) | x <- [0 .. 19], y <- [0 .. 16]] -- 斜めの積の集合 line3 :: [Int] line3 = [calc (minus x) (plus y) | x <- [3 .. 19], y <- [0 .. 16]] -- 斜めの積の集合 line4 :: [Int] line4 = [calc (plus x) (plus x) | x <- [0 .. 16], y <- [0 .. 16]] problem011 :: Int problem011 = maximum (line1 ++ line2 ++ line3 ++ line4) main = print problem011
 わかりやすさに重点を置いたため、少々長いコードになりました。

 リストより配列の方がデータにアクセスする時間は短いだろうと思い、同じアルゴリズムで「配列版」も作ってみました。

import Data.Array nums :: Array (Int, Int) Int nums = listArray ((0, 0), (19, 19)) [08,02,22,97,38,15,00,40,00,75,04,05,07,78,52,12,50,77,91,08, 49,49,99,40,17,81,18,57,60,87,17,40,98,43,69,48,04,56,62,00, 81,49,31,73,55,79,14,29,93,71,40,67,53,88,30,03,49,13,36,65, 52,70,95,23,04,60,11,42,69,24,68,56,01,32,56,71,37,02,36,91, 22,31,16,71,51,67,63,89,41,92,36,54,22,40,40,28,66,33,13,80, 24,47,32,60,99,03,45,02,44,75,33,53,78,36,84,20,35,17,12,50, 32,98,81,28,64,23,67,10,26,38,40,67,59,54,70,66,18,38,64,70, 67,26,20,68,02,62,12,20,95,63,94,39,63,08,40,91,66,49,94,21, 24,55,58,05,66,73,99,26,97,17,78,78,96,83,14,88,34,89,63,72, 21,36,23,09,75,00,76,44,20,45,35,14,00,61,33,97,34,31,33,95, 78,17,53,28,22,75,31,67,15,94,03,80,04,62,16,14,09,53,56,92, 16,39,05,42,96,35,31,47,55,58,88,24,00,17,54,24,36,29,85,57, 86,56,00,48,35,71,89,07,05,44,44,37,44,60,21,58,51,54,17,58, 19,80,81,68,05,94,47,69,28,73,92,13,86,52,17,77,04,89,55,40, 04,52,08,83,97,35,99,16,07,97,57,32,16,26,26,79,33,27,98,66, 88,36,68,87,57,62,20,72,03,46,33,67,46,55,12,32,63,93,53,69, 04,42,16,73,38,25,39,11,24,94,72,18,08,46,29,32,40,62,76,36, 20,69,36,41,72,30,23,88,34,62,99,69,82,67,59,85,74,04,36,16, 20,73,35,29,78,31,90,01,74,31,49,71,48,86,81,16,23,57,05,54, 01,70,54,71,83,51,54,69,16,92,33,48,61,43,52,01,89,19,67,48] -- 数列の生成 same :: Int -> [Int] same i = [i, i, i, i] plus :: Int -> [Int] plus i = [i + n | n <- [0 .. 3]] minus :: Int -> [Int] minus i = [i - n | n <- [0 .. 3]] -- 4 個の数の積 calc :: [Int] -> [Int] -> Int calc xs ys = product $ map (nums !) $ zip xs ys -- 横の積の集合 line1 :: [Int] line1 = [calc (plus x) (same y) | x <- [0 .. 16], y <- [0 .. 19]] -- 縦の積の集合 line2 :: [Int] line2 = [calc (same x) (plus y) | x <- [0 .. 19], y <- [0 .. 16]] -- 斜めの積の集合 line3 :: [Int] line3 = [calc (minus x) (plus y) | x <- [3 .. 19], y <- [0 .. 16]] -- 斜めの積の集合 line4 :: [Int] line4 = [calc (plus x) (plus x) | x <- [0 .. 16], y <- [0 .. 16]] problem011 :: Int problem011 = maximum (line1 ++ line2 ++ line3 ++ line4) main = print problem011
 結果はと言うと……「リスト版」と「配列版」の実効速度はほとんど一緒でした。扱うデータの量がもっと多くなれば、結果も違ってくるのでしょうか?

 自分なりの「Haskell らしい書式」は、まだまだ見つかっていませんが、「実効速度や行数よりも、わかりやすさやメンテナンスのしやすさを優先していくほうがいいのかな?」という気はしてきています。

« 2009年12月 | トップページ | 2010年2月 »

2016年7月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト

最近のトラックバック

無料ブログはココログ