## Prime Sieves in Agda

Prime numbers in Agda are *slow*. First, they’re Peano-based, so a huge chunk of optimizations we might make in other languages are out of the window. Second, we really often want to *prove* that they’re prime, so the generation code has to carry verification logic with it (I won’t do that today, though). And third, as always in Agda, you have to convince the compiler of termination. With all of that in mind, let’s try and write a (very slow, very basic) prime sieve in Agda.

First, we can make an “array” of numbers that we cross off as we go.

```
: ∀ n → List (Fin n)
primes = []
primes zero (suc zero) = []
primes (suc (suc zero)) = []
primes (suc (suc (suc m))) = sieve (tabulate (just ∘ Fin.suc))
primes where
: Fin _ → List (Maybe (Fin _)) → List (Maybe (Fin _))
cross-off
: List (Maybe (Fin _)) → List (Fin _)
sieve = []
sieve [] (nothing ∷ xs) = sieve xs
sieve (just x ∷ xs) = suc x ∷ sieve (cross-off x xs)
sieve
= foldr f (const []) fs p
cross-off p fs where
= ∀ {i} → Fin i → List (Maybe (Fin (2 + m)))
B
: Maybe (Fin (2 + m)) → B → B
f _ xs zero = nothing ∷ xs p
f (suc y) = x ∷ xs y f x xs
```

Very simple so far: we run through the list, filtering out the multiples of each prime as we see it. Unfortunately, this won’t pass the termination checker. This recursive call to `sieve`

is the problem:

`(just x ∷ xs) = suc x ∷ sieve (cross-off x xs) sieve `

Agda finds if a function is terminating by checking that at least one argument gets (structurally) smaller on every recursive call. `sieve`

only takes one argument (the input list), so that’s the one that needs to get smaller. In the line above, if we replaced it with the following:

`(just x ∷ xs) = suc x ∷ sieve xs sieve `

We’d be good to go: `xs`

is definitely smaller than `(just x ∷ xs)`

. `cross-off x xs`

, though? The thing is, `cross-off`

returns a list of the same length that it’s given. But the function call is opaque: Agda can’t automatically see the fact that the length stays the same. Reaching for a proof here is the wrong move, though: you can get all of the same benefit by switching out the list for a length-indexed vector.

```
: ∀ n → List (Fin n)
primes = []
primes zero (suc zero) = []
primes (suc (suc zero)) = []
primes (suc (suc (suc m))) = sieve (tabulate (just ∘ Fin.suc))
primes where
: ∀ {n} → Fin _ → Vec (Maybe _) n → Vec (Maybe _) n
cross-off
: ∀ {n} → Vec (Maybe (Fin (2 + m))) n → List (Fin (3 + m))
sieve = []
sieve [] (nothing ∷ xs) = sieve xs
sieve (just x ∷ xs) = suc x ∷ sieve (cross-off x xs)
sieve
= foldr B f (const []) fs p
cross-off p fs where
= λ n → ∀ {i} → Fin i → Vec (Maybe (Fin (2 + m))) n
B
: ∀ {n} → Maybe (Fin (2 + m)) → B n → B (suc n)
f _ xs zero = nothing ∷ xs p
f (suc y) = x ∷ xs y f x xs
```

Actually, my explanation above is a little bit of a lie. Often, the way I think about dependently-typed programs has a lot to do with my intuition for “proofs” and so on. But this leads you down the wrong path (and it’s why writing a proof that `cross-off`

returns a list of the same length is the wrong move).

The actual termination checking algorithm is very simple, albeit strict: the argument passed recursively must be *structurally* smaller. That’s it. Basically, the recursive argument has to be contained in one of the arguments passed. It has nothing to do with Agda “seeing” inside the function `cross-off`

or anything like that. What we’ve done above (to make it terminate) is add another argument to the function: the length of the vector. The argument is implicit, but if we were to make it explicit in the recursive call:

`{suc n} (just x ∷ xs) = suc x ∷ sieve {n} (cross-off x xs) sieve `

We can see that it does indeed get structurally smaller.

# Adding the Squaring Optimization

A simple improvement we should be able to make is stopping once we hit the square root of the limit. Since we don’t want to be squaring as we go, we’ll use the following identity:

$(n + 1)^2 = n^2 + 2n + 1$

to figure out the square of the next number from the previous. In fact, we’ll just pass in the limit, and reduce it by $2n + 1$ each time, until it reaches zero:

```
: ∀ n → List (Fin n)
primes = []
primes zero (suc zero) = []
primes (suc (suc zero)) = []
primes (suc (suc (suc m))) = sieve 1 m (Vec.tabulate (just ∘ Fin.suc ∘ Fin.suc))
primes where
: ∀ {n} → ℕ → Vec (Maybe _) n → Vec (Maybe _) n
cross-off
: ∀ {n} → ℕ → ℕ → Vec (Maybe (Fin (3 + m))) n → List (Fin (3 + m))
sieve _ zero = List.mapMaybe id ∘ Vec.toList
sieve _ (suc _) [] = []
sieve (suc l) (nothing ∷ xs) = sieve (suc i) (l ∸ i ∸ i) xs
sieve i (suc l) (just x ∷ xs) = x ∷ sieve (suc i) (l ∸ i ∸ i) (cross-off i xs)
sieve i
= Vec.foldr B f (const []) fs p
cross-off p fs where
= λ n → ℕ → Vec (Maybe (Fin (3 + m))) n
B
: ∀ {i} → Maybe (Fin (3 + m)) → B i → B (suc i)
f _ xs zero = nothing ∷ xs p
f (suc y) = x ∷ xs y f x xs
```

