-- Copyright (c) David Amos, 2008. All rights reserved.

module Math.Combinatorics.FiniteGeometry where

import Data.List as L
import qualified Data.Set as S
-- import qualified Data.Map as M -- not really required
import Math.Algebra.Field.Base
import Math.Algebra.Field.Extension hiding ( (<+>) )
import Math.Algebra.LinearAlgebra -- hiding ( det )

-- import PermutationGroup
-- import SchreierSims as SS


-- !! This should really live somewhere else
-- subsets of size k
combinationsOf 0 _ = [[]]
combinationsOf _ [] = []
combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs



ptsAG 0 fq = [[]]
ptsAG n fq = [x:xs | x <- fq, xs <- ptsAG (n-1) fq]

ptsPG 0 _ = [[1]]
ptsPG n fq = map (0:) (ptsPG (n-1) fq) ++ map (1:) (ptsAG n fq)


-- "projective normal form"
pnf (0:xs) = 0 : pnf xs
pnf (1:xs) = 1 : xs
pnf (x:xs) = 1 : map (* x') xs where x' = recip x

ispnf (0:xs) = ispnf xs
ispnf (1:xs) = True
ispnf _ = False

closureAG ps =
    let multipliers = [ (1 - sum xs) : xs | xs <- ptsAG (k-1) fq ] -- k-vectors over fq whose sum is 1
    in S.toList $ S.fromList [foldl1 (<+>) $ zipWith (*>) m ps | m <- multipliers]
    where n = length $ head ps -- the dimension of the space we're working in
          k = length ps        -- the dimension of the flat
          fq = eltsFq undefined


-- closure of points in PG(n,Fq)
-- take all linear combinations of the points (ie the subspace generated by the points, considered as points in Fq ^(n+1) )
-- then discard all which aren't in PNF (thus dropping back into PG(n,Fq))
closurePG ps = L.sort $ filter ispnf $ map (<*>> ps) $ ptsAG k fq where
    k = length ps
    fq = eltsFq undefined



-- van Lint & Wilson, p325, 332
qtorial n q | n >= 0 = product [(q^i - 1) `div` (q-1) | i <- [1..n]]

-- van Lint & Wilson, p326
qnomial n k q = (n `qtorial` q) `div` ( (k `qtorial` q) * ((n-k) `qtorial` q) )

-- Cameron, p129
numFlatsPG n q k = qnomial (n+1) (k+1) q -- because it's the number of subspaces in AG n+1 

-- Cameron, p137
numFlatsAG n q k = q^(n-k) * qnomial n k q


qtorials q = scanl (*) 1 [(q^i - 1) `div` (q-1) | i <- [1..]]

qnomials q = iterate succ [1] where
    succ xs = L.zipWith3 (\l r c -> l+c*r) (0:xs) (xs++[0]) (iterate (*q) 1)
    -- succ xs = zipWith (+) (0:xs) $ zipWith (*) (xs++[0]) $ iterate (*q) 1 
-- This amounts to saying
-- [n+1,k]_q = [n,k-1]_q + q^k [n,k]_q
-- Cameron, Combinatorics, p126



-- FLATS VIA REDUCED ROW ECHELON FORMS
-- Suggested by Cameron p125


data ZeroOneStar = Zero | One | Star deriving (Eq)

instance Show ZeroOneStar where
    show Zero = "0"
    show One  = "1"
    show Star = "*"

-- reduced row echelon forms
rrefs n k = map (rref 1 1) (combinationsOf k [1..n]) where
    rref r c (x:xs) =
        if c == x
        then zipWith (:) (oneColumn r) (rref (r+1) (c+1) xs)
        else zipWith (:) (starColumn r) (rref r (c+1) (x:xs))
    rref _ c [] = replicate k (replicate (n+1-c) Star)
    oneColumn r = replicate (r-1) Zero ++ One : replicate (k-r) Zero
    starColumn r = replicate (r-1) Star ++ replicate (k+1-r) Zero

flatsPG n fq k = concatMap substStars $ rrefs (n+1) (k+1) where
    substStars (r:rs) = [r':rs' | r' <- substStars' r, rs' <- substStars rs]
    substStars [] = [[]]
    substStars' (Star:xs) = [x':xs' | x' <- fq, xs' <- substStars' xs]
    substStars' (Zero:xs) = map (0:) $ substStars' xs
    substStars' (One:xs) = map (1:) $ substStars' xs
    substStars' [] = [[]]


-- Flats in AG(n,Fq) are just the flats in PG(n,Fq) which are not "at infinity"
flatsAG n fq k = [map tail (r : map (r <+>) rs) | r:rs <- flatsPG n fq k, head r == 1]
-- The head r == 1 condition is saying that we want points which are in the "finite" part of PG(n,Fq), not points at infinity
-- The reason we add r to each of the rs is to bring them into the "finite" part
-- (If you don't do this, it can lead to incorrect results, for example some of the flats having the same closure)