Haskellで最小の自然数を求めるアルゴリズム

このサイトは学部では早稲田で物理を, 修士では東大で数学を専攻し, 今も非アカデミックの立場で数学や物理と向き合っている一市民の奮闘の記録です. 運営者情報および運営理念についてはこちらをご覧ください.

中高の数学の復習から専門的な数学・物理までいろいろな情報を発信しています.
中高数学に関しては自然を再現しようや役に立つ中高数学中高数学お散歩コース
大学数学に関しては現代数学観光ツアーなどの無料の通信講座があります.
その他にも無料の通信講座はこちらのページにまとまっています.
ご興味のある方はぜひお気軽にご登録ください!

1 はじめに
- 1.1 追記
2 1. The smallest free number (最小自然数) 関数プログラミング珠玉のアルゴリズムデザイン

はじめに

数学や物理, 特に物理を学ぶのに関数プログラミング,
特に Haskell が非常にうまく使えるという話がある.
それを見て Haskell の勉強が急務だと思い色々な形で Haskell への知見を深めている.

例えば別企画で次のような記事 (群) もある.

そして何となくアルゴリズムをきちんと勉強しないといけないのかな,
とずっと思って放置したままになっている.
とりあえずいい加減何かやってみようと思い,
何となく買ってみたのが Richard Bird の pearls of functional algorithm design だ.

英語版を手に入れたのだが翻訳があったらしいと後で気付いた.

ウルトラ難しく 1 章を読んだ時点で既に心が折れている.
ただせっかく読んだ分は記録に残しておきたい.
他の数学、物理系プログラミングを進めていて,
うんざりしてきたら現実逃避でこちらを進めるというのもあっていいだろうとは思っている.

この本, Graham の本と同じスタイル,
つまり正確なコードを書かずに一部疑似コード (疑似記号?) を使っていて,
慣れていないとすぐに翻訳できないのが最高に鬱陶しい.
例えば第 1 章で言うと本の P.1 で出てくる $\notin$ は
not . flip elem のように表現しないといけないよう (本当によくわかっていない) だし,
本の P.2 で出てくる $\lor$ は本来の Haskell コードでは || だ.

先に書いたように最近また精力的に数学や物理, プログラミングのコンテンツを作っているけれども,
こういうのはよくないな, と痛感させられるコンテンツだ.
そういう意味でも厳しいコンテンツを見てみることは役に立つと思って,
無理せずゆるくやっていこう.

追記

翻訳者の山下伸夫さんからコメントを頂いてしまった.

Functional Pearls of Algorithm Designのプログラム表記と、ASCII記号によるHaskellプログラム表記との対応については、 http://pfad.sampou.org/ に追記しましたので、御笑覧くださいませ。

引用しておこう.

演算子記号とASCII文字列との対応

本書の記号: HaskellのASCII文字

≤: <=

≥: >=

∨: ||

∧: &&

⋅: .

≠: /=

∈: elem

∉: notElem

⊑: isPrefixOf

div: div

mod: mod

min: min

max: max

knows: knows

演算子以外の構文記号の一部については，GHCの言語拡張UnicodeSyntaxを有効にするとソースコード中に記述可能です．詳細については，GHCユーザーガイド 9.3.1 Unicode syntaxをご覧ください．

演算子に関しては，Unicodeの記号が使えますので，たとえば，1章については，以下のような定義モジュールをインポートすれば，そのままコードで表現できます．

