# Theory Pratt_Certificate

theory Pratt_Certificate
imports Lehmer Code_Target_Numeral
theory Pratt_Certificate
imports
Complex_Main
"../Lehmer/Lehmer"
"~~/src/HOL/Library/Code_Target_Numeral"
begin

(* TODO: Move? *)
function mod_exp :: "nat ⇒ nat ⇒ nat ⇒ nat" where
"mod_exp b 0 m = 1 mod m"
| "e > 0 ⟹ even e ⟹ mod_exp b e m = mod_exp ((b^2) mod m) (e div 2) m"
| "odd e ⟹ mod_exp b e m = (b * mod_exp ((b^2) mod m) (e div 2) m) mod m"
by fast simp_all
termination
by (relation "measure (λ(_,e,_) ⇒ e)") (auto intro: div_less_dividend odd_pos)

lemma mod_exp_correct: "mod_exp b e m = (b ^ e) mod m"
proof (induction b e m rule: mod_exp.induct, goal_cases)
case (2 e b m)
from 2 have "mod_exp b e m = b ^ (2 * (e div 2)) mod m"
by (subst power_mult) (simp add: power_mod)
also from 2 have "2 * (e div 2) = e" by simp
finally show ?case .
next
case (3 e b m)
from 3 have "mod_exp b e m = b ^ (Suc (2 * (e div 2))) mod m"
by (simp only: power_mult power_Suc) (simp add: power_mod mod_mult_right_eq [symmetric])
also from 3 have "Suc (2 * (e div 2)) = e" by simp
finally show ?case .
qed simp

declare mod_exp.simps [simp del]

lemma eval_mod_exp [simp]:
"mod_exp b 0 (Suc 0) = 0"
"m ≠ 1 ⟹ mod_exp b 0 m = 1"
"mod_exp b 1 m = b mod m"
"mod_exp b (Suc 0) m = b mod m"
"mod_exp b (numeral (num.Bit0 n)) m = mod_exp (b⇧2 mod m) (numeral n) m"
"mod_exp b (numeral (num.Bit1 n)) m = b * mod_exp (b⇧2 mod m) (numeral n) m mod m"
proof -
have "1 mod m = 1" if "m ≠ 1" using that
by (cases m) auto
show "m ≠ 1 ⟹ mod_exp b 0 m = 1" "mod_exp b 0 (Suc 0) = 0"
by (cases m) (simp_all add: mod_exp.simps)
show "mod_exp b 1 m = b mod m" by (simp add: mod_exp_correct)
show "mod_exp b (Suc 0) m = b mod m" by (simp add: mod_exp_correct)
have "numeral (num.Bit0 n) = (2 * numeral n :: nat)"
by (subst numeral.numeral_Bit0) (simp del: arith_simps)
also have "mod_exp b … m = mod_exp (b⇧2 mod m) (numeral n) m"
by (subst mod_exp.simps(2)) (simp_all del: arith_simps)
finally show "mod_exp b (numeral (num.Bit0 n)) m = mod_exp (b⇧2 mod m) (numeral n) m" by simp
have "numeral (num.Bit1 n) = Suc (2 * numeral n :: nat)"
by (subst numeral.numeral_Bit1) (simp del: arith_simps)
also have "mod_exp b … m = b * mod_exp (b⇧2 mod m) (numeral n) m mod m"
by (subst mod_exp.simps(3)) (simp_all del: arith_simps)
finally show "mod_exp b (numeral (num.Bit1 n)) m = b * mod_exp (b⇧2 mod m) (numeral n) m mod m" .
qed

(*
Need to use @{term "b*b::nat"} here instead of squaring, otherwise code_unfold produces
infinite recursion
FIXME: Better performance using divMod or bit tests/bit shifts. Tail-recursive would also be good.
*)
lemma mod_exp_code [code]:
"mod_exp b e m =
(if e = 0 then if m = 1 then 0 else 1
else if even e then mod_exp ((b*b) mod m) (e div 2) m
else (b * mod_exp ((b*b) mod m) (e div 2) m) mod m)"
by (simp_all add: mod_exp.simps(2,3) power2_eq_square)

lemmas mod_exp_code_unfold [code_unfold] = mod_exp_correct [symmetric]

(* END TODO *)

