前回

前回は挿入ソートの正当性を証明しました．
今回は同様の方法でヒープソートの正当性を証明します．
ただし，ヒープソートは難易度が桁違いなので前編ではソートされることだけを示します．

リポジトリ

リポジトリはここにあります．
今回作成したソースコードは以下の3つです．

Heap.agda ヒープソートの実装
HeapProp/HeapMerge.agda ヒープ木のマージにおける諸性質，補題
HeapProp/HeapSorted.agda ヒープソートによってソートされることの証明

ヒープソートのアルゴリズム

ヒープソートはヒープ木を作って分解することでソートする手法です．
C言語などではヒープ木の実装に配列を用いますが，Agdaでは配列が無いので別の手段が必要です(配列があったとしても証明が面倒そうなので別の手法が欲しい)．
といっても別に複雑な仕組みではなく代数的データ型で2分木を作ってそれをヒープ木に見たてればいいです．

今回は以下のように実装しました．

Heap.agda

data Heap (A : Set) : Set where
    E : Heap A
    ∧ : A → Heap A → Heap A → Heap A

mergeHeap : Heap ℕ → Heap ℕ → Heap ℕ
mergeHeap E y = y
mergeHeap x E = x
mergeHeap x@(∧ n rx lx) y@(∧ m ry ly) with n ≤? m
... | yes n≤m = ∧ n (mergeHeap y lx) rx
... | no n≰m  = ∧ m (mergeHeap x ly) ry

insertHeap : ℕ → Heap ℕ → Heap ℕ
insertHeap n h = mergeHeap (∧ n E E) h

makeHeap : List ℕ → Heap ℕ
makeHeap = foldr insertHeap E

fromHeap : ℕ → Heap ℕ → List ℕ
fromHeap zero x = []
fromHeap (suc n) E = []
fromHeap (suc n) (∧ x l r) = x ∷ fromHeap n (mergeHeap l r)

heapSize : {A : Set} → Heap A → ℕ
heapSize E = zero
heapSize (∧ x l r) = suc (heapSize l + heapSize r)

hsort : List ℕ → List ℕ
hsort xs = fromHeap (heapSize h) h
    where h = makeHeap xs

step-insert≥ : {n x : ℕ}{rx lx : Heap ℕ} → ¬ (n ≤ x) → insertHeap n (∧ x rx lx) ≡ ∧ x (mergeHeap (∧ n E E) lx) rx
step-insert≥ {n}{x}{rx}{lx} n≰x with n ≤? x
... | yes n≤x = ⊥-elim (n≰x n≤x)
... | no n≰x′ = refl

step-insert≤ : {n x : ℕ}{rx lx : Heap ℕ} → n ≤ x → insertHeap n (∧ x rx lx) ≡ ∧ n (∧ x rx lx) E
step-insert≤ {n}{x}{rx}{lx} n≤x with n ≤? x
... | yes n≤x′ = refl
... | no n≰x  = ⊥-elim (n≰x n≤x)

step-merge≤ : {x y : ℕ}{lx rx ly ry : Heap ℕ} → x ≤ y
        → mergeHeap (∧ x lx rx) (∧ y ly ry) ≡ (∧ x (mergeHeap (∧ y ly ry) rx) lx)
step-merge≤ {x}{y} x≤y with x ≤? y
... | yes x≤y′ = refl
... | no x≰y  = ⊥-elim (x≰y x≤y)

step-merge≥ : {x y : ℕ}{lx rx ly ry : Heap ℕ} → ¬ (x ≤ y)
        → mergeHeap (∧ x lx rx) (∧ y ly ry) ≡ (∧ y (mergeHeap (∧ x lx rx) ry) ly)
step-merge≥ {x}{y} x≰y with x ≤? y
... | yes x≤y = ⊥-elim (x≰y x≤y)
... | no x≰y′  = refl

証明

前回と同様に適当な補題を示して，目的の命題を証明します．前回と異なりヒープソートは中間にヒープ木を挟むので2つの補題が必要になります．

makeHeapによって正しくヒープ木が作られること
fromHeapに正しいヒープ木を与えられるとソート列が得られること

新たに"正しいヒープ木である"という述語が必要になるので，それを定義します．

HeapMerge.agda

data IsHeap : Heap ℕ → Set where
    h[] : IsHeap E
    h[n] : (n : ℕ) → IsHeap (∧ n E E)
    hE← : {x y : ℕ} {rx lx : Heap ℕ} → IsHeap (∧ x rx lx)
        → y ≤ x → IsHeap (∧ y E (∧ x rx lx))
    hE→ : {x y : ℕ} {rx lx : Heap ℕ} → IsHeap (∧ x rx lx)
        → y ≤ x → IsHeap (∧ y (∧ x rx lx) E)
    h∧ : {x y z : ℕ} {rx lx ry ly : Heap ℕ}
        → IsHeap (∧ x rx lx) → IsHeap (∧ y ry ly)
        → z ≤ x × z ≤ y → IsHeap (∧ z (∧ x rx lx) (∧ y ry ly))

