Palindromic prime pyramids (aka truncatable palindromic primes)

PalindromicPrimes.fst

(**
  A truncatable palindromic prime is either a prime digit, or a palindromic
  prime number that contains no 0 and such that removing its leading and
  trailing digits yields a truncatable palindromic prime.

  Here we prove that there are no truncatable palindromic primes with more
  than 9 digits, and that the largest is 733929337.

  See http://www.utm.edu/staff/caldwell/preprints/JRM_prime_pyramids.pdf

  Tested with FStar@dcaba9436

  @summary Palindromic prime pyramids (aka truncatable palindromic primes)
  @author  Santiago Zanella-Beguelin
  @date    27 October 2017
 *)
module PalindromicPrimes

open FStar.Mul
open FStar.List.Tot

(** b divides a when a % b = 0 *)
val divides: pos -> nat -> bool
let divides b a = a % b = 0

(** Computes the greatest divisor of [a] by trial division up to its square root *)
val greatest_divisor: a:pos -> b:pos -> d:pos{d < b} -> Tot pos (decreases (a - b))
let rec greatest_divisor a b d =
  if a < b * b then d
  else greatest_divisor a (b + 1) (if b `divides` a then b else d)

(** A prime is a number greater than 1 whose greatest divisor is 1 *)
val is_prime: nat -> bool
let is_prime a = 1 < a && greatest_divisor a 2 1 = 1

(** Digit reversal *)
val reverse_acc: nat -> nat -> nat
let rec reverse_acc n m = 
  if n <= 9 then m + n else reverse_acc (n / 10) (10 * (m + n % 10))

val reverse: nat -> nat
let reverse n = reverse_acc n 0

(** A palindrome is a number which equals its reversal *)
val is_palindrome: nat -> bool
let is_palindrome n = (n = reverse n)

type palindrome = n:nat{is_palindrome n}

(** Remove the digits at both ends of a palindrome *)
val truncate: palindrome -> nat
let truncate n = reverse (n / 10) / 10

(** The length of a number in digits *)
val length: n:nat -> pos
let rec length n = if n <= 9 then 1 else 1 + length (n / 10)

(** A digression to prove that [length a < length b ==> a < b] *)

val pow10: nat -> pos
let rec pow10 = function
  | 0 -> 1
  | n -> 10 * pow10 (n-1)

(** [10^(length n - 1) <= n < 10^(length n)] *)
val length_interval: n:pos -> Lemma
  (pow10 (length n - 1) <= n /\ n < pow10 (length n))
let rec length_interval n =
  if 9 < n then length_interval (n / 10)

val pow10_monotone: n:nat -> m:nat{n <= m} -> Lemma (pow10 n <= pow10 m)
let rec pow10_monotone n m =
  if n <> m then pow10_monotone n (m - 1)

(** Shorter length implies smaller magnitude *)
val length_lt: a:pos -> b:pos{length a < length b} -> Lemma (a < b)
let length_lt a b =
  length_interval a; length_interval b;
  pow10_monotone (length a) (length b - 1)

val length_reverse_acc: n:nat -> m:pos -> Lemma
  (length (reverse_acc n (10 * m)) = length n + length m)
let rec length_reverse_acc n m =
  if 9 < n then length_reverse_acc (n / 10) (10 * m + n % 10)

(** [reverse] preserves length if the trailing digit is not 0 *)
val length_reverse: n:nat{0 < n % 10} -> Lemma (length (reverse n) = length n)
let length_reverse n =
  if 9 < n then length_reverse_acc (n / 10) (n % 10)

(** Check if a number contains 0 digits *)
val no_zero: nat -> bool
let rec no_zero n =
  if n <= 9 then 0 < n % 10
  else 0 < n % 10 && no_zero (n / 10)

val no_zero_reverse_acc: n:nat{no_zero n} -> m:nat{no_zero m} -> Lemma
  (no_zero (reverse_acc n (10 * m)))
let rec no_zero_reverse_acc n m =
  if 9 < n then no_zero_reverse_acc (n / 10) (10 * m + n % 10)

val no_zero_reverse: n:nat{no_zero n} -> Lemma (no_zero (reverse n))
let no_zero_reverse n =
  if 9 < n then no_zero_reverse_acc (n / 10) (n % 10)

val no_zero_truncate: n:palindrome{2 < length n /\ no_zero n} -> Lemma
  (no_zero (truncate n))
