1(*
2 * Copyright 2020, Data61, CSIRO (ABN 41 687 119 230)
3 *
4 * SPDX-License-Identifier: BSD-2-Clause
5 *)
6
7theory IsPrime
8imports
9  "AutoCorres.AutoCorres"
10  "HOL-Computational_Algebra.Primes"
11begin
12
13external_file "is_prime.c"
14
15(* Parse the input file. *)
16install_C_file "is_prime.c"
17
18(* Abstract the input file. *)
19autocorres [ ts_rules = nondet, unsigned_word_abs = is_prime_linear is_prime ] "is_prime.c"
20
21context is_prime begin
22
23definition
24  "partial_prime p (n :: nat) \<equiv>
25        (1 < p \<and> (\<forall>i \<in> {2 ..< min p n}. \<not> i dvd p))"
26
27lemma partial_prime_ge [simp]:
28     "\<lbrakk> p' \<ge> p \<rbrakk> \<Longrightarrow> partial_prime p p' = prime p"
29  by (clarsimp simp: partial_prime_def prime_nat_iff' min_def)
30
31lemma divide_self_plus_one [simp]: "(x dvd Suc x) = (x = 1)"
32  apply (case_tac "x = 0", simp)
33  apply (case_tac "x = 1", simp)
34  apply (clarsimp simp: dvd_def)
35  apply (case_tac "k = 0", simp)
36  apply (case_tac "k = 1", simp)
37  apply (subgoal_tac "x * 2 \<le> x * k")
38   apply (subgoal_tac "Suc x < x * 2")
39    apply simp
40   apply clarsimp
41  apply clarsimp
42  done
43
44lemma partial_prime_Suc [simp]:
45  "partial_prime p (Suc n)
46              = (partial_prime p n \<and> (1 < n \<and> Suc n < p \<longrightarrow> \<not> n dvd p))"
47   apply (clarsimp simp: partial_prime_def min_def atLeastLessThanSuc not_le)
48  apply (case_tac "p = Suc n")
49   apply (clarsimp simp: atLeastLessThanSuc)
50  apply fastforce
51  done
52
53lemma mod_to_dvd:
54   "(n mod i \<noteq> 0) = (\<not> i dvd (n :: nat))"
55  by (clarsimp simp: dvd_eq_mod_eq_0)
56
57lemma prime_of_product [simp]: "prime ((a::nat) * b) = ((a = 1 \<and> prime b) \<or> (prime a \<and> b = 1))"
58  by (metis mult.commute mult.right_neutral prime_product)
59
60lemma partial_prime_2 [simp]: "(partial_prime a 2) = (a > 1)"
61  by (clarsimp simp: partial_prime_def)
62
63definition [simp]:
64  "is_prime_linear_inv n i s \<equiv> (1 < i \<and> 1 < n \<and> i \<le> n \<and> partial_prime n i)"
65
66theorem is_prime_correct:
67    "\<lbrace> \<lambda>s. n \<le> UINT_MAX \<rbrace> is_prime_linear' n \<lbrace> \<lambda>r s. (r \<noteq> 0) \<longleftrightarrow> prime n \<rbrace>!"
68  apply (rule validNF_assume_pre)
69  apply (case_tac "n = 0")
70   apply (clarsimp simp: is_prime_linear'_def, wp, simp)[1]
71  apply (case_tac "n = 1")
72   apply (clarsimp simp: is_prime_linear'_def, wp, simp)[1]
73  apply (unfold is_prime_linear'_def)
74  apply (subst whileLoopE_add_inv [
75      where I="\<lambda>r s. is_prime_linear_inv n r s"
76                  and M="(\<lambda>(r, s). n - r)"])
77  apply (wp, auto simp: mod_to_dvd [simplified])
78  done
79
80lemma not_prime:
81    "\<lbrakk> \<not> prime (a :: nat); a > 1 \<rbrakk> \<Longrightarrow> \<exists>x y. x * y = a \<and> 1 < x \<and> 1 < y \<and> x * x \<le> a"
82  apply (clarsimp simp: prime_nat_iff dvd_def)
83  apply (case_tac "m > k")
84   apply (metis Suc_lessD Suc_lessI less_imp_le_nat mult.commute nat_0_less_mult_iff nat_mult_less_cancel_disj)
85  apply fastforce
86  done
87
88lemma sqrt_prime:
89  "\<lbrakk> a * a > n; \<forall>x<a. (x dvd n) = (x = Suc 0 \<or> x = n); 1 < n \<rbrakk> \<Longrightarrow> prime n"
90  apply (rule ccontr)
91  apply (drule not_prime)
92   apply clarsimp
93  apply (metis dvd_triv_right less_le_trans mult.commute mult_le_cancel2
94           One_nat_def less_eq_nat.simps(1) less_not_refl2
95           mult_eq_self_implies_10 not_less)
96  done
97
98lemma partial_prime_sqr:
99     "\<lbrakk> n * n > p \<rbrakk> \<Longrightarrow> partial_prime p n = prime p"
100  apply (case_tac "n \<ge> p")
101   apply clarsimp
102  apply (case_tac "partial_prime p n")
103   apply clarsimp
104   apply (erule sqrt_prime)
105    apply (clarsimp simp: partial_prime_def)
106    apply (case_tac "x = 0", simp)
107    apply (case_tac "x = 1", simp)
108    apply (frule_tac x=x in bspec)
109     apply (clarsimp simp: min_def)
110    apply clarsimp
111   apply (clarsimp simp: not_le partial_prime_def)
112  apply (case_tac "p = 0", simp)
113  apply (case_tac "p = 1", simp)
114  apply (auto simp: not_le partial_prime_def min_def prime_nat_iff')
115  done
116
117definition "SQRT_UINT_MAX \<equiv> 65536 :: nat"
118
119lemma uint_max_factor [simp]:
120  "UINT_MAX = SQRT_UINT_MAX * SQRT_UINT_MAX - 1"
121  by (clarsimp simp: UINT_MAX_def SQRT_UINT_MAX_def)
122
123lemma prime_dvd:
124    "\<lbrakk> prime (p::nat) \<rbrakk> \<Longrightarrow> (r dvd p) = (r = 1 \<or> r = p)"
125  by (fastforce simp: prime_nat_iff)
126
127definition is_prime_inv
128  where [simp]: "is_prime_inv n i s \<equiv> (1 < i \<and> i \<le> n \<and> i \<le> SQRT_UINT_MAX \<and> i * i \<le> SQRT_UINT_MAX * SQRT_UINT_MAX \<and> partial_prime n i)"
129
130lemma nat_leE: "\<lbrakk> (a::nat) \<le> b; a < b \<Longrightarrow> R; a = b \<Longrightarrow> R \<rbrakk> \<Longrightarrow> R"
131  apply atomize_elim
132    apply clarsimp
133  done
134
135lemma sqr_less_mono [simp]:
136    "((i::nat) * i < j * j) = (i < j)"
137  by (metis (full_types) le0 less_not_refl3 linorder_neqE_nat
138        mult_strict_mono' order.strict_trans)
139
140lemma [simp]: "(a - b < a - c) = ((b::nat) > c \<and> c < a)"
141  by arith
142
143declare mult_le_mono [intro]
144
145lemma sqr_le_sqr_minus_1 [simp]:
146    "\<lbrakk> b \<noteq> 0 \<rbrakk> \<Longrightarrow> (a * a \<le> b * b - Suc 0) = (a < b)"
147  by (metis gr0I is_prime.sqr_less_mono nat_0_less_mult_iff nat_le_Suc_less)
148
149theorem is_prime_faster_correct:
150  notes times_nat.simps(2) [simp del] mult_Suc_right [simp del]
151  shows "\<lbrace> \<lambda>s. n \<le> UINT_MAX \<rbrace> is_prime' n \<lbrace> \<lambda>r s. (r \<noteq> 0) \<longleftrightarrow> prime n \<rbrace>!"
152  apply (rule validNF_assume_pre)
153  apply (case_tac "n = 0")
154   apply (clarsimp simp: is_prime'_def, wp, simp)[1]
155  apply (case_tac "n = 1")
156   apply (clarsimp simp: is_prime'_def, wp, simp)[1]
157  apply (unfold is_prime'_def dvd_eq_mod_eq_0 [symmetric] SQRT_UINT_MAX_def [symmetric])
158  apply (subst whileLoopE_add_inv [
159      where I="\<lambda>r s. is_prime_inv n r s"
160      and M="(\<lambda>(r, s). (Suc n) * (Suc n) - r * r)"])
161   apply wp
162    apply clarsimp
163    apply (metis One_nat_def Suc_leI Suc_lessD nat_leE prime_dvd leD mult_le_mono n_less_n_mult_m)
164   apply (fastforce elim: nat_leE simp: partial_prime_sqr)
165  apply (clarsimp simp: SQRT_UINT_MAX_def)
166  done
167
168end
169
170end
171