IsHeapはSortedによく似てますが，ヒープ木は2分木であるためコンストラクタが2つ増えました．
そして単純な方法ではこのコンストラクタによって場合分けをして補題を証明することになるので非常に面倒です．しかし，他の方法を思いつかなかったので力技で示しました．
実際に証明した補題が以下になります．

HeapMerge.agda

lem-insertHeap : (n : ℕ) → (x : Heap ℕ) → IsHeap x → IsHeap (insertHeap n x)
lem-insertHeap n .E h[] = h[n] n
lem-insertHeap n .(∧ x E E) (h[n] x) with n ≤? x
... | yes n≤x = hE→ (h[n] x) n≤x
... | no n≰x  = hE→ (h[n] n) (lem-≤ n≰x)
lem-insertHeap n (∧ x E (∧ y ry ly)) (hE← hy x≤y) with n ≤? x
... | yes n≤x = hE→ (hE← hy x≤y) n≤x
... | no n≰x with n ≤? y
...     | yes n≤y = hE→ (subst IsHeap (step-insert≤ n≤y) (lem-insertHeap n (∧ y ry ly) hy)) (lem-≤ n≰x)
...     | no n≰y  = hE→ (subst IsHeap (step-insert≥ n≰y) (lem-insertHeap n (∧ y ry ly) hy)) x≤y
lem-insertHeap n (∧ x (∧ y ry ly) E) (hE→ hy x≤y) with n ≤? x
... | yes n≤x = hE→ (hE→ hy x≤y) n≤x
... | no n≰x  = h∧ (h[n] n) hy (lem-≤ n≰x , x≤y)
lem-insertHeap n (∧ x (∧ y ry ly) (∧ z rz lz)) (h∧ hy hz (x≤y , x≤z)) with n ≤? x
... | yes n≤x = hE→ (h∧ hy hz (x≤y , x≤z)) n≤x
... | no n≰x with n ≤? z
...     | yes n≤z = h∧ (subst IsHeap (step-insert≤ n≤z) (lem-insertHeap n (∧ z rz lz) hz)) hy (lem-≤ n≰x , x≤y)
...     | no n≰z  = h∧ (subst IsHeap (step-insert≥ n≰z) (lem-insertHeap n (∧ z rz lz) hz)) hy (x≤z , x≤y)


lem-mergeHeap : (x y : Heap ℕ) → IsHeap x → IsHeap y → IsHeap (mergeHeap x y)
-- 省略

lem-makeHeap : (xs : List ℕ) → IsHeap (makeHeap xs)
lem-makeHeap [] = h[]
lem-makeHeap (x ∷ xs) = lem-insertHeap x (makeHeap xs) (lem-makeHeap xs)

lem-mergreHeapは引数にIsHeapが2つあり，それぞれ場合分けして証明するととても長くなりますが，方針としてはlem-insertHeapと全く同じなので省略しました．
次にヒープ木をfromHeapで分解するとソート列が得られることを示すのですが，fromHeapにおける再帰呼出の引数はmergeHeap x yの形なのでコンストラクタが減っているかが自明ではありません．
なのでヒープ木の大きさに関する帰納法を明示的に書く必要があります．
実際に書くと以下の様になります．

HeapSorted.agda

mergeSize : {x y : Heap ℕ} → heapSize (mergeHeap x y) ≡ heapSize x + heapSize y
mergeSize {E} {y} = refl
mergeSize {∧ x l₁ r₁} {E} = cong suc (sym (+-identityʳ (heapSize l₁ + heapSize r₁)))
mergeSize {∧ x l₁ r₁} {∧ y l₂ r₂} with x ≤? y
... | yes x≤y = cong suc (
            heapSize (mergeHeap (∧ y l₂ r₂) r₁) + heapSize l₁
        ≡⟨ cong (λ e → e + heapSize l₁) (mergeSize {∧ y l₂ r₂} {r₁}) ⟩
            suc (heapSize l₂ + heapSize r₂ + heapSize r₁ + heapSize l₁)
        ≡⟨ cong suc (+-assoc (heapSize l₂ + heapSize r₂) (heapSize r₁) (heapSize l₁)) ⟩
            suc (heapSize l₂ + heapSize r₂ + (heapSize r₁ + heapSize l₁))
        ≡⟨ cong suc (+-comm (heapSize l₂ + heapSize r₂) (heapSize r₁ + heapSize l₁)) ⟩
            suc (heapSize r₁ + heapSize l₁ + (heapSize l₂ + heapSize r₂))
        ≡⟨ cong (λ e → suc (e + (heapSize l₂ + heapSize r₂))) (+-comm (heapSize r₁) (heapSize l₁)) ⟩
            suc (heapSize l₁ + heapSize r₁ + (heapSize l₂ + heapSize r₂))
        ≡⟨ sym (+-suc (heapSize l₁ + heapSize r₁) (heapSize l₂ + heapSize r₂)) ⟩
            heapSize l₁ + heapSize r₁ + suc (heapSize l₂ + heapSize r₂)
        ∎)
