-- | CSM classes of the (open) strata in the set of /ordered/ @n@-tuples,
-- that is, @Q^n = P^1 x P^1 x ... x P^1@
--
-- Of special interest is the open stratum of distinct points, 
-- since any other stratum can be computed from that stratum 
-- by a simple push-forward.
-- 
-- The open stratum of distinct points can be computed recursively, 
-- since the full space @Q^n@ is the disjoint union of all stratums 
-- (indexed by /set partitions/).
-- 
-- But we also have a recursive formula, which makes the computation 
-- significantly faster.
--

{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances,
             ScopedTypeVariables, Rank2Types, GADTs
  #-}

module Math.RootLoci.CSM.Equivariant.Ordered 
  ( -- * The product of projective lines @P^1 x ... x P^1@
    tangentChernClass
    -- * Diagonal embedding
  , j_star 
  , smallDiagonal
    -- * Recursive computation of the CSM of the strata
  , computeOpenStratumCSM     
  , computeAnyStratumCSM
  , computeClosureOfAnyStratumCSM
    -- * The structure lemma
  , QPow(..)
  , umbralDistinctFormula
  , umbralSubstQPow
  , computeQPolys
    -- * The recursive formula for the @Q_k(a,b)@ polynomials
  , formulaQPoly 
    -- * Formula for the CSM class of the stratum of distinct points
  , formulaDistinctCSM
  , formulaAnyStratumCSM
  ) 
  where

--------------------------------------------------------------------------------

-- Semigroup became a superclass of Monoid
#if MIN_VERSION_base(4,11,0)     
import Data.Foldable
import Data.Semigroup
#endif

import Math.Combinat.Classes
import Math.Combinat.Numbers
import Math.Combinat.Sign
import Math.Combinat.Partitions.Integer ( Partition(..) )
import Math.Combinat.Partitions.Set
import Math.Combinat.Sets

import Math.RootLoci.Algebra
import Math.RootLoci.Geometry
import Math.RootLoci.Misc

import qualified Data.Set as Set ; import Data.Set (Set)

import qualified Math.Algebra.Polynomial.FreeModule as ZMod

import Math.RootLoci.CSM.Equivariant.PushForward

--------------------------------------------------------------------------------
-- * The product of projective lines @P^1 x ... x P^1@

-- | Chern class of the tangent bundle of a product of projective lines.
--
-- The formula is:
--
-- > c(T(P^1 x P^1 ... x P^1)) = prod_i (1 + alpha + beta + 2*omega_i)
--
-- because
--
-- > c(T(PV)) = \prod_k (1 + w_i + omega)  `mod`  prod_k (w_i + omega) 
--
-- and
-- 
-- > (1+alpha+omega) * (1+beta+omega) = 1 + alpha + beta + 2*omega 
--
-- since the quadratic term is c_2 of a line bundle which is zero
--
tangentChernClass
  :: ChernBase base 
  => Int                  -- ^ the number of projective lines
  -> ZMod (Omega base)    -- ^ the tangent chern class of their product
tangentChernClass :: Int -> ZMod (Omega base)
tangentChernClass Int
n = (FreeMod Integer (Omega AB), FreeMod Integer (Omega Chern))
-> ChernBase base => ZMod (Omega base)
forall (f :: * -> *) (g :: * -> *) base.
(f (g AB), f (g Chern)) -> ChernBase base => f (g base)
select2 
  ( Int -> FreeMod Integer (Omega AB)
tangentChernClassAB    Int
n
  , Int -> FreeMod Integer (Omega Chern)
tangentChernClassChern Int
n
  )

tangentChernClassAB
  :: Int                  -- ^ The number of projective lines
  -> ZMod (Omega AB)
tangentChernClassAB :: Int -> FreeMod Integer (Omega AB)
tangentChernClassAB Int
d = [FreeMod Integer (Omega AB)] -> FreeMod Integer (Omega AB)
forall b c.
(Ord b, Monoid b, Eq c, Num c) =>
[FreeMod c b] -> FreeMod c b
ZMod.product [ Int -> FreeMod Integer (Omega AB)
forall c. (Eq c, Num c) => Int -> FreeMod c (Omega AB)
entry Int
i | Int
i<-[Int
1..Int
d] ] where
  entry :: Int -> FreeMod c (Omega AB)
entry Int
i = [(Omega AB, c)] -> FreeMod c (Omega AB)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList
    [ ([Int] -> AB -> Omega AB
forall ab. [Int] -> ab -> Omega ab
Omega []  (Int -> Int -> AB
AB Int
0 Int
0) , c
1)
    , ([Int] -> AB -> Omega AB
forall ab. [Int] -> ab -> Omega ab
Omega []  (Int -> Int -> AB
AB Int
1 Int
0) , c
1)
    , ([Int] -> AB -> Omega AB
forall ab. [Int] -> ab -> Omega ab
Omega []  (Int -> Int -> AB
AB Int
0 Int
1) , c
1)
    , ([Int] -> AB -> Omega AB
forall ab. [Int] -> ab -> Omega ab
Omega [Int
i] (Int -> Int -> AB
AB Int
0 Int
0) , c
2)      -- 2x !
    ]

tangentChernClassChern
  :: Int                  -- ^ The number of projective lines
  -> ZMod (Omega Chern)
tangentChernClassChern :: Int -> FreeMod Integer (Omega Chern)
tangentChernClassChern Int
d = [FreeMod Integer (Omega Chern)] -> FreeMod Integer (Omega Chern)
forall b c.
(Ord b, Monoid b, Eq c, Num c) =>
[FreeMod c b] -> FreeMod c b
ZMod.product [ Int -> FreeMod Integer (Omega Chern)
forall c. (Eq c, Num c) => Int -> FreeMod c (Omega Chern)
entry Int
i | Int
i<-[Int
1..Int
d] ] where
  entry :: Int -> FreeMod c (Omega Chern)
entry Int
i = [(Omega Chern, c)] -> FreeMod c (Omega Chern)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList
    [ ([Int] -> Chern -> Omega Chern
forall ab. [Int] -> ab -> Omega ab
Omega []  (Int -> Int -> Chern
Chern Int
0 Int
0) , c
1)
    , ([Int] -> Chern -> Omega Chern
forall ab. [Int] -> ab -> Omega ab
Omega []  (Int -> Int -> Chern
Chern Int
1 Int
0) , c
1)
    , ([Int] -> Chern -> Omega Chern
forall ab. [Int] -> ab -> Omega ab
Omega [Int
i] (Int -> Int -> Chern
Chern Int
0 Int
0) , c
2)      -- 2x !
    ]

--------------------------------------------------------------------------------
-- * Diagonal embedding

-- | Diagonal embeddings of ordered products of P^1-s
j_star :: ChernBase base => [[Int]] -> ZMod (Omega base) -> ZMod (Omega base)
j_star :: [[Int]] -> ZMod (Omega base) -> ZMod (Omega base)
j_star [[Int]]
indices = FreeMod Integer (Eta base) -> ZMod (Omega base)
forall coeff ab.
(Eq coeff, Num coeff, Ord ab) =>
FreeMod coeff (Eta ab) -> FreeMod coeff (Omega ab)
unsafeEtaToOmega (FreeMod Integer (Eta base) -> ZMod (Omega base))
-> (ZMod (Omega base) -> FreeMod Integer (Eta base))
-> ZMod (Omega base)
-> ZMod (Omega base)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Int]] -> ZMod (Omega base) -> FreeMod Integer (Eta base)
forall base.
ChernBase base =>
[[Int]] -> ZMod (Omega base) -> ZMod (Eta base)
delta_star' [[Int]]
indices where

