From 35ca2704765a7276dccc96d722512abccc628408 Mon Sep 17 00:00:00 2001
From: Nathan van Doorn <nvd1234@gmail.com>
Date: Thu, 11 Jul 2024 14:34:19 +0200
Subject: [PATCH] Implement (bounded) sieve of Eratosthenes

---
 eratosthenes.agda-lib |  3 ++
 src/Demo.agda         | 12 ++++++
 src/Eratosthenes.agda | 98 +++++++++++++++++++++++++++++++++++++++++++
 src/SplayHeap.agda    | 90 +++++++++++++++++++++++++++++++++++++++
 4 files changed, 203 insertions(+)
 create mode 100644 eratosthenes.agda-lib
 create mode 100644 src/Demo.agda
 create mode 100644 src/Eratosthenes.agda
 create mode 100644 src/SplayHeap.agda

diff --git a/eratosthenes.agda-lib b/eratosthenes.agda-lib
new file mode 100644
index 0000000..a835a3c
--- /dev/null
+++ b/eratosthenes.agda-lib
@@ -0,0 +1,3 @@
+name: eratosthenes
+include: src
+depend: standard-library-2.1
diff --git a/src/Demo.agda b/src/Demo.agda
new file mode 100644
index 0000000..866207d
--- /dev/null
+++ b/src/Demo.agda
@@ -0,0 +1,12 @@
+{-# OPTIONS --guardedness #-}
+
+module Demo where
+
+open import Data.Nat.Show
+open import Function.Base
+open import IO
+
+open import Eratosthenes
+
+main : Main
+main = run $ putStrLn ∘ show ⟨ List.mapM′ ⟩ primes 1000000
diff --git a/src/Eratosthenes.agda b/src/Eratosthenes.agda
new file mode 100644
index 0000000..0ef1dc3
--- /dev/null
+++ b/src/Eratosthenes.agda
@@ -0,0 +1,98 @@
+------------------------------------------------------------------------
+-- Implementation of Melissa O'Neill's variant of the prime sieve.
+-- I only implement the bounded case and I don't implement wheels yet.
+-- This is also not proven correct, you'll have to take my word for it
+-- that it does what I claim.
+------------------------------------------------------------------------
+{-# OPTIONS --safe --without-K #-}
+
+module Eratosthenes where
+
+open import Data.Nat.Base
+open import Data.Nat.Induction using (<-wellFounded-fast)
+open import Data.Nat.Properties hiding (≤-total; ≤-isTotalOrder; ≤-totalOrder)
+open import Data.List.Base hiding (upTo)
+open import Data.Product.Base
+open import Data.Sum.Base using (inj₁; inj₂)
+open import Function.Base
+open import Induction.WellFounded
+import Relation.Binary.Construct.On as On
+open import Relation.Binary.Bundles using (TotalOrder)
+open import Relation.Binary.Definitions using (Total)
+open import Relation.Binary.Structures using (IsTotalOrder)
+open import Relation.Binary.PropositionalEquality
+open import Relation.Nullary.Construct.Add.Extrema
+open import Relation.Nullary.Decidable
+
+-- Reimplementations of a couple of things from stdlib because
+-- their existing definitions at time of writing are slow
+
+-- ≤-total is currently defined in stdlib using unary arithmetic. This makes it
+-- terrible to use as a conditional. Our heap implementation is generic over a
+-- total order, so we redefine this and the bundle we care ultimately care
+-- about.
+≤-total : Total _≤_
+≤-total m n with m ≤? n
+... | yes m≤n = inj₁ m≤n
+... | no m≰n = inj₂ (≰⇒≥ m≰n)
+
+≤-isTotalOrder : IsTotalOrder _≡_ _≤_
+≤-isTotalOrder = record
+  { isPartialOrder = ≤-isPartialOrder
+  ; total = ≤-total
+  }
+
+≤-totalOrder : TotalOrder _ _ _
+≤-totalOrder = record { isTotalOrder = ≤-isTotalOrder }
+
+-- upTo in stdlib creates larger and larger closures. This causes some slowdown.
+-- We implement it in a tail recursive manner instead.
+upFromThen : ℕ → ℕ → List ℕ
+upFromThen from zero = []
+upFromThen from (suc then) = from ∷ upFromThen (suc from) then
+
+upTo : ℕ → List ℕ
+upTo = upFromThen 0
+
+open import SplayHeap (On.totalOrder ≤-totalOrder (proj₂ {A = ℕ})) public
+
+insertPrime : ℕ → Heap → Heap
+insertPrime p table = insert (p , p * p) table
+
+sieve′ : List ℕ → (t : Heap) → .{{NonEmpty t}} → List ℕ
+sieve′ [] table = []
+sieve′ (x ∷ xs) table with (p , nextComposite) , table′ ← findDeleteMin table with nextComposite ≤? x
+... | no _ = x ∷ sieve′ xs (insertPrime x table)
+... | yes _ = sieve′ xs adjusted
+  where
+    P : Heap → Set
+    P _ = List (ℕ × ℕ) × Heap
+
+    adjust-rec : ∀ (t : Heap) → (∀ {t′ : Heap} → # t′ < # t → P t′) → P t
+    adjust-rec t rec with nonEmpty? t
+    ... | no _ = [] , empty
+    ... | yes ne with (p , c) , t′ ← findDeleteMin t {{ne}} in eq with c ≤? x
+    ... | no _ = [] , t
+    ... | yes _ = map₁ ((p , c + p) ∷_) ∘ rec $ begin-strict
+      # t′                            ≡˘⟨ cong (#_ ∘ proj₂) eq ⟩
+      # proj₂ (findDeleteMin t)       <⟨ n<1+n (# proj₂ (findDeleteMin t)) ⟩
+      suc (# proj₂ (findDeleteMin t)) ≡⟨ #-findDeleteMin t ⟩
+      # t                             ∎
+      where
+        -- I wish I could have instance patterns!
+        instance _ = ne
+        open ≤-Reasoning
+
+    adjust : Heap → List (ℕ × ℕ) × Heap
+    adjust = All.wfRec (On.wellFounded #_ <-wellFounded-fast) _ P adjust-rec
+    
+    adjusted : Heap
+    adjusted with next , t′ ← adjust table′ = insert (p , nextComposite + p) (foldr insert t′ next)
+
+sieve : List ℕ → List ℕ
+sieve [] = []
+sieve (x ∷ xs) = x ∷ sieve′ xs (insertPrime x empty)
+
+primes : ℕ → List ℕ
+primes n = sieve (drop 2 (upTo n))
+
diff --git a/src/SplayHeap.agda b/src/SplayHeap.agda
new file mode 100644
index 0000000..3536350
--- /dev/null
+++ b/src/SplayHeap.agda
@@ -0,0 +1,90 @@
+------------------------------------------------------------------------
+-- Implementation of a splay heap.
+--
+-- Based on Okasaki's Purely Functional Data Structures and McBride's
+-- How to Keep Your Neighbours in Order.
+------------------------------------------------------------------------
+
+{-# OPTIONS --safe --without-K #-}
+
+open import Relation.Binary.Bundles
+
+module SplayHeap {c ℓ₁ ℓ₂} (T : TotalOrder c ℓ₁ ℓ₂) where
+
+open TotalOrder T hiding (refl)
+
+open import Data.Nat.Base using (ℕ; zero; suc; _+_)
+open import Data.Nat.Properties using (n<1+n; +-suc; +-assoc)
+open import Data.Product.Base
+open import Data.Sum.Base using (inj₁; inj₂)
+open import Function.Base
+open import Level using (_⊔_)
+open import Relation.Binary.Construct.Add.Extrema.NonStrict _≤_
+open import Relation.Binary.PropositionalEquality using (_≡_; refl; cong; module ≡-Reasoning)
+open import Relation.Nullary.Construct.Add.Extrema
+open import Relation.Nullary.Decidable.Core
+
+data Tree (l u : Carrier ±) : Set (c ⊔ ℓ₂) where
+  leaf : .(l ≤± u) → Tree l u
+  node : (x : Carrier) → Tree l [ x ] → Tree [ x ] u → Tree l u
+
+Heap : Set (c ⊔ ℓ₂)
+Heap = Tree ⊥± ⊤±
+
+data NonEmpty {l u} : Tree l u → Set (c ⊔ ℓ₂) where
+  instance nonEmpty : ∀ {x l r} → NonEmpty (node x l r)
+
+nonEmpty? : ∀ {l u} (t : Tree l u) → Dec (NonEmpty t)
+nonEmpty? (leaf _) = no λ ()
+nonEmpty? (node _ _ _) = yes nonEmpty
+
+private
+  -- linear in length of spine :(
+  -- It seems to compile away fine, though.
+  slacken⊥ : ∀ {l l′ u} → Tree l u → .(l′ ≤± l) → Tree l′ u
+  slacken⊥ (leaf l≤u) l′≤l = leaf (≤±-trans trans l′≤l l≤u)
+  slacken⊥ (node x a b) l≤l′ = node x (slacken⊥ a l≤l′) b
+
+  partition : ∀ {l u} (p : Carrier) → .(l ≤± [ p ]) → .([ p ] ≤± u) → Tree l u → Tree l [ p ] × Tree [ p ] u
+  partition p l≤p p≤u (leaf _) = leaf l≤p , leaf p≤u
+  partition p l≤p p≤u (node x a b) with total x p
+  partition p l≤p p≤u (node x a (leaf _)) | inj₁ x≤p = node x a (leaf [ x≤p ]) , leaf p≤u
+  partition p l≤p p≤u (node x a (node y b c)) | inj₁ x≤p with total y p
+  ... | inj₁ y≤p with small , big ← partition p [ y≤p ] p≤u c = node y (node x a b) small , big
+  ... | inj₂ p≤y with small , big ← partition p [ x≤p ] [ p≤y ] b = node x a small , node y big c
+  partition p l≤p p≤u (node x (leaf _) a) | inj₂ p≤x = leaf l≤p , node x (leaf [ p≤x ]) a
+  partition p l≤p p≤u (node x (node y a b) c) | inj₂ p≤x with total y p
+  ... | inj₁ y≤p with small , big ← partition p [ y≤p ] [ p≤x ] b = node y a small , node x big c
+  ... | inj₂ p≤y with small , big ← partition p l≤p [ p≤y ] a = small , (node y big (node x b c))
+
+insert : Carrier → Heap → Heap
+insert x t with l , r ← partition x ⊥±≤[ x ] [ x ]≤⊤± t = node x l r
+
+findDeleteMin : ∀ {l u} (t : Tree l u) → .{{NonEmpty t}} → Carrier × Tree l u
+findDeleteMin (node x (leaf l≤x) a) = x , slacken⊥ a l≤x
+findDeleteMin (node x (node y (leaf l≤y) a) b) = y , node x (slacken⊥ a l≤y) b
+findDeleteMin (node x (node y a@(node _ _ _) b) c) with z , a′ ← findDeleteMin a = z , node y a′ (node x b c)
+
+empty : Heap
+empty = leaf ⊥±≤⊤±
+
+-- number of elements in the heap
+-- useful for induction proofs hopefully
+#_ : ∀ {l u} → Tree l u → ℕ
+# leaf _ = 0
+# node x a b = suc (# a + # b)
+
+private
+  #-slacken⊥ : ∀ {l l′ u} (t : Tree l u) .(l′≤l : l′ ≤± l) → # slacken⊥ t l′≤l ≡ # t
+  #-slacken⊥ (leaf x) l′≤l = refl
+  #-slacken⊥ (node x a b) l′≤l = cong (suc ∘ (_+ # b)) (#-slacken⊥ a l′≤l)
+
+#-findDeleteMin : ∀ {l u} (t : Tree l u) → .{{_ : NonEmpty t}} → suc (# proj₂ (findDeleteMin t)) ≡ # t
+#-findDeleteMin (node x (leaf l≤x) b) = cong suc (#-slacken⊥ b l≤x)
+#-findDeleteMin (node x (node y (leaf l≤y) b) c) = cong (suc ∘ suc ∘ (_+ # c)) (#-slacken⊥ b l≤y)
+#-findDeleteMin (node x (node y a@(node _ s t) b) c) = cong (suc ∘ suc) $ begin
+  # proj₂ (findDeleteMin a) + suc (# b + # c)   ≡⟨ +-suc (# proj₂ (findDeleteMin a)) (# b + # c) ⟩
+  suc (# proj₂ (findDeleteMin a)) + (# b + # c) ≡⟨ cong (_+ (# b + # c)) (#-findDeleteMin a) ⟩
+  suc (# s + # t) + (# b + # c)                 ≡˘⟨ cong suc (+-assoc (# s + # t) (# b) (# c)) ⟩
+  suc (# s + # t + # b + # c)                   ∎
+  where open ≡-Reasoning