... | no x≰y  = cong suc (
            heapSize (mergeHeap (∧ x l₁ r₁) r₂) + heapSize l₂
        ≡⟨ cong (λ e → e + heapSize l₂) (mergeSize {∧ x l₁ r₁} {r₂}) ⟩
            suc (heapSize l₁ + heapSize r₁ + heapSize r₂ + heapSize l₂)
        ≡⟨ cong suc (+-assoc (heapSize l₁ + heapSize r₁) (heapSize r₂) (heapSize l₂)) ⟩
            suc (heapSize l₁ + heapSize r₁ + (heapSize r₂ + heapSize l₂))
        ≡⟨ cong (λ e → suc (heapSize l₁ + heapSize r₁ + e)) (+-comm (heapSize r₂) (heapSize l₂)) ⟩
            suc (heapSize l₁ + heapSize r₁ + (heapSize l₂ + heapSize r₂))
        ≡⟨ sym (+-suc (heapSize l₁ + heapSize r₁) (heapSize l₂ + heapSize r₂)) ⟩
            heapSize l₁ + heapSize r₁ + suc (heapSize l₂ + heapSize r₂)
        ∎)

mergeSize< : {x : ℕ}{l r : Heap ℕ} → heapSize (mergeHeap l r) < heapSize (∧ x l r)
mergeSize< {x}{l}{r} rewrite mergeSize {l}{r} = s≤s ≤-refl

_⟪_ : {A : Set} → Rel (Heap A) Level.zero
_⟪_ = _<_ on heapSize

⟪-well-founded : {A : Set} → Well-founded (_⟪_ {A})
⟪-well-founded = Inverse-image.well-founded heapSize <-well-founded

module _ {ℓ} {A : Set} where
  open WF.All (⟪-well-founded {A}) ℓ public
    renaming ( wfRec-builder to ⟪-rec-builder
             ; wfRec to ⟪-rec
             )

P : ∀ x → Set
P x = IsHeap x → Sorted (fromHeap (heapSize x) x)

step : ∀ x → (∀ y → y ⟪ x → P y) → P x
step .E rec h[] = s[]
step .(∧ n E E) rec (h[n] n) = s[n]
step (∧ x E (∧ y l r)) rec (hE← h x≤y) = s∷ x≤y (rec (∧ y l r) (s≤s (s≤s ≤-refl)) h)
step (∧ x (∧ y l r) E) rec (hE→ h x≤y) = s∷ x≤y (subst Sorted (cong
    (λ e → y ∷ fromHeap e (mergeHeap l r)) (sym (+-identityʳ (heapSize l + heapSize r))))
    (rec (∧ y l r) (s≤s (s≤s (m≤m+n (heapSize l + heapSize r) zero))) h))
step (∧ x t₁@(∧ x₁ l₁ r₁) t₂@(∧ x₂ l₂ r₂)) rec (h∧ h₁ h₂ (x≤x₁ , x≤x₂))
    with x₁ ≤? x₂