-- | The CSM of the small diagonal in @P^1 x ... x P^1@
smallDiagonal :: forall base. ChernBase base => Int -> ZMod (Omega base)
smallDiagonal :: Int -> ZMod (Omega base)
smallDiagonal Int
n = [Int] -> ZMod (Omega base)
smallDiagonal' [Int
1..Int
n] where

  smallDiagonal' :: [Int] -> ZMod (Omega base)
  smallDiagonal' :: [Int] -> ZMod (Omega base)
smallDiagonal' [Int]
indices = [[Int]] -> ZMod (Omega base) -> ZMod (Omega base)
forall base.
ChernBase base =>
[[Int]] -> ZMod (Omega base) -> ZMod (Omega base)
j_star [[Int]
indices] (Int -> ZMod (Omega base)
forall base. ChernBase base => Int -> ZMod (Omega base)
tangentChernClass Int
1)

--------------------------------------------------------------------------------
-- * CSM of the strata

-- | Recursively compute the CSM of the Zariski-open set @U^n@ of distinct ordered points
-- in @Q^d = P^1 x ... x P^1@. We can compute this by we can subtract all the distinct 
-- fat diagonals from the Chern class of @Q^d@, and the diagonals are just pushforwards 
-- of the same thing for smaller @d@-s.
--
-- NOTE: We also have a more explicit formula for the result (which is /much/ faster to compute)
-- and we can compare the two.
--
-- Note: Forgetting the alpha\/beta part, this should equal to
--
-- > (1-h1-h2-...-hd)^(d-3)
--
-- But, remember that in this formula, @h_i^2 = 0@ for all i!
--
-- Including also @alpha@ and @beta@ we have instead the umbral formula
--
-- > (q-h1-h2-...-hd)^(d-3)
-- 
-- where we also have to do the umbral substitution @q^k -> Q_k@, and the polynomials @Q_k(alpha,beta)@ 
-- are defined recursively, and are defined for @k >= -3@.
--
computeOpenStratumCSM :: ChernBase base => Int -> ZMod (Omega base)
computeOpenStratumCSM :: Int -> ZMod (Omega base)
computeOpenStratumCSM = (forall base. ChernBase base => Int -> ZMod (Omega base))
-> forall base. ChernBase base => Int -> ZMod (Omega base)
forall key (f :: * -> *) (g :: * -> *).
CacheKey key =>
(forall base. ChernBase base => key -> f (g base))
-> forall base. ChernBase base => key -> f (g base)
polyCache2 forall base. ChernBase base => Int -> ZMod (Omega base)
calcOpenStratumCSM  where
             
  calcOpenStratumCSM :: forall b. ChernBase b => Int -> ZMod (Omega b)
  calcOpenStratumCSM :: Int -> ZMod (Omega b)
