examples/euler_sieve.mlw

   1 (** Euler's Sieve
   2
   3     This is a variant of Eratosthenes's sieve with complexity O(N),
   4     where N is the upper limit of the sieve.
   5     Cf https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Euler's_sieve
   6
   7     Euler's Sieve makes use of a linked list of all integers,
   8     and marks composite numbers as in Eratosthenes's sieve.
   9
  10     In the program below, the linked list is implemented with an array
  11     (where each integer is mapped to the next one) and negative values
  12     are used to represent marked integers.
  13     In addition, the list only contains the odd integers (to save half
  14     of space).
  15
  16     Author: Josué Moreau (Université Paris-Saclay)
  17 *)
  18
  19 module ArithmeticResults
  20
  21   use int.Int
  22   use number.Divisibility
  23   use number.Prime
  24   use int.EuclideanDivision
  25
  26   let lemma mult_croissance_locale (n a: int)
  27     requires { n > 0 /\ 0 <= a }
  28     ensures  { n * a < n * (a + 1) } = ()
  29
  30   let rec lemma mult_croissance (n a b: int)
  31     requires { n > 0 /\ 0 <= a < b }
  32     ensures  { n * a < n * b }
  33     variant  { b - a }
  34     = if a + 1 = b then ()
  35       else mult_croissance n (a + 1) b
  36
  37   let lemma comp_mult_2 (n k: int)
  38     requires { n > 0 /\ k >= 2 }
  39     ensures  { n * k >= n * 2 } = ()
  40
  41   let lemma div_croissance_locale1 (i n: int)
  42     requires { 0 <= i /\ n > 0 }
  43     ensures  { div i n <= div (i + 1) n } = ()
  44
  45   let rec lemma div_croissance1 (i j n: int)
  46     requires { 0 <= i <= j /\ n > 0 }
  47     ensures  { div i n <= div j n }
  48     variant  { j - i }
  49     = if i < j then
  50         div_croissance1 i (j - 1) n
  51
  52   let rec lemma div_croissance_locale2 (m i: int)
  53     requires { i > 0 /\ m >= 0 }
  54     ensures  { div m (i + 1) <= div m i }
  55     variant  { m }
  56     = if m > 0 then div_croissance_locale2 (m - 1) i
  57
  58   let rec lemma div_croissance2 (m i j: int)
  59     requires { 0 < i <= j /\ m >= 0 }
  60     ensures  { div m i >= div m j }
  61     variant  { j - i }
  62     = if i < j then div_croissance2 m i (j - 1)
  63
  64   let lemma div_mult_1 (n k max: int)
  65     requires { n > 0 /\ max >= n /\ n * k <= max }
  66     ensures  { k = div (n * k) n <= div max n } = ()
  67
  68   let lemma mult_borne_sous_exp (n a b: int)
  69     requires { a >= 1 /\ b >= 1 /\ n >= 1 /\ a * b < n }
  70     ensures  { a < n /\ b < n } = ()
  71
  72   let rec lemma sq_ineq (a b: int)
  73     requires { a >= 0 /\ b >= 0 }
  74     requires { a * a < b * b }
  75     ensures  { a < b }
  76     variant  { b * b - a * a }
  77     = if (b - 1) * (b - 1) > a * a then
  78         sq_ineq a (b - 1)
  79
  80 end
  81
  82 module DivisibilityResults
  83
  84   use int.Int
  85   use int.EuclideanDivision
  86   use number.Divisibility
  87   use number.Prime
  88
  89   let lemma divides_div (n m k: int)
  90     requires { 2 <= n /\ 2 <= m /\ 2 <= k < n /\ divides n m /\ not divides k m }
  91     ensures  { not divides k (div m n) } = ()
  92
  93   let lemma divides_inf (n m: int)
  94     requires { 2 <= n /\ 2 <= m /\ divides n m /\
  95                forall k. 2 <= k < n -> not divides k m }
  96     ensures  { forall k. 2 <= k < n -> not divides k (div m n) } = ()
  97
  98   let lemma not_prime_divider_limits (n: int)
  99     requires { not prime n /\ 2 <= n }
 100     ensures  { exists i. 2 <= i /\ i * i <= n /\ divides i n } = ()
 101
 102   let lemma no_prod_impl_no_divider (n: int)
 103     requires { n >= 0 }
 104     ensures  { forall i.
 105                  2 <= i < n * n ->
 106                  not (exists k l. 2 <= k < n /\ 2 <= l < i /\ k * l = i ) ->
 107                  not (exists k. 2 <= k < n /\ k <> i /\ divides k i) } = ()
 108
 109   use ArithmeticResults
 110
 111   let lemma not_prime_impl_divisor_under_sqrt (n: int)
 112     requires { n >= 2 }
 113     ensures  { forall i.
 114                  2 <= i < n * n ->
 115                  not prime i ->
 116                  exists j. 2 <= j < i /\ ((j < n)
 117                            by j * j < n * n /\ j >= 0) /\ divides j i } = ()
 118
 119 end
 120
 121 module EulerSieveSpec
 122
 123   use int.Int
 124   use number.Divisibility
 125   use number.Prime
 126   use seq.Seq
 127   use int.EuclideanDivision
 128
 129   use ArithmeticResults
 130   use DivisibilityResults
 131
 132
 133 (*******************************************************************************)
 134 (*                                                                             *)
 135 (*         INVARIANTS SUR LES RELATIONS ENTRE STRUCTURES DE DONNEES            *)
 136 (*                                                                             *)
 137 (*******************************************************************************)
 138
 139   (* Borne sur le suivant de chaque élément de la liste chaînée *)
 140   predicate inv_nexts (nexts: seq int) (n: int) =
 141     forall i. 0 <= i < n ->
 142               i < nexts[i] <= n
 143
 144   (* Tous les éléments éliminés de la liste chaînée sont marqués *)
 145   predicate all_eliminated_marked (marked: seq bool)
 146                                   (nexts: seq int) =
 147     marked.length = nexts.length /\
 148     forall i. 0 <= i < marked.length ->
 149               forall j. i < j < nexts[i] ->
 150                         marked[j]
 151
 152   (* Tous les éléments éliminés, à partir d'un certain indice, sont marqués *)
 153   predicate all_eliminated_marked_partial (marked: seq bool)
 154                                           (nexts: seq int)
 155                                           (min: int) =
 156     marked.length = nexts.length /\
 157     forall i. min <= i < marked.length ->
 158               forall j. i < j < nexts[i] ->
 159                         marked[j]
 160
 161   (* L'élément suivant d'un élément non marqué inférieur à max / n est
 162      non marqué lorsqu'il est lui-même inférieur à max / n *)
 163   predicate not_marked_impl_next_not_marked (marked_old: seq bool)
 164                                             (nexts: seq int)
 165                                             (n: int) =
 166     marked_old.length = nexts.length /\
 167     marked_old.length >= 2 /\
 168     n >= 2 /\
 169     forall i. 2 <= i <= div (marked_old.length - 1) n ->
 170               nexts[i] <= div (marked_old.length - 1) n ->
 171               not marked_old[i] ->
 172               not marked_old[nexts[i]]
 173
 174   (* Les deux tableaux marked sont identiques *)
 175   predicate is_copy (marked: seq bool)
 176                     (marked_old: seq bool) =
 177     marked.length = marked_old.length /\
 178     forall i. 0 <= i < marked.length ->
 179               not marked[i] ->
 180               not marked_old[i]
 181
 182   (* L'élément suivant d'un élément non marqué inférieur à p <= max / n est
 183      non marqué lorsqu'il est inférieur à max / n *)
 184   predicate not_marked_impl_next_not_marked_partial (marked: seq bool)
 185                                                     (nexts: seq int)
 186                                                     (n: int)
 187                                                     (p: int) =
 188     marked.length = nexts.length /\
 189     marked.length >= 2 /\
 190     n >= 2 /\
 191     p <= div (marked.length - 1) n /\
 192     forall i. 2 <= i < p ->
 193               nexts[i] <= div (marked.length - 1) n ->
 194               not marked[i] ->
 195               not marked[nexts[i]]
 196
 197
 198 (*******************************************************************************)
 199 (*                                                                             *)
 200 (*                   PROPRIETES LIEES AUX NOMBRES PREMIERS                     *)
 201 (*                                                                             *)
 202 (*******************************************************************************)
 203
 204   (* Tous les nombres jusqu'à n sont non marqués si et seulement si ils sont
 205      premiers *)
 206   predicate all_primes (marked: seq bool) (n: int) =
 207     forall i. 0 <= i < n -> not marked[i] <-> prime i
 208
 209   (* Tous les multiples de i sont marqués *)
 210   predicate all_multiples_marked (marked: seq bool) (i max: int) =
 211     2 <= i < marked.length /\
 212     forall k. 2 <= k < max ->
 213               i * k < marked.length ->
 214               marked[i * k]
 215
 216   (* Tous les multiples des 2 <= i < n sont marqués *)
 217   predicate previously_marked_multiples (marked: seq bool) (n: int) =
 218     forall i. 2 <= i < n ->
 219               all_multiples_marked marked i marked.length
 220
 221   (* Les multiples des 2 <= i < n sont les uniques nombres marqués *)
 222   predicate only_multiples_marked (marked: seq bool) (n: int) =
 223     forall k. 2 <= k < marked.length ->
 224               marked[k] ->
 225               exists i j. 2 <= i < n /\ 2 <= j < marked.length /\ i * j = k
 226
 227   (* Les multiples de n avec les éléments non marqués précédement sont
 228      marqués *)
 229   predicate prime_multiples_marked (marked_old: seq bool)
 230                                    (marked: seq bool)
 231                                    (n max: int) =
 232     marked_old.length = marked.length /\
 233     n < max <= marked.length /\
 234     forall i. n <= i < max ->
 235               not marked_old[i] ->
 236               n * i < marked_old.length ->
 237               marked[n * i]
 238
 239   (* Invariants de base de remove_products *)
 240   predicate inv_remove_products (nexts: seq int)
 241                                 (marked: seq bool)
 242                                 (n: int) =
 243     nexts.length = marked.length /\
 244     not marked[2] /\
 245     all_primes marked n /\
 246     prime n /\
 247     not marked[n] /\
 248     inv_nexts nexts nexts.length
 249
 250
 251 (*******************************************************************************)
 252 (*                                                                             *)
 253 (*             QUELQUES LEMMES DE CONSERVATION DES STRUCTURES                  *)
 254 (*                                                                             *)
 255 (*******************************************************************************)
 256
 257   let lemma conservation_all_eliminated_marked_on_marked_change
 258       (marked: seq bool)
 259       (nexts: seq int)
 260       (i: int)
 261     requires { marked.length = nexts.length }
 262     requires { inv_nexts nexts nexts.length }
 263     requires { all_eliminated_marked marked nexts }
 264     requires { 0 <= i < marked.length }
 265     ensures  { all_eliminated_marked marked[i <- true] nexts } = ()
 266
 267   let lemma conservation_all_eliminated_marked_on_nexts_change
 268       (marked: seq bool)
 269       (nexts: seq int)
 270       (i v: int)
 271     requires { marked.length = nexts.length }
 272     requires { all_eliminated_marked marked nexts }
 273     requires { inv_nexts nexts marked.length }
 274     requires { 0 <= i < marked.length }
 275     requires { i < v <= marked.length }
 276     requires { forall j. i < j < v -> marked[j] }
 277     ensures  { all_eliminated_marked marked nexts[i <- v] } = ()
 278
 279 (*******************************************************************************)
 280 (*                                                                             *)
 281 (*                         QUELQUES PREDICATS UTILES                           *)
 282 (*                                                                             *)
 283 (*******************************************************************************)
 284
 285   predicate ordered (a: seq int) (n: int) =
 286     forall i j. 0 <= i < j < n -> a[i] < a[j]
 287
 288   predicate all_inf_or_eq (a: seq int) (n k: int) =
 289     forall i. 0 <= i < n -> a[i] <= k
 290
 291 end
 292
 293 module EulerSieve
 294
 295   use int.Int
 296   use number.Divisibility
 297   use number.Prime
 298   use seq.Seq
 299   use int.EuclideanDivision
 300
 301   use ArithmeticResults
 302   use DivisibilityResults
 303
 304   use EulerSieveSpec
 305
 306 (*******************************************************************************)
 307 (*                                                                             *)
 308 (*                   PROPRIETES LIEES AUX NOMBRES PREMIERS                     *)
 309 (*                                                                             *)
 310 (*******************************************************************************)
 311
 312   (* Si tous les multiples des nombres i < n sont marqués,
 313      alors les mutliples de ces multiples sont marqués *)
 314   let lemma multiples_of_marked_are_marked (marked: seq bool) (n: int)
 315     requires { 2 <= n <= marked.length }
 316     requires { previously_marked_multiples marked n }
 317     ensures  { forall k i.
 318                 2 <= k < n ->
 319                 2 <= i < marked.length ->
 320                 k * i < marked.length ->
 321                 forall j.
 322                   1 <= j < marked.length ->
 323                   k * i * j < marked.length ->
 324                   (marked[k * i * j]
 325                    by 1 <= i * j < marked.length
 326                    by k >= 2 /\ i >= 2 /\ j >= 1 /\ k * (i * j) < marked.length) }
 327     = ()
 328
 329   (* Lemme essentiel : si tous les multiples des nombres i < n sont marqués
 330                        et qu'on a marqué tous les produits de n avec les nombres
 331                        non marqués
 332                        alors tous les multiples de n sont marqués *)
 333   let lemma prev_and_new_impl_all_multiples_marked
 334       (marked_old marked: seq bool)
 335       (n max: int)
 336     requires { 2 <= n < marked.length /\ 2 <= max < marked.length /\
 337                marked_old.length = marked.length }
 338     requires { is_copy marked marked_old }
 339     requires { previously_marked_multiples marked_old n }
 340     requires { previously_marked_multiples marked n }
 341     requires { only_multiples_marked marked_old n }
 342     requires { prime_multiples_marked marked_old marked n marked.length }
 343     ensures  { all_multiples_marked marked n marked.length } =
 344     assert { forall i. 2 <= i < n /\ 2 <= i < marked.length ->
 345                        not marked_old[i] ->
 346                        i * n < marked.length ->
 347                        marked[i * n] };
 348     assert { forall k. 2 <= k < marked_old.length ->
 349                        marked_old[k] ->
 350                        n * k < marked_old.length ->
 351                        exists i j. 2 <= i < n /\ 2 <= j < marked_old.length /\
 352                                    i * j = k /\ marked_old[i * j] };
 353     multiples_of_marked_are_marked marked_old n;
 354     multiples_of_marked_are_marked marked n;
 355     assert { marked.length = marked_old.length };
 356     assert { forall k. 2 <= k < marked_old.length ->
 357                        marked_old[k] ->
 358                        n * k < marked_old.length ->
 359                        marked[k * n]
 360              by forall k. 2 <= k < marked_old.length ->
 361                        marked_old[k] ->
 362                        n * k < marked_old.length ->
 363                        marked_old[k * n] };
 364     assert { prime_multiples_marked marked_old marked n marked.length };
 365     assert { (forall k. 2 <= k < marked.length ->
 366                        not marked_old[k] -> n * k < marked.length ->
 367                        marked[n * k]) }
 368
 369 (*******************************************************************************)
 370 (*                                                                             *)
 371 (*             LEMMES DE CONSERVATION LIES AUX NOMBRES PREMIERS                *)
 372 (*                                                                             *)
 373 (*******************************************************************************)
 374
 375   let lemma conservation_only_multiples_marked (marked: seq bool) (n i j: int)
 376     requires { 2 <= i < n /\ 2 <= j < marked.length /\ i * j < marked.length }
 377     requires { only_multiples_marked marked n }
 378     ensures  { only_multiples_marked marked[(i * j) <- true] n } =
 379     assert { forall k. 0 <= k < marked.length ->
 380                        k <> i * j ->
 381                        marked[(i * j) <- true][k] = marked[k] };
 382     assert { forall k.
 383                2 <= k < marked.length ->
 384                k <> i * j ->
 385                marked[(i * j) <- true][k] ->
 386                exists x y. 2 <= x < n /\ 2 <= y < marked.length /\ k = x * y };
 387     assert { marked[i * j] -> 2 <= i < n -> 2 <= j < marked.length ->
 388              exists x y. 2 <= x < n /\ 2 <= y < marked.length /\ x * y = i * j }
 389
 390   let lemma conservation_previously_marked_multiples (marked: seq bool)
 391                                                      (nexts: seq int)
 392                                                      (n: int)
 393     requires { 2 <= n < marked.length /\ marked.length = nexts.length /\
 394                nexts[n] <= marked.length }
 395     requires { previously_marked_multiples marked n }
 396     requires { only_multiples_marked marked (n + 1) }
 397     requires { all_eliminated_marked marked nexts }
 398     requires { not_marked_impl_next_not_marked marked nexts nexts[n] }
 399     requires { all_multiples_marked marked n marked.length }
 400     ensures  { previously_marked_multiples marked nexts[n] } =
 401     assert { previously_marked_multiples marked (n + 1) };
 402     assert { forall i. n < i < nexts[n] -> marked[i] };
 403     assert { forall i. n < i < nexts[n] ->
 404               ((exists k l. 2 <= k <= n /\ 2 <= l < marked.length /\ k * l = i)
 405                by marked[i]) };
 406     multiples_of_marked_are_marked marked (n + 1);
 407     assert { forall i j.
 408               n < i < nexts[n] ->
 409               2 <= j < marked.length ->
 410               i * j < marked.length ->
 411               (marked[i * j]
 412                by (exists k l. 2 <= k <= n /\ 2 <= l < marked.length /\ k * l = i)) };
 413     assert { forall i.
 414               n < i < nexts[n] ->
 415               all_multiples_marked marked i marked.length }
 416
 417   lemma conservation_previously_marked_multiples_on_marked_change:
 418     forall marked n.
 419       previously_marked_multiples marked n ->
 420       forall i. 0 <= i < marked.length ->
 421                 previously_marked_multiples marked[i <- true] n
 422
 423 (*******************************************************************************)
 424 (*                                                                             *)
 425 (*                      INVARIANTS SIMPLES DE FONCTIONS                        *)
 426 (*                                                                             *)
 427 (*******************************************************************************)
 428
 429   let lemma conservation_not_marked_impl_next_not_marked
 430       (marked: seq bool)
 431       (nexts: seq int)
 432       (max n p: int)
 433     requires { not_marked_impl_next_not_marked marked nexts n /\
 434                  nexts[n] > n > 0 /\
 435                  div max nexts[n] <= div max n
 436                  by forall i. p < i < nexts[p] -> marked[i] /\
 437                     not_marked_impl_next_not_marked_partial marked nexts n p }
 438     ensures  { not_marked_impl_next_not_marked marked nexts nexts[n] } = ()
 439
 440   let lemma unchanged_other_elements (s1: seq 'a) (s2: seq 'a) (i: int) (v: 'a)
 441     requires { 0 <= i < length s1 /\ length s1 = length s2 }
 442     requires { s1 = s2[i <- v] }
 443     ensures  { forall j. 0 <= j < length s1 -> i <> j -> s1[j] = s2[j] } = ()
 444
 445 (*******************************************************************************)
 446 (*                                                                             *)
 447 (*              DEFINITION DE LA STRUCTURE DE DONNEE ABSTRAITE                 *)
 448 (*                                                                             *)
 449 (*******************************************************************************)
 450
 451   use mach.int.Int63
 452
 453   type t = private {
 454     mutable ghost nexts: seq int;
 455     mutable ghost marked: seq bool;
 456     max: int63;
 457   }
 458     invariant { max_int > max >= 3 }
 459     invariant { nexts.length = marked.length = max + 1 }
 460     invariant { inv_nexts nexts nexts.length }
 461     invariant { all_eliminated_marked marked nexts }
 462     invariant { forall i. 3 <= i <= max -> mod i 2 = 0 -> marked[i] }
 463     invariant { forall i. 3 <= i < max - 1 ->
 464                           mod i 2 = 1 ->
 465                           mod (Seq.get nexts i) 2 = 1 \/ Seq.get nexts i = max + 1 }
 466     invariant { Seq.get nexts max = max + 1 /\
 467                 (mod (max - 1) 2 = 0 -> Seq.get nexts (max - 1) = max) /\
 468                 (mod (max - 1) 2 = 1 -> Seq.get nexts (max - 1) = max + 1) }
 469
 470     by { nexts = Seq.create 4 (fun i -> i + 1);
 471          marked = Seq.create 4 (fun i -> i < 2);
 472          max = 3 }
 473
 474   val create (max: int63) : t
 475     requires { max_int > max >= 3 }
 476     ensures  { result.max = max }
 477     ensures  { Seq.get result.marked 0 = Seq.get result.marked 1 = true /\
 478                not Seq.get result.marked 2 }
 479     ensures  { forall i. 3 <= i <= max ->
 480                          mod i 2 = 0 <-> Seq.get result.marked i }
 481     ensures  { forall i. 3 <= i < max - 1 ->
 482                          mod i 2 = 0 -> Seq.get result.nexts i = i + 1 }
 483     ensures  { forall i. 3 <= i < max - 1 ->
 484                          mod i 2 = 1 -> Seq.get result.nexts i = i + 2 }
 485     ensures  { forall i. 0 <= i <= max ->
 486                          Seq.get result.marked i -> i < 2 \/ divides 2 i }
 487
 488   val set_next (t: t) (i v: int63) : unit
 489     requires { 0 <= i <= t.max /\ i < v <= t.max + 1 }
 490     requires { mod i 2 = 1 }
 491     requires { forall j. i < j < v -> Seq.get t.marked j }
 492     requires { not Seq.get t.marked i }
 493     requires { mod v 2 = 1 \/ v = t.max + 1 }
 494     writes   { t.nexts }
 495     ensures  { t.nexts = (old t.nexts)[i <- v] }
 496
 497   val get_next (t: t) (i: int63) : int63
 498     requires { 3 <= i <= t.max }
 499     requires { mod i 2 = 1 }
 500     ensures  { 3 <= result <= t.max + 1 }
 501     ensures  { result = t.nexts[i] }
 502     ensures  { mod result 2 = 1 \/ result = t.max + 1 }
 503
 504   val set_mark (t: t) (i: int63) : unit
 505     requires { 0 <= i <= t.max }
 506     requires { mod i 2 = 1 }
 507     writes   { t.marked }
 508     ensures  { t.marked = (old t.marked)[i <- true] }
 509
 510   val get_mark (t: t) (i: int63) : bool
 511     requires { 0 <= i <= t.max }
 512     requires { mod i 2 = 1 }
 513     ensures { result = t.marked[i] }
 514
 515   val get_max (t: t) : int63
 516     ensures { result = t.max /\ result >= 2 }
 517
 518 (*******************************************************************************)
 519 (*                                                                             *)
 520 (*     PREUVE DE REMOVE_PRODUCTS DANS LA STRUCTURE DE DONNEE ABSTRAITE         *)
 521 (*                                                                             *)
 522 (*******************************************************************************)
 523
 524   let remove_products (t: t) (n: int63) : unit
 525     requires { 3 <= n <= t.max /\ n * n <= t.max }
 526     requires { inv_remove_products t.nexts t.marked n }
 527     requires { previously_marked_multiples t.marked n }
 528     requires { only_multiples_marked t.marked n }
 529     requires { not_marked_impl_next_not_marked t.marked t.nexts n }
 530     ensures  { inv_remove_products t.nexts t.marked n }
 531     ensures  { not_marked_impl_next_not_marked t.marked t.nexts t.nexts[n] }
 532     ensures  { previously_marked_multiples t.marked t.nexts[n] }
 533     ensures  { only_multiples_marked t.marked t.nexts[n] }
 534     = let d = get_max t / n in
 535       let ghost max = t.max in
 536       let ghost marked_old = t.marked in
 537       let rec loop (p: int63) (ghost x: int) =
 538         requires { n <= p <= max /\ 3 <= n <= max /\ mod p 2 = 1 /\
 539                    p <= x < t.nexts[p] /\ t.nexts[x] = t.nexts[p] /\
 540                    t.marked[n * n] }
 541         requires { inv_remove_products t.nexts t.marked n }
 542         requires { previously_marked_multiples t.marked n }
 543         requires { not t.marked[p] }
 544         requires { is_copy t.marked marked_old }
 545         requires { all_eliminated_marked_partial marked_old t.nexts x }
 546         requires { not_marked_impl_next_not_marked marked_old t.nexts n }
 547         requires { prime_multiples_marked marked_old t.marked n t.nexts[x] }
 548         requires { not_marked_impl_next_not_marked_partial t.marked t.nexts n p }
 549         requires { only_multiples_marked t.marked (n + 1) }
 550         ensures  { inv_remove_products t.nexts t.marked n }
 551         ensures  { not_marked_impl_next_not_marked t.marked t.nexts t.nexts[n] }
 552         ensures  { is_copy t.marked marked_old }
 553         ensures  { previously_marked_multiples t.marked n }
 554         ensures  { prime_multiples_marked marked_old t.marked n t.marked.length }
 555         ensures  { only_multiples_marked t.marked (n + 1) }
 556         variant  { max - t.nexts[p] }
 557         let next = get_next t p in
 558         if 0 <= next <= get_max t then
 559           if next <= d then begin
 560             assert { n * next <= max by n * next <= n * d by next <= d };
 561             ghost (conservation_only_multiples_marked t.marked (to_int n + 1)
 562                               (to_int n) t.nexts[to_int p]);
 563             let ghost marked_copy = t.marked in
 564             set_mark t (n * next);
 565             unchanged_other_elements t.marked marked_copy
 566                                      (to_int n * to_int next) true;
 567             assert { p < 2 * p < 2 * t.nexts[p] <= n * t.nexts[p] = n * next };
 568             assert { not t.marked[p] by n * next > p /\ p <= length t.marked };
 569             assert { not marked_old[t.nexts[p]] by 2 <= p < next <= div max n };
 570             assert { prime_multiples_marked marked_old t.marked n t.nexts[next]
 571                      by prime_multiples_marked marked_old t.marked n t.nexts[x] };
 572             assert { forall i. 0 <= i < p -> t.nexts[i] <= p
 573                      by forall i. 0 <= i < p ->
 574                                   forall j. i < j < t.nexts[i] -> t.marked[j] };
 575             if get_mark t next then begin
 576               assert { t.nexts[p] <> t.marked.length };
 577               let ghost nexts_copy = t.nexts in
 578               set_next t p (get_next t next);
 579               unchanged_other_elements t.nexts nexts_copy (to_int p)
 580                                        nexts_copy[to_int next];
 581               assert { all_eliminated_marked_partial marked_old t.nexts next
 582                        by p <= x < next /\
 583                           forall j. next <= j < length t.nexts ->
 584                                     t.nexts[j] = nexts_copy[j] };
 585               loop p (to_int next)
 586             end else begin
 587               assert { forall i. p < i < t.nexts[p] -> t.marked[i] };
 588               loop next (to_int next)
 589             end
 590           end else begin
 591             assert { n * next > max
 592                      by n * (div max n + 1) > max >= n * div max n /\
 593                         n * next >= n * (d + 1) by next >= d + 1 };
 594             ghost (conservation_not_marked_impl_next_not_marked t.marked t.nexts
 595                                   (to_int max) (to_int n) (to_int p))
 596         end else
 597           ghost (conservation_not_marked_impl_next_not_marked t.marked t.nexts
 598                                   (to_int max) (to_int n) (to_int p)) in
 599         ghost (conservation_only_multiples_marked t.marked (to_int n + 1)
 600                                                   (to_int n) (to_int n));
 601         let ghost marked_copy = t.marked in
 602         set_mark t (n * n);
 603         unchanged_other_elements t.marked marked_copy (to_int n * to_int n) true;
 604         assert { forall i. 0 <= i < n -> t.nexts[i] <> n * n
 605                  by forall i. 0 <= i < n -> t.nexts[i] <= n < n * n
 606                  by not t.marked[n] /\
 607                     forall i. 0 <= i < n ->
 608                               forall j. i < j < t.nexts[i] -> t.marked[j] };
 609         loop n (to_int n);
 610         ghost (prev_and_new_impl_all_multiples_marked marked_old t.marked
 611                                                       (to_int n) (to_int max));
 612         ghost (conservation_previously_marked_multiples t.marked t.nexts
 613                                                         (to_int n))
 614
 615 (*******************************************************************************)
 616 (*                                                                             *)
 617 (*           QUELQUES LEMMES NECESSAIRES A LA SUITE DE LA PREUVE               *)
 618 (*                                                                             *)
 619 (*******************************************************************************)
 620
 621   let lemma previously_marked_multiples_impl_prime (marked: seq bool) (n: int)
 622     requires { 2 <= n < marked.length /\ not marked[n] }
 623     requires { previously_marked_multiples marked n }
 624     ensures  { prime n } =
 625     assert { forall k.
 626                2 <= k < n -> not divides k n
 627              by forall k.
 628                   2 <= k < n ->
 629                   not (exists i. 2 <= i < marked.length /\ n = k * i) }
 630
 631   let lemma only_multiples_marked_impl_not_marked (marked: seq bool)
 632                                                  (nexts: seq int)
 633                                                  (n: int)
 634     requires { 2 <= n < marked.length }
 635     requires { only_multiples_marked marked nexts[n] }
 636     requires { prime n }
 637     ensures  { not marked[n] } =
 638     assert { forall i j. 2 <= i < n -> 2 <= j < marked.length -> i * j <> n }
 639
 640 end
 641
 642 module EulerSieveImpl
 643
 644   use int.Int
 645   use seq.Seq
 646   use mach.int.Int63
 647   use mach.array.ArrayInt63
 648
 649   use int.Abs
 650   use int.EuclideanDivision
 651   use number.Divisibility
 652   use number.Prime
 653
 654   use DivisibilityResults
 655   use EulerSieveSpec
 656
 657 (*******************************************************************************)
 658 (*                                                                             *)
 659 (*              LEMMES DE CONSERVATION LIES A L'IMPLEMENTATION                 *)
 660 (*                                                                             *)
 661 (*******************************************************************************)
 662
 663   let lemma conservation_inv_arr_on_mark (arr: seq int) (i: int)
 664     requires { forall j k. 0 <= j < length arr ->
 665                            j < k < div (abs arr[j]) 2 -> arr[k] < 0 }
 666     requires { forall i. 0 <= i < length arr ->
 667                          i < div (abs arr[i]) 2 <= length arr }
 668     requires { 0 <= i < length arr }
 669     requires { arr[i] >= 0 }
 670     ensures  { forall j k. 0 <= j < length arr ->
 671                            j < k < div (abs arr[i <- - arr[i]][j]) 2 ->
 672                            arr[i <- - arr[i]][k] < 0 } =
 673     assert { forall j. 0 <= j < length arr ->
 674                        arr[j] < i -> arr[j] = arr[i <- - arr[i]][j] };
 675     assert { forall j. 0 <= j < length arr ->
 676                        arr[i] < j -> arr[j] = arr[i <- - arr[i]][j] }
 677
 678   let lemma conservation_inv_arr_on_jump (arr: seq int) (min i: int)
 679     requires { min >= 0 }
 680     requires { forall j k. min <= j < length arr ->
 681                            j < k < div (abs arr[j]) 2 -> arr[k] < 0 }
 682     requires { forall i. min <= i < length arr ->
 683                          i < div (abs arr[i]) 2 <= length arr }
 684     requires { min <= i < length arr }
 685     requires { 0 <= div arr[i] 2 < length arr }
 686     requires { arr[div arr[i] 2] < 0 }
 687     ensures  { forall j k. min <= j < length arr ->
 688                            j < k < div (abs arr[i <- - arr[div arr[i] 2]][j]) 2 ->
 689                            arr[i <- - arr[div arr[i] 2]][k] < 0 } =
 690     let ghost next = div (Seq.get arr i) 2 in
 691     let ghost arr1 = arr[i <- - Seq.get arr next] in
 692     assert { forall j. min <= j < i -> arr[j] = arr1[j] };
 693     assert { forall j. i < j < length arr -> arr[j] = arr1[j] }
 694
 695 (*******************************************************************************)
 696 (*                                                                             *)
 697 (*         DEFINITION DE L'IMPLEMENTATION DE LA STRUCTURE DE DONNEE            *)
 698 (*                                                                             *)
 699 (*******************************************************************************)
 700
 701   type t = {
 702     mutable ghost nexts: seq int;
 703     mutable ghost marked: seq bool;
 704     arr: array63;
 705     max: int63;
 706     max_arr: int63
 707   }
 708     invariant { max_int > max >= 3 }
 709     invariant { Seq.length nexts = Seq.length marked = max + 1 }
 710     invariant { div (max - 1) 2 = max_arr }
 711     invariant { length arr = max_arr + 1 }
 712     invariant { inv_nexts nexts (Seq.length nexts) }
 713     invariant { all_eliminated_marked marked nexts }
 714     invariant { forall i. 3 <= i <= max -> mod i 2 = 0 -> Seq.get marked i }
 715     invariant { forall i. 3 <= i < max - 1 ->
 716                           mod i 2 = 1 ->
 717                           mod (Seq.get nexts i) 2 = 1 \/ Seq.get nexts i = max + 1 }
 718     invariant { Seq.get nexts max = max + 1 /\
 719                 (mod (max - 1) 2 = 0 -> Seq.get nexts (max - 1) = max) /\
 720                 (mod (max - 1) 2 = 1 -> Seq.get nexts (max - 1) = max + 1) }
 721     (* glueing invariant *)
 722     invariant { forall i. 0 <= i <= max_arr -> -(max + 1) <= arr[i] <= max + 1 }
 723     invariant { forall i. 0 <= i <= max_arr ->
 724                           Seq.get marked (2 * i + 1) <-> arr[i] < 0 }
 725     invariant { forall i. 0 <= i <= max_arr ->
 726                           not Seq.get marked (2 * i + 1) ->
 727                           arr[i] = Seq.get nexts (2 * i + 1) }
 728     invariant { forall i. 0 <= i <= max_arr ->
 729                           Seq.get marked (2 * i + 1) ->
 730                           arr[i] = - Seq.get nexts (2 * i + 1) }
 731     invariant { forall i. 0 <= i <= max_arr ->
 732                           i < div (abs arr[i]) 2 <= max_arr + 1 /\
 733                           abs arr[i] <= max + 1 }
 734     invariant { forall i j. 0 <= i <= max_arr ->
 735                             i < j < div (abs arr[i]) 2 -> arr[j] < 0 }
 736     by { nexts = Seq.create 4 (fun i -> i + 1);
 737          marked = Seq.create 4 (fun i -> i < 2);
 738          arr = ArrayInt63.set (ArrayInt63.make 2 (-2)) 1 4;
 739          max = 3; max_arr = 1 }
 740
 741   let create (max: int63) : t
 742     requires { max_int > max >= 3 }
 743     ensures  { result.max = max }
 744     ensures  { Seq.get result.marked 0 = Seq.get result.marked 1 = true /\
 745                not Seq.get result.marked 2 }
 746     ensures  { forall i. 1 <= i <= div (max - 1) 2 ->
 747                          not Seq.get result.marked (2 * i + 1) }
 748     ensures  { forall i. 2 <= i <= div (max + 1) 2 ->
 749                          2 * i <= max ->
 750                          Seq.get result.marked (2 * i) }
 751     ensures  { forall i. 1 <= i <= div (max - 1) 2 ->
 752                          2 * i + 1 < max - 1 ->
 753                          Seq.get result.nexts (2 * i + 1) = 2 * i + 3 }
 754     ensures  { forall i. 2 <= i <= div (max - 1) 2 ->
 755                          2 * i < max - 1 ->
 756                          Seq.get result.nexts (2 * i) = 2 * i + 1 }
 757     ensures  { forall i. 0 <= i <= max ->
 758                          Seq.get result.marked i -> i < 2 \/ divides 2 i }
 759   = let ghost len = pure { max + 1 } in
 760     let len_arr = (max - 1) / 2 + 1 in
 761     let ghost nexts = Seq.create len (pure { fun (i: int) ->
 762           if i = max then max + 1
 763           else if i = max - 1 then if mod i 2 = 0 then max else max + 1
 764           else if i < 3 || mod i 2 = 0 then i + 1
 765           else i + 2}) in
 766     let ghost marked = Seq.create len (fun i ->
 767           i = 0 || i = 1 || (i > 2 && mod i 2 = 0)) in
 768     let ghost f = pure { fun i -> if i = len_arr - 1 then max + 1
 769           else if i = 0 then -2
 770           else 2 * i + 3 } in
 771     let arr = ArrayInt63.make len_arr (-2) in
 772     for i = 1 to len_arr - 1 do
 773       invariant { forall j. 0 <= j < i -> arr[j] = f j }
 774       arr[i] <- if i = len_arr - 1 then max + 1 else 2 * i + 3
 775     done;
 776     { nexts = nexts;
 777       marked = marked;
 778       arr = arr;
 779       max = max;
 780       max_arr = (max - 1) / 2 }
 781
 782   let set_next (t: t) (i v: int63) : unit
 783     requires { 0 <= i <= t.max /\ i < v <= t.max + 1 }
 784     requires { mod i 2 = 1 }
 785     requires { forall j. i < j < v -> Seq.get t.marked j }
 786     requires { not Seq.get t.marked i }
 787     requires { mod v 2 = 1 \/ v = t.max + 1 }
 788     writes   { t.nexts, t.arr }
 789     ensures  { t.nexts = Seq.set (old t.nexts) i v }
 790   = ghost ( conservation_all_eliminated_marked_on_nexts_change
 791                       t.marked t.nexts (to_int i) (to_int v) );
 792     t.arr[i / 2] <- v;
 793     ghost (t.nexts <- Seq.set t.nexts (to_int i) (to_int v))
 794
 795   let get_next (t: t) (i: int63) : int63
 796     requires { 3 <= i <= t.max }
 797     requires { mod i 2 = 1 }
 798     ensures  { 3 <= result <= t.max + 1 }
 799     ensures  { result = Seq.get t.nexts i }
 800     ensures  { mod result 2 = 1 \/ result = t.max + 1 }
 801   = let x = i / 2 in
 802     if t.arr[x] < 0 then - t.arr[x] else t.arr[x]
 803
 804   let set_mark (t: t) (i: int63) : unit
 805     requires { 0 <= i <= t.max }
 806     requires { mod i 2 = 1 }
 807     writes   { t.marked, t.arr }
 808     ensures  { t.marked = Seq.set (old t.marked) i true }
 809   = let x = i / 2 in
 810     ghost ( conservation_all_eliminated_marked_on_marked_change
 811                               t.marked t.nexts (to_int i));
 812     if t.arr[x] >= 0 then begin
 813       ghost ( conservation_inv_arr_on_mark t.arr.elts (to_int x) );
 814       t.arr[x] <- - t.arr[x]
 815     end;
 816     ghost (t.marked <- Seq.set t.marked (to_int i) true)
 817
 818   let get_mark (t: t) (i: int63) : bool
 819     requires { 0 <= i <= t.max }
 820     requires { mod i 2 = 1 }
 821     ensures  { result = Seq.get t.marked i }
 822   = t.arr[i / 2] < 0
 823
 824   let get_max (t: t) : int63
 825     ensures { result = t.max /\ result >= 2 }
 826   = t.max
 827
 828   clone EulerSieve with type t = t,
 829     val create = create,
 830     val set_next = set_next, val get_next = get_next,
 831     val set_mark = set_mark, val get_mark = get_mark,
 832     val get_max = get_max
 833
 834   predicate inv_count (arr: seq int) (min: int) =
 835     forall i. min <= i < arr.length ->
 836               (i < div (abs arr[i]) 2 <= arr.length /\
 837                - max_int <= arr[i] <= max_int /\
 838                (forall j. i < j < div (abs arr[i]) 2 -> arr[j] < 0) /\
 839                (forall j. 2 * i + 1 < j < abs arr[i] -> not prime j))
 840
 841 (*******************************************************************************)
 842 (*                                                                             *)
 843 (*                         PREUVE DU CRIBLE D'EULER                            *)
 844 (*                                                                             *)
 845 (*******************************************************************************)
 846
 847   let euler_sieve (max: int63) : array63
 848     requires { max_int > max >= 3 }
 849     ensures  { forall i j. 0 <= i < j < result.length -> result[i] < result[j] }
 850     ensures  { forall i. 0 <= i < result.length -> 2 <= result[i] <= max }
 851     ensures  { forall i. 0 <= i < result.length -> prime result[i] }
 852     ensures  { forall i. 2 <= i <= max -> prime i ->
 853                          exists j. 0 <= j < result.length /\ result[j] = i }
 854   = [@vc:do_not_keep_trace] (* traceability info breaks the proof with Vampire on sub goal 13.2 *)
 855     let t = create max in
 856     let rec loop (n: int63) =
 857       requires { 3 <= n <= max }
 858       requires { n * n <= max }
 859       requires { previously_marked_multiples t.marked n }
 860       requires { only_multiples_marked t.marked n }
 861       requires { not_marked_impl_next_not_marked t.marked t.nexts n }
 862       requires { inv_remove_products t.nexts t.marked n }
 863       ensures { forall i. 0 <= i < Seq.length t.marked ->
 864                           not Seq.get t.marked i <-> prime i }
 865       variant { max + 1 - n }
 866         remove_products t n;
 867         let nn = get_next t n in
 868         if nn <= max / nn then begin
 869           assert { nn * nn <= max };
 870           ghost (previously_marked_multiples_impl_prime t.marked (to_int nn));
 871           assert { (forall i. n < i < nn -> not prime i)
 872                    by (forall i.
 873                       n < i < nn ->
 874                       (exists j k.
 875                          2 <= j < nn /\
 876                          2 <= k < Seq.length t.marked /\
 877                          j * k = i)
 878                          by Seq.get t.marked i) };
 879           loop (nn)
 880         end else begin
 881           ghost (not_prime_impl_divisor_under_sqrt (to_int nn));
 882           assert { forall i.
 883                      2 <= i < Seq.length t.marked ->
 884                      not Seq.get t.marked i ->
 885                      not (exists k l. 2 <= k < nn /\ 2 <= l < i /\
 886                                       k * l = i) };
 887           ghost (no_prod_impl_no_divider (to_int nn));
 888           assert { Seq.length t.marked <= nn * nn
 889                    by nn * nn > max
 890                    by nn * (div max nn + 1) > max >= nn * div max nn
 891                       /\ nn * nn >= nn * (div max nn + 1)
 892                    by nn >= div max nn + 1 };
 893           assert { forall i. 2 <= i < Seq.length t.marked ->
 894                              not Seq.get t.marked i ->
 895                              prime i
 896                    by forall i.
 897                         2 <= i < Seq.length t.marked ->
 898                         not Seq.get t.marked i ->
 899                         (not (exists k. 2 <= k < nn /\ k <> i /\
 900                                         divides k i)
 901                          by forall k.
 902                               2 <= k ->
 903                               k <> i ->
 904                               (divides k i <-> exists l. 2 <= l < i /\ k * l = i)
 905                               by divides k i -> 2 <= div i k < i) }
 906         end in
 907     if max >= 9 then loop 3;
 908     assert { forall i. 1 <= i <= t.max_arr -> t.arr[i] >= 0 <-> prime (2 * i + 1) };
 909     assert { forall p i. 1 <= p <= t.max_arr ->
 910                          2 * p + 1 < i < abs t.arr[p] -> not prime i
 911              by mod i 2 = 0 -> (Seq.get t.marked i
 912              by Seq.get t.nexts (2 * i + 1) = t.arr[i]) };
 913     assert { forall j. 1 <= j <= t.max_arr -> t.arr[0 <- 2][j] = t.arr[j] };
 914     let ref cnt = 1:int63 in
 915     let ref p = 1 in
 916     t.arr[0] <- 2;
 917     while 2 * p + 1 <= max do
 918       invariant { 1 <= p <= t.max_arr + 1 }
 919       invariant { 1 <= cnt <= p }
 920       invariant { 2 * p + 1 <= max -> t.arr[p] >= 0 }
 921       invariant { 2 * p + 1 <= max -> prime (2 * p + 1) }
 922       invariant { inv_count t.arr cnt }
 923       invariant { ordered t.arr cnt }
 924       invariant { all_inf_or_eq t.arr cnt (2 * p) }
 925       invariant { forall i. cnt <= i <= t.max_arr ->
 926                             t.arr[i] >= 0 <-> prime (2 * i + 1) }
 927       invariant { forall i. 0 <= i < cnt -> prime t.arr[i] }
 928       invariant { forall i. 2 <= i < 2 * p + 1 -> prime i ->
 929                             exists j. 0 <= j < cnt /\ t.arr[j] = i }
 930       variant { t.max_arr + 1 - p , max + 1 - t.arr[p] }
 931       let next = t.arr[p] / 2 in
 932       if next <= t.max_arr then begin
 933         assert { 1 <= p <= t.max_arr /\ t.arr[p] >= 0 /\ prime (2 * p + 1) };
 934         if t.arr[next] < 0 then begin
 935           label BeforeAssign in
 936           t.arr[p] <- - t.arr[next];
 937           unchanged_other_elements t.arr.elts (pure { t.arr.elts at BeforeAssign })
 938                                    (to_int p) (- to_int t.arr[next]);
 939           assert { forall i. p < i < div t.arr[p] 2 -> t.arr[i] < 0 };
 940           assert { t.arr[p] = abs (t.arr at BeforeAssign)[next]
 941                             > (t.arr at BeforeAssign)[p] > p };
 942           conservation_inv_arr_on_jump (pure { t.arr.elts at BeforeAssign })
 943                                        (to_int cnt) (to_int p);
 944           assert { cnt <= next < length t.arr };
 945           assert { (forall j. 2 * p + 1 < j < (t.arr at BeforeAssign)[p] ->
 946                               not prime j) /\
 947                    (forall j. 2 * next + 1 < j < abs (t.arr at BeforeAssign)[next] ->
 948                               not prime j) }
 949         end else begin
 950           label BeforeAssign in
 951           t.arr[cnt] <- 2 * p + 1;
 952           unchanged_other_elements t.arr.elts (pure { t.arr.elts at BeforeAssign })
 953                                               (to_int cnt) (2 * to_int p + 1);
 954           cnt <- cnt + 1;
 955           p <- next;
 956           assert { forall k. cnt <= k < length (t.arr at BeforeAssign) ->
 957                              t.arr[k] = (t.arr at BeforeAssign)[k] }
 958         end
 959       end else begin
 960         label BeforeAssign in
 961         t.arr[cnt] <- 2 * p + 1;
 962         unchanged_other_elements t.arr.elts (pure { t.arr.elts at BeforeAssign })
 963                                  (to_int cnt) (2 * to_int p + 1);
 964         cnt <- cnt + 1;
 965         p <- t.max_arr + 1;
 966         assert { forall k. cnt <= k < length (t.arr at BeforeAssign) ->
 967                            t.arr[k] = (t.arr at BeforeAssign)[k] }
 968       end
 969     done;
 970     sub t.arr 0 cnt
 971
 972 end