section {* Pratt's Primality Certificates *}
text_raw {* \label{sec:pratt} *}

text {*
This work formalizes Pratt's proof system as described in his article
Every Prime has a Succinct Certificate''\cite{pratt1975certificate}.
The proof system makes use of two types of predicates:
\begin{itemize}
\item $\text{Prime}(p)$: $p$ is a prime number
\item $(p, a, x)$: @{text "∀q ∈ prime_factors(x). [a^((p - 1) div q) ≠ 1] (mod p)"}
\end{itemize}
We represent these predicates with the following datatype:
*}

datatype pratt = Prime nat | Triple nat nat nat

text {*
Pratt describes an inference system consisting of the axiom $(p, a, 1)$
and the following inference rules:
\begin{itemize}
\item R1: If we know that $(p, a, x)$ and @{text "[a^((p - 1) div q) ≠ 1] (mod p)"} hold for some
prime number $q$ we can conclude $(p, a, qx)$ from that.
\item R2: If we know that $(p, a, p - 1)$ and  @{text "[a^(p - 1) = 1] (mod p)"} hold, we can
infer $\text{Prime}(p)$.
\end{itemize}
Both rules follow from Lehmer's theorem as we will show later on.

A list of predicates (i.e., values of type @{type pratt}) is a \emph{certificate}, if it is
built according to the inference system described above. I.e., a list @{term "x # xs :: pratt list"}
is a certificate if @{term "xs :: pratt list"} is a certificate and @{term "x :: pratt"} is
either an axiom or all preconditions of @{term "x :: pratt"} occur in @{term "xs :: pratt list"}.

We call a certificate @{term "xs :: pratt list"} a \emph{certificate for @{term p}},
if @{term "Prime p"} occurs in @{term "xs :: pratt list"}.

The function @{text valid_cert} checks whether a list is a certificate.
*}

fun valid_cert :: "pratt list ⇒ bool" where
"valid_cert [] = True"
| R2: "valid_cert (Prime p#xs) ⟷ 1 < p ∧ valid_cert xs
∧ (∃ a . [a^(p - 1) = 1] (mod p) ∧ Triple p a (p - 1) ∈ set xs)"
| R1: "valid_cert (Triple p a x # xs) ⟷ p > 1 ∧ 0 < x  ∧ valid_cert xs ∧ (x=1 ∨
(∃q y. x = q * y ∧ Prime q ∈ set xs ∧ Triple p a y ∈ set xs
∧ [a^((p - 1) div q) ≠ 1] (mod p)))"

text {*
We define a function @{term size_cert} to measure the size of a certificate, assuming
a binary encoding of numbers. We will use this to show that there is a certificate for a
prime number $p$ such that the size of the certificate is polynomially bounded in the size
of the binary representation of $p$.
*}
fun size_pratt :: "pratt ⇒ real" where
"size_pratt (Prime p) = log 2 p" |
"size_pratt (Triple p a x) = log 2 p + log 2 a + log 2 x"

fun size_cert :: "pratt list ⇒ real" where
"size_cert [] = 0" |
"size_cert (x # xs) = 1 + size_pratt x + size_cert xs"

subsection {* Soundness *}

text {*
In Section \ref{sec:pratt} we introduced the predicates $\text{Prime}(p)$ and $(p, a, x)$.
In this section we show that for a certificate every predicate occuring in this certificate
holds. In particular, if $\text{Prime}(p)$ occurs in a certificate, $p$ is prime.
*}