calcOpenStratumCSM Int
d
    | Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0     =  ZMod (Omega b)
forall b c. (Monoid b, Num c) => FreeMod c b
ZMod.one 
    | Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1     =  Int -> ZMod (Omega b)
forall base. ChernBase base => Int -> ZMod (Omega base)
tangentChernClass Int
1
    | Bool
otherwise  = (Int -> ZMod (Omega b)
forall base. ChernBase base => Int -> ZMod (Omega base)
tangentChernClass Int
d) ZMod (Omega b) -> ZMod (Omega b) -> ZMod (Omega b)
forall b c.
(Ord b, Eq c, Num c) =>
FreeMod c b -> FreeMod c b -> FreeMod c b
`ZMod.sub` ([ZMod (Omega b)] -> ZMod (Omega b)
forall b c. (Ord b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b
ZMod.sum [ZMod (Omega b)]
diagonals)
    where
      diagonals :: [ZMod (Omega b)]
diagonals = 
        [ SetPartition -> ZMod (Omega b)
forall base. ChernBase base => SetPartition -> ZMod (Omega base)
computeAnyStratumCSM SetPartition
setp
        | SetPartition
setp <- Int -> [SetPartition]
setPartitions Int
d 
        , let k :: Int
k = SetPartition -> Int
forall a. HasNumberOfParts a => a -> Int
numberOfParts SetPartition
setp
        , Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
d
        ]


-- | Simply the pushforward of the CSM of the open stratum along the
-- diagonal map corresponding to the given set partition 
--
computeAnyStratumCSM :: ChernBase base => SetPartition -> ZMod (Omega base)
computeAnyStratumCSM :: SetPartition -> ZMod (Omega base)
computeAnyStratumCSM (SetPartition [[Int]]
pps) = ([[Int]] -> ZMod (Omega base) -> ZMod (Omega base)
forall base.
ChernBase base =>
[[Int]] -> ZMod (Omega base) -> ZMod (Omega base)
j_star [[Int]]
pps (ZMod (Omega base) -> ZMod (Omega base))
-> ZMod (Omega base) -> ZMod (Omega base)
forall a b. (a -> b) -> a -> b
$ Int -> ZMod (Omega base)
forall base. ChernBase base => Int -> ZMod (Omega base)
computeOpenStratumCSM (Int -> ZMod (Omega base)) -> Int -> ZMod (Omega base)
forall a b. (a -> b) -> a -> b
$ [[Int]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Int]]
pps)

