|
| 1 | +---------------------------------------------------- |
| 2 | +-- Example 3: Multi-Modal Proof |
| 3 | +---------------------------------------------------- |
| 4 | + |
| 5 | +-- Spec: This program checks whether a given natural number is non-prime. |
| 6 | +-- |
| 7 | +-- Comment (not part of the spec): The time complexity of this algorithm is |
| 8 | +-- O(√n): we check for divisors up to the square root of n. The program quite |
| 9 | +-- simple, but to verify its correctness, we need to prove a property: |
| 10 | +-- |
| 11 | +-- n is prime ↔ n > 1 ∧ ∀ d, 2 ≤ d ≤ √n → n % d ≠ 0 |
| 12 | +-- |
| 13 | +-- It is obvious for people who know number theory, but is painful to formalise |
| 14 | +-- it in a proof assistant. |
| 15 | + |
| 16 | +import Mathlib.Tactic |
| 17 | + |
| 18 | +import Loom.MonadAlgebras.NonDetT.Extract |
| 19 | +import Loom.MonadAlgebras.WP.Tactic |
| 20 | +import Loom.MonadAlgebras.WP.DoNames' |
| 21 | + |
| 22 | +import CaseStudies.Velvet.Std |
| 23 | +import CaseStudies.TestingUtil |
| 24 | + |
| 25 | +open PartialCorrectness DemonicChoice Lean.Elab.Term.DoNames |
| 26 | + |
| 27 | +set_option loom.semantics.termination "partial" |
| 28 | +set_option loom.semantics.choice "demonic" |
| 29 | + |
| 30 | +--------------------------------------------------------- |
| 31 | +-- The program and the spec |
| 32 | +--------------------------------------------------------- |
| 33 | + |
| 34 | +-- Helper definition to count divisors of a natural number |
| 35 | +-- This counts all positive divisors from 1 to n |
| 36 | +def countDivisors (n: Nat) : Nat := |
| 37 | + (List.range (n + 1)).filter (fun d => d > 0 ∧ n % d = 0) |>.length |
| 38 | + |
| 39 | +-- Helper definition for prime numbers |
| 40 | +-- A number is prime if and only if: |
| 41 | +-- 1. It is greater than 1 |
| 42 | +-- 2. It has exactly 2 positive divisors (1 and itself) |
| 43 | +def isPrime (n: Nat) : Prop := |
| 44 | + n > 1 ∧ countDivisors n = 2 |
| 45 | + |
| 46 | +method IsNonPrime (n: Nat) |
| 47 | + return (result: Bool) |
| 48 | + ensures result ↔ ¬isPrime n |
| 49 | + do |
| 50 | + if n ≤ 1 then |
| 51 | + return true |
| 52 | + let mut i: Nat := 2 |
| 53 | + let mut ret: Bool := false |
| 54 | + while i * i ≤ n |
| 55 | + invariant 2 ≤ i |
| 56 | + invariant (ret = false ↔ ∀ d, 2 ≤ d ∧ d < i → n % d ≠ 0) |
| 57 | + invariant (i - 1) * (i - 1) ≤ n |
| 58 | + do |
| 59 | + if n % i = 0 then |
| 60 | + ret := true |
| 61 | + i := i + 1 |
| 62 | + return ret |
| 63 | + |
| 64 | +------------------------------------------------ |
| 65 | +-- Program verification |
| 66 | +------------------------------------------------ |
| 67 | + |
| 68 | +theorem goal2 |
| 69 | +(n : ℕ) |
| 70 | +(i : ℕ) |
| 71 | +(ret : Bool) |
| 72 | +(i_1 : ℕ) |
| 73 | +(ret_1 : Bool) |
| 74 | +(if_neg : 1 < n) |
| 75 | +(invariant_1 : 2 ≤ i_1) |
| 76 | +(invariant_2 : ret_1 = false ↔ ∀ (d : ℕ), 2 ≤ d → d < i_1 → ¬n % d = 0) |
| 77 | +(invariant_3 : (i_1 - 1) * (i_1 - 1) ≤ n) |
| 78 | +(done_1 : n < i_1 * i_1) |
| 79 | +(i_2 : i = i_1 ∧ ret = ret_1) |
| 80 | +: ret_1 = true ↔ ¬isPrime n := |
| 81 | + by |
| 82 | + -- If ret_1 is true, then there exists a divisor d between 2 and i_1-1, making n not prime. |
| 83 | + have h_true : ret_1 = Bool.true → ¬isPrime n := by |
| 84 | + aesop; |
| 85 | + unfold isPrime at a; aesop; |
| 86 | + unfold countDivisors at right; contrapose! right; aesop; |
| 87 | + -- Since $w$ divides $n$, and $1$ and $n$ are also divisors of $n$, the list of divisors must contain at least these three elements. |
| 88 | + have h_divisors : 1 ∈ List.filter (fun d => 0 < d ∧ n % d = 0) (List.range (n + 1)) ∧ w ∈ List.filter (fun d => 0 < d ∧ n % d = 0) (List.range (n + 1)) ∧ n ∈ List.filter (fun d => 0 < d ∧ n % d = 0) (List.range (n + 1)) := by |
| 89 | + simp_all +decide [ List.mem_filter, List.mem_range ]; |
| 90 | + exact ⟨ ⟨ by linarith, Nat.mod_one _ ⟩, ⟨ by nlinarith only [ left, left_1, left_2, invariant_3, Nat.sub_add_cancel ( by linarith : 1 ≤ i ) ], by linarith ⟩, by linarith ⟩; |
| 91 | + have h_divisors_card : List.toFinset (List.filter (fun d => 0 < d ∧ n % d = 0) (List.range (n + 1))) ⊇ {1, w, n} := by |
| 92 | + aesop_cat; |
| 93 | + have h_divisors_card : Finset.card (List.toFinset (List.filter (fun d => 0 < d ∧ n % d = 0) (List.range (n + 1)))) ≥ 3 := by |
| 94 | + refine' le_trans _ ( Finset.card_mono h_divisors_card ); |
| 95 | + rw [ Finset.card_insert_of_notMem, Finset.card_insert_of_notMem ] <;> norm_num; |
| 96 | + · nlinarith [ Nat.sub_add_cancel ( by linarith : 1 ≤ i ) ]; |
| 97 | + · grind; |
| 98 | + exact h_divisors_card.not_lt ( lt_of_le_of_lt ( List.toFinset_card_le _ ) ( by aesop ) ); |
| 99 | + -- If ret_1 is false, then by invariant_2, there are no divisors of n in the range 2 to i_1-1. Since n is less than i_1 squared, and i_1 is at least 2, this implies that n has no divisors other than 1 and itself. Hence, n is prime. |
| 100 | + have h_false : ret_1 = Bool.false → isPrime n := by |
| 101 | + intro h |
| 102 | + have h_no_divisors : ∀ d, 1 < d → d < n → ¬n % d = 0 := by |
| 103 | + intros d hd1 hd2 hd3; |
| 104 | + have h_div : d ≤ i_1 - 1 ∨ n / d ≤ i_1 - 1 := by |
| 105 | + exact Classical.or_iff_not_imp_left.2 fun h => by nlinarith [ Nat.div_mul_cancel ( Nat.dvd_of_mod_eq_zero hd3 ), Nat.sub_add_cancel ( by linarith : 1 ≤ i_1 ) ] ; |
| 106 | + bound; |
| 107 | + · exact invariant_2.mp rfl d hd1 ( by linarith [ Nat.sub_add_cancel ( by linarith : 1 ≤ i ) ] ) hd3; |
| 108 | + · exact invariant_2.mp rfl ( n / d ) ( by nlinarith [ Nat.div_mul_cancel ( Nat.dvd_of_mod_eq_zero hd3 ) ] ) ( by omega ) ( Nat.mod_eq_zero_of_dvd <| Nat.div_dvd_of_dvd <| Nat.dvd_of_mod_eq_zero hd3 ) |
| 109 | + have h_prime : Nat.Prime n := by |
| 110 | + exact Nat.prime_def_lt'.mpr ⟨ if_neg, fun d hd₁ hd₂ hd₃ => h_no_divisors d hd₁ hd₂ <| Nat.mod_eq_zero_of_dvd hd₃ ⟩ |
| 111 | + exact (by |
| 112 | + constructor <;> aesop; |
| 113 | + -- Since n is prime, its only divisors are 1 and itself. |
| 114 | + have h_divisors : List.toFinset (List.filter (fun d => d > 0 ∧ n % d = 0) (List.range (n + 1))) = {1, n} := by |
| 115 | + ext d; aesop; |
| 116 | + · exact Classical.or_iff_not_imp_left.2 fun h => by have := Nat.dvd_of_mod_eq_zero right; rw [ Nat.dvd_prime h_prime ] at this; aesop; |
| 117 | + · linarith; |
| 118 | + · rw [ Nat.mod_one ]; |
| 119 | + · grind; |
| 120 | + -- Since the Finset {1, n} has cardinality 2, we can conclude that the list's length is also 2. |
| 121 | + have h_card : List.length (List.filter (fun d => d > 0 ∧ n % d = 0) (List.range (n + 1))) = Finset.card (List.toFinset (List.filter (fun d => d > 0 ∧ n % d = 0) (List.range (n + 1)))) := by |
| 122 | + rw [ List.toFinset_card_of_nodup ]; |
| 123 | + refine' List.Nodup.filter _ _; |
| 124 | + grind; |
| 125 | + exact h_card.trans ( h_divisors.symm ▸ by rw [ Finset.card_insert_of_notMem, Finset.card_singleton ] ; aesop )); |
| 126 | + grind |
| 127 | + |
| 128 | +---------------------------------------------------------------- |
| 129 | +-- Putting it all together |
| 130 | +---------------------------------------------------------------- |
| 131 | + |
| 132 | +prove_correct IsNonPrime by |
| 133 | + loom_solve <;> try simp_all |
| 134 | + · { rw [isPrime]; aesop } |
| 135 | + · apply (goal2 n i ret i_1 ret_1 if_neg invariant_1 invariant_2 invariant_3 done_1 i_2) |
0 commit comments