# Finding Prime Factors

A slight variation on the code above (the first version) will give us the prime factors of a number:

```
: ∀ n → List (Fin n)
primeFactors = []
primeFactors zero (suc zero) = []
primeFactors (suc (suc zero)) = []
primeFactors (suc (suc (suc m))) = sieve (Vec.tabulate (just ∘ Fin.suc))
primeFactors where
: ∀ {n} → Vec (Maybe (Fin (2 + m))) n → List (Fin (3 + m))
sieve = []
sieve [] (nothing ∷ xs) = sieve xs
sieve (just x ∷ xs) = Vec.foldr B remove b xs sieve x
sieve where
= λ n → ∀ {i}
B → (Vec (Maybe (Fin (2 + m))) n
→ List (Fin (3 + m)))
→ Fin i
→ List (Fin (3 + m))
: B 0
b = suc x ∷ k []
b k zero (suc _) = k []
b k
: ∀ {n} → Maybe (Fin (2 + m)) → B n → B (suc n)
remove = ys (k ∘ (nothing ∷_)) x
remove y ys k zero (suc j) = ys (k ∘ (y ∷_)) j remove y ys k
```

Adding the squaring optimization complicates things significantly:

```
: ∀ n → List (Fin n)
primeFactors = []
primeFactors zero (suc zero) = []
primeFactors (suc (suc zero)) = []
primeFactors (suc (suc (suc m))) = sqr (suc m) m suc sieve
primeFactors where
_2F-_ : ∀ {n} → ℕ → Fin n → ℕ
= x
x 2F- zero = zero
zero 2F- suc y = zero
suc zero 2F- suc y (suc x) 2F- suc y = x 2F- y
suc
: ∀ n
sqr → ℕ
→ (Fin n → Fin (2 + m))
→ (∀ {i} → Vec (Maybe (Fin (2 + m))) i → ℕ → List (Fin (3 + m)))
→ List (Fin (3 + m))
= k [] n
sqr n zero f k (suc l) f k = k [] zero
sqr zero (suc n) (suc l) f k =
sqr let x = f zero
in sqr n (l 2F- x) (f ∘ suc) (k ∘ (just x ∷_))
: ∀ {n} → Vec (Maybe (Fin (2 + m))) n → ℕ → List (Fin (3 + m))
sieve = go xs′
sieve xs′ i where
: ∀ {n} → Vec (Maybe (Fin (2 + m))) n → List (Fin (3 + m))
go = []
go [] (nothing ∷ xs) = go xs
go (just x ∷ xs) = Vec.foldr B remove (b i) xs x go
go where
= λ n → ∀ {i}
B → Fin i
→ (Vec (Maybe (Fin (2 + m))) n → List (Fin (3 + m)))
→ List (Fin (3 + m))
: ℕ → B 0
b = suc x ∷ k []
b zero zero k (suc y) k = k []
b zero (suc n) zero k = b n x k
b (suc n) (suc y) k = b n y k
b
: ∀ {n} → Maybe (Fin (2 + m)) → B n → B (suc n)
remove = ys x (k ∘ (nothing ∷_))
remove y ys zero k (suc j) k = ys j (k ∘ (y ∷_)) remove y ys
```

# Infinitude

The above sieve aren’t “true” in that each `remove`

is linear, so the performance is $\mathcal{O}(n^2)$ overall. This is the same problem we ran into with the naive infinite sieve in Haskell.

Since it bears such a similarity to the infinite sieve, we have to ask: can *this* sieve be infinite? Agda supports a notion of infinite data, so it would seem like it:

```
infixr 5 _◂_
record Stream (A : Set) : Set where
constructor _◂_
coinductive
field
: A
head : Stream A
tail open Stream
: Stream ℕ
primes = sieve 1 nats
primes where
: Stream ℕ
nats = 0
head nats = nats
tail nats
: ℕ → Stream ℕ → Stream ℕ
sieve (sieve i xs) = suc i
head (sieve i xs) = remove i (head xs) (tail xs) (sieve ∘ suc ∘ (_+ i))
tail where
: ℕ → ℕ → Stream ℕ → (ℕ → Stream ℕ → Stream ℕ) → Stream ℕ
remove = remove i (head zs) (tail zs) (k ∘ suc)
remove zero zero zs k (suc z) zs k = remove i z zs (k ∘ suc)
remove zero (suc y) zero zs k = k zero (remove y (head zs) (tail zs) _◂_)
remove (suc y) (suc z) zs k = remove y z zs (k ∘ suc) remove
```

But this won’t pass the termination checker. What we actually need to prove to do so is that there are infinitely many primes: a nontrivial task in Agda.