{-# LANGUAGE UnicodeSyntax #-}
module Operators where

import qualified Data.List

infixr 5 \

(\) ∷ Eq a ⇒ [a] → [a] → [a]
(\) = (Data.List.\)

infix 4 ∈, ∉

(∈) ∷ Eq a ⇒ a → [a] → Bool
(∈) = elem

(∉) ∷ Eq a ⇒ a → [a] → Bool
(∉) = notElem

infixr 2 ∨

(∨) ∷ Bool → Bool → Bool
(∨) = (||)

infixr 3 ∧

(∧) ∷ Bool → Bool → Bool
(∧) = (&&)

infix 4 ≤, ≥

(≤) ∷ Ord a ⇒ a → a → Bool
(≤) = (<=)

(≥) ∷ Ord a ⇒ a → a → Bool
(≥) = (>=)

数学関係だとほとんどコメント来ないのに (素人くさい) プログラミングの話だと
やたら本質的に私が助かるコメントが来るのでありがたさしかない.

1. The smallest free number (最小自然数) 関数プログラミング珠玉のアルゴリズムデザイン

集合 $X \subset \mathbb{N}$ の最小値を求めようという問題を考えよう.
もちろん適当に派生的な問題もあるし現実的な応用もある.
解は $X$ がどのように表現されているかによる.

任意の数のリストを線型時間でソートするアルゴリズムを考えてみる.
次の minfree のような挙動をする関数を作りたい.


-- Int を Nat にしたい: Nat じたいは例えば Data.Nat にはある. 適当にインポートしたりする必要あり?
minfree :: [Int] -> Int
--minfree xs = head $ [0..] Data.List.\\ xs
minfree xs = head $ [0..] Main.\\ xs

-- Data.List.\\ にある関数: 本の記述では notElem が TeX でいう $\notin$ の記号で書かれていた.
(\\) :: Eq a => [a] -> [a] -> [a]
us \\ vs = filter (not . flip elem vs) us

main = do
  let a = [8, 23, 9, 0, 12, 11, 1, 10, 13, 7, 41, 4, 14, 21, 5, 17, 3, 19, 2, 6]
--  print $ Data.List.sort a
  print $ minfree a

RESULT

この関数 minfree は $O(n^2)$ のオーダーなので困る.
これをどうにかしたい.

Haskell の配列による解法

[0 .. length xs] の全ての要素が xs に入っているわけではないことに注意する.
xs に入っていない最小の数は filter (<= (length xs)) xs に入っていない最小の数だ.

結果から書くと次のようになる.
この実装は xs の要素が全て違う必要はないものの自然数であることは要請する.


import Data.Array

minfree' = search . checklist

search :: Array Int Bool -> Int
search = length . takeWhile id . elems

checklist :: [Int] -> Array Int Bool
checklist xs = accumArray (||) False (0, n) $ zip (filter (<= n) xs) (repeat True)
    where n = length xs

main = do
  print $ minfree' [8, 23, 9, 0, 12, 11, 1, 10, 13, 7, 41, 4, 14, 21, 5, 17, 3, 19, 2, 6]

RESULT

こんなさっくり書かれて意味がわかるわけもない.
当然 search と checklist をそれぞれじっくり追うしかない.
まずは checklist を見てみる.

checklist

たぶん慣れていないと accumArray がめんどい.
というかいまだにきちんと理解しきれていない.
次の節で別立てで見ることにして,
まずは accumArray に食わせるリストを作るところからはじめよう.
つまり zip (filter ( <= n) xs) (repeat True) だ.

repeat True は全要素がブーリアン True の無限リストを作る.
zip ではめて有限リストを切り出している.

zip にかけるもう一方のリストは filter (<= n) xs で作る.
これは特に言うことない.


main = do
  let xs = [8, 23, 9, 0, 12, 11, 1, 10, 13, 7, 41, 4, 14, 21, 5, 17, 3, 19, 2, 6]
  let n = length xs
  let a = zip (filter (<= n) xs) (repeat True)
  -- print $ repeat True
  print $ filter (<= n) xs
  print a

RESULT

[8,9,0,12,11,1,10,13,7,4,14,5,17,3,19,2,6]
[(8,True),(9,True),(0,True),(12,True),(11,True),(1,True),(10,True),(13,True),(7,True),(4,True),(14,True),(5,True),(17,True),(3,True),(19,True),(2,True),(6,True)]

これでタプルのリストができる.
このタプルのリストを accumArray にくわせる.

`accumArray`

多分これが一番めんどい.
泥臭くコード片を追いかけないとわからない.
ふだん (手続き型で?) 書くときはループで書くところを関数でスパっと書いている感じ.
逐次 print デバッグみたいなのがやりづらいので何をやっているかまだ追いづらい.
多分慣れればもっとサクサク書けるのだろうなとは思う.

とにかく泥臭く追いかけよう.
foldr などと同じくコレクションに関数を順次適用させて新たなコレクションを作る系の処理だ.
実コード例


import Data.Array

main = do
  -- 配列の添字は [1..3]
  -- 配列の各要素は第 2 引数 0 で初期化される
  let a = [(1, 2), (3, 4), (1, 5)]
  print $ accumArray (+) 0 (1,3) a

RESULT

array (1,3) [(1,7),(2,0),(3,4)]

accumArray の第 1 引数は当てる関数,
第 2 引数は初期値, 第 3 引数は新たな配列の添字の範囲,
第 4 引数はなめていく配列だ.
Python あたりの標準的な記法を使って配列、リストを表すことにしよう.
accumArray の結果の配列を b とする.

b[1] からはじまるので b[1] がどうなるかを調べよう.
結果から行くと a[0] と a[2] の第 1 要素の和,
つまり b[1] = a[0][1] + a[2][1] になる.
これを見ていく.

基本は a をなめていくのだ.

a[0][0] を調べると 1: これと b[1] の値を関連づける.
初期値は 0 で + を使うことになる.
+ のもう一方の引数として a[0][1] を使う.
結果として b[1] = 0 + a[0][1] になる.
次に a[1][0] を調べると 3: これと b[3] の値を関連づける.
初期値は 0 で + を使うことになる.
+ のもう一方の引数として a[1][1] を使う.
結果として b[3] = 0 + a[1][1] になる.
次に a[2][0] を調べると 1: これと b[1] の値を関連づける.
初期値は 0 で + を使うことになる.
+ のもう一方の引数として a[2][1] を使う.
結果として b[1] = 0 + a[0][1] + a[2][1] になる.
一度も出てこない b[2] は初期値そのまま b[2] = 0.

これが上の accumArray の結果.

checklist は 0 から $n$ 個までの $n+1$ スロットある.
最初は全て False だ.

search

改めて書くと次のように定義される.


minfree' = search . checklist

search :: Array Int Bool -> Int
search = length . takeWhile id . elems

checklist :: [Int] -> Array Int Bool
checklist xs = accumArray (||) False (0, n) $ zip (filter (<= n) xs) (repeat True)
    where n = length xs

Boolean の配列 (Data.Array) から整数を取る.
この整数はもとのリストの最小値にあたる.
この他に関数の合成 (.) も使っている.

elems から見ていこう.
これは配列の要素をリストに格納して返す.


import Data.Array

main = do
  let a = listArray (0,3) [1,2,3,4]
  print $ a
  print $ elems a

RESULT

array (0,3) [(0,1),(1,2),(2,3),(3,4)]
[1,2,3,4]

これも特に問題ないだろう.
よくわからないのは id.
なぜこれを挟む必要があるのだろうか?
外したらエラーになった.
takeWhile に predicate を食わせないといけないからその役割なのだろうか.
そう思えば必要な理由は理解できる.
今の私のレベルでは自分で書こうとすると多分解決できずそのまま死にそう.

takeWhile は predicate を食わせてそれが True を返すうちはリストを返し続けてくれる関数.
False が来たら打ち切ってそこまでの長さを測れば目的の値が取れる.

分割統治法 (divide & conquer)

xs の要素は全て異なることを前提にしている.
基本的なアイデアはリストを分けてそのリストのそれぞれに対処すること.

最初のコードは次の通り.
徐々に改善していく.


import Data.List
minfreeDivideAndConquer b xs =
  if null ([0..b-1] Data.List.\\ us)
  then head ([b..] Data.List.\\ vs)
  else head ([0..] Data.List.\\ us)
  where (us, vs) = partition (< b) xs

main = do
  let a = [8, 23, 9, 0, 12, 11, 1, 10, 13, 7, 41, 4, 14, 21, 5, 17, 3, 19, 2, 6]
  print $ Data.List.sort a
  print $ minfreeDivideAndConquer 4 a

RESULT

[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,17,19,21,23,41]
15

いくつか関数の挙動を調べておこう.
まずは null.


main = do
  print $ null []
  print $ null [1]

RESULT

True
False

次に partition.


import Data.List
main = do
  print $ Data.List.partition (< 2) [1,2,3,4]

RESULT

([1],[2,3,4])

最初の改良点は null ([0..b-1] Data.List.\\ us) で,
$O(n^2)$ かかるこの処理を軽くしたい.
リストに重複がなければ null ([0..b-1] \\ us) == (length us == b) であることを使う.
ここでリストに重複がないことを使っている.
これと次の minfrom を使う.


minfrom :: Int -> [Int] -> Int
minfrom a xs = head ([a..] \\ xs)

RESULT

そこで次のコードが出てくる.


import Data.List
b = 5
minfreeDivideAndConquer' xs = minfrom 0 xs
minfrom a xs | null xs            = a
             | length us == b - a = minfrom b vs
             | otherwise          = minfrom a us
  where (us, vs) = partition (< b) xs

main = do
  print $ minfreeDivideAndConquer' [8, 23, 9, 0, 12, 11, 1, 10, 13, 7, 41, 4, 14, 21, 5, 17, 3, 19, 2, 6]

このコードを runhaskell で動かそうと思ったら動かない.
エラーになるわけではなく処理が回りっぱなしという感じ.
何なのだろう.

何はともあれ, 上のコードは b=5 のハードコードが気に入らない.
これを消した以外にも length us を何度も計算しないように
where に押し込めたりと工夫したのが次のコード.


import Data.List
minfreeDivideAndConquer'' xs = minfrom 0 (length xs, xs)
minfrom a (n, xs) | n == 0     = a
                  | m == b - a = minfrom b (n - m, vs)
                  | otherwise  = minfrom a (m, us)
  where (us, vs) = partition (< b) xs
        b        = a + 1 + n `div` 2
        m        = length us

main = do
  print $ minfreeDivideAndConquer'' [8, 23, 9, 0, 12, 11, 1, 10, 13, 7, 41, 4, 14, 21, 5, 17, 3, 19, 2, 6]

RESULT

今回はこれでおしまい.

1. The smallest free number (最小自然数) 関数プログラミング珠玉のアルゴリズムデザイン

はじめに

追記

1. The smallest free number (最小自然数) 関数プログラミング珠玉のアルゴリズムデザイン

Haskell の配列による解法

checklist

`accumArray`

search

分割統治法 (divide & conquer)

関連記事

『学部生のための Sage』に対抗して『よくわからない Sage』を作ろう

【フランスは、会社のデスクで飯食うのは法的に厳禁です】あと CNRS 情報

漸化式から微分へ: 微分方程式のシミュレーションの観点から／中高数学駆け込み寺第 8 回

手書きストロークをフーリエ級数展開し円運動として可視化する Char2Fourier

数列とは何か？プログラム解説中高数学駆け込み寺

TeX: 式の中で文章を入れて改行させたいときは parbox か minipage

このサイトについて

YouTube チャンネル登録

メルマガ登録

役に立つ・面白い記事があればクリックを！

最近の投稿

最近のコメント

アーカイブ

カテゴリー

メタ情報

1. The smallest free number (最小自然数) 関数プログラミング 珠玉のアルゴリズムデザイン

はじめに

追記

1. The smallest free number (最小自然数) 関数プログラミング 珠玉のアルゴリズムデザイン

Haskell の配列による解法

checklist

accumArray

search

分割統治法 (divide & conquer)

共有:

関連記事

『学部生のための Sage』に対抗して『よくわからない Sage』を作ろう

【フランスは、会社のデスクで飯食うのは法的に厳禁です】あと CNRS 情報

漸化式から微分へ: 微分方程式のシミュレーションの観点から／中高数学駆け込み寺 第 8 回

手書きストロークをフーリエ級数展開し円運動として可視化する Char2Fourier

数列とは何か？プログラム解説 中高数学駆け込み寺

TeX: 式の中で文章を入れて改行させたいときは parbox か minipage

このサイトについて

YouTube チャンネル登録

メルマガ登録

役に立つ・面白い記事があればクリックを！

最近の投稿

最近のコメント

アーカイブ

カテゴリー

メタ情報

1. The smallest free number (最小自然数) 関数プログラミング珠玉のアルゴリズムデザイン

1. The smallest free number (最小自然数) 関数プログラミング珠玉のアルゴリズムデザイン

`accumArray`

漸化式から微分へ: 微分方程式のシミュレーションの観点から／中高数学駆け込み寺第 8 回

数列とは何か？プログラム解説中高数学駆け込み寺