lemma prime_factors_one [simp]: shows "prime_factors (Suc 0) = {}"
using prime_factorization_1 [where ?'a = nat] by simp

lemma prime_factors_of_prime: fixes p :: nat assumes "prime p" shows "prime_factors p = {p}"
using assms by (fact prime_prime_factors)

theorem pratt_sound:
assumes 1: "valid_cert c"
assumes 2: "t ∈ set c"
shows "(t = Prime p ⟶ prime p) ∧
(t = Triple p a x ⟶ ((∀ q ∈ prime_factors x . [a^((p - 1) div q) ≠ 1] (mod p)) ∧ 0<x))"
using assms
proof (induction c arbitrary: p a x t)
case Nil then show ?case by force
next
case (Cons y ys)
{ assume "y=Triple p a x" "x=1"
then have "(∀ q ∈ prime_factors x . [a^((p - 1) div q) ≠ 1] (mod p)) ∧ 0<x" by simp
}
moreover
{ assume x_y: "y=Triple p a x" "x~=1"
hence "x>0" using Cons.prems by auto
obtain q z where "x=q*z" "Prime q ∈ set ys ∧ Triple p a z ∈ set ys"
and cong:"[a^((p - 1) div q) ≠ 1] (mod p)" using Cons.prems x_y by auto
then have factors_IH:"(∀ r ∈ prime_factors z . [a^((p - 1) div r) ≠ 1] (mod p))" "prime q" "z>0"
using Cons.IH Cons.prems x>0 y=Triple p a x
by force+
then have "prime_factors x = prime_factors z ∪ {q}"  using x =q*z x>0
by (simp add: prime_factors_product prime_factors_of_prime)
then have "(∀ q ∈ prime_factors x . [a^((p - 1) div q) ≠ 1] (mod p)) ∧ 0 < x"
using factors_IH cong by (simp add: x>0)
}
ultimately have y_Triple:"y=Triple p a x ⟹ (∀ q ∈ prime_factors x .
[a^((p - 1) div q) ≠ 1] (mod p)) ∧ 0<x" by linarith
{ assume y: "y=Prime p" "p>2" then
obtain a where a:"[a^(p - 1) = 1] (mod p)" "Triple p a (p - 1) ∈ set ys"
using Cons.prems by auto
then have Bier:"(∀q∈prime_factors (p - 1). [a^((p - 1) div q) ≠ 1] (mod p))"
using Cons.IH Cons.prems(1) by (simp add:y(1))
then have "prime p" using lehmers_theorem[OF _ _a(1)] p>2 by fastforce
}
moreover
{ assume "y=Prime p" "p=2" hence "prime p" by simp }
moreover
{ assume "y=Prime p" then have "p>1"  using Cons.prems  by simp }
ultimately have y_Prime:"y = Prime p ⟹ prime p" by linarith

show ?case
proof (cases "t ∈ set ys")
case True
show ?thesis using Cons.IH[OF _ True] Cons.prems(1) by (cases y) auto
next
case False
thus ?thesis using Cons.prems(2) y_Prime y_Triple by force
qed
qed

subsection {* Completeness *}

text {*
In this section we show completeness of Pratt's proof system, i.e., we show that for
every prime number $p$ there exists a certificate for $p$. We also give an upper
bound for the size of a minimal certificate

The prove we give is constructive. We assume that we have certificates for all prime
factors of $p - 1$ and use these to build a certificate for $p$ from that. It is
important to note that certificates can be concatenated.
*}

lemma valid_cert_appendI:
assumes "valid_cert r"
assumes "valid_cert s"
shows "valid_cert (r @ s)"
using assms
proof (induction r)
case (Cons y ys) then show ?case by (cases y) auto
qed simp

lemma valid_cert_concatI: "(∀x ∈ set xs . valid_cert x) ⟹ valid_cert (concat xs)"
by (induction xs) (auto simp add: valid_cert_appendI)

lemma size_pratt_le:
fixes d::real
assumes "∀ x ∈ set c. size_pratt x ≤ d"
shows "size_cert c ≤ length c * (1 + d)" using assms
by (induction c) (simp_all add: algebra_simps)

fun build_fpc :: "nat ⇒ nat ⇒ nat ⇒ nat list ⇒ pratt list" where
"build_fpc p a r [] = [Triple p a r]" |
"build_fpc p a r (y # ys) = Triple p a r # build_fpc p a (r div y) ys"

text {*
The function @{term build_fpc} helps us to construct a certificate for $p$ from
the certificates for the prime factors of $p - 1$. Called as
@{term "build_fpc p a (p - 1) qs"} where $@{term "qs"} = q_1 \ldots q_n$
is prime decomposition of $p - 1$ such that $q_1 \cdot \dotsb \cdot q_n = @{term "p - 1 :: nat"}$,
it returns the following list of predicates:
$(p,a,p-1), (p,a,\frac{p - 1}{q_1}), (p,a,\frac{p - 1}{q_1 q_2}), \ldots, (p,a,\frac{p-1}{q_1 \ldots q_n}) = (p,a,1)$

I.e., if there is an appropriate $a$ and and a certificate @{term rs} for all
prime factors of $p$, then we can construct a certificate for $p$ as
@{term [display] "Prime p # build_fpc p a (p - 1) qs @ rs"}
*}

text {*
The following lemma shows that @{text "build_fpc"} extends a certificate that
satisfies the preconditions described before to a correct certificate.
*}

lemma correct_fpc:
assumes "valid_cert xs" "p > 1"
assumes "prod_list qs = r" "r ≠ 0"
assumes "∀ q ∈ set qs . Prime q ∈ set xs"
assumes "∀ q ∈ set qs . [a^((p - 1) div q) ≠ 1] (mod p)"
shows "valid_cert (build_fpc p a r qs @ xs)"
using assms
proof (induction qs arbitrary: r)
case Nil thus ?case by auto
next
case (Cons y ys)
have "prod_list ys = r div y" using Cons.prems by auto
then have T_in: "Triple p a (prod_list ys) ∈ set (build_fpc p a (r div y) ys @ xs)"
by (cases ys) auto