-- | We sum over the closure
computeClosureOfAnyStratumCSM :: ChernBase base => SetPartition -> ZMod (Omega base)
computeClosureOfAnyStratumCSM :: SetPartition -> ZMod (Omega base)
computeClosureOfAnyStratumCSM SetPartition
setp = [ZMod (Omega base)] -> ZMod (Omega base)
forall b c. (Ord b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b
ZMod.sum
  [ SetPartition -> ZMod (Omega base)
forall base. ChernBase base => SetPartition -> ZMod (Omega base)
computeAnyStratumCSM SetPartition
p | SetPartition
p <- Set SetPartition -> [SetPartition]
forall a. Set a -> [a]
Set.toList (SetPartition -> Set SetPartition
closureSetOfSetPartition SetPartition
setp) ] 

--------------------------------------------------------------------------------
-- * The structure lemma

-- | A formal monomial @q^k@
newtype QPow = QPow Int deriving (QPow -> QPow -> Bool
(QPow -> QPow -> Bool) -> (QPow -> QPow -> Bool) -> Eq QPow
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: QPow -> QPow -> Bool
$c/= :: QPow -> QPow -> Bool
== :: QPow -> QPow -> Bool
$c== :: QPow -> QPow -> Bool
Eq,Eq QPow
Eq QPow
-> (QPow -> QPow -> Ordering)
-> (QPow -> QPow -> Bool)
-> (QPow -> QPow -> Bool)
-> (QPow -> QPow -> Bool)
-> (QPow -> QPow -> Bool)
-> (QPow -> QPow -> QPow)
-> (QPow -> QPow -> QPow)
-> Ord QPow
QPow -> QPow -> Bool
QPow -> QPow -> Ordering
QPow -> QPow -> QPow
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: QPow -> QPow -> QPow
$cmin :: QPow -> QPow -> QPow
max :: QPow -> QPow -> QPow
$cmax :: QPow -> QPow -> QPow
>= :: QPow -> QPow -> Bool
$c>= :: QPow -> QPow -> Bool
> :: QPow -> QPow -> Bool
$c> :: QPow -> QPow -> Bool
<= :: QPow -> QPow -> Bool
$c<= :: QPow -> QPow -> Bool
< :: QPow -> QPow -> Bool
$c< :: QPow -> QPow -> Bool
compare :: QPow -> QPow -> Ordering
$ccompare :: QPow -> QPow -> Ordering
$cp1Ord :: Eq QPow
Ord,Int -> QPow -> ShowS
[QPow] -> ShowS
QPow -> String
(Int -> QPow -> ShowS)
-> (QPow -> String) -> ([QPow] -> ShowS) -> Show QPow
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [QPow] -> ShowS
$cshowList :: [QPow] -> ShowS
show :: QPow -> String
$cshow :: QPow -> String
showsPrec :: Int -> QPow -> ShowS
$cshowsPrec :: Int -> QPow -> ShowS
Show)

#if MIN_VERSION_base(4,11,0)     

instance Semigroup QPow where
  <> :: QPow -> QPow -> QPow
(<>) (QPow Int
e) (QPow Int
f) = Int -> QPow
QPow (Int
eInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
f)

instance Monoid QPow where
  mempty :: QPow
mempty = Int -> QPow
QPow Int
0

#else

instance Monoid QPow where
  mempty = QPow 0
  mappend (QPow e) (QPow f) = QPow (e+f)

#endif

instance Pretty QPow where
  pretty :: QPow -> String
pretty (QPow Int
k) = String -> Int -> String
showVarPower String
"q" Int
k

--------------------------------------------------------------------------------

-- | The umbral formula for the open stratum of the CSM of distinct ordered point:
--
-- > (q - u1 - u2 - ... - un)^(n-3)
--
-- where @u_i^2 = 1@. This also works @n = 0,1,2,3@
-- For these we have the expansion:
--
-- > (q - u1 - u2 - u3)^0   =  q^0
-- > (q - u1 - u2     )^-1  =  1/q + u1/q^2 + u2/q^2 + (2*u1*u2)/q^3
-- > (q - u1          )^-2  =  1/q^2 + (2*u1)/q^3
-- > (q               )^-3  =  1/q^3
--
umbralDistinctFormula :: Int -> ZMod (Omega QPow)
umbralDistinctFormula :: Int -> ZMod (Omega QPow)
umbralDistinctFormula Int
n
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
0  = String -> ZMod (Omega QPow)
forall a. HasCallStack => String -> a
error String
"umbralDistinct: n should be nonnegative"
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0  = Omega QPow -> ZMod (Omega QPow)
forall c b. Num c => b -> FreeMod c b
ZMod.generator (Omega QPow -> ZMod (Omega QPow))
-> Omega QPow -> ZMod (Omega QPow)
forall a b. (a -> b) -> a -> b
$ [Int] -> Int -> Omega QPow
monom [] (-Int
3)
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1  = [(Omega QPow, Integer)] -> ZMod (Omega QPow)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList  
                [ ([Int] -> Int -> Omega QPow
monom []    (-Int
2) , Integer
1) 
                , ([Int] -> Int -> Omega QPow
monom [Int
1]   (-Int
3) , Integer
2)
                ]
  | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2  = [(Omega QPow, Integer)] -> ZMod (Omega QPow)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList  
                [ ([Int] -> Int -> Omega QPow
monom []    (-Int
1) , Integer
1)
                , ([Int] -> Int -> Omega QPow
monom [Int
1]   (-Int
2) , Integer
1)
                , ([Int] -> Int -> Omega QPow
monom [Int
2]   (-Int
2) , Integer
1)
                , ([Int] -> Int -> Omega QPow
monom [Int
1,Int
2] (-Int
3) , Integer
2)
                ]
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
3  = [ZMod (Omega QPow)] -> ZMod (Omega QPow)
forall b c. (Ord b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b
ZMod.sum
                [ Integer -> ZMod (Omega QPow) -> ZMod (Omega QPow)
forall b c. (Ord b, Eq c, Num c) => c -> FreeMod c b -> FreeMod c b
ZMod.scale Integer
coeff (ZMod (Omega QPow) -> ZMod (Omega QPow))
-> ZMod (Omega QPow) -> ZMod (Omega QPow)
forall a b. (a -> b) -> a -> b
$ (Int -> [Omega QPow] -> ZMod (Omega QPow)
forall a. (Ord a, Monoid a) => Int -> [a] -> ZMod a
symPoly (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k) [Omega QPow]
us) ZMod (Omega QPow) -> ZMod (Omega QPow) -> ZMod (Omega QPow)
forall a. Num a => a -> a -> a
* (Omega QPow -> ZMod (Omega QPow)
forall c b. Num c => b -> FreeMod c b
ZMod.generator (Omega QPow -> ZMod (Omega QPow))
-> Omega QPow -> ZMod (Omega QPow)
forall a b. (a -> b) -> a -> b
$ [Int] -> Int -> Omega QPow
monom [] Int
k)
                | Int
k<-[Int
0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3]
                , let coeff :: Integer
coeff = Int -> Integer -> Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k) (Int -> Integer
forall a. Integral a => a -> Integer
factorial (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Int -> Integer
forall a. Integral a => a -> Integer
factorial Int
k)
                ]

  where
    monom :: [Int] -> Int -> Omega QPow
monom [Int]
xs Int
k = [Int] -> QPow -> Omega QPow
forall ab. [Int] -> ab -> Omega ab
Omega [Int]
xs (Int -> QPow
QPow Int
k)
    us :: [Omega QPow]
us = [ [Int] -> Int -> Omega QPow
monom [Int
i] Int
0 | Int
i<-[Int
1..Int
n] ]

-- | Given a function specifying what to substitute in the place of @q^k@, we do the substitution.
umbralSubstQPow :: (ChernBase base) => (QPow -> ZMod base) -> ZMod (Omega QPow) -> ZMod (Omega base)
umbralSubstQPow :: (QPow -> ZMod base) -> ZMod (Omega QPow) -> ZMod (Omega base)
umbralSubstQPow QPow -> ZMod base
subst1 ZMod (Omega QPow)
input = [ZMod (Omega base)] -> ZMod (Omega base)
forall b c. (Ord b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b
ZMod.sum 
  [ [(Omega base, Integer)] -> ZMod (Omega base)
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList 
      [ ([Int] -> base -> Omega base
forall ab. [Int] -> ab -> Omega ab
Omega [Int]
us base
ab , Integer
cInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
coeff) 
      | (base
ab,Integer
c) <- ZMod base -> [(base, Integer)]
forall c b. FreeMod c b -> [(b, c)]
ZMod.toList (QPow -> ZMod base
subst1 QPow
qpow) 
      ] 
  | (Omega [Int]
us QPow
qpow , Integer
coeff) <- ZMod (Omega QPow) -> [(Omega QPow, Integer)]
forall c b. FreeMod c b -> [(b, c)]
ZMod.toList ZMod (Omega QPow)
input  
  ]

--------------------------------------------------------------------------------

-- | It is not hard to prove (by considering the pushforward along
-- the map forgetting one of the points), that the CSM of the locus
-- @U^n@ of the distinct points has the following form (for @n>=3@):
--
-- > csm(U^n) = sum_{k=0}^{n-3} \frac{(n-3)!}{k!} (-1)^{n-3-k} \sigma_{n-3-k}(u) Q_k(a,b)
-- 
-- We can already compute all CSM-s recursively, and from that information we can
-- determine these polynomials.
--
-- Which then we can compare with the recursive formula for the
-- polynomials itself (which is /much/ faster to evaluate)
--
computeQPolys :: Int -> ZMod AB
computeQPolys :: Int -> ZMod AB
computeQPolys = ZMod AB -> Int -> (Int -> ZMod AB) -> Int -> ZMod AB
forall a. a -> Int -> (Int -> a) -> Int -> a
icache' ZMod AB
forall c b. FreeMod c b
ZMod.zero (-Int
3) Int -> ZMod AB
calcComputeQPolys where

  calcComputeQPolys :: Int -> ZMod AB
  calcComputeQPolys :: Int -> ZMod AB
calcComputeQPolys Int
n 
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  -Int
3    = String -> ZMod AB
forall a. HasCallStack => String -> a
error String
"computeQPolys: n >= -3 is required"
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
3    = ZMod AB
forall b c. (Monoid b, Num c) => FreeMod c b
ZMod.one
    | Bool
otherwise  = (Omega AB -> AB) -> FreeMod Integer (Omega AB) -> ZMod AB
forall a b c.
(Ord a, Ord b, Eq c, Num c) =>
(a -> b) -> FreeMod c a -> FreeMod c b
ZMod.mapBase Omega AB -> AB
project FreeMod Integer (Omega AB)
almost
    where

      almost :: FreeMod Integer (Omega AB)
almost = FreeMod Integer (Omega AB)
open FreeMod Integer (Omega AB)
-> FreeMod Integer (Omega AB) -> FreeMod Integer (Omega AB)
forall a. Num a => a -> a -> a
- FreeMod Integer (Omega AB)
smaller
      open :: FreeMod Integer (Omega AB)
open   = Int -> FreeMod Integer (Omega AB)
forall base. ChernBase base => Int -> ZMod (Omega base)
computeOpenStratumCSM (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3)     -- we should use this as the basis of the computation, unfortunately it's rather slow
    
      umbSmaller :: ZMod (Omega QPow)
umbSmaller = Int -> ZMod (Omega QPow)
umbralDistinctFormula (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
3) ZMod (Omega QPow) -> ZMod (Omega QPow) -> ZMod (Omega QPow)
forall a. Num a => a -> a -> a
- ZMod (Omega QPow)
umbHighest
      umbHighest :: ZMod (Omega QPow)
umbHighest = Omega QPow -> ZMod (Omega QPow)
forall c b. Num c => b -> FreeMod c b
ZMod.generator ([Int] -> QPow -> Omega QPow
forall ab. [Int] -> ab -> Omega ab
Omega [] (Int -> QPow
QPow Int
n))        -- q^n
      smaller :: FreeMod Integer (Omega AB)
smaller     = (QPow -> ZMod AB)
-> ZMod (Omega QPow) -> FreeMod Integer (Omega AB)
forall base.
ChernBase base =>
(QPow -> ZMod base) -> ZMod (Omega QPow) -> ZMod (Omega base)
umbralSubstQPow (\(QPow Int
k) -> Int -> ZMod AB
computeQPolys Int
k) ZMod (Omega QPow)
umbSmaller

{-
      smaller = ZMod.sum 
        [ ZMod.scale coeff $ 
            (symPoly (n-k) us) * (embed $ computeQPolys k)
        | k<-[0..n-1]
        , let coeff = negateIfOdd (n+k) (factorial n `div` factorial k)
        ]
      us = [ Omega [i] (AB 0 0) | i<-[1..n+3] ]
      embed = ZMod.mapBase $ \ab -> Omega [] ab
-}

      project :: Omega AB -> AB
project (Omega [Int]
us AB
ab) = case [Int]
us of
        [] -> AB
ab
        [Int]
_  -> String -> AB
forall a. HasCallStack => String -> a
error (String -> AB) -> String -> AB
forall a b. (a -> b) -> a -> b
$ String
"computeQPolys: cannot project u terms:\n  " String -> ShowS
forall a. [a] -> [a] -> [a]
++ FreeMod Integer (Omega AB) -> String
forall a. Pretty a => a -> String
pretty FreeMod Integer (Omega AB)
almost

--------------------------------------------------------------------------------
-- * The recursive formula for the @Q_k(a,b)@ polynomials

-- | The Fibonacci-type recursive formula for the @Q_k(a,b)@ polynomials
--
-- > Q_{-3} = 1
-- > Q_k    = Q_{k-1} * (1 - (k+1)*(a+b)) - Q_{k-2} * a*b * (k-1)*(k+2)
-- >        = Q_{k-1} * (1 - (k+1)* c_1 ) - Q_{k-2} * c_2 * (k-1)*(k+2)
--
-- We provide both the Chern root and the Chern class version in a uniform
-- way for convenience.
formulaQPoly :: ChernBase base => Int -> ZMod base
formulaQPoly :: Int -> ZMod base
formulaQPoly Int
n = (ZMod AB, FreeMod Integer Chern) -> ChernBase base => ZMod base
forall (f :: * -> *) base.
(f AB, f Chern) -> ChernBase base => f base
select1 
  ( Int -> ZMod AB
formulaQPolyAB   Int
n 
  , Int -> FreeMod Integer Chern
formulaQPolyChern Int
n
  )

formulaQPolyAB :: Int -> ZMod AB
formulaQPolyAB :: Int -> ZMod AB
formulaQPolyAB = ZMod AB -> Int -> (Int -> ZMod AB) -> Int -> ZMod AB
forall a. a -> Int -> (Int -> a) -> Int -> a
icache' ZMod AB
forall c b. FreeMod c b
ZMod.zero (-Int
3) Int -> ZMod AB
calcQPoly where
  
  calcQPoly :: Int -> ZMod AB
  calcQPoly :: Int -> ZMod AB
calcQPoly Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  -Int
3   = ZMod AB
forall c b. FreeMod c b
ZMod.zero
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
3   = Integer -> ZMod AB
forall b c. (Monoid b, Eq c, Num c) => c -> FreeMod c b
ZMod.konst Integer
1
    | Bool
otherwise = ZMod AB
mult1 ZMod AB -> ZMod AB -> ZMod AB
forall a. Num a => a -> a -> a
* ZMod AB
prev1 ZMod AB -> ZMod AB -> ZMod AB
forall a. Num a => a -> a -> a
+ ZMod AB
mult2 ZMod AB -> ZMod AB -> ZMod AB
forall a. Num a => a -> a -> a
* ZMod AB
prev2
    where
      prev1 :: ZMod AB
prev1 = Int -> ZMod AB
formulaQPolyAB (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
      prev2 :: ZMod AB
prev2 = Int -> ZMod AB
formulaQPolyAB (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2)

      Pair ZMod AB
mult1 ZMod AB
mult2 = Int -> Pair (ZMod AB)
forall base. ChernBase base => Int -> Pair (ZMod base)
qpolyRecursionCoeffs Int
n

-- | Chern class version of the @Q_k@ formula (should be faster then the Chern root version, because the are less terms).
formulaQPolyChern :: Int -> ZMod Chern
formulaQPolyChern :: Int -> FreeMod Integer Chern
formulaQPolyChern = FreeMod Integer Chern
-> Int
-> (Int -> FreeMod Integer Chern)
-> Int
-> FreeMod Integer Chern
forall a. a -> Int -> (Int -> a) -> Int -> a
icache' FreeMod Integer Chern
forall c b. FreeMod c b
ZMod.zero (-Int
3) Int -> FreeMod Integer Chern
calcQPoly where
  
  calcQPoly :: Int -> ZMod Chern
  calcQPoly :: Int -> FreeMod Integer Chern
calcQPoly Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  -Int
3   = FreeMod Integer Chern
forall c b. FreeMod c b
ZMod.zero
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
3   = Integer -> FreeMod Integer Chern
forall b c. (Monoid b, Eq c, Num c) => c -> FreeMod c b
ZMod.konst Integer
1
    | Bool
otherwise = FreeMod Integer Chern
mult1 FreeMod Integer Chern
-> FreeMod Integer Chern -> FreeMod Integer Chern
forall a. Num a => a -> a -> a
* FreeMod Integer Chern
prev1 FreeMod Integer Chern
-> FreeMod Integer Chern -> FreeMod Integer Chern
forall a. Num a => a -> a -> a
+ FreeMod Integer Chern
mult2 FreeMod Integer Chern
-> FreeMod Integer Chern -> FreeMod Integer Chern
forall a. Num a => a -> a -> a
* FreeMod Integer Chern
prev2
    where
      nn :: Integer
nn = Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n :: Integer

      prev1 :: FreeMod Integer Chern
prev1 = Int -> FreeMod Integer Chern
formulaQPolyChern (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
      prev2 :: FreeMod Integer Chern
prev2 = Int -> FreeMod Integer Chern
formulaQPolyChern (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2)

      Pair FreeMod Integer Chern
mult1 FreeMod Integer Chern
mult2 = Int -> Pair (FreeMod Integer Chern)
forall base. ChernBase base => Int -> Pair (ZMod base)
qpolyRecursionCoeffs Int
n

qpolyRecursionCoeffs :: ChernBase base => Int -> Pair (ZMod base)
qpolyRecursionCoeffs :: Int -> Pair (ZMod base)
qpolyRecursionCoeffs Int
n = (Pair (ZMod AB), Pair (FreeMod Integer Chern))
-> ChernBase base => Pair (ZMod base)
forall (f :: * -> *) (g :: * -> *) base.
(f (g AB), f (g Chern)) -> ChernBase base => f (g base)
select2 
  (  ZMod AB -> ZMod AB -> Pair (ZMod AB)
forall a. a -> a -> Pair a
Pair  ZMod AB
mult1_AB    ZMod AB
mult2_AB 
  ,  FreeMod Integer Chern
-> FreeMod Integer Chern -> Pair (FreeMod Integer Chern)
forall a. a -> a -> Pair a
Pair  FreeMod Integer Chern
mult1_Chern FreeMod Integer Chern
mult2_Chern
  )
  where

    mult1_AB :: ZMod AB
mult1_AB = [(AB, Integer)] -> ZMod AB
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList 
      [ ( Int -> Int -> AB
AB Int
0 Int
0 ,     Integer
1 )
      , ( Int -> Int -> AB
AB Int
1 Int
0 , -Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1 )
      , ( Int -> Int -> AB
AB Int
0 Int
1 , -Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1 )
      ]
    mult2_AB :: ZMod AB
mult2_AB = AB -> Integer -> ZMod AB
forall b c. (Ord b, Num c, Eq c) => b -> c -> FreeMod c b
ZMod.singleton (Int -> Int -> AB
AB Int
1 Int
1) (-(Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
2)) 
  
    mult1_Chern :: FreeMod Integer Chern
mult1_Chern = [(Chern, Integer)] -> FreeMod Integer Chern
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList 
      [ ( Int -> Int -> Chern
Chern Int
0 Int
0 ,     Integer
1 )
      , ( Int -> Int -> Chern
Chern Int
1 Int
0 , -Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1 )
      ]
    mult2_Chern :: FreeMod Integer Chern
mult2_Chern = Chern -> Integer -> FreeMod Integer Chern
forall b c. (Ord b, Num c, Eq c) => b -> c -> FreeMod c b
ZMod.singleton (Int -> Int -> Chern
Chern Int
0 Int
1) (-(Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
nnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
2))

    nn :: Integer
nn = Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n :: Integer

--------------------------------------------------------------------------------
-- small @Q_k@ polynomials

{-
polyZMod :: ZMod AB -> (forall base. ChernBase base => ZMod base)
polyZMod ab = select1 (ab, abToChern ab)

-- | @Q_0 = ( 1 - a + b) ( 1 + a - b) = 1 - a^2 - b^2 + 2ab = 1 - c_1^2 + 4c_2@
konstQ0 :: ChernBase base => ZMod base
konstQ0 = polyZMod q0 where 
  q0 = ZMod.fromList [ ( AB 0 0 ,  1 )  , ( AB 2 0 , -1 )  , ( AB 0 2 , -1 )  , ( AB 1 1 ,  2 )  ]  

-- | @Q_-1 = 1 + a + b + 2 a*b = 1 + c_1 + 2c_2@
konstQminus1 :: ChernBase base => ZMod base
konstQminus1 = polyZMod qminus1 where
  qminus1 = ZMod.fromList [ ( AB 0 0 ,  1 ) ,  ( AB 1 0 ,  1 )  , ( AB 0 1 ,  1 )  , ( AB 1 1 ,  2 ) ]

-- | @Q_-2 = 1 + a + b = 1 + c_1@
konstQminus2 :: ChernBase base => ZMod base
konstQminus2 = polyZMod qminus2 where
  qminus2 = ZMod.fromList [ ( AB 0 0 ,  1 ) , ( AB 1 0 ,  1 ) , ( AB 0 1 ,  1 ) ]

-- | @Q_-3 = 1@
konstQminus3 :: ChernBase base => ZMod base
konstQminus3 = ZMod.konst 1
-}

--------------------------------------------------------------------------------
-- * Formula for the CSM class of the stratum of distinct points

-- | The formula for the CSM of the set of distinct ordered points
-- using the formula for the Q_k(a,b) polynomials above
--
formulaDistinctCSM :: ChernBase base => Int -> ZMod (Omega base)
formulaDistinctCSM :: Int -> ZMod (Omega base)
formulaDistinctCSM Int
n 
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0     = String -> ZMod (Omega base)
forall a. HasCallStack => String -> a
error String
"formulaDistinctCSM: dimension should be nonnegative"
  | Bool
otherwise = (QPow -> ZMod base) -> ZMod (Omega QPow) -> ZMod (Omega base)
forall base.
ChernBase base =>
(QPow -> ZMod base) -> ZMod (Omega QPow) -> ZMod (Omega base)
umbralSubstQPow QPow -> ZMod base
forall base. ChernBase base => QPow -> ZMod base
fun 
              (ZMod (Omega QPow) -> ZMod (Omega base))
-> ZMod (Omega QPow) -> ZMod (Omega base)
forall a b. (a -> b) -> a -> b
$ Int -> ZMod (Omega QPow)
umbralDistinctFormula Int
n
  where
    fun :: QPow -> ZMod base
fun (QPow Int
k) = Int -> ZMod base
forall base. ChernBase base => Int -> ZMod base
formulaQPoly Int
k
{-
  | n < 3     = computeOpenStratumCSM n
  | otherwise = ZMod.sum 
      [ ZMod.scale coeff poly
      | k <- [0..n-3] 
      , let coeff = paritySignValue (n-3-k) * div (factorial (n-3)) (factorial k)
      , let qk    = formulaQPoly k
      , let sym   = choose (n-3-k) [1..n]
      , let poly  = ZMod.fromList [ (Omega xs ab, k) | xs <- sym, (ab,k) <- ZMod.toList qk ]
      ]
-}

-- | Just the pushforward of the previous along @Delta_mu@
formulaAnyStratumCSM :: ChernBase base => SetPartition -> ZMod (Omega base)
formulaAnyStratumCSM :: SetPartition -> ZMod (Omega base)
formulaAnyStratumCSM SetPartition
setp = FreeMod Integer (Eta base) -> ZMod (Omega base)
forall coeff ab.
(Eq coeff, Num coeff, Ord ab) =>
FreeMod coeff (Eta ab) -> FreeMod coeff (Omega ab)
unsafeEtaToOmega (FreeMod Integer (Eta base) -> ZMod (Omega base))
-> FreeMod Integer (Eta base) -> ZMod (Omega base)
forall a b. (a -> b) -> a -> b
$ SetPartition -> ZMod (Omega base) -> FreeMod Integer (Eta base)
forall base.
ChernBase base =>
SetPartition -> ZMod (Omega base) -> ZMod (Eta base)
delta_star SetPartition
setp (Int -> ZMod (Omega base)
forall base. ChernBase base => Int -> ZMod (Omega base)
formulaDistinctCSM Int
k) where
  k :: Int
k = SetPartition -> Int
forall a. HasNumberOfParts a => a -> Int
numberOfParts SetPartition
setp
  
--------------------------------------------------------------------------------