let no_zero_truncate n =
  no_zero_reverse (n / 10);
  length_reverse (n / 10)

(** Truncating a number with more than 2 digits and no 0 removes exactly 2 digits *)
val length_truncate: n:palindrome{2 < length n && no_zero n} -> Lemma
  (length (truncate n) = length n - 2)
let length_truncate n =
  length_reverse (n / 10)

val divides_le_greatest_divisor: a:pos -> b:pos
  -> d:pos{d < b /\ d * d <= a /\ d `divides` a /\
          (forall d'.{:pattern (d' `divides` a)}
            d < d' /\ d' < b ==> ~(d' `divides` a))} ->
  Lemma (ensures d <= greatest_divisor a b d) (decreases (a - b))
let rec divides_le_greatest_divisor a b d =
  if a < b * b then ()
  else divides_le_greatest_divisor a (b + 1) (if b `divides` a then b else d)

val even_composite: a:pos{2 < a} -> Lemma
  (requires a % 10 = 0 \/ a % 10 = 2 \/ a % 10 = 4 \/ a % 10 = 6 \/ a % 10 = 8)
  (ensures  ~(is_prime a))
let even_composite a = divides_le_greatest_divisor a 3 2

val five_composite: a:pos{5 < a} -> Lemma
  (requires a % 10 = 5)
  (ensures  ~(is_prime a))
let five_composite a =
  if 2 `divides` a || (a < 3 * 3 && 3 `divides` a) || a < 5 * 5 then ()
  else divides_le_greatest_divisor a 6 5

(** All primes other than 2 and 5 end in 1, 3, 7, or 9 *)
val prime_last_digit: n:nat{is_prime n /\ n <> 2 /\ n <> 5} -> Lemma
  (let d = n % 10 in d = 1 \/ d = 3 \/ d = 7 \/ d = 9)
let prime_last_digit n =
  let d = n % 10 in
  if d = 0 || d = 2 || d = 4 || d = 6 || d = 8 then even_composite n
  else if d = 5 then five_composite n

(** [truncate n] is smaller than [n] if [n] has no zeros *)
val truncate_lt: n:palindrome{no_zero n} -> Lemma (truncate n < n) [SMTPat (truncate n << n)]
let truncate_lt n =
  if 2 < length n then (length_truncate n; length_lt (truncate n) n)

(**
 * The base of a prime pyramid is either a single prime digit, or
 * a palindromic prime with no 0 that when truncated is also the
 * base of a prime pyramid.
 *)
val is_prime_pyramid: nat -> bool
let rec is_prime_pyramid n =
  if n <= 9 then is_prime n
  else is_prime n && is_palindrome n && no_zero n && is_prime_pyramid (truncate n)

type prime_pyramid = n:nat{is_prime_pyramid n}

(** Add a digit at both ends of a number *)
val add: n:nat{n <= 9} -> nat -> nat
let add d n = 10 * reverse (10 * n + d) + d

(**
 * Adding the truncated digits back to a prime pyramid gives the original number.
 * This is the only lemma we assume.
 *)
assume val add_truncate: n:prime_pyramid{9 < n} ->
  Lemma (n = add (n % 10) (truncate n))

(** Odd numbers *)
let is_odd (n:nat) = n % 2 = 1

type odd = n:nat{is_odd n}

val odd_minus: n:odd{1 < n} -> Lemma (is_odd (n - 2))
let odd_minus n = ()

val odd_plus: n:odd -> Lemma (is_odd (n + 2))
let odd_plus n = ()

#set-options "--initial_fuel 4 --max_fuel 4 --z3rlimit 50"

(** The length of the base of a prime pyramid is odd *)
val length_prime_pyramid: n:prime_pyramid -> Lemma (is_odd (length n))
let rec length_prime_pyramid n =
  if length n <= 3 then ()
  else if length n = 4 then length_truncate n
  else
    begin
    length_prime_pyramid (truncate n);
    length_truncate n;
    odd_plus (length (truncate n))
    end

(** Primes that can be obtained by adding a digit at both ends of a number *)
val next_prime_palindromes: nat -> list nat
let next_prime_palindromes n =
  filter is_prime [ add 1 n; add 3 n; add 7 n; add 9 n ]

(** All prime pyramids in increasing length order *)
let p1 : list nat = [ 2; 3; 5; 7 ]
let p3 : list nat = [ 727; 929; 131; 151; 353; 757; 373 ]
let p5 : list nat = [ 37273; 39293; 11311; 71317; 31513; 33533; 37573; 97579; 93739 ]
let p7 : list nat = [ 3392933; 7392937; 3315133; 1335331; 9375739 ]
let p9 : list nat = [ 733929337; 373929373 ]
let p11: list nat = [ ]

(** List of prime pyramids of odd length [i] *)
val prime_pyramids: i:odd -> list nat
let rec prime_pyramids i =
  if i = 1 then p1
  else flatten (map next_prime_palindromes (prime_pyramids (i - 2)))

val prime_pyramids9: unit -> Lemma (prime_pyramids 9 = p9)
let prime_pyramids9 _ = assert_norm (prime_pyramids 9 = p9)

val prime_pyramids11: unit -> Lemma (prime_pyramids 11 = [])
let prime_pyramids11 _ = assert_norm (prime_pyramids 11 = [])

(** There are no prime pyramids with 11 or more digits *)
val prime_pyramids_i: i:odd{11 <= i} -> Lemma (prime_pyramids i = [])
let rec prime_pyramids_i i = 
  if i = 11 then prime_pyramids11 () else prime_pyramids_i (i - 2)

(** [p1] is an exhaustive list of all prime pyramids of length 1 *)
val p1_complete: n:prime_pyramid{length n = 1} -> Lemma (mem n p1)
let p1_complete n = ()

(** Property of [collect] missing in the standard library *)
val mem_collect: #a:eqtype -> #b:eqtype -> f:(a -> list b) -> l:list a -> x:a -> y:b ->
  Lemma
  (requires mem x l /\ mem y (f x))
  (ensures  mem y (flatten (map f l)))
let rec mem_collect #a #b f l x y = 
  match l with
  | [] -> ()
  | hd :: tl -> 
    append_mem (f hd) (flatten (map f tl)) y;
    if mem y (f hd) then () else mem_collect #a #b f tl x y

(** Property of [filter] missing in the standard library *)
val mem_filter: #a:eqtype -> f:(a -> bool) -> l:list a -> x:a ->
  Lemma
  (requires f x /\ mem x l)
  (ensures  mem x (filter f l))
let rec mem_filter #a f l x =
  match l with
  | [] -> ()
  | hd :: tl -> if x = hd then () else mem_filter #a f tl x 

(** A prime pyramid is directly obtainable from its truncation *)
val mem_next_prime_palindromes_truncate: n:prime_pyramid{length n >= 3} -> 
  Lemma (mem n (next_prime_palindromes (truncate n)))
let mem_next_prime_palindromes_truncate n =
  prime_last_digit n;
  add_truncate n

(** [prime_pyramids i] is an exhaustive list of prime pyramids of length [i] *)
val pi_complete: i:odd -> n:prime_pyramid{length n = i} -> 
  Lemma (mem n (prime_pyramids i))
let rec pi_complete i n =
  if i = 1 then p1_complete n
  else
    begin
    length_truncate n;
    odd_minus i;
    pi_complete (i - 2) (truncate n);
    mem_next_prime_palindromes_truncate n;
    mem_filter is_prime (next_prime_palindromes (truncate n)) n;
    mem_collect next_prime_palindromes (prime_pyramids (i - 2)) (truncate n) n
    end

(** A prime pyramid with 11 digits or more is absurd *)
val prime_pyramids_max: n:prime_pyramid{11 <= length n} -> Lemma False
let prime_pyramids_max n =
  length_prime_pyramid n;
  pi_complete (length n) n;
  prime_pyramids_i (length n)

(** A prime pyramid has at most 9 digits *)
val length_bound: n:prime_pyramid -> Lemma (length n <= 9)
let length_bound n =
  if 9 < length n then (length_prime_pyramid n; prime_pyramids_max n)

(** The largest prime pyramid base *)
let largest: prime_pyramid =
  assert_norm (no_zero 733929337);
  assert_norm (is_prime_pyramid 733929337);
  733929337

#reset-options

(** 733929337 is the largest prime pyramid base *)
val largest_prime_pyramid: n:prime_pyramid -> Lemma (n <= 733929337)
let largest_prime_pyramid n =
  assert_norm (length largest = 9);
  if length n = 9 then
    begin
    pi_complete 9 n;
    prime_pyramids9 ();
    assert_norm (mem n p9 ==> n = largest \/ n = 373929373)
    end
  else if length n < 9 then length_lt n largest
  else length_bound n