have "valid_cert (build_fpc p a (r div y) ys @ xs)"
using Cons.prems by (intro Cons.IH) auto
then have "valid_cert (Triple p a r # build_fpc p a (r div y) ys @ xs)"
using r ≠ 0 T_in Cons.prems by auto
then show ?case by simp
qed

lemma length_fpc:
"length (build_fpc p a r qs) = length qs + 1" by (induction qs arbitrary: r) auto

lemma div_gt_0:
fixes m n :: nat assumes "m ≤ n" "0 < m" shows "0 < n div m"
proof -
have "0 < m div m" using 0 < m div_self by auto
also have "m div m ≤ n div m" using m ≤ n by (rule div_le_mono)
finally show ?thesis .
qed

lemma size_pratt_fpc:
assumes "a ≤ p" "r ≤ p" "0 < a" "0 < r" "0 < p" "prod_list qs = r"
shows "∀x ∈ set (build_fpc p a r qs) . size_pratt x ≤ 3 * log 2 p" using assms
proof (induction qs arbitrary: r)
case Nil
then have "log 2 a ≤ log 2 p" "log 2 r ≤ log 2 p" by auto
then show ?case by simp
next
case (Cons q qs)
then have "log 2 a ≤ log 2 p" "log 2 r ≤ log 2 p" by auto
then have  "log 2 a + log 2 r ≤ 2 * log 2 p" by arith
moreover have "r div q > 0" using Cons.prems by (fastforce intro: div_gt_0)
moreover hence "prod_list qs = r div q" using Cons.prems(6) by auto
moreover have "r div q ≤ p" using r≤p div_le_dividend[of r q] by linarith
ultimately show ?case using Cons by simp
qed

lemma concat_set:
assumes "∀ q ∈ qs . ∃ c ∈ set cs . Prime q ∈ set c"
shows "∀ q ∈ qs . Prime q ∈ set (concat cs)"
using assms by (induction cs) auto

lemma p_in_prime_factorsE:
fixes n :: nat
assumes "p ∈ prime_factors n" "0 < n"
obtains "2 ≤ p" "p ≤ n" "p dvd n" "prime p"
proof
from assms show "prime p" by auto
then show "2 ≤ p" by (auto dest: prime_gt_1_nat)

from assms show "p dvd n" by auto
then show "p ≤ n" using  0 < n by (rule dvd_imp_le)
qed

lemma prime_factors_list_prime:
fixes n :: nat
assumes "prime n"
shows "∃ qs. prime_factors n = set qs ∧ prod_list qs = n ∧ length qs = 1"
using assms by (auto simp add: prime_factorization_prime intro: exI [of _ "[n]"])

lemma prime_factors_list:
fixes n :: nat assumes "3 < n" "¬ prime n"
shows "∃ qs. prime_factors n = set qs ∧ prod_list qs = n ∧ length qs ≥ 2"
using assms
proof (induction n rule: less_induct)
case (less n)
obtain p where "p ∈ prime_factors n" using n > 3 prime_factors_elem by force
then have p':"2 ≤ p" "p ≤ n" "p dvd n" "prime p"
using 3 < n by (auto elim: p_in_prime_factorsE)
{ assume "n div p > 3" "¬ prime (n div p)"
then obtain qs
where "prime_factors (n div p) = set qs" "prod_list qs = (n div p)" "length qs ≥ 2"
using p' by atomize_elim (auto intro: less simp: div_gt_0)
moreover
have "prime_factors (p * (n div p)) = insert p (prime_factors (n div p))"
using 3 < n 2 ≤ p p ≤ n prime p
by (auto simp: prime_factors_product div_gt_0 prime_factors_of_prime)
ultimately
have "prime_factors n = set (p # qs)" "prod_list (p # qs) = n" "length (p#qs) ≥ 2"
using p dvd n by simp_all
hence ?case by blast
}
moreover
{ assume "prime (n div p)"
then obtain qs
where "prime_factors (n div p) = set qs" "prod_list qs = (n div p)" "length qs = 1"
using prime_factors_list_prime by blast
moreover
have "prime_factors (p * (n div p)) = insert p (prime_factors (n div p))"
using 3 < n 2 ≤ p p ≤ n prime p
by (auto simp: prime_factors_product div_gt_0 prime_factors_of_prime)
ultimately
have "prime_factors n = set (p # qs)" "prod_list (p # qs) = n" "length (p#qs) ≥ 2"
using p dvd n by simp_all
hence ?case by blast
} note case_prime = this
moreover
{ assume "n div p = 1"
hence "n = p" using n>3  using One_leq_div[OF p dvd n] p'(2) by force
hence ?case using prime p ¬ prime n by auto
}
moreover
{ assume "n div p = 2"
hence ?case using case_prime by force
}
moreover
{ assume "n div p = 3"
hence ?case using p' case_prime by force
}
ultimately show ?case using p' div_gt_0[of p n] case_prime by fastforce

