examples/conjugate.mlw

   1 (*
   2   Conjugate of a partition.
   3
   4   Author: Jean-Christophe Filliatre (CNRS)
   5
   6
   7   A partition of an integer n is a way of writing n as a sum of positive
   8   integers. For instance 7 = 3 + 2 + 2 + 1. Such a partition can be
   9   displayed as a Ferrer diagram:
  10
  11             3 * * *
  12             2 * *
  13             2 * *
  14             1 *
  15
  16   The conjugate of that a partition is another partition of 7, obtained
  17   by flipping the diagram above along its main diagonal. We get
  18
  19             4 * * * *
  20             3 * * *
  21             1 *
  22
  23   Equivalently, this is the number of cells in each column of the original
  24   Ferrer diagram:
  25
  26               4 3 1
  27             3 * * *
  28             2 * *
  29             2 * *
  30             1 *
  31
  32   The following program computes the conjugate of a partition.
  33   A partition is represented as an array of integers, sorted in
  34   non-increasing order, with a least a 0 at the end.
  35
  36   This was inspired by this article:
  37
  38     Franck Butelle, Florent Hivert, Micaela Mayero, and Frédéric
  39     Toumazet. Formal Proof of SCHUR Conjugate Function. In Proceedings
  40     of Calculemus 2010, pages 158-171. Springer-Verlag LNAI, 2010.
  41 *)
  42
  43 module Conjugate
  44
  45   use int.Int
  46   use ref.Refint
  47   use array.Array
  48
  49   predicate is_partition (a: array int) =
  50     (* at least one element *)
  51     a.length > 0 &&
  52     (* sorted in non-increasing order *)
  53     (forall i j. 0 <= i <= j < a.length -> a[i] >= a[j]) &&
  54     (* at least one 0 sentinel *)
  55     a[a.length - 1] = 0
  56
  57   (* values in a[0..n[ are > v, and a[n] <= v *)
  58   predicate numofgt (a: array int) (n: int) (v: int) =
  59     0 <= n < a.length &&
  60     (forall j. 0 <= j < n -> v < a[j]) &&
  61     v >= a[n]
  62
  63   predicate is_conjugate (a b: array int) =
  64     b.length > a[0] &&
  65     forall j. 0 <= j < b.length -> numofgt a b[j] j
  66
  67   let conjugate (a: array int) : array int
  68     requires { is_partition a }
  69     ensures  { is_conjugate a result }
  70   =
  71     let b = Array.make (a[0] + 1) 0 in
  72     let ref partc = 0 in
  73     while a[partc] <> 0 do
  74       invariant { 0 <= partc < a.length }
  75       invariant { forall j. a[partc] <= j < b.length -> numofgt a b[j] j }
  76       variant   { a.length - partc }
  77       label L in
  78       let ghost start = partc in
  79       let edge = a[partc] in
  80       incr partc;
  81       while a[partc] = edge do
  82         invariant { start <= partc < a.length }
  83         invariant { forall j. start <= j < partc -> a[j] = edge }
  84         variant   { a.length - partc }
  85         incr partc
  86       done;
  87       for i = a[partc] to edge - 1 do
  88         invariant { forall j. edge <= j < b.length -> b[j] = (b at L)[j] }
  89         invariant { forall j. a[partc] <= j < i -> b[j] = partc }
  90         b[i] <- partc
  91       done
  92     done;
  93     assert { a[partc] = 0 };
  94     b
  95
  96 end
  97
  98 module Conjugate32
  99
 100   use int.Int
 101   use ref.Ref
 102   use mach.int.Int32
 103   use mach.array.Array32
 104
 105   predicate is_partition (a: array int32) =
 106     (* at least one element *)
 107     a.length > 0 &&
 108     (* sorted in non-increasing order *)
 109     (forall i j. 0 <= i <= j < a.length -> a[i] >= a[j]) &&
 110     (* at least one 0 sentinel *)
 111     a[a.length - 1] = 0
 112
 113   (* values in a[0..n[ are > v, and a[n] <= v *)
 114   predicate numofgt (a: array int32) (n: int) (v: int) =
 115     0 <= n < a.length &&
 116     (forall j. 0 <= j < n -> v < a[j]) &&
 117     v >= a[n]
 118
 119   predicate is_conjugate (a b: array int32) =
 120     b.length > a[0] &&
 121     forall j. 0 <= j < b.length -> numofgt a b[j] j
 122
 123   let conjugate (a: array int32) (b: array int32) : unit
 124     requires { is_partition a }
 125     requires { length b = a[0] + 1 }
 126     requires { forall i. 0 <= i < b.length -> b[i] = 0 }
 127     ensures  { is_conjugate a b }
 128   =
 129     let ref partc = 0 in
 130     while a[partc] <> 0 do
 131       invariant { 0 <= partc < a.length }
 132       invariant { forall j. a[partc] <= j < b.length -> numofgt a b[j] j }
 133       variant   { a.length - partc }
 134       label L in
 135       let ghost start = partc in
 136       let edge = a[partc] in
 137       partc := partc + 1;
 138       while a[partc] = edge do
 139         invariant { start <= partc < a.length }
 140         invariant { forall j. start <= j < partc -> a[j] = edge }
 141         variant   { a.length - partc }
 142         partc := partc + 1
 143       done;
 144       for i = a[partc] to edge - 1 do
 145         invariant { forall j. edge <= j < b.length -> b[j] = (b at L)[j] }
 146         invariant { forall j. a[partc] <= j < i -> b[j] = partc }
 147         b[i] <- partc
 148       done
 149     done;
 150     assert { a[partc] = 0 }
 151
 152 end
 153
 154 module Test
 155
 156   use int.Int
 157   use array.Array
 158   use Conjugate
 159
 160   let test () ensures { result.length >= 4 } =
 161     let a = make 5 0 in
 162     a[0] <- 3; a[1] <- 2; a[2] <- 2; a[3] <- 1;
 163     conjugate a
 164
 165   exception BenchFailure
 166
 167   let bench () raises { BenchFailure -> true } =
 168     let a = test () in
 169     if a[0] <> 4 then raise BenchFailure;
 170     if a[1] <> 3 then raise BenchFailure;
 171     if a[2] <> 1 then raise BenchFailure;
 172     if a[3] <> 0 then raise BenchFailure;
 173     a
 174
 175 end
 176
 177
 178 (*
 179   Original C code from SCHUR
 180
 181   Note that arrays are one-based
 182   (that code was translated from Pascal code where arrays were one-based)
 183
 184   #define MAX 100
 185
 186   void conjgte (int A[MAX], int B[MAX]) {
 187     int i, partc = 1, edge = 0;
 188     while (A[partc] != 0) {
 189       edge = A[partc];
 190       do
 191         partc = partc + 1;
 192       while (A[partc] == edge);
 193       for (i = A[partc] + 1; i <= edge; i++)
 194         B[i] = partc - 1;
 195     }
 196   }
 197 *)