---------------------------------------------------------------------
--
-- Permut.hs -- Haskell functions for generating permutations
--              in linear time (unrankA060117 & unrankA060118)
--              plus a few shoelace filtering functions.
--
-- Coded by Antti Karttunen, December 2002
--
-- See http://www.haskell.org/
--     http://www.iki.fi/~kartturi/
--
-- Load as :load karttu/Permut.hs
--
--
---------------------------------------------------------------------


module Permut(pos0, perminv) where

import Array

pos0 _ [] = error "pos0{Permut}: pos0 []"

pos0 y (x:xs) = if x == y then 0 else 1+(pos0 y xs)

naturals = 0:map (+1) naturals

ints  :: [Int]
ints = 0:map (+1) ints

dbl :: [Int] -> [Int]
dbl (il) = map (\i -> 2*i) il

oddnumbers = [(n+1) | n <- dbl [0..]]


factorials :: [Integer]
factorials = 1:1:[b*((div b a)+1) | (a,b) <- zip factorials (tail factorials)]

facts :: [Int]
facts = 1:1:[b*((div b a)+1) | (a,b) <- zip facts (tail facts)]

nth :: Int -> [Int] -> Int
nth n il = head (reverse (take (n+1) il))

fact :: Int -> Int
fact n = nth n facts

-- fact 0 = 1
-- fact n = n * (fact n-1)

perminv p = take (length p) (map (\x -> pos0 x p) naturals)

-- factexp 23 --> [3,2,1]

factexp      :: Int -> [Int]
factexp 0    = [0]
factexp n    = reverse (factexpX n 2)

factexpX     :: Int -> Int -> [Int]
factexpX 0 i = []
factexpX n i = (mod n i) : factexpX (div n i) (i+1)

-- Note that (factlen n) == (length (factexp n)) for all n >= 1. See:
-- take 721 (map factlen ints)
-- take 721 (map (length . factexp) ints)


factlen      :: Int -> Int
factlen 0    = 0
factlen n    = factlenX n 2

factlenX     :: Int -> Int -> Int
factlenX 0 i = i-2
factlenX n i = factlenX (div n i) (i+1)


-- Perm0Vec = zero-based permutation vectors (lists)
-- implemented here as arrays.

newtype Perm0Vec = MakePerm0Vec (Array Int Int)

toPerm0Vec  :: (Array Int Int) -> Perm0Vec
toPerm0Vec x = MakePerm0Vec x

fromPerm0Vec :: Perm0Vec -> (Array Int Int)
fromPerm0Vec (MakePerm0Vec x) = x

perm0VecToPerm1List :: Perm0Vec -> [Int]
perm0VecToPerm1List p = [(a!i)+1 | i <- range (bounds a)] where a = (fromPerm0Vec p)


instance Show Perm0Vec where
    showsPrec _ p = showList (perm0VecToPerm1List p)



idperm   :: Int -> Perm0Vec
idperm n = toPerm0Vec (array (0,n-1) [(i,i) | i <- [0..n-1]])

--
-- fromPerm0Vec (swapels (idperm 6) 3 0)
-- --> array (0,5) [(0,3),(1,1),(2,2),(3,0),(4,4),(5,5)]
--