qed

lemma prod_list_ge:
fixes xs::"nat list"
assumes "∀ x ∈ set xs . x ≥ 1"
shows "prod_list xs ≥ 1" using assms by (induction xs) auto

lemma sum_list_log:
fixes b::real
fixes xs::"nat list"
assumes b: "b > 0" "b ≠ 1"
assumes xs:"∀ x ∈ set xs . x ≥ b"
shows "(∑x←xs. log b x) = log b (prod_list xs)"
using assms
proof (induction xs)
case Nil
thus ?case by simp
next
case (Cons y ys)
have "real (prod_list ys) > 0" using prod_list_ge Cons.prems by fastforce
thus ?case using log_mult[OF Cons.prems(1-2)] Cons by force
qed

lemma concat_length_le:
fixes g :: "nat ⇒ real"
assumes "∀ x ∈ set xs . real (length (f x)) ≤ g x"
shows "length (concat (map f xs)) ≤ (∑x←xs. g x)" using assms
by (induction xs) force+

lemma prime_gt_3_impl_p_minus_one_not_prime:
fixes p::nat
assumes "prime p" "p>3"
shows "¬ prime (p - 1)"
proof
assume "prime (p - 1)"
have "¬ even p" using assms by (simp add: prime_odd_nat)
hence "2 dvd (p - 1)" by presburger
then obtain q where "p - 1 = 2 * q" ..
then have "2 ∈ prime_factors (p - 1)" using p>3
by (auto simp: prime_factorization_times_prime)
thus False using prime_factors_of_prime p>3 prime (p - 1) by auto
qed

text {*
We now prove that Pratt's proof system is complete and derive upper bounds for
the length and the size of the entries of a minimal certificate.
*}

theorem pratt_complete':
assumes "prime p"
shows "∃c. Prime p ∈ set c ∧ valid_cert c ∧ length c ≤ 6*log 2 p - 4 ∧ (∀ x ∈ set c. size_pratt x ≤ 3 * log 2 p)" using assms
proof (induction p rule: less_induct)
case (less p)
{ assume [simp]: "p = 2"
have "Prime p ∈ set [Prime 2, Triple 2 1 1]" by simp
then have ?case by fastforce }
moreover
{ assume [simp]: "p = 3"
let ?cert = "[Prime 3, Triple 3 2 2, Triple 3 2 1, Prime 2, Triple 2 1 1]"

have "length ?cert ≤ 6*log 2 p - 4
⟷ 2 powr 9 ≤ 2 powr (log 2 p * 6)" by auto
also have "… ⟷ True"
by (simp add: powr_powr[symmetric] powr_numeral)
finally have ?case
by (intro exI[where x="?cert"]) (simp add: cong_nat_def)
}
moreover
{ assume "p > 3"
have qlp: "∀q ∈ prime_factors (p - 1) . q < p" using prime p
by (metis One_nat_def Suc_pred le_imp_less_Suc lessI less_trans p_in_prime_factorsE prime_gt_1_nat zero_less_diff)
hence factor_certs:"∀q ∈ prime_factors (p - 1) . (∃c . ((Prime q ∈ set c) ∧ (valid_cert c)
∧ length c ≤ 6*log 2 q - 4) ∧ (∀ x ∈ set c. size_pratt x ≤ 3 * log 2 q))"
by (auto intro: less.IH)
obtain a where a:"[a^(p - 1) = 1] (mod p) ∧ (∀ q. q ∈ prime_factors (p - 1)
⟶ [a^((p - 1) div q) ≠ 1] (mod p))" and a_size: "a > 0" "a < p"
using converse_lehmer[OF prime p] by blast

have "¬ prime (p - 1)" using p>3 prime_gt_3_impl_p_minus_one_not_prime prime p by auto
have "p ≠ 4" using prime p by auto
hence "p - 1 > 3" using p > 3 by auto

