largest_prime_factor.mlw 2.59 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
(*
Euler project Problem 3: Largest prime factor

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?

*)

module PrimeFactor

12 13 14 15
  use mach.int.Int
  use number.Divisibility
  use number.Prime
  use number.Coprime
16

17 18
  (** returns the smallest divisor of `n` greater than or equal to `d`,
      assuming no divisor between 2 and `d`. *)
19 20 21 22 23 24 25 26 27 28 29 30 31
  let rec smallest_divisor (d n:int) : int
    requires { 2 <= n }
    requires { 2 <= d <= n }
    requires { forall i:int. 2 <= i < d -> not divides i n }
    ensures { d <= result <= n }
    ensures { divides result n }
    ensures { forall i:int. 2 <= i < result -> not divides i n }
    variant { n - d }
  = if d * d > n then begin
      assert { forall i:int. 2 <= i < n /\ divides i n ->
        i >= d && let u = div n i in u * i = n && divides u n &&
        u * i = n && (u >= d -> n >= d * i >= d * d && false)
        && u >= 2 && u < n && false } ; n
Mário Pereira's avatar
Mário Pereira committed
32
    end else if d >= 2 && n % d = 0 then d else
33 34
  smallest_divisor (d+1) n

35 36
  use ref.Ref
  use list.List
37 38 39 40 41 42 43 44

  let largest_prime_factor (n:int) : int
    requires { 2 <= n }
    ensures { prime result }
    ensures { divides result n }
    ensures { forall i:int. result < i <= n -> not (prime i /\ divides i n) }
  = let d = smallest_divisor 2 n in
    let factor = ref d in
Mário Pereira's avatar
Mário Pereira committed
45
    let target = ref (n / d) in
46 47 48 49 50 51 52 53 54 55 56 57 58
    assert { !target * d = n && divides !target n } ;
    assert { forall i:int. prime i /\ divides i n /\ i > d ->
      coprime d i && divides i !target };
    while !target >= 2 do
      invariant { 1 <= !target <= n }
      invariant { 2 <= !factor <= n }
      invariant { divides !factor n }
      invariant { prime !factor }
      invariant { forall i:int. divides i !target /\ i >= 2 ->
        i >= !factor /\ divides i n }
      invariant { forall i:int. prime i /\ divides i n /\ i > !factor ->
        divides i !target }
      variant { !target }
Mário Pereira's avatar
Mário Pereira committed
59
      let oldt = ghost !target in
60 61 62 63 64
      let ghost oldf = !factor in
      assert { divides !target !target && !target >= 2 && !target >= !factor };
      let d = smallest_divisor !factor !target in
      assert { prime d };
      factor := d;
Mário Pereira's avatar
Mário Pereira committed
65
      target := !target / d;
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
      assert { !target * d = oldt && divides !target oldt } ;
      assert { forall i:int. prime i /\ divides i n /\ i > d ->
        i > oldf && divides i oldt && 1 <= d < i
        && coprime d i && divides i !target }
    done;
    !factor

  let test () =
    largest_prime_factor 13195 (* should be 29 *)

  let solve () =
    largest_prime_factor 600851475143

end


(***
Local Variables:
compile-command: "why3ide largest_prime_factor.mlw"
End:
*)