swapels  :: Perm0Vec -> Int -> Int -> Perm0Vec
swapels a i j = if i == j then a else toPerm0Vec(b // ([(i,b!j),(j,b!i)]))
                where b = (fromPerm0Vec a)



-- fromPerm0Vec (perm0VecInv (swapels (swapels (idperm 4) 0 1) 1 3))
-- --> array (0,3) [(0,3),(1,0),(2,2),(3,1)]  

perm0VecInv :: Perm0Vec -> Perm0Vec
perm0VecInv a = toPerm0Vec(array (bounds b) (map trinds (assocs b)))
                where { b = (fromPerm0Vec a); trinds (a,b) = (b,a); }


-- The following two algorithms are slight modifications of unrank1
-- algorithm as presented by W. Myrvold and F. Ruskey, in
-- Ranking and Unranking Permutations in Linear Time,
-- Inform. Process. Lett. 79 (2001), no. 6, 281-284.
-- Available on-line: http://www.cs.uvic.ca/~fruskey/Publications/RankPerm.html


unrankA060117  :: Int -> Perm0Vec
unrankA060117 n = unrankA060117x n 1 (idperm (1 + (factlen n)))

unrankA060117x :: Int -> Int -> Perm0Vec -> Perm0Vec
unrankA060117x 0 i p = p
unrankA060117x r i p = swapels (unrankA060117x (div r (i+1)) (i+1) p) i (i-(mod r (i+1)))



unrankA060118  :: Int -> Perm0Vec
unrankA060118 n = unrankA060118x n 1 (idperm (1 + (factlen n)))

-- fixed length unrank:
flunrank  :: Int -> Int -> Perm0Vec
flunrank s n = unrankA060118x n 1 (idperm s)

unrankA060118x :: Int -> Int -> Perm0Vec -> Perm0Vec
unrankA060118x 0 i p = p
unrankA060118x r i p = unrankA060118x (div r (i+1)) (i+1) (swapels p i (i-(mod r (i+1))))


-- take 25 (map (unrankA060117) ints)
-- [[1],[2,1],[1,3,2],[3,1,2],[3,2,1],[2,3,1],
--  [1,2,4,3],[2,1,4,3],[1,4,2,3],[4,1,2,3],[4,2,1,3],[2,4,1,3],
--  [1,4,3,2],[4,1,3,2],[1,3,4,2],[3,1,4,2],[3,4,1,2],[4,3,1,2],
--  [4,2,3,1],[2,4,3,1],[4,3,2,1],[3,4,2,1],[3,2,4,1],[2,3,4,1],[1,2,3,5,4]]
--


-- take 25 (map (unrankA060118) ints)
-- [[1],[2,1],[1,3,2],[2,3,1],[3,2,1],[3,1,2],
--  [1,2,4,3],[2,1,4,3],[1,3,4,2],[2,3,4,1],[3,2,4,1],[3,1,4,2],
--  [1,4,3,2],[2,4,3,1],[1,4,2,3],[2,4,1,3],[3,4,1,2],[3,4,2,1],
--  [4,2,3,1],[4,1,3,2],[4,3,2,1],[4,3,1,2],[4,2,1,3],[4,1,2,3],[1,2,3,5,4]]


a060118permut = [unrankA060118(i) | i <- [0..]]

nTrues :: [Bool] -> Int
nTrues (bl) = nTruesAux bl 0

nTruesAux :: [Bool] -> Int -> Int
nTruesAux [] i = i
nTruesAux (b:bs) i
            | (b == True) = nTruesAux bs i+1
            | otherwise = nTruesAux bs i


------------------------------------------------------------------------
--
-- Filtering functions for various simple necklaces & other subclasses
-- of permutations.
--
------------------------------------------------------------------------


--
-- Note: instead of the external eyelet-ordering
--
--   n     n+1
--  n-1    n+2
--  n-2    n+3
--  ...    ...
--   3     2n-2
--   2     2n-1
--   1     2n
--
-- used in OEIS and the Nature article, we will internally
-- use the eyelet-ordering:
--
--  2n-2    1              n-1
--  2n-4    3              n-2
--  2n-6    5              n-3
--  ...    ...             ...
--   4     2n-5             2
--   2     2n-3             1
--   0     2n-1             0
--                        level
--
-- which means that the eyelets on the left side
-- are all tagged with even number, and those on
-- the right side with odd number respectively.



first_and_last_fixed  :: Perm0Vec -> Bool
first_and_last_fixed p = (a!first == first) && (a!last == last)
                         where a = (fromPerm0Vec p)
                               (first,last) = (bounds a)


eyeletLevel :: Int -> Int -> Int
eyeletLevel maxodd eye
            | even(eye) = (div eye 2)
            | odd(eye)  = (div (maxodd - eye) 2)


eyeletLevels  :: Perm0Vec -> [Int]
eyeletLevels p = [(eyeletLevel maxodd (a!i)) | i <- range (bounds a)]
                 where a = (fromPerm0Vec p)
                       second (a,b) = b
                       maxodd = second(bounds a)


permDeltas  :: Perm0Vec -> [Int]
permDeltas p = [(a!i)-(a!(i-1)) | i <- tail (range (bounds a))]
               where a = (fromPerm0Vec p)


deltas :: [Int] -> [Int]
deltas [] = []
deltas (a:[]) = []
deltas whole@(a:b:_) = b-a : deltas(tail whole)


waxing_et_waning :: [Int] -> Bool
waxing_et_waning [] = False
waxing_et_waning (d:ds)
                 | d < 0     = (waning_only ds)
                 | otherwise = (waxing_et_waning ds)


waning_only :: [Int] -> Bool
waning_only [] = True
waning_only (d:ds)
            | d > 0     = False
            | otherwise = (waning_only ds)

monotoneLacing  :: Perm0Vec -> Bool
monotoneLacing p = waxing_et_waning (deltas (eyeletLevels p))


alternating_deltas :: [Int] -> Bool
alternating_deltas []       = True
alternating_deltas (a:[])   = True
alternating_deltas (a:b:[]) = (a*b < 0)
alternating_deltas whole@(a:b:_)
                   | (a*b > 0) = False
                   | otherwise = alternating_deltas (tail whole)


alternatingPermutation  :: Perm0Vec -> Bool
alternatingPermutation p = (ds == []) || (((head ds) > 0) && (alternating_deltas ds))
                           where ds = (permDeltas p)


-- When the lacing crosses sides every time, it means that
-- the array a (which implements the permutation vector p)
-- contains only even values at even indices, and odd
-- values at odd indices, thus if we sum all (index,value)
-- pairs (obtained with assocs a) in modulo 2, we should
-- get zero if the crossing condition holds.

crossesEveryTime  :: Perm0Vec -> Bool
crossesEveryTime p = 0 == sum (map (\x -> mod x 2) (map add (assocs a)))
                     where { a = (fromPerm0Vec p); add (a,b) = a+b; }


-- We suppose here that elems gives the elements of one-dimensional
-- array in order (0,1,2,3,...)

-- This will return False if there are three even or three odd
-- elements in succession, otherwise True.
-- That is, the lace must cross at least every second time.

noThriceOnTheSameSide  :: Perm0Vec -> Bool
noThriceOnTheSameSide p = noThriceOnTheSameSideX (elems (fromPerm0Vec p))

noThriceOnTheSameSideX :: [Int] -> Bool
noThriceOnTheSameSideX [] = True
noThriceOnTheSameSideX (_:[]) = True
noThriceOnTheSameSideX (_:_:[]) = True
noThriceOnTheSameSideX whole@(a:b:c:_)
                              | (even(a+b) && even(b+c)) = False
                              | otherwise = noThriceOnTheSameSideX (tail whole)



no3ConsecutiveEylets  :: Perm0Vec -> Bool
no3ConsecutiveEylets p = no3ConsecutiveEyletsX (elems (fromPerm0Vec p))

no3ConsecutiveEyletsX :: [Int] -> Bool
no3ConsecutiveEyletsX [] = True
no3ConsecutiveEyletsX (_:[]) = True
no3ConsecutiveEyletsX (_:_:[]) = True
no3ConsecutiveEyletsX whole@(a:b:c:_)
                              | ((b == a+2) && (c == b+2)) = False
                              | ((a == b+2) && (b == c+2)) = False
                              | otherwise = no3ConsecutiveEyletsX (tail whole)


-- 1,2,24,720, etc.
a0xxx00 = [(nTrues (take (fact n) (map first_and_last_fixed (map (\x -> flunrank n x) ints)))) | n <- dbl [1..]]

-- 0,8,64,416
a0xxx01 = [(nTrues (take (fact n) (map monotoneLacing (map (\x -> flunrank n x) ints)))) | n <- dbl [1..]]

-- 1,5,61,1385
a0xxx02 = [(nTrues (take (fact n) (map alternatingPermutation (map (\x -> flunrank n x) ints)))) | n <- dbl [1..]]

-- 1,4,36,576
a0xxx03 = [(nTrues (take (fact n) (map crossesEveryTime (map (\x -> flunrank n x) ints)))) | n <- dbl [1..]]

-- 2,24,504,
a0xxx04 = [(nTrues (take (fact n) (map noThriceOnTheSameSide (map (\x -> flunrank n x) ints)))) | n <- dbl [1..]]

-- 2,24,632,
a0xxx05 = [(nTrues (take (fact n) (map no3ConsecutiveEylets (map (\x -> flunrank n x) ints)))) | n <- dbl [1..]]