then obtain qs where prod_qs_eq:"prod_list qs = p - 1"
and qs_eq:"set qs = prime_factors (p - 1)" and qs_length_eq: "length qs ≥ 2"
using prime_factors_list[OF _ ¬ prime (p - 1)] by auto
obtain f where f:"∀q ∈ prime_factors (p - 1) . ∃ c. f q = c
∧ ((Prime q ∈ set c) ∧ (valid_cert c) ∧ length c ≤ 6*log 2 q - 4)
∧ (∀ x ∈ set c. size_pratt x ≤ 3 * log 2 q)"
using factor_certs by metis
let ?cs = "map f qs"
have cs: "∀q ∈ prime_factors (p - 1) . (∃c ∈ set ?cs . (Prime q ∈ set c) ∧ (valid_cert c)
∧ length c ≤ 6*log 2 q - 4
∧ (∀ x ∈ set c. size_pratt x ≤ 3 * log 2 q))"
using f qs_eq by auto

have cs_cert_size: "∀c ∈ set ?cs . ∀ x ∈ set c. size_pratt x ≤ 3 * log 2 p"
proof
fix c assume "c ∈ set (map f qs)"
then obtain q where "c = f q" and "q ∈ set qs" by auto
hence *:"∀ x ∈ set c. size_pratt x ≤ 3 * log 2 q" using f qs_eq by blast
have "q < p" "q > 0" using qlp q ∈ set qs qs_eq prime_factors_gt_0_nat by auto
show "∀ x ∈ set c. size_pratt x ≤ 3 * log 2 p"
proof
fix x assume "x ∈ set c"
hence "size_pratt x ≤ 3 * log 2 q" using * by fastforce
also have "… ≤ 3 * log 2 p" using q < p q > 0 p > 3 by simp
finally show "size_pratt x ≤ 3 * log 2 p" .
qed
qed

have cs_valid_all: "∀c ∈ set ?cs . valid_cert c"
using f qs_eq by fastforce

have "∀x ∈ set (build_fpc p a (p - 1) qs). size_pratt x ≤ 3 * log 2 p"
using cs_cert_size a_size p > 3 prod_qs_eq by (intro size_pratt_fpc) auto
hence "∀x ∈ set (build_fpc p a (p - 1) qs @ concat ?cs) . size_pratt x ≤ 3 * log 2 p"
using cs_cert_size by auto
moreover
have "Triple p a (p - 1) ∈ set (build_fpc p a (p - 1) qs @ concat ?cs)" by (cases qs) auto
moreover
have "valid_cert ((build_fpc p a (p - 1) qs)@ concat ?cs)"
proof (rule correct_fpc)
show "valid_cert (concat ?cs)"
using cs_valid_all by (auto simp: valid_cert_concatI)
show "prod_list qs = p - 1" by (rule prod_qs_eq)
show "p - 1 ≠ 0" using prime_gt_1_nat[OF prime p] by arith
show "∀ q ∈ set qs . Prime q ∈ set (concat ?cs)"
using concat_set[of "prime_factors (p - 1)"] cs qs_eq by blast
show "∀ q ∈ set qs . [a^((p - 1) div q) ≠ 1] (mod p)" using qs_eq a by auto
qed (insert ‹p > 3›, simp_all)
moreover
{ let ?k = "length qs"

have qs_ge_2:"∀q ∈ set qs . q ≥ 2" using qs_eq
by (auto intro: prime_ge_2_nat)

have "∀x∈set qs. real (length (f x)) ≤ 6 * log 2 (real x) - 4" using f qs_eq by blast
hence "length (concat ?cs) ≤ (∑q←qs. 6*log 2 q - 4)" using concat_length_le
by fast
hence "length (Prime p # ((build_fpc p a (p - 1) qs)@ concat ?cs))
≤ ((∑q←(map real qs). 6*log 2 q - 4) + ?k + 2)"
by (simp add: o_def length_fpc)
also have "… = (6*(∑q←(map real qs). log 2 q) + (-4 * real ?k) + ?k + 2)"
by (simp add: o_def sum_list_subtractf sum_list_triv sum_list_const_mult)
also have "… ≤ 6*log 2 (p - 1) - 4" using ?k≥2 prod_qs_eq sum_list_log[of 2 qs] qs_ge_2
by force
also have "… ≤ 6*log 2 p - 4" using log_le_cancel_iff[of 2 "p - 1" p] p>3 by force
ultimately have "length (Prime p # ((build_fpc p a (p - 1) qs)@ concat ?cs))
≤ 6*log 2 p - 4" by linarith }
ultimately obtain c where c:"Triple p a (p - 1) ∈ set c" "valid_cert c"
"length (Prime p #c) ≤ 6*log 2 p - 4"
"(∀ x ∈ set c. size_pratt x ≤ 3 * log 2 p)" by blast
hence "Prime p ∈ set (Prime p # c)" "valid_cert (Prime p # c)"
"(∀ x ∈ set (Prime p # c). size_pratt x ≤ 3 * log 2 p)"
using a prime p by (auto simp: Primes.prime_gt_Suc_0_nat)
hence ?case using c by blast
}
moreover have "p≥2" using less by (simp add: prime_ge_2_nat)
ultimately show ?case using less by fastforce
qed

