Implement (bounded) sieve of Eratosthenes
This commit is contained in:
commit
35ca270476
|
|
@ -0,0 +1,3 @@
|
|||
name: eratosthenes
|
||||
include: src
|
||||
depend: standard-library-2.1
|
||||
|
|
@ -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
|
||||
|
|
@ -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))
|
||||
|
||||
|
|
@ -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
|
||||
Loading…
Reference in New Issue