... | yes x₁≤x₂ = s∷ x≤x₁ (subst Sorted (
            fromHeap
            (heapSize (mergeHeap (∧ x₁ l₁ r₁) (∧ x₂ l₂ r₂)))
            (mergeHeap (∧ x₁ l₁ r₁) (∧ x₂ l₂ r₂))
        ≡⟨ cong₂ (λ e₁ e₂ → fromHeap (heapSize e₁) e₂) (step-merge≤ x₁≤x₂) (step-merge≤ x₁≤x₂) ⟩
            x₁ ∷ fromHeap
            (heapSize (mergeHeap (∧ x₂ l₂ r₂) r₁) + heapSize l₁)
            (mergeHeap (mergeHeap (∧ x₂ l₂ r₂) r₁) l₁)
        ≡⟨ cong
            (λ e → x₁ ∷ fromHeap e (mergeHeap (mergeHeap (∧ x₂ l₂ r₂) r₁) l₁))
            (   (heapSize (mergeHeap (∧ x₂ l₂ r₂) r₁) + heapSize l₁)
            ≡⟨ cong (λ e → e + heapSize l₁) (mergeSize {t₂}{r₁}) ⟩
                suc (heapSize l₂ + heapSize r₂ + heapSize r₁ + heapSize l₁)
            ≡⟨ cong suc (+-assoc (heapSize l₂ + heapSize r₂) (heapSize r₁) (heapSize l₁)) ⟩
                suc (heapSize l₂ + heapSize r₂ + (heapSize r₁ + heapSize l₁))
            ≡⟨ cong (λ e → suc (heapSize l₂ + heapSize r₂ + e)) (+-comm (heapSize r₁) (heapSize l₁)) ⟩
                suc (heapSize l₂ + heapSize r₂ + (heapSize l₁ + heapSize r₁))
            ≡⟨ cong suc (+-comm (heapSize l₂ + heapSize r₂) (heapSize l₁ + heapSize r₁)) ⟩
                suc ((heapSize l₁ + heapSize r₁) + (heapSize l₂ + heapSize r₂))
            ≡⟨ sym (+-suc (heapSize l₁ + heapSize r₁) (heapSize l₂ + heapSize r₂)) ⟩
                heapSize l₁ + heapSize r₁ + suc (heapSize l₂ + heapSize r₂)
            ∎) ⟩
            x₁ ∷ fromHeap (heapSize l₁ + heapSize r₁ + suc (heapSize l₂ + heapSize r₂))
            (mergeHeap (mergeHeap (∧ x₂ l₂ r₂) r₁) l₁)
        ∎) (rec (mergeHeap t₁ t₂) (mergeSize< {x}{t₁}{t₂}) (lem-mergeHeap t₁ t₂ h₁ h₂)))
... | no x₁≰x₂  = s∷ x≤x₂ (subst Sorted (
            fromHeap
            (heapSize (mergeHeap (∧ x₁ l₁ r₁) (∧ x₂ l₂ r₂)))
            (mergeHeap (∧ x₁ l₁ r₁) (∧ x₂ l₂ r₂))
        ≡⟨ cong₂ (λ e₁ e₂ → fromHeap (heapSize e₁) e₂) (step-merge≥ x₁≰x₂) (step-merge≥ x₁≰x₂) ⟩
            x₂ ∷ fromHeap
            (heapSize (mergeHeap t₁ r₂) + heapSize l₂)
            (mergeHeap (mergeHeap t₁ r₂) l₂)
        ≡⟨ cong (λ e → x₂ ∷ fromHeap e (mergeHeap (mergeHeap t₁ r₂) l₂)) (
                heapSize (mergeHeap t₁ r₂) + heapSize l₂
            ≡⟨ cong (λ e → e + heapSize l₂) (mergeSize {t₁}{r₂}) ⟩
                suc (heapSize l₁ + heapSize r₁ + heapSize r₂ + heapSize l₂)
            ≡⟨ cong suc (+-assoc (heapSize l₁ + heapSize r₁) (heapSize r₂) (heapSize l₂)) ⟩
                suc (heapSize l₁ + heapSize r₁ + (heapSize r₂ + heapSize l₂))
            ≡⟨ cong (λ e → suc (heapSize l₁ + heapSize r₁ + e)) (+-comm (heapSize r₂) (heapSize l₂)) ⟩
                suc (heapSize l₁ + heapSize r₁ + (heapSize l₂ + heapSize r₂))
            ≡⟨ sym (+-suc (heapSize l₁ + heapSize r₁) (heapSize l₂ + heapSize r₂)) ⟩
                heapSize l₁ + heapSize r₁ + suc (heapSize l₂ + heapSize r₂)
            ∎) ⟩
            x₂ ∷ fromHeap
            (heapSize l₁ + heapSize r₁ + suc (heapSize l₂ + heapSize r₂))
            (mergeHeap (mergeHeap (∧ x₁ l₁ r₁) r₂) l₂)
        ∎)(rec (mergeHeap t₁ t₂) (mergeSize< {x}{t₁}{t₂}) (lem-mergeHeap t₁ t₂ h₁ h₂)))

lem-fromHeap : (x : Heap ℕ) → IsHeap x → Sorted (fromHeap (heapSize x) x)
lem-fromHeap x h = ⟪-rec P step x h

ここまで示せたら目的である"ヒープソートによってソートされる"ことは簡単に示せます．

HeapSorted.agda

hSorted : (xs : List ℕ) → Sorted (hsort xs)
hSorted xs = lem-fromHeap (makeHeap xs) (lem-makeHeap xs)

感想

明示的に帰納法を回す方法がわからなかったが，今回を通じてなんとなくわかってきた．
記事を書いている途中で気付いたが，fromHeapも帰納法を用いた方がよかった．その内書き直したい．
後編ではヒープソートによって内容が変わらないことを示します．

参考文献

agda-stdlib
お気楽 Haskell プログラミング入門
ヒープソートの実装法を参考にしました．
TPPmark2017のAgda版解説
Agdaで帰納法を使うやり方を参考にしました．