text {*
We now recapitulate our results. A number $p$ is prime if and only if there
is a certificate for $p$. Moreover, for a prime $p$ there always is a certificate
whose size is polynomially bounded in the logarithm of $p$.
*}

corollary pratt:
"prime p ⟷ (∃c. Prime p ∈ set c ∧ valid_cert c)"
using pratt_complete' pratt_sound(1) by blast

corollary pratt_size:
assumes "prime p"
shows "∃c. Prime p ∈ set c ∧ valid_cert c ∧ size_cert c ≤ (6 * log 2 p - 4) * (1 + 3 * log 2 p)"
proof -
obtain c where c: "Prime p ∈ set c" "valid_cert c"
and len: "length c ≤ 6*log 2 p - 4" and "(∀ x ∈ set c. size_pratt x ≤ 3 * log 2 p)"
using pratt_complete' assms by blast
hence "size_cert c ≤ length c * (1 + 3 * log 2 p)" by (simp add: size_pratt_le)
also have "… ≤ (6*log 2 p - 4) * (1 + 3 * log 2 p)" using len by simp
finally show ?thesis using c by blast
qed

subsection ‹Proof Method Setup›

text ‹
Using the following lazy disjunction', we can force the simplifier to fully simplify the
first operand of the disjunction first and not evaluate the second one at all if the first
one simplifies to @{term True}. This improved performance three-fold in our simple test.
›

definition LAZY_DISJ where
"LAZY_DISJ a b ⟷ a ∨ b"

lemma eval_LAZY_DISJ [simp]: "LAZY_DISJ False b = b" "LAZY_DISJ True b = True"
by (simp_all add: LAZY_DISJ_def)

lemma LAZY_DISJ_cong [cong]: "a = a' ⟹ LAZY_DISJ a b = LAZY_DISJ a' b"
by simp

text ‹
The following alternative definitions of valid certificates are better suited for
evaluation by the simplifier than the original ones.
›

context
begin

lemma valid_cert_Cons1:
"valid_cert (Prime p#xs) ⟷
p > 1 ∧ (∃t∈set xs. case t of Prime _ ⇒ False |
Triple p' a x ⇒ p' = p ∧ x = p - 1 ∧ mod_exp a (p-1) p = 1 ) ∧ valid_cert xs"
(is "?lhs = ?rhs")
proof
assume ?lhs thus ?rhs by (auto simp: mod_exp_correct cong_nat_def split: pratt.splits)
next
assume ?rhs
hence "p > 1" "valid_cert xs" by blast+
moreover from ‹?rhs› obtain t where "t ∈ set xs" "case t of Prime _ ⇒ False |
Triple p' a x ⇒ p' = p ∧ x = p - 1 ∧ [a^(p-1) = 1] (mod p)"
by (auto simp: cong_nat_def mod_exp_correct cong: pratt.case_cong)
ultimately show ?lhs by (cases t) auto
qed

private lemma Suc_0_mod_eq_Suc_0_iff:
"Suc 0 mod n = Suc 0 ⟷ n ≠ Suc 0"
proof -
consider "n = 0" | "n = Suc 0" | "n > 1" by (cases n) auto
thus ?thesis by cases auto
qed

private lemma Suc_0_eq_Suc_0_mod_iff:
"Suc 0 = Suc 0 mod n ⟷ n ≠ Suc 0"
using Suc_0_mod_eq_Suc_0_iff by (simp add: eq_commute)

