-- | Calculs de combinaisons.
module Htirage.Combin where

-- | @n`nCk`k@ retourne le nombre de combinaisons
-- de longueur 'k' d’un ensemble de longueur 'n'.
nCk :: Integral i => i -> i -> i
n`nCk`k | n<0||k<0||n<k = undefined
        | k > n `div` 2 = go n (n - k) -- more efficient and safe with smaller numbers
        | otherwise     = go n k
        where go i j | j == 0    = 1
                     | otherwise = go i (j-1) * (i-j+1) `div` j

-- | @combinOfRank n k r@ retourne les indices de permutation
-- de la combinaison de 'k' entiers parmi @[1..n]@
-- au rang lexicographique 'r' dans @[0..(n`nCk`k)-1]@.
-- 
-- Construit chaque choix de la combinaison en prenant le prochain plus grand
-- dont le successeur engendre un nombre de combinaisons
-- qui dépasse le rang restant à atteindre.
--
-- DOC: <http://www.site.uottawa.ca/~lucia/courses/5165-09/GenCombObj.pdf>, p.26
combinOfRank :: Integral i => i -> i -> i -> [i]
combinOfRank n k rk | rk<0||n`nCk`k<rk = undefined
                    | otherwise = for1K 1 1 rk
	where
	for1K i j r | i <  k    = uptoRank i j r
	            | i == k    = [j+r] -- because when i == k, nbCombs is always 1
	            | otherwise = []
	uptoRank i j r | nbCombs <- (n-j)`nCk`(k-i)
	               , nbCombs <= r = uptoRank i (j+1) (r-nbCombs)
	               | otherwise    = j : for1K (i+1) (j+1) r

-- | @rankOfCombin n ns@ retourne le rang lexicographique dans @[0..(n`nCk`length ns)-1]@
-- de la combinaison 'ns' d’entiers parmi @[1..n]@.
--
-- WARNING: 'ns' doit être triée de manière ascendante.
--
-- Compte le nombre de combinaisons précédant celle de rang 'r'.
--
-- DOC: <http://www.site.uottawa.ca/~lucia/courses/5165-09/GenCombObj.pdf>, pp.24-25
--
-- @
-- 'rankOfCombin' n ('combinOfRank' n k r) == r
-- 'combinOfRank' n ('length' ns) ('rankOfCombin' n ns) == ns
-- @
rankOfCombin :: Integral i => i -> [i] -> i
rankOfCombin n ns | any (\x -> x<1||n<x) ns || n<k = undefined
                  | otherwise = for1K 1 0 0 ns
	where
	k = fromInteger (toInteger (length ns))
	for1K _ r _ []      = r
	for1K i r x1 (x:xs) = for1K (i+1) r' x xs
		where r' = r + sum [ (n-j)`nCk`(k-i)
		                   | j <- [x1+1..x-1]
		                   ]

-- | @permute ps xs@ remplace chaque élément de 'ps'
-- par l’élement qu’il indexe dans 'xs' entre @[1..'length' xs]@.
permute :: [Int] -> [a] -> [a]
permute ps xs = [xs !! pred p | p <- ps]