conjugate.mlw 3.59 KB
Newer Older
1 2 3 4 5 6
(*
  Conjugate of a partition.

  Author: Jean-Christophe Filliatre (CNRS)


7
  A partition of an integer n is a way of writing n as a sum of positive
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
  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 import ref.Refint
  use import array.Array

  predicate is_partition (a: array int) =
    (* at least one element *)
    a.length > 0 &&
    (* sorted in non-increasing order *)
    (forall i j: int. 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: int. 0 <= j < n -> v < a[j]) &&
    v >= a[n]

  predicate is_conjugate (a b: array int) =
    b.length > a[0] &&
    forall j: int. 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 partc = ref 0 in
    while a[!partc] <> 0 do
      invariant { 0 <= !partc < a.length }
      invariant { forall j: int. a[!partc] <= j < b.length -> numofgt a b[j] j }
      variant   { a.length - !partc }
      'L:
      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: int. start <= j < !partc -> a[j] = edge }
        variant   { a.length - !partc }
        incr partc
      done;
      for i = a[!partc] to edge - 1 do
        invariant { forall j: int. edge <= j < b.length -> b[j] = (at b 'L)[j] }
        invariant { forall j: int. a[!partc] <= j < i -> b[j] = !partc }
        b[i] <- !partc
      done
    done;
    assert { a[!partc] = 0 };
    b

end

97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
module Test

  use import int.Int
  use import array.Array
  use import 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;
115 116
    if a[3] <> 0 then raise BenchFailure;
    a
117 118 119 120

end


121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
(*
  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;
    }
  }
*)