why3/examples/conjugate.mlw

(*
  Conjugate of a partition.

  Author: Jean-Christophe Filliatre (CNRS)


  A partition of an integer n is a way of writing n as a sum of positive
  integers. For instance 7 = 3 + 2 + 2 + 1. Such a partition can be
  displayed as a Ferrer diagram:

            3 * * *
            2 * *
            2 * *
            1 *

  The conjugate of that a partition is another partition of 7, obtained
  by flipping the diagram above along its main diagonal. We get

            4 * * * *
            3 * * *
            1 *

  Equivalently, this is the number of cells in each column of the original
  Ferrer diagram:

              4 3 1
            3 * * *
            2 * *
            2 * *
            1 *

  The following program computes the conjugate of a partition.
  A partition is represented as an array of integers, sorted in
  non-increasing order, with a least a 0 at the end.

  This was inspired by this article:

    Franck Butelle, Florent Hivert, Micaela Mayero, and Frédéric
    Toumazet. Formal Proof of SCHUR Conjugate Function. In Proceedings
    of Calculemus 2010, pages 158-171. Springer-Verlag LNAI, 2010.
*)

module Conjugate

  use int.Int
  use ref.Refint
  use array.Array

  predicate is_partition (a: array int) =
    (* at least one element *)
    a.length > 0 &&
    (* sorted in non-increasing order *)
    (forall i j. 0 <= i <= j < a.length -> a[i] >= a[j]) &&
    (* at least one 0 sentinel *)
    a[a.length - 1] = 0

  (* values in a[0..n[ are > v, and a[n] <= v *)
  predicate numofgt (a: array int) (n: int) (v: int) =
    0 <= n < a.length &&
    (forall j. 0 <= j < n -> v < a[j]) &&
    v >= a[n]

  predicate is_conjugate (a b: array int) =
    b.length > a[0] &&
    forall j. 0 <= j < b.length -> numofgt a b[j] j

  let conjugate (a: array int) : array int
    requires { is_partition a }
    ensures  { is_conjugate a result }
  =
    let b = Array.make (a[0] + 1) 0 in
    let ref partc = 0 in
    while a[partc] <> 0 do
      invariant { 0 <= partc < a.length }
      invariant { forall j. a[partc] <= j < b.length -> numofgt a b[j] j }
      variant   { a.length - partc }
      label L in
      let ghost start = partc in
      let edge = a[partc] in
      incr partc;
      while a[partc] = edge do
        invariant { start <= partc < a.length }
        invariant { forall j. start <= j < partc -> a[j] = edge }
        variant   { a.length - partc }
        incr partc
      done;
      for i = a[partc] to edge - 1 do
        invariant { forall j. edge <= j < b.length -> b[j] = (b at L)[j] }
        invariant { forall j. a[partc] <= j < i -> b[j] = partc }
        b[i] <- partc
      done
    done;
    assert { a[partc] = 0 };
    b

end

module Conjugate32

  use int.Int
  use ref.Ref
  use mach.int.Int32
  use mach.array.Array32

  predicate is_partition (a: array int32) =
    (* at least one element *)
    a.length > 0 &&
    (* sorted in non-increasing order *)
    (forall i j. 0 <= i <= j < a.length -> a[i] >= a[j]) &&
    (* at least one 0 sentinel *)
    a[a.length - 1] = 0

  (* values in a[0..n[ are > v, and a[n] <= v *)
  predicate numofgt (a: array int32) (n: int) (v: int) =
    0 <= n < a.length &&
    (forall j. 0 <= j < n -> v < a[j]) &&
    v >= a[n]

  predicate is_conjugate (a b: array int32) =
    b.length > a[0] &&
    forall j. 0 <= j < b.length -> numofgt a b[j] j

  let conjugate (a: array int32) (b: array int32) : unit
    requires { is_partition a }
    requires { length b = a[0] + 1 }
    requires { forall i. 0 <= i < b.length -> b[i] = 0 }
    ensures  { is_conjugate a b }
  =
    let ref partc = 0 in
    while a[partc] <> 0 do
      invariant { 0 <= partc < a.length }
      invariant { forall j. a[partc] <= j < b.length -> numofgt a b[j] j }
      variant   { a.length - partc }
      label L in
      let ghost start = partc in
      let edge = a[partc] in
      partc := partc + 1;
      while a[partc] = edge do
        invariant { start <= partc < a.length }
        invariant { forall j. start <= j < partc -> a[j] = edge }
        variant   { a.length - partc }
        partc := partc + 1
      done;
      for i = a[partc] to edge - 1 do
        invariant { forall j. edge <= j < b.length -> b[j] = (b at L)[j] }
        invariant { forall j. a[partc] <= j < i -> b[j] = partc }
        b[i] <- partc
      done
    done;
    assert { a[partc] = 0 }

end

module Test

  use int.Int
  use array.Array
  use Conjugate

  let test () ensures { result.length >= 4 } =
    let a = make 5 0 in
    a[0] <- 3; a[1] <- 2; a[2] <- 2; a[3] <- 1;
    conjugate a

  exception BenchFailure

  let bench () raises { BenchFailure -> true } =
    let a = test () in
    if a[0] <> 4 then raise BenchFailure;
    if a[1] <> 3 then raise BenchFailure;
    if a[2] <> 1 then raise BenchFailure;
    if a[3] <> 0 then raise BenchFailure;
    a

end


(*
  Original C code from SCHUR

  Note that arrays are one-based
  (that code was translated from Pascal code where arrays were one-based)

  #define MAX 100

  void conjgte (int A[MAX], int B[MAX]) {
    int i, partc = 1, edge = 0;
    while (A[partc] != 0) {
      edge = A[partc];
      do
        partc = partc + 1;
      while (A[partc] == edge);
      for (i = A[partc] + 1; i <= edge; i++)
        B[i] = partc - 1;
    }
  }
*)