lemma valid_cert_Cons2:
"valid_cert (Triple p a x # xs) ⟷ x > 0 ∧ p > 1 ∧ (LAZY_DISJ (x = 1) (
(∃t∈set xs. case t of Prime _ ⇒ False |
Triple p' a' y ⇒ p' = p ∧ a' = a ∧ y dvd x ∧
(let q = x div y in Prime q ∈ set xs ∧ mod_exp a ((p-1) div q) p ≠ 1)))) ∧ valid_cert xs"
(is "?lhs = ?rhs")
proof
assume ?lhs
from ‹?lhs› have pos: "x > 0" and gt_1: "p > 1" and valid: "valid_cert xs" by simp_all
show ?rhs
proof (cases "x = 1")
case True
with ‹?lhs› show ?thesis by auto
next
case False
with ‹?lhs› have "(∃q y. x = q * y ∧ Prime q ∈ set xs ∧ Triple p a y ∈ set xs
∧ [a^((p - 1) div q) ≠ 1] (mod p))" by auto
then guess q y by (elim exE conjE) note qy = this
hence "(∃t∈set xs. case t of Prime _ ⇒ False |
Triple p' a' y ⇒ p' = p ∧ a' = a ∧ y dvd x ∧
(let q = x div y in Prime q ∈ set xs ∧ mod_exp a ((p-1) div q) p ≠ 1))"
using pos by (intro bexI[of _ "Triple p a y"])
(auto simp: Suc_0_mod_eq_Suc_0_iff Suc_0_eq_Suc_0_mod_iff cong_nat_def mod_exp_correct)
with pos gt_1 valid show ?thesis unfolding LAZY_DISJ_def by blast
qed
next
assume ?rhs
hence pos: "x > 0" and gt_1: "p > 1" and valid: "valid_cert xs" by simp_all
show ?lhs
proof (cases "x = 1")
case True
with ‹?rhs› show ?thesis by auto
next
case False
with ‹?rhs› obtain t where t: "t ∈ set xs" "case t of Prime x ⇒ False
| Triple p' a' y ⇒ p' = p ∧ a' = a ∧ y dvd x ∧ (let q = x div y
in Prime q ∈ set xs ∧ mod_exp a ((p - 1) div q) p ≠ 1)" by auto
then obtain y where y: "t = Triple p a y" "y dvd x" "let q = x div y in Prime q ∈ set xs ∧
mod_exp a ((p - 1) div q) p ≠ 1"
by (cases t rule: pratt.exhaust) auto
with gt_1 have y': "let q = x div y in Prime q ∈ set xs ∧ [a^((p - 1) div q) ≠ 1] (mod p)"
by (auto simp: cong_nat_def Let_def mod_exp_correct Suc_0_mod_eq_Suc_0_iff Suc_0_eq_Suc_0_mod_iff)
define q where "q = x div y"
have "∃q y. x = q * y ∧ Prime q ∈ set xs ∧ Triple p a y ∈ set xs
∧ [a^((p - 1) div q) ≠ 1] (mod p)"
by (rule exI[of _ q], rule exI[of _ y]) (insert t y y', auto simp: Let_def q_def)
with pos gt_1 valid show ?thesis unfolding LAZY_DISJ_def by simp
qed
qed

declare valid_cert.simps(2,3) [simp del]

lemmas eval_valid_cert = valid_cert.simps(1) valid_cert_Cons1 valid_cert_Cons2

private lemma bex_simps_lazy: "Bex {} P = False" "Bex (insert x A) P = LAZY_DISJ (P x) (Bex A P)"
unfolding LAZY_DISJ_def by simp_all

lemma pratt_primeI:
assumes "valid_cert xs" "Prime p ∈ set xs"
shows   "prime p"
using pratt_sound[OF assms] by simp

ML_file "pratt.ML"

end

method_setup pratt = ‹
Scan.option (Scan.lift (Args.bracks Pratt.parse_cert)) >>
(SIMPLE_METHOD o HEADGOAL oo Pratt.pratt_tac)
› "Prove primality of natural numbers using Pratt certificates."

text ‹
The following two theorems serve as regression tests -- the first one computes the
certificate in ML automatically, whereas the second one uses a pre-computed certificate.
The first example should not take more than a few seconds; the second one no more than
30 seconds or so.
›

lemma "prime (2503 :: nat)"
by pratt

lemma "prime (131059 :: nat)"
by (pratt [131059, (131059, 2, 131058), (131059, 2, 65529), (131059, 2, 21843),
(131059, 2, 7281), (131059, 2, 2427), (131059, 2, 809), (131059, 2, 1), 809,
(809, 3, 808), (809, 3, 404), (809, 3, 202), (809, 3, 101), (809, 3, 1), 101,
(101, 2, 100), (101, 2, 50), (101, 2, 25), (101, 2, 5), (101, 2, 1), 5,
(5, 2, 4), (5, 2, 2), (5, 2, 1), 3, (3, 2, 2), (3, 2, 1), 2, (2, 1, 1)])

end
`