examples/largest_prime_factor.mlw

   1 (*
   2 Euler project Problem 3: Largest prime factor
   3
   4 The prime factors of 13195 are 5, 7, 13 and 29.
   5
   6 What is the largest prime factor of the number 600851475143 ?
   7
   8 *)
   9
  10 module PrimeFactor
  11
  12   use mach.int.Int
  13   use number.Divisibility
  14   use number.Prime
  15   use number.Coprime
  16
  17   (** returns the smallest divisor of `n` greater than or equal to `d`,
  18       assuming no divisor between 2 and `d`. *)
  19   let rec smallest_divisor (d n:int) : int
  20     requires { 2 <= n }
  21     requires { 2 <= d <= n }
  22     requires { forall i:int. 2 <= i < d -> not divides i n }
  23     ensures { d <= result <= n }
  24     ensures { divides result n }
  25     ensures { forall i:int. 2 <= i < result -> not divides i n }
  26     variant { n - d }
  27   = if d * d > n then begin
  28       assert { forall i:int. 2 <= i < n /\ divides i n ->
  29         i >= d && let u = div n i in u * i = n && divides u n &&
  30         u * i = n && (u >= d -> n >= d * i >= d * d && false)
  31         && u >= 2 && u < n && false } ; n
  32     end else if d >= 2 && n % d = 0 then d else
  33   smallest_divisor (d+1) n
  34
  35   use ref.Ref
  36   use list.List
  37
  38   let largest_prime_factor (n:int) : int
  39     requires { 2 <= n }
  40     ensures { prime result }
  41     ensures { divides result n }
  42     ensures { forall i:int. result < i <= n -> not (prime i /\ divides i n) }
  43   = let d = smallest_divisor 2 n in
  44     let factor = ref d in
  45     let target = ref (n / d) in
  46     assert { !target * d = n && divides !target n } ;
  47     assert { forall i:int. prime i /\ divides i n /\ i > d ->
  48       coprime d i && divides i !target };
  49     while !target >= 2 do
  50       invariant { 1 <= !target <= n }
  51       invariant { 2 <= !factor <= n }
  52       invariant { divides !factor n }
  53       invariant { prime !factor }
  54       invariant { forall i:int. divides i !target /\ i >= 2 ->
  55         i >= !factor /\ divides i n }
  56       invariant { forall i:int. prime i /\ divides i n /\ i > !factor ->
  57         divides i !target }
  58       variant { !target }
  59       let oldt = ghost !target in
  60       let ghost oldf = !factor in
  61       assert { divides !target !target && !target >= 2 && !target >= !factor };
  62       let d = smallest_divisor !factor !target in
  63       assert { prime d };
  64       factor := d;
  65       target := !target / d;
  66       assert { !target * d = oldt && divides !target oldt } ;
  67       assert { forall i:int. prime i /\ divides i n /\ i > d ->
  68         i > oldf && divides i oldt && 1 <= d < i
  69         && coprime d i && divides i !target }
  70     done;
  71     !factor
  72
  73   let test () =
  74     largest_prime_factor 13195 (* should be 29 *)
  75
  76   let solve () =
  77     largest_prime_factor 600851475143
  78
  79 end
  80
  81
  82 (***
  83 Local Variables:
  84 compile-command: "why3ide largest_prime_factor.mlw"
  85 End:
  86 *)