Htirage/Combin.hs

   1 -- | Calculs de combinaisons.
   2 module Htirage.Combin where
   3
   4 import Data.Bool
   5 import Data.Eq (Eq(..))
   6 import Data.Foldable (any, sum)
   7 import Data.Int (Int)
   8 import Data.List ((!!), length)
   9 import Data.Ord (Ord(..))
  10 import Prelude (Integral(..), Num(..), pred, undefined)
  11
  12 -- | @'nCk' n k@ retourne le nombre de combinaisons
  13 -- de longueur 'k' d’un ensemble de longueur 'n'.
  14 --
  15 -- Computed using the formula:
  16 -- @'nCk' n (k+1) == 'nCk' n (k-1) * (n-k+1) / k@
  17 nCk :: Integral i => i -> i -> i
  18 n`nCk`k | n<0||k<0||n<k = undefined
  19         | otherwise     = go 1 1
  20         where
  21         go i acc = if k' < i then acc else go (i+1) (acc * (n-i+1) `div` i)
  22         -- Use a symmetry to compute over smaller numbers,
  23         -- which is more efficient and safer
  24         k' = if n`div`2 < k then n-k else k
  25
  26 -- | @combinOfRank n k r@ retourne les indices de permutation
  27 -- de la combinaison de 'k' entiers parmi @[1..n]@
  28 -- au rang lexicographique 'r' dans @[0..'nCk' n k - 1]@.
  29 --
  30 -- Construit chaque choix de la combinaison en prenant le prochain plus grand
  31 -- dont le successeur engendre un nombre de combinaisons
  32 -- qui dépasse le rang restant à atteindre.
  33 --
  34 -- DOC: <http://www.site.uottawa.ca/~lucia/courses/5165-09/GenCombObj.pdf>, p.26
  35 combinOfRank :: Integral i => i -> i -> i -> [i]
  36 combinOfRank n k rk | rk<0||n`nCk`k<rk = undefined
  37                     | otherwise = for1K 1 1 rk
  38         where
  39         for1K i j r | i <  k    = uptoRank i j r
  40                     | i == k    = [j+r] -- because when i == k, nbCombs is always 1
  41                     | otherwise = []
  42         uptoRank i j r | nbCombs <- (n-j)`nCk`(k-i)
  43                        , nbCombs <= r = uptoRank i (j+1) (r-nbCombs)
  44                        | otherwise    = j : for1K (i+1) (j+1) r
  45
  46 -- | @rankOfCombin n ns@ retourne le rang lexicographique dans @[0..'nCk' n (length ns) - 1]@
  47 -- de la combinaison 'ns' d’entiers parmi @[1..n]@.
  48 --
  49 -- WARNING: 'ns' doit être triée de manière ascendante.
  50 --
  51 -- Compte le nombre de combinaisons précédant celle de rang 'r'.
  52 --
  53 -- DOC: <http://www.site.uottawa.ca/~lucia/courses/5165-09/GenCombObj.pdf>, pp.24-25
  54 --
  55 -- @
  56 -- 'rankOfCombin' n ('combinOfRank' n k r) == r
  57 -- 'combinOfRank' n ('length' ns) ('rankOfCombin' n ns) == ns
  58 -- @
  59 rankOfCombin :: Integral i => i -> [i] -> i
  60 rankOfCombin n ns | any (\x -> x<1||n<x) ns || n<k = undefined
  61                   | otherwise = for1K 1 0 0 ns
  62         where
  63         k = fromInteger (toInteger (length ns))
  64         for1K _ r _ []      = r
  65         for1K i r x1 (x:xs) = for1K (i+1) r' x xs
  66                 where r' = r + sum [ (n-j)`nCk`(k-i)
  67                                    | j <- [x1+1..x-1]
  68                                    ]
  69
  70 -- | @permute ps xs@ remplace chaque élément de 'ps'
  71 -- par l’élement qu’il indexe dans 'xs' entre @[1..'length' xs]@.
  72 permute :: [Int] -> [a] -> [a]
  73 permute ps xs = [xs !! pred p | p <- ps]