packages feed

coincident-root-loci (empty) → 0.2

raw patch · 37 files changed

+4354/−0 lines, 37 filesdep +arraydep +basedep +coincident-root-locisetup-changed

Dependencies added: array, base, coincident-root-loci, combinat, containers, random, tasty, tasty-hunit, transformers

Files

+ LICENSE view
@@ -0,0 +1,29 @@+Copyright (c) 2015-2017, Balazs Komuves+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++- Redistributions of source code must retain the above copyright notice,+this list of conditions and the following disclaimer.+ +- Redistributions in binary form must reproduce the above copyright notice,+this list of conditions and the following disclaimer in the documentation+and/or other materials provided with the distribution.+ +- Neither names of the copyright holders nor the names of the contributors+may be used to endorse or promote products derived from this software without+specific prior written permission. ++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ coincident-root-loci.cabal view
@@ -0,0 +1,100 @@+Name:                coincident-root-loci+Version:             0.2+Synopsis:            Equivariant CSM classes of coincident root loci++Description:         This library contians a set of function to compute, among+                     others, the @GL(2)@-equivariant Chern-Schwartz-MacPherson+                     classes of coincident root loci, which are subvarieties+                     of the space of unordered @n@-tuples of points in the complex+                     projective line. To such an @n@-tuples we can associate +                     a partition of @n@ given by the multiplicities of the distinct+                     points; this stratifies the set of all @n@-tuples, and we+                     call these strata \"coincident root loci\".++                     This package is supplementary software for a forthcoming paper.++License:             BSD3+License-file:        LICENSE+Author:              Balazs Komuves+Copyright:           (c) 2015-2017 Balazs Komuves+Maintainer:          bkomuves (plus) hackage (at) gmail (dot) com+Homepage:            http://code.haskell.org/~bkomuves/+Stability:           Experimental+Category:            Math+Tested-With:         GHC == 8.0.2+Cabal-Version:       >= 1.18+Build-Type:          Simple++--------------------------------------------------------------------------------++Library++  Build-Depends:       base >= 4 && < 5, +                       array >= 0.5, containers, random, transformers,+                       combinat >= 0.2.8.2++  Exposed-Modules:     +                       -- Math.RootLoci+                       Math.RootLoci.Classic+                       -- Math.RootLoci.Dual+                       Math.RootLoci.Dual.Restriction+                       Math.RootLoci.Dual.Localization+                       -- Math.RootLoci.CSM+                       -- Math.RootLoci.CSM.Equivariant+                       Math.RootLoci.CSM.Equivariant.Direct+                       Math.RootLoci.CSM.Equivariant.Recursive+                       Math.RootLoci.CSM.Equivariant.Ordered+                       Math.RootLoci.CSM.Equivariant.PushForward+                       Math.RootLoci.CSM.Equivariant.Umbral+                       Math.RootLoci.CSM.Aluffi+                       Math.RootLoci.CSM.Projective+                       Math.RootLoci.Geometry+                       Math.RootLoci.Geometry.Forget+                       Math.RootLoci.Geometry.Cohomology+                       Math.RootLoci.Geometry.Mobius+                       -- Math.RootLoci.Applications+                       -- Math.RootLoci.Applications.FlexLines+                       Math.RootLoci.Algebra+                       Math.RootLoci.Algebra.FreeMod+                       Math.RootLoci.Algebra.Polynomial+                       Math.RootLoci.Algebra.SymmPoly+                       Math.RootLoci.Misc+                       Math.RootLoci.Misc.Pretty+                       Math.RootLoci.Misc.PTable+                       Math.RootLoci.Misc.Common++  Default-Extensions:  CPP, BangPatterns+  Other-Extensions:    MultiParamTypeClasses, ScopedTypeVariables, +                       GeneralizedNewtypeDeriving++  Default-Language:    Haskell2010++  Hs-Source-Dirs:      src++  ghc-options:         -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports+    +--------------------------------------------------------------------------------+    +test-suite test++  default-language:    Haskell2010+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             testSuite.hs+  +  build-depends:       base >= 4 && < 5, containers >= 0.4, array >= 0.5,+                       tasty >= 0.11, tasty-hunit >= 0.9,+                       combinat >= 0.2.8.2,+                       coincident-root-loci >= 0.2+                       +  other-modules:       Tests.Common+                       Tests.Dual+                       Tests.Pushforward                  +                       Tests.CSM.Equivariant                  +                       Tests.CSM.Projective                  +                       Tests.RootVsClass.Check                  +                       Tests.RootVsClass.Direct                  +                       Tests.RootVsClass.Ordered                  +                       Tests.RootVsClass.PushForward                  +                       Tests.RootVsClass.Recursive                  +                       Tests.RootVsClass.Umbral
+ src/Math/RootLoci/Algebra.hs view
@@ -0,0 +1,23 @@++-- | Re-exporting the Algebra.* modules.+-- +-- Because of limitations of the import-export mechanism, you still have+-- to do+--+-- > import Math.RootLoci.Algebra+-- > import qualified Math.RootLoci.Algebra.FreeMod as ZMod+--++module Math.RootLoci.Algebra+  ( ZMod , QMod , FreeMod +  , module Math.RootLoci.Algebra.Polynomial+  , module Math.RootLoci.Algebra.SymmPoly+--  , module ZMod                               -- apparently this does not work+  )+  where++import Math.RootLoci.Algebra.FreeMod ( ZMod , QMod , FreeMod )+import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra.Polynomial+import Math.RootLoci.Algebra.SymmPoly
+ src/Math/RootLoci/Algebra/FreeMod.hs view
@@ -0,0 +1,214 @@++-- | Free modules. +--+-- This module should be imported qualified++{-# LANGUAGE BangPatterns, FlexibleInstances, TypeSynonymInstances #-}+module Math.RootLoci.Algebra.FreeMod where++--------------------------------------------------------------------------------++import Prelude   hiding ( sum , product )+import Data.List hiding ( sum , product )++import Data.Monoid+import Data.Ratio+import Data.Maybe++import Math.Combinat.Sets ( choose )++import qualified Data.Map.Strict as Map+import Data.Map.Strict (Map)++--------------------------------------------------------------------------------++-- | Free module over a coefficient ring with the given base. Internally a map+-- storing the coefficients. We maintain the invariant that the coefficients+-- are never zero.+newtype FreeMod coeff base = FreeMod { unFreeMod :: Map base coeff } deriving (Eq,Show)++-- | Free module with integer coefficients+type ZMod base = FreeMod Integer  base++-- | Free module with rational coefficients+type QMod base = FreeMod Rational base++--------------------------------------------------------------------------------++instance (Monoid b, Ord b, Eq c, Num c) => Num (FreeMod c b) where+  (+)    = add+  (-)    = sub+  negate = neg+  (*)    = mul+  fromInteger = konst . fromInteger+  abs    = error "FreeMod/abs"+  signum = error "FreeMod/signum"++--------------------------------------------------------------------------------+-- * Sanity checking++-- | Should be the identity function+normalize :: (Ord b, Eq c, Num c) => FreeMod c b -> FreeMod c b+normalize = FreeMod . Map.filter (/=0) . unFreeMod++-- | Safe equality testing (should be identical to @==@)+safeEq :: (Ord b, Eq b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> Bool+safeEq x y = normalize x == normalize y++--------------------------------------------------------------------------------+-- * Constructing and deconstructing++-- | The additive unit+zero :: FreeMod c b+zero = FreeMod $ Map.empty++-- | A module generator+generator :: Num c => b -> FreeMod c b +generator x = FreeMod $ Map.singleton x 1++-- | A single generator with a coefficient+singleton :: (Ord b) => b -> c -> FreeMod c b+singleton b c = FreeMod $ Map.singleton b c++-- | Conversion from list. +-- Note that we assume here that each generator appears at most once!+fromList :: (Eq c, Num c, Ord b) => [(b,c)] -> FreeMod c b+fromList = FreeMod . Map.fromList . filter cond where+  cond (b,x) = (x/=0)++-- | Conversion to list +toList :: FreeMod c b -> [(b,c)]+toList = Map.toList . unFreeMod++-- | Extract the coefficient of a generator+coeffOf :: (Ord b, Num c) => b -> FreeMod c b -> c+coeffOf b (FreeMod x) = case Map.lookup b x of+  Just c  -> c+  Nothing -> 0++-- | Finds the term with the largest generator (in the natural ordering of the generators)+findMaxTerm :: (Ord b) => FreeMod c b -> Maybe (b,c)+findMaxTerm (FreeMod m) = if Map.null m+  then Nothing+  else Just (Map.findMax m)++-- | Finds the term with the smallest generator (in the natural ordering of the generators)+findMinTerm :: (Ord b) => FreeMod c b -> Maybe (b,c)+findMinTerm (FreeMod m) = if Map.null m+  then Nothing+  else Just (Map.findMin m)++--------------------------------------------------------------------------------+-- * Basic operations++-- | Negation+neg :: Num c => FreeMod c b -> FreeMod c b +neg (FreeMod m) = FreeMod (Map.map negate m)++-- | Additions+add :: (Ord b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+add (FreeMod m1) (FreeMod m2) = FreeMod (Map.mergeWithKey f id id m1 m2) where+  f _ x y = case x+y of { 0 -> Nothing ; z -> Just z }++-- | Subtraction+sub :: (Ord b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+sub (FreeMod m1) (FreeMod m2) = FreeMod (Map.mergeWithKey f id (Map.map negate) m1 m2) where+  f _ x y = case x-y of { 0 -> Nothing ; z -> Just z }++-- | Scaling by a number+scale :: (Ord b, Eq c, Num c) => c -> FreeMod c b -> FreeMod c b+scale 0 _           = zero+scale x (FreeMod m) = FreeMod (Map.mapMaybe f m) where+  f y = case x*y of { 0 -> Nothing ; z -> Just z }++-- | Dividing by a number (assuming that the coefficient ring is integral, and each coefficient+-- is divisible by the given number)+invScale :: (Ord b, Eq c, Integral c, Show c) => c -> FreeMod c b -> FreeMod c b+invScale d (FreeMod m) = FreeMod (Map.mapMaybe f m) where+  f a = case divMod a d of+    (b,0) -> case b of { 0 -> Nothing ; z -> Just z }+    _     -> error $ "FreeMod/invScale: not divisible by " ++ show d++--------------------------------------------------------------------------------++-- | Summation+sum :: (Ord b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b+sum []  = zero+sum zms = (foldl1' add) zms++-- | Linear combination+linComb :: (Ord b, Eq c, Num c) => [(c, FreeMod c b)] -> FreeMod c b+linComb = sumWith where++   sumWith :: (Ord b, Eq c, Num c) => [(c, FreeMod c b)] -> FreeMod c b+   sumWith []  = zero+   sumWith zms = sum [ scale c zm | (c,zm) <- zms ]++-- | Expand each generator into a term in another module and then sum the results+flatMap :: (Ord b1, Ord b2, Eq c, Num c) => (b1 -> FreeMod c b2) -> FreeMod c b1 -> FreeMod c b2+flatMap f = sum . map g . Map.assocs . unFreeMod where+  g (x,c) = scale c (f x)++flatMap' :: (Ord b1, Ord b2, Eq c2, Num c2) => (c1 -> c2) -> (b1 -> FreeMod c2 b2) -> FreeMod c1 b1 -> FreeMod c2 b2+flatMap' embed f = sum . map g . Map.assocs . unFreeMod where+  g (x,c) = scale (embed c) (f x)++-- | The histogram of a multiset of generators is naturally an element of the given Z-module.+{-# SPECIALIZE histogram :: Ord b => [b] -> ZMod b #-} +histogram :: (Ord b, Num c) => [b] -> FreeMod c b+histogram xs = FreeMod $ foldl' f Map.empty xs where+  f old x = Map.insertWith (+) x 1 old+  +--------------------------------------------------------------------------------+-- * Rings++-- | The multiplicative unit+one :: (Monoid b, Num c) => FreeMod c b+one = konst 1++-- | A constant+konst :: (Monoid b) => c -> FreeMod c b+konst c = FreeMod (Map.singleton mempty c)++-- | Multiplying two ring elements+mul :: (Ord b, Monoid b, Eq c, Num c) => FreeMod c b -> FreeMod c b -> FreeMod c b+mul xx yy = sum [ (f x c) | (x,c) <- toList xx ] where+  f x c = FreeMod $ Map.fromList [ (x<>y, cd) | (y,d) <- ylist , let cd = c*d , cd /= 0 ]+  ylist = toList yy++-- | Product+product :: (Ord b, Monoid b, Eq c, Num c) => [FreeMod c b] -> FreeMod c b+product []  = generator mempty+product xs  = foldl1' mul xs++-- | Multiplies by a monomial+mulMonom :: (Ord b, Monoid b) => b -> FreeMod c b -> FreeMod c b+mulMonom monom = FreeMod . Map.mapKeys (mappend monom) . unFreeMod++--------------------------------------------------------------------------------+-- * Misc++-- | A symmetric polynomial of some generators+symPoly :: (Ord a, Monoid a) => Int -> [a] -> ZMod a+symPoly k xs = fromList $ map (\x -> (x,1)) $ (map mconcat $ choose k xs) ++-- | Changing the base set+mapBase :: (Ord a, Ord b) => (a -> b) -> FreeMod c a -> FreeMod c b+mapBase f = onFreeMod (Map.mapKeys f)++-- | Changing the coefficient ring+mapCoeff :: (Ord b) => (c1 -> c2) -> FreeMod c1 b -> FreeMod c2 b+mapCoeff f = onFreeMod' (Map.map f)++-- | Extract a subset of terms+filterBase :: (Ord a, Ord b) => (a -> Maybe b) -> FreeMod c a -> FreeMod c b+filterBase f = onFreeMod (Map.fromList . mapMaybe g . Map.toList) where+  g (k,x) = case f k of { Just k' -> Just (k',x) ; Nothing -> Nothing }++onFreeMod :: (Ord a, Ord b) => (Map a c -> Map b c) -> FreeMod c a -> FreeMod c b+onFreeMod f = FreeMod . f . unFreeMod++onFreeMod' :: (Ord a, Ord b) => (Map a c -> Map b d) -> FreeMod c a -> FreeMod d b+onFreeMod' f = FreeMod . f . unFreeMod++--------------------------------------------------------------------------------
+ src/Math/RootLoci/Algebra/Polynomial.hs view
@@ -0,0 +1,102 @@++-- | Univariate polynomials++{-# LANGUAGE GeneralizedNewtypeDeriving #-}+module Math.RootLoci.Algebra.Polynomial where++--------------------------------------------------------------------------------++import Data.Array ( assocs ) ++import Math.Combinat.Numbers++import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod+import Math.RootLoci.Algebra.FreeMod ( FreeMod , ZMod , QMod )++--------------------------------------------------------------------------------+-- * Polynomials++-- | Standard univariate polynomials+newtype Poly coeff = Poly { fromPoly :: FreeMod coeff X } deriving (Eq,Num,Show)++-- | Univariate polynomials using /rising factorials/ as a basis function+newtype RisingPoly  coeff = RisingPoly  { fromRisingPoly  :: FreeMod coeff RisingF }  deriving (Eq,Show)++-- | Univariate polynomials using /falling factorials/ as a basis function+newtype FallingPoly coeff = FallingPoly { fromFallingPoly :: FreeMod coeff FallingF } deriving (Eq,Show)++instance (Num c, Show c, Eq c, IsSigned c) => Pretty (Poly        c) where pretty (Poly        p) = pretty p +instance (Num c, Show c, Eq c, IsSigned c) => Pretty (RisingPoly  c) where pretty (RisingPoly  p) = pretty p +instance (Num c, Show c, Eq c, IsSigned c) => Pretty (FallingPoly c) where pretty (FallingPoly p) = pretty p ++--------------------------------------------------------------------------------+-- * Monomials ++-- | A power of @x@ (that is, a monomial of the form @x^i@)+newtype X = X Int deriving (Eq,Ord,Show)++instance Monoid X where+  mempty = X 0+  mappend (X e) (X f) = X (e+f)++instance Pretty X where+  pretty (X e) = case e of+    0 -> "1"+    1 -> "x"+    _ -> "x^" ++ show e++--------------------------------------------------------------------------------+-- * Rising and falling factorials ++-- | Rising factorial @x^(k) = x(x+1)(x+2)...(x+k-1)@+newtype RisingF = RF Int deriving (Eq,Ord,Show)++-- | Falling factorial @x_(k) = x(x-1)(x-2)...(x-k+1)@+newtype FallingF = FF Int deriving (Eq,Ord,Show)++instance Pretty RisingF where+  pretty (RF k) = case k of+    0 -> "1"+    1 -> "x"+    _ -> "x^(" ++ show k ++ ")"++instance Pretty FallingF where+  pretty (FF k) = case k of+    0 -> "1"+    1 -> "x"+    _ -> "x_(" ++ show k ++ ")"++risingPoly :: RisingF -> Poly Integer+risingPoly (RF k) = Poly $ ZMod.fromList+  [ (X p, abs c) | (p,c) <- assocs (signedStirling1stArray k) ]++fallingPoly :: FallingF -> Poly Integer+fallingPoly (FF k) = Poly $ ZMod.fromList+  [ (X p,     c) | (p,c) <- assocs (signedStirling1stArray k) ]++--------------------------------------------------------------------------------+-- * Lagrange interpolation++lagrangeInterp :: [(Rational,Rational)] -> Poly Rational+lagrangeInterp = Poly . lagrangeInterp'++lagrangeInterp' :: [(Rational,Rational)] -> QMod X+lagrangeInterp' xys = final where+  final = ZMod.sum [ ZMod.scale (ys!!j) (lagrangePoly' xs j) | j<-[0..m-1] ]  where+  m = length xys+  (xs,ys) = unzip xys++lagrangePoly' :: [Rational] -> Int -> QMod X+lagrangePoly' xs j = ZMod.scale (1/denom) numer where+  numer  = ZMod.product [ term i    | i<-[0..m-1] , i /= j ]+  denom  = product      [ x j - x i | i<-[0..m-1] , i /= j ]+  m      = length xs+  x i    = xs !! i+  term i = ZMod.fromList +    [ (X 1 ,     1 )+    , (X 0 , - x i )+    ]++--------------------------------------------------------------------------------
+ src/Math/RootLoci/Algebra/SymmPoly.hs view
@@ -0,0 +1,366 @@++-- | Symmetric polynomials in two variables @alpha@ and @beta@.+--+-- We provide three representation:+--+-- * symmetric polynomials in @alpha@ and @beta@ (Chern roots)+--+-- * polynomials in the elementary symmetric polynomials @c1=alpha+beta@ and @c2=alpha*beta@ (Chern classes)+--+-- * Schur polynomials @s[i,j]@+--+-- The monomials of the first two of these form monoids (the product of +-- monomials is again a monomial), and can be used uniformly with the+-- help of some type-level hackery.+--+-- How to use the unified interface?+-- Suppose you have a function like this:+--+-- > tau :: ChernBase base => Int -> ZMod base+--+-- When calling it, you want to specify the output type (either @ZMod AB@ or @ZMod Chern@).+-- You can do that three ways:+--+-- > x = tau @AB 10                  -- this needs -XTypeApplications+-- > x = (tau 10 :: ZMod AB)+-- > x = spec1' ChernRoot $ tau 10+--+-- The first one is the most convenient, but it only works with GHC 8 and later.+-- The other two work with older GHC versions, too.+--+++{-# LANGUAGE DataKinds, TypeFamilies, Rank2Types, GADTs, StandaloneDeriving #-}+module Math.RootLoci.Algebra.SymmPoly where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Sign+import Math.Combinat.Numbers++import qualified Data.Map.Strict as Map++import Control.Monad+import System.Random++import Math.RootLoci.Algebra.FreeMod (ZMod)+import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Misc.Pretty++import Unsafe.Coerce as Unsafe++--------------------------------------------------------------------------------+-- * Base monomials++-- | Chern roots: @alpha^i * beta^j@, monomial base of @Z[alpha,beta]@+data AB = AB !Int !Int deriving (Eq,Ord,Show)++-- | Chern classes: @c1^i * c2^j@, monomial base of @Z[c1,c2]@+data Chern = Chern !Int !Int deriving (Eq,Ord,Show)++-- | Schur basis function: @S[i,j]@+data Schur = Schur !Int !Int deriving (Eq,Ord,Show) ++alpha, beta :: AB+alpha = AB 1 0 +beta  = AB 0 1    ++--------------------------------------------------------------------------------++-- | @alpha * beta = c2@+alphaBeta :: AB+alphaBeta = AB 1 1    ++-- | @c1 = alpha + beta@+c1 :: Chern+c1 = Chern 1 0     ++-- | @c2 = alpha * beta@+c2 :: Chern+c2 = Chern 0 1     ++--------------------------------------------------------------------------------+-- * Unified interface++-- | A singleton for distinguishing the two cases +data Sing base where+  ChernRoot  :: Sing AB+  ChernClass :: Sing Chern++deriving instance Eq  (Sing base)+deriving instance Ord (Sing base)++-- | Common interface to work with Chern classes and Chern roots uniformly+class (Eq base, Ord base, Monoid base, Graded base, Pretty base) => ChernBase base where+  chernTag  :: base       -> Sing base+  chernTag1 :: f base     -> Sing base+  chernTag2 :: f (g base)     -> Sing base+  chernTag3 :: f (g (h base)) -> Sing base+  fromAB    :: ZMod AB    -> ZMod base  +  fromChern :: ZMod Chern -> ZMod base  +  fromSchur :: ZMod Schur -> ZMod base+  toAB      :: ZMod base  -> ZMod AB  +  toChern   :: ZMod base  -> ZMod Chern+  toSchur   :: ZMod base  -> ZMod Schur++instance ChernBase AB where+  chernTag  _ = ChernRoot+  chernTag1 _ = ChernRoot+  chernTag2 _ = ChernRoot+  chernTag3 _ = ChernRoot+  fromAB     = id+  fromChern  = chernToAB+  fromSchur  = schurToAB+  toAB       = id+  toChern    = abToChern+  toSchur    = abToSchur++instance ChernBase Chern where+  chernTag  _ = ChernClass+  chernTag1 _ = ChernClass+  chernTag2 _ = ChernClass+  chernTag3 _ = ChernClass+  fromAB     = abToChern+  fromChern  = id+  fromSchur  = schurToChern+  toAB       = chernToAB+  toChern    = id+  toSchur    = chernToSchur++--------------------------------------------------------------------------------+-- * Helper functions for constructing and specializing uniform things++-- | Constructing uniform things+select0 :: (AB, Chern) -> (ChernBase base => base)+select0 what = let final = select0' what (chernTag final) in final++select1 :: (f AB, f Chern) -> (ChernBase base => f base)+select1 what = let final = select1' what (chernTag1 final) in final++select2 :: (f (g AB), f (g Chern)) -> (ChernBase base => f (g base))+select2 what = let final = select2' what (chernTag2 final) in final++select3 :: (f (g (h AB)), f (g (h Chern))) -> (ChernBase base => f (g (h base)))+select3 what = let final = select3' what (chernTag3 final) in final++-- | Constructing unifom things using a tag+select0' :: (AB, Chern) -> (ChernBase base => Sing base -> base)+select0' (ab,ch) = \sing -> case sing of { ChernRoot -> ab ; ChernClass -> ch }++select1' :: (f AB, f Chern) -> (ChernBase base => Sing base -> f base)+select1' (ab,ch) = \sing -> case sing of { ChernRoot -> ab ; ChernClass -> ch }++select2' :: (f (g AB), f (g Chern)) -> (ChernBase base => Sing base -> f (g base))+select2' (ab,ch) = \sing -> case sing of { ChernRoot -> ab ; ChernClass -> ch }++select3' :: (f (g (h AB)), f (g (h Chern))) -> (ChernBase base => Sing base -> f (g (h base)))+select3' (ab,ch) = \sing -> case sing of { ChernRoot -> ab ; ChernClass -> ch }++-- | Specializing uniform things+spec0' :: ChernBase base => Sing base -> (forall b. ChernBase b => b) -> base+spec0' _ x = x++spec1' :: ChernBase base => Sing base -> (forall b. ChernBase b => f b) -> f base+spec1' _ x = x++spec2' :: ChernBase base => Sing base -> (forall b. ChernBase b => f (g b)) -> f (g base)+spec2' _ x = x++spec3' :: ChernBase base => Sing base -> (forall b. ChernBase b => f (g (h b))) -> f (g (h base))+spec3' _ x = x++{-+proxyOf :: a -> Proxy a+proxyOf _ = Proxy++proxyOf1 :: f a -> Proxy a+proxyOf1 _ = Proxy++proxyOf2 :: g (f a) -> Proxy a+proxyOf2 _ = Proxy+-}++--------------------------------------------------------------------------------++instance Monoid AB where+  mempty = AB 0 0 +  (AB a1 b1) `mappend` (AB a2 b2) = AB (a1+a2) (b1+b2)++instance Monoid Chern where+  mempty = Chern 0 0 +  (Chern e1 f1) `mappend` (Chern e2 f2) = Chern (e1+e2) (f1+f2)++instance Monoid Schur where+  mempty  = Schur 0 0+  mappend = error "Schur/mappend: not a monoid"++--------------------------------------------------------------------------------++instance Pretty AB where+  pretty ab = case ab of+    AB 0 0 -> "" +    AB e 0 -> showVarPower "a" e+    AB 0 f -> showVarPower "b" f+    AB e f -> showVarPower "a" e ++ "*" ++ showVarPower "b" f+ +instance Pretty Chern where+  pretty (Chern 0 0) = ""+  pretty (Chern e 0) = showVarPower "c1" e+  pretty (Chern 0 f) = showVarPower "c2" f+  pretty (Chern e f) = showVarPower "c1" e ++ "*" ++ showVarPower "c2" f++instance Pretty Schur where+  pretty (Schur a b) +    | b == 0     = "s[" ++ show a ++ "]"+    | otherwise  = "s[" ++ show a ++ "," ++ show b ++ "]"++--------------------------------------------------------------------------------+-- * Grading++class Graded a where+  grade :: a -> Int++instance Graded AB    where grade (AB    a b) = a + b+instance Graded Chern where grade (Chern e f) = e + 2*f+instance Graded Schur where grade (Schur i j) = i + j++filterGrade :: (Ord b, Graded b) => Int -> ZMod b -> ZMod b+filterGrade g = ZMod.onFreeMod filt where+  filt = Map.filterWithKey $ \x _ -> (grade x == g)++--------------------------------------------------------------------------------+-- * Conversions++chernToAB :: ZMod Chern -> ZMod AB +chernToAB = ZMod.flatMap expandToAlphaBeta_1 where++  -- | c1^k * c2^n = (alpha+beta)^k * (alpha*beta)^n+  expandToAlphaBeta_1 :: Chern -> ZMod AB +  expandToAlphaBeta_1 (Chern k n) = ZMod.fromList [ (AB (n+i) (n+k-i) , binomial k i) | i<-[0..k] ]++--------------------------------------------------------------------------------++-- | Converts a symmetric polynomial in the AB base (Chern roots) +-- to the Chern base (elementary symmetric polynomials or Chern classes)+abToChern :: ZMod AB -> ZMod Chern+abToChern ab = case symmetricReduction ab of+  Right c -> c+  Left _  -> error "abToChern: input was not symmetric"++-- | @Left@ means there is a non-symmetric remainder; @Right@ means+-- that input was symmetric.+symmetricReduction :: ZMod AB -> Either (ZMod Chern, ZMod AB) (ZMod Chern)+symmetricReduction = go [] where++  go sofar zmod = case ZMod.findMaxTerm zmod of+    Nothing          -> Right q+    Just (AB n m, k) -> if n < m+      then Left (q,zmod)+      else go ((ch,k):sofar) (zmod - this) where+        ch   = Chern (n-m) m+        this = ZMod.scale k $ expandToAlphaBeta_1 ch+    where+      q = ZMod.fromList sofar++  -- | c1^k * c2^n = (alpha+beta)^k * (alpha*beta)^n+  expandToAlphaBeta_1 :: Chern -> ZMod AB +  expandToAlphaBeta_1 (Chern k n) = ZMod.fromList [ (AB (n+i) (n+k-i) , binomial k i) | i<-[0..k] ]+            +--------------------------------------------------------------------------------++-- | Convert Schur to Chern roots+schurToAB :: ZMod Schur -> ZMod AB+schurToAB = ZMod.flatMap schurExpandAB_1 where++  schurExpandAB_1 :: Schur -> ZMod AB+  schurExpandAB_1 (Schur a b)+    | b > a     = error "schurExpandAB"+    | b < 0     = error "schurExpandAB"+    | otherwise = ZMod.fromList [ ( AB (a-j) (b+j) , 1 ) | j <- [0..a-b] ]++  {-+    schurab[i_, j_] := +     Expand[Factor[ Det[{{a^(i + 1), b^(i + 1)}, {a^j, b^j}}]] / +       Det[{{a, b}, {1, 1}}] ]+  -}++--------------------------------------------------------------------------------++-- | Convert Schur to Chern classes (elementary symmetric polynomials)+schurToChern :: ZMod Schur -> ZMod Chern+schurToChern = ZMod.flatMap schurExpandChern_1 where++  schurExpandChern_1 :: Schur -> ZMod Chern+  schurExpandChern_1 (Schur a b) +    | b > a     = error "schurExpandChern_1"+    | b < 0     = error "schurExpandChern_1"+    | otherwise = ZMod.fromList [ ( Chern (a-b-2*j) (b+j) , paritySignValue j * binomial (a-b-j) j ) | j <- [0..(div (a-b) 2)] ]++  --  schurcd[i_, j_] := SymmetricReduction[schurab[i, j], {a, b}, {c1, c2}][[1]]++--------------------------------------------------------------------------------++chernToSchur :: ZMod Chern -> ZMod Schur+chernToSchur = ZMod.flatMap chernExpandSchur_1 where++  chernExpandSchur_1 :: Chern -> ZMod Schur+  chernExpandSchur_1 (Chern e f)+    | e < 0 || f < 0 = error "chernExpandSchur"+    | otherwise      = ZMod.fromList [ ( Schur (e+f-i) (f+i) , catalanTriangle (e-i) i) | i<-[0..(div e 2)] ]++--------------------------------------------------------------------------------++abToSchur :: ZMod AB -> ZMod Schur+abToSchur = chernToSchur . abToChern++chernToSchurNaive :: ZMod Chern -> ZMod Schur+chernToSchurNaive = ZMod.fromList . go where++  go zmod = case ZMod.findMaxTerm zmod of +    Nothing             ->  []+    Just (Chern a b, k) -> ( s , k ) : go (zmod - this) where+      this = ZMod.scale k $ schurExpandChern_1 s+      s    = Schur (a+b) b++  schurExpandChern_1 :: Schur -> ZMod Chern+  schurExpandChern_1 (Schur a b) +    | b > a     = error "schurExpandChern_1"+    | b < 0     = error "schurExpandChern_1"+    | otherwise = ZMod.fromList [ ( Chern (a-b-2*j) (b+j) , paritySignValue j * binomial (a-b-j) j ) | j <- [0..(div (a-b) 2)] ]++--------------------------------------------------------------------------------+-- * random polynomials for testing++randomChernMonom :: IO Chern+randomChernMonom = do+  a <- randomRIO (0,30)+  b <- randomRIO (0,15)+  return (Chern a b)++randomSchurMonom :: IO Schur+randomSchurMonom = do+  a <- randomRIO (0,30)+  b <- randomRIO (0,30)+  return (Schur (a+b) b)++withRandomCoeff :: IO a -> IO (a,Integer)+withRandomCoeff rnd = do+  k <- randomRIO (-100,100)+  x <- rnd+  return (x,k)++randomChernPoly :: IO (ZMod Chern)   +randomChernPoly = do+  n <- randomRIO (0,100)+  ZMod.fromList <$> replicateM n (withRandomCoeff randomChernMonom)++randomSchurPoly :: IO (ZMod Schur)   +randomSchurPoly = do+  n <- randomRIO (0,100)+  ZMod.fromList <$> replicateM n (withRandomCoeff randomSchurMonom)++--------------------------------------------------------------------------------+
+ src/Math/RootLoci/CSM/Aluffi.hs view
@@ -0,0 +1,105 @@++-- | Aluffi's computation of the non-equivariant CSM in @P^n@+--+-- See: Paolo Aluffi: Characteristic classes of discriminants and enumerative geometry, Comm. in Algebra 26(10), 3165-3193 (1998).+--+{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}++module Math.RootLoci.CSM.Aluffi where++--------------------------------------------------------------------------------++import Data.List++import Control.Monad++import Math.Combinat.Classes+import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer+import Math.Combinat.Sets++import qualified Data.Map as Map ; import Data.Map (Map)+import qualified Data.Set as Set ; import Data.Set (Set)++import Data.Array (Array)+import Data.Array.IArray++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Classic+import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++--------------------------------------------------------------------------------+-- * CSM computation++-- | Paolo Aluffi's explicit formula for the (non-equivariant) CSM of open coincident root loci+aluffiOpenCSM :: Partition -> ZMod G+aluffiOpenCSM part@(Partition ps) = ZMod.invScale (aut part) xsum where+  n = sum ps+  d = length ps+  xsum = ZMod.fromList [ ( G (n-d+k) , coeff k ) | k<-[0..d] ] +  coeff k = negateIfOdd k +          $ signedBinomial (d-3) k * factorial k * factorial (d-k) * symPolyNum (d-k) (map fromIntegral ps)++-- | Summing together the open loci CSMs, we got the CSMs of the closures of the strata+aluffiClosedCSM :: Partition -> ZMod G+aluffiClosedCSM part@(Partition ps) = ZMod.sum opens where+  opens = [ aluffiOpenCSM q | q <- Set.toList (closureSet part) ]++--------------------------------------------------------------------------------+-- * Euler characteristics++-- | Euler characteristic, computed form 'aluffiOpenCSM'+aluffiOpenEuler :: Partition -> Integer+aluffiOpenEuler p = ZMod.coeffOf (G n) (aluffiOpenCSM p) where+  n = partitionWeight p++-- | Euler characteristic, computed form 'aluffiClosedCSM'+aluffiClosedEuler :: Partition -> Integer+aluffiClosedEuler p = ZMod.coeffOf (G n) (aluffiClosedCSM p) where+  n = partitionWeight p++--------------------------------------------------------------------------------++-- | It is easy to see from Aluffi\'s formula that only dimensions 1 and 2 has nonzero Euler characteristic.+-- This function implements the resulting rather trivial formula:+--+-- > chi( X_{n}   ) = 2+-- > chi( X_{p,q} ) = if p==q then 1 else 2+-- > chi( X_{...} ) = 0+--+openEulerChar :: Partition -> Integer+openEulerChar (Partition ps) = case ps of+  [n]   -> 2+  [a,b] -> if a==b then 1 else 2+  _     -> 0++--------------------------------------------------------------------------------+-- * General linear sections++-- | Converts the CSM class of a (locally closed?) projective variety Z to the Euler characteristics+-- of general linear sections of Z (so the first number will be @chi(Z)@, the second will be+-- @chi(Z cap H1)@, the third @chi(Z cap H1 cap H2)@ with @H1@, @H2@... being generic hyperplanes.+-- Finally the codim-th number will be the degree.+--+-- See: Paolo Aluffi: EULER CHARACTERISTICS OF GENERAL LINEAR SECTIONS AND POLYNOMIAL CHERN CLASSES,+-- Proposition 2.6+-- +csmToEulerOfLinearSections +  :: Int             -- ^ the dimension of the ambient projective space @P^n@+  -> ZMod G          -- ^ the CSM class+  -> [Integer]       -- ^ the resulting sequence of Euler characteristics+csmToEulerOfLinearSections n csm = [ euler i | i<-[0..n] ] where+  csmArr  = accumArray (flip const) 0 (0,n) [ (i,c) | (G i, c) <- ZMod.toList csm ] :: Array Int Integer+  euler k = foldl' (+) 0 [ signedBinomial (-k) i * csmArr ! (n-k-i) | i<-[0..n-k] ]++-- | We can compute the degree of the closures of the strata by intersection them+-- with @dim(X)@ generic hiperplanes.+aluffiDegree :: Partition -> Integer+aluffiDegree part = list !! dimension part where+  list = csmToEulerOfLinearSections (weight part) (aluffiClosedCSM part)++--------------------------------------------------------------------------------
+ src/Math/RootLoci/CSM/Equivariant/Direct.hs view
@@ -0,0 +1,65 @@++-- | We compute the open CSM classes directly, generalizing Aluffi's argument+-- to the equivariant case:+--+-- First we compute the CSM of set of the distinct /ordered/ points, then+-- push that forward first with @delta_*@ then with @pi_*@ to get the+-- CSM of the distinct unordered points with given multiplicities.+--+-- After that, we can get the closed CSM classes by summing over the+-- strata in the closure.+--+-- This is faster, especially since we have a (recursive) formula for the +-- CSM of the distinct ordered points.++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}+module Math.RootLoci.CSM.Equivariant.Direct +  ( directOpenCSM+  , directClosedCSM+  )+  where++--------------------------------------------------------------------------------++import Math.Combinat.Partitions.Integer++import qualified Data.Set as Set++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.CSM.Equivariant.PushForward+import qualified Math.RootLoci.CSM.Equivariant.Ordered as Ordered++--------------------------------------------------------------------------------++-- | CSM class of the open strata.+--  +-- We just push-forward first with Delta then down with Pi the conjectured +-- (recursive) formula for the CSM of the set of distinct ordered points+-- +directOpenCSM :: ChernBase base => Partition -> ZMod (Gam base)+directOpenCSM = polyCache2 directCalcOpenCSM where++  directCalcOpenCSM :: ChernBase base => Partition -> ZMod (Gam base)+  directCalcOpenCSM part@(Partition xs) = result where+    m = partitionWeight part+    result   = ZMod.invScale (aut part) $ pi_star m middle+    middle   = delta_star_ part distinct+    distinct = Ordered.formulaDistinctCSM (length xs)++--------------------------------------------------------------------------------++-- | To compute the CSM of the closed loci, we just some over the open strata+-- in the closure.+directClosedCSM :: ChernBase base => Partition -> ZMod (Gam base)+directClosedCSM = polyCache2 calc where+  +  calc :: ChernBase base => Partition -> ZMod (Gam base)+  calc part = ZMod.sum [ directOpenCSM q | q <- Set.toList (closureSet part) ] ++--------------------------------------------------------------------------------+
+ src/Math/RootLoci/CSM/Equivariant/Ordered.hs view
@@ -0,0 +1,412 @@++-- | 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++--------------------------------------------------------------------------------++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.RootLoci.Algebra.FreeMod 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 n = select2 +  ( tangentChernClassAB    n+  , tangentChernClassChern n+  )++tangentChernClassAB+  :: Int                  -- ^ The number of projective lines+  -> ZMod (Omega AB)+tangentChernClassAB d = ZMod.product [ entry i | i<-[1..d] ] where+  entry i = ZMod.fromList+    [ (Omega []  (AB 0 0) , 1)+    , (Omega []  (AB 1 0) , 1)+    , (Omega []  (AB 0 1) , 1)+    , (Omega [i] (AB 0 0) , 2)      -- 2x !+    ]++tangentChernClassChern+  :: Int                  -- ^ The number of projective lines+  -> ZMod (Omega Chern)+tangentChernClassChern d = ZMod.product [ entry i | i<-[1..d] ] where+  entry i = ZMod.fromList+    [ (Omega []  (Chern 0 0) , 1)+    , (Omega []  (Chern 1 0) , 1)+    , (Omega [i] (Chern 0 0) , 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 indices = unsafeEtaToOmega . delta_star' 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 n = smallDiagonal' [1..n] where++  smallDiagonal' :: [Int] -> ZMod (Omega base)+  smallDiagonal' indices = j_star [indices] (tangentChernClass 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 = polyCache2 calcOpenStratumCSM  where+             +  calcOpenStratumCSM :: forall b. ChernBase b => Int -> ZMod (Omega b)+  calcOpenStratumCSM d+    | d == 0     =  ZMod.one +    | d == 1     =  tangentChernClass 1+    | otherwise  = (tangentChernClass d) `ZMod.sub` (ZMod.sum diagonals)+    where+      diagonals = +        [ computeAnyStratumCSM setp+        | setp <- setPartitions d +        , let k = numberOfParts setp+        , k < 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 pps) = (j_star pps $ computeOpenStratumCSM $ length pps)++-- | We sum over the closure+computeClosureOfAnyStratumCSM :: ChernBase base => SetPartition -> ZMod (Omega base)+computeClosureOfAnyStratumCSM setp = ZMod.sum+  [ computeAnyStratumCSM p | p <- Set.toList (closureSetOfSetPartition setp) ] ++--------------------------------------------------------------------------------+-- * The structure lemma++-- | A formal monomial @q^k@+newtype QPow = QPow Int deriving (Eq,Ord,Show)++instance Monoid QPow where+  mempty = QPow 0+  mappend (QPow e) (QPow f) = QPow (e+f)++instance Pretty QPow where+  pretty (QPow k) = showVarPower "q" 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 n+  | n <  0  = error "umbralDistinct: n should be nonnegative"+  | n == 0  = ZMod.generator $ monom [] (-3)+  | n == 1  = ZMod.fromList  +                [ (monom []    (-2) , 1) +                , (monom [1]   (-3) , 2)+                ]+  | n == 2  = ZMod.fromList  +                [ (monom []    (-1) , 1)+                , (monom [1]   (-2) , 1)+                , (monom [2]   (-2) , 1)+                , (monom [1,2] (-3) , 2)+                ]+  | n >= 3  = ZMod.sum+                [ ZMod.scale coeff $ (ZMod.symPoly (n-3-k) us) * (ZMod.generator $ monom [] k)+                | k<-[0..n-3]+                , let coeff = negateIfOdd (n-3+k) (factorial (n-3) `div` factorial k)+                ]++  where+    monom xs k = Omega xs (QPow k)+    us = [ monom [i] 0 | i<-[1..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 subst1 input = ZMod.sum +  [ ZMod.fromList +      [ (Omega us ab , c*coeff) +      | (ab,c) <- ZMod.toList (subst1 qpow) +      ] +  | (Omega us qpow , coeff) <- ZMod.toList 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 = icache' ZMod.zero (-3) calcComputeQPolys where++  calcComputeQPolys :: Int -> ZMod AB+  calcComputeQPolys n +    | n <  -3    = error "computeQPolys: n >= -3 is required"+    | n == -3    = ZMod.one+    | otherwise  = ZMod.mapBase project almost+    where++      almost = open - smaller+      open   = computeOpenStratumCSM (n+3)     -- we should use this as the basis of the computation, unfortunately it's rather slow+    +      umbSmaller = umbralDistinctFormula (n+3) - umbHighest+      umbHighest = ZMod.generator (Omega [] (QPow n))        -- q^n+      smaller     = umbralSubstQPow (\(QPow k) -> computeQPolys k) umbSmaller++{-+      smaller = ZMod.sum +        [ ZMod.scale coeff $ +            (ZMod.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 us ab) = case us of+        [] -> ab+        _  -> error $ "computeQPolys: cannot project u terms:\n  " ++ pretty 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 n = select1 +  ( formulaQPolyAB   n +  , formulaQPolyChern n+  )++formulaQPolyAB :: Int -> ZMod AB+formulaQPolyAB = icache' ZMod.zero (-3) calcQPoly where+  +  calcQPoly :: Int -> ZMod AB+  calcQPoly n+    | n <  -3   = ZMod.zero+    | n == -3   = ZMod.konst 1+    | otherwise = mult1 * prev1 + mult2 * prev2+    where+      prev1 = formulaQPolyAB (n-1)+      prev2 = formulaQPolyAB (n-2)++      Pair mult1 mult2 = qpolyRecursionCoeffs 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 = icache' ZMod.zero (-3) calcQPoly where+  +  calcQPoly :: Int -> ZMod Chern+  calcQPoly n+    | n <  -3   = ZMod.zero+    | n == -3   = ZMod.konst 1+    | otherwise = mult1 * prev1 + mult2 * prev2+    where+      nn = fromIntegral n :: Integer++      prev1 = formulaQPolyChern (n-1)+      prev2 = formulaQPolyChern (n-2)++      Pair mult1 mult2 = qpolyRecursionCoeffs n++qpolyRecursionCoeffs :: ChernBase base => Int -> Pair (ZMod base)+qpolyRecursionCoeffs n = select2 +  (  Pair  mult1_AB    mult2_AB +  ,  Pair  mult1_Chern mult2_Chern+  )+  where++    mult1_AB = ZMod.fromList +      [ ( AB 0 0 ,     1 )+      , ( AB 1 0 , -nn-1 )+      , ( AB 0 1 , -nn-1 )+      ]+    mult2_AB = ZMod.singleton (AB 1 1) (-(nn-1)*(nn+2)) +  +    mult1_Chern = ZMod.fromList +      [ ( Chern 0 0 ,     1 )+      , ( Chern 1 0 , -nn-1 )+      ]+    mult2_Chern = ZMod.singleton (Chern 0 1) (-(nn-1)*(nn+2))++    nn = fromIntegral 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 n +  | n < 0     = error "formulaDistinctCSM: dimension should be nonnegative"+  | otherwise = umbralSubstQPow fun +              $ umbralDistinctFormula n+  where+    fun (QPow k) = formulaQPoly 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 setp = unsafeEtaToOmega $ delta_star setp (formulaDistinctCSM k) where+  k = numberOfParts setp+  +--------------------------------------------------------------------------------
+ src/Math/RootLoci/CSM/Equivariant/PushForward.hs view
@@ -0,0 +1,305 @@++-- | Compute the pushforward maps @pi_*@ and @delta_*@ between the+-- @GL2@-equivariant cohomology rings+--+-- Recall that:+--+-- * @Delta_nu : Q^d -> Q^n@+--+-- * @pi : Q^n -> P^n@+--+-- and @Q^n = P^1 x P^1 x ... x P^1@.+--  ++{-# LANGUAGE +      BangPatterns, TypeSynonymInstances, FlexibleInstances, FlexibleContexts,+      ScopedTypeVariables, TypeFamilies +  #-}++module Math.RootLoci.CSM.Equivariant.PushForward +  ( -- * The function tau+    tau , tauEta+    -- * pushforward along the diagonal map @Delta_{nu} : Q^d -> Q^n@+  , delta_star_ , delta_star , delta_star' +    -- * pushforward along the order-forgetting map @pi : Q^n -> P^n@+  , pi_star_table+  , compute_pi_star+  , pi_star+    -- * Fibonacci-type recursion formula for @pi_*@+  , piStarTableAff +  , piStarTableProj+  )+  where++--------------------------------------------------------------------------------++import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set+import Math.Combinat.Sets+import Math.Combinat.Tuples++import Data.Array (Array)+import Data.Array.IArray++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++--------------------------------------------------------------------------------+-- * The function tau++-- | @tau_k := ( a^(k+1) - b^(k+1) ) / ( a - b )@+tau :: ChernBase base => Int -> ZMod base+tau k = select1 ( tauAB k , tauChern k ) ++-- | In chern classes, the coefficients of tau are (signed) binomial coefficients; cf. A011973+tauChern :: Int -> ZMod Chern+tauChern k +  | k <  -1    = error "tau: negative index is not implemented"+  | k == -1    = ZMod.zero +  | otherwise  = ZMod.fromList [ ( Chern (k - 2*j) j , negateIfOdd j $ binomial (k-j) j )  | j<-[0..div k 2] ]++tauChernUnsafe :: Int -> ZMod Chern+tauChernUnsafe = icache $ \k -> abToChern (tauAB k)++tauAB :: Int -> ZMod AB+tauAB k +  | k <  -1   = error "tau: negative index is not implemented"+  | k == -1   = ZMod.zero+  | otherwise = ZMod.fromList [ (AB j (k-j) , 1) | j <- [0..k] ]++tauEta :: ChernBase base => Int -> ZMod (Eta base)+tauEta k = injectZMod (tau k)++--------------------------------------------------------------------------------+-- * @Delta_{\nu} : Q^d -> Q^n@++-- | Input: diagonal eta indices, and whether we are pushing forward 1 or the generator u/xi+delta_star_single :: ChernBase base => [Int] -> Bool -> ZMod (Eta base)+delta_star_single ks xi = +  if xi+    then bbb +    else aaa +  where+    n = length ks++    aaa = ZMod.sum [ sigma (n-1-i) * (tauEta i) | i<-[0..n-1] ]+    bbb = full - rest++    ab   = ZMod.generator $ Eta [] $ select0 (alphaBeta, c2)+    full = ZMod.generator (Eta ks mempty)      -- == sigma_n(eta)+    rest = ZMod.sum [ sigma (n-i) * tauEta (i-2) * ab | i<-[2..n] ]++    sigma k = ZMod.symPoly k [ Eta [k] mempty | k<-ks ]+  +-- | a group generator on the left is a subset (=product) of U-s, which+-- we map to a linear combinaton of H-s+delta_star_1 :: ChernBase base => Partition -> Omega base -> ZMod (Eta base)+delta_star_1 part = delta_star_1' (linearIndices part)++-- | a group generator on the left is a subset (=product) of U-s, which+-- we map to a linear combinaton of H-s+delta_star_1' :: forall base. ChernBase base => [[Int]] -> Omega base -> ZMod (Eta base)+delta_star_1' idxtable (Omega us ab) = final where+  +  final = mulInjMonom ab $ ZMod.product $ go 1 idxtable+          +  go :: Int -> [[Int]] -> [ZMod (Eta base)]+  go _ []       = []+  go k (is:iss) = this : go (k+1) iss where+    this = delta_star_single is (k `elem` us)++delta_star_ :: ChernBase base => Partition -> ZMod (Omega base) -> ZMod (Eta base)+delta_star_ part = ZMod.flatMap (delta_star_1 part)++delta_star :: ChernBase base => SetPartition -> ZMod (Omega base) -> ZMod (Eta base)+delta_star setp = ZMod.flatMap (delta_star_1' (fromSetPartition setp))++-- | We can give an explicit indexing scheme (set partition), instead of the linear indexing+-- used above. This will be useful when computing the \"open\" part+delta_star' :: ChernBase base => [[Int]] -> ZMod (Omega base) -> ZMod (Eta base)+delta_star' indices = ZMod.flatMap (delta_star_1' indices)++--------------------------------------------------------------------------------+-- * @pi : Q^n -> P^n@++-- | This is upside the class where @[0:1]@ is a root with multiplicity @k@ and @[1:0]@ is a root with multiplicity l+up_root_xy :: Int -> (Int,Int) -> ZMod (Eta AB)+up_root_xy n (k,l) = as * bs where++  as = ZMod.product [ abh      i  1 0 | i<-[1..k] ]+  bs = ZMod.product [ abh (n+1-j) 0 1 | j<-[1..l] ]++  -- (eta_i + na*alpha + nb*beta)+  abh i na nb = ZMod.fromList +    [ (Eta [i] (AB 0 0) , 1 )  +    , (Eta []  (AB 1 0) , na)  +    , (Eta []  (AB 0 1) , nb) +    ]++-- | This is downside the class where @[0:1]@ is a root with multiplicity @k@ and @[1:0]@ is a root with multiplicity l.+-- It should be true that @pi_* up_root_xy = down_root_xy@+down_root_xy :: Int -> (Int,Int) -> ZMod (Gam AB)+down_root_xy n (k,l) = as * bs where++  as = ZMod.product [ abg (n-i) (i) | i<-[0..k-1] ]+  bs = ZMod.product [ abg (j) (n-j) | j<-[0..l-1] ]++  -- (na*alpha + nb*beta + gamma)+  abg na nb = ZMod.fromList +    [ (Gam 1 (AB 0 0) , 1 ) +    , (Gam 0 (AB 1 0) , fromIntegral na) +    , (Gam 0 (AB 0 1) , fromIntegral nb) +    ]++pi_star_0 :: Int -> Int -> ZMod (Gam AB)+pi_star_0 n k = ZMod.sum+  [ ZMod.scale +      (negateIfOdd i $ binomial k i * factorial (n-k+i)) +      (mulAB (AB i 0) $ down_root_xy n (k-i,0)) +  | i<-[0..k] ]++-- | Table of @pi_*( eta_1*eta_2*...*eta_k )@, computed by breaking the symmetry.+pi_star_table :: Int -> Array Int (ZMod (Gam AB))+pi_star_table = monoCache calc where+  calc n = listArray (0,n) [ pi_star_0 n k | k<-[0..n] ]++-- | Slow implementation of @pi_star@, using @pi_star_table@+compute_pi_star +  :: Int               -- ^ the number of points @m@ (recall the pi : @Q^m -> P^m@)+  -> ZMod (Eta AB) +  -> ZMod (Gam AB)+compute_pi_star m = ZMod.flatMap f where +  table = pi_star_table m+  f (Eta hs ab) = mulAB ab (table ! length hs)++--------------------------------------------------------------------------------+-- * Fibonacci-type recursion formula for @pi_*@++-- | However it should faster to just use the recursion for the @P_j(m)@ polynomials,+-- which this function does.+pi_star +  :: forall base. (ChernBase base) +  => Int                      -- ^ the number of points @m@ (recall the pi : @Q^m -> P^m@)+  -> ZMod (Eta base) +  -> ZMod (Gam base)+pi_star m = ZMod.flatMap f where +  table = piStarTableProj m :: Array Int (ZMod (Gam base))+  f (Eta hs ab) = mulInjMonom ab (table ! length hs)++piStarTableAff :: ChernBase base => Int -> Array Int (ZMod base)+piStarTableAff = polyCache2 calc where+  calc n = select2 ( aff_fibPiStar_AB n , aff_fibPiStar_Chern n )++piStarTableProj :: ChernBase base => Int -> Array Int (ZMod (Gam base))+piStarTableProj = polyCache3 calc where+  calc n = select3 ( proj_fibPiStar_AB n , proj_fibPiStar_Chern n )++{-++class ChernBase (PiStarBase tgtmonom) => PiStar tgtmonom where+  type PiStarBase tgtmonom :: *+  piStarTable :: Int -> Array Int (ZMod tgtmonom)++instance PiStar (Gam Chern) where { piStarTable = proj_fibPiStar_Chern ; type PiStarBase (Gam Chern) = Chern }+instance PiStar (Gam AB   ) where { piStarTable = proj_fibPiStar_AB    ; type PiStarBase (Gam AB   ) = AB    }+instance PiStar      Chern  where { piStarTable = aff_fibPiStar_Chern  ; type PiStarBase (Chern    ) = Chern }+instance PiStar      AB     where { piStarTable = aff_fibPiStar_AB     ; type PiStarBase (AB       ) = AB    }++-- instance PiStar (Gam Schur) where { piStarTable = proj_fibPiStar_Schur ; type PiStarBase = Gam Schur }+-- instance PiStar      Schur  where { piStarTable = aff_fibPiStar_Schur  ; type PiStarBase = Schur     }+-}++proj_fibPiStar_Chern :: Int -> Array Int (ZMod (Gam Chern))+proj_fibPiStar_Chern m = listArray (0,m) $ take (m+1) fib where++  fib :: [ZMod (Gam Chern)]+  fib = ZMod.konst                    (factorial  m   )+      : ZMod.singleton (Gam 1 mempty) (factorial (m-1))+      : zipWith3 g [1..] (tail fib) fib ++  g :: Integer -> ZMod (Gam Chern) -> ZMod (Gam Chern) -> ZMod (Gam Chern)+  g k prev1 prev2 +    = ZMod.invScale (mm-k)+    $ mulGam prev1 + ZMod.scale k (mulInjMonom c1 prev1) +                   + ZMod.scale k (mulInjMonom c2 prev2) ++  mm = fromIntegral m :: Integer++--  c1 = Chern 1 0+--  c2 = Chern 0 1++----------------------------------------++aff_fibPiStar_Chern :: Int -> Array Int (ZMod Chern)+aff_fibPiStar_Chern m = listArray (0,m) $ take (m+1) fib where++  fib :: [ZMod Chern]+  fib = ZMod.konst (factorial m)+      : ZMod.zero+      : zipWith3 g [1..] (tail fib) fib ++  g :: Integer -> ZMod Chern -> ZMod Chern -> ZMod Chern+  g k prev1 prev2 +    = ZMod.invScale (mm-k)+    $ ZMod.scale    (   k)+    $ (ZMod.mulMonom c1 prev1 + ZMod.mulMonom c2 prev2) ++  mm = fromIntegral m :: Integer++--  c1 = Chern 1 0+--  c2 = Chern 0 1++----------------------------------------++proj_fibPiStar_AB :: Int -> Array Int (ZMod (Gam AB))+proj_fibPiStar_AB m = fmap (convertGam chernToAB) (proj_fibPiStar_Chern m)++proj_fibPiStar_Schur :: Int -> Array Int (ZMod (Gam Schur))+proj_fibPiStar_Schur m = fmap (convertGam chernToSchur) (proj_fibPiStar_Chern m)++aff_fibPiStar_AB :: Int -> Array Int (ZMod AB)+aff_fibPiStar_AB m =  fmap chernToAB (aff_fibPiStar_Chern m)++aff_fibPiStar_Schur :: Int -> Array Int (ZMod Schur)+aff_fibPiStar_Schur m =  fmap chernToSchur (aff_fibPiStar_Chern m)++--------------------------------------------------------------------------------+-- * helpers++-- | Multiplies by an injected monomial +mulInjMonom :: (Functor f, Monoid ab, Ord (f ab)) => ab -> ZMod (f ab) -> ZMod (f ab)+mulInjMonom monom = ZMod.mapBase f where+  f = fmap (mappend monom)++-- | Multiplies by @(alpha^i * beta^j)@+mulAB :: (Functor f, Ord (f AB)) => AB -> ZMod (f AB) -> ZMod (f AB)+mulAB = mulInjMonom++-- | Multiplies with @gamma@+mulGam :: Ord ab => ZMod (Gam ab) -> ZMod (Gam ab)+mulGam = ZMod.mapBase f where +  f (Gam k x) = Gam (k+1) x++{-+-- | Multiplies by alpha^i beta^j+omegaMulAB :: AB -> ZMod (Omega AB) -> ZMod (Omega AB)+omegaMulAB (AB i j) = Map.mapKeys f where+  f (Omega us (AB a b)) = Omega us (AB (a+i) (b+j))++-- | Multiplies by alpha^i beta^j+etaMulAB :: AB -> ZMod (Eta AB)-> ZMod (Eta AB)+etaMulAB (AB i j) = Map.mapKeys f where+  f (Eta hs (AB a b)) = Eta hs (AB (a+i) (b+j))++-- | Multiplies by alpha^i beta^j+gamMulAB :: AB -> ZMod (Gam AB) -> ZMod (Gam AB)+gamMulAB (AB i j) = Map.mapKeys f where+  f (Gam g (AB a b)) = Gam g (AB (a+i) (b+j))+-}++--------------------------------------------------------------------------------
+ src/Math/RootLoci/CSM/Equivariant/Recursive.hs view
@@ -0,0 +1,97 @@++-- | We compute the @GL2@-equivariant open and closed CSM classes recursively,+-- starting from smallest strata. +--+-- The idea is that we have a smooth resolution of the /closure/ of the strata @X_mu@, +-- namely, the set of @n=length(mu)@ ordered points: @Q^n = P^1 x ... x P^1@+--+-- We can pushforward this to @Q^m@, and get a linear combination of the strata of+-- the CSM-s we want to compute. Since the smallest strata is actually closed,+-- we know that, and can work upward from that.+--+-- This is rather slow, however as it's a very different algorithm copmared to+-- the direct approach, it's useful for checking if the two agrees.+--++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}+module Math.RootLoci.CSM.Equivariant.Recursive where++--------------------------------------------------------------------------------++import qualified Data.Set as Set ; import Data.Set (Set)+import qualified Data.Map as Map ; import Data.Map (Map)++import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set++import qualified Math.RootLoci.CSM.Equivariant.Ordered as Ordered+import Math.RootLoci.CSM.Equivariant.PushForward++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++--------------------------------------------------------------------------------+-- * CSM calculation++-- | This is just the pushforward along @Delta_nu@ of the tangent Chern class.+--+-- As @Delta@ is injective, the resulting class is just the CSM class of the+-- closed /ordered/ strata corresponding to one of the set partitions which+-- matches the given partition +----+upperClass :: ChernBase base => SetPartition -> ZMod (Eta base)+upperClass = polyCache2 calcUpper where++  calcUpper :: ChernBase base => SetPartition -> ZMod (Eta base)+  calcUpper part@(SetPartition ps) = delta_star part (Ordered.tangentChernClass d) where+    d = length ps++-- | pushforward of `upperCSM` to the space of unordered points+lowerClass :: ChernBase base => Partition -> ZMod (Gam base)+lowerClass = polyCache2 calcLower where++  calcLower :: ChernBase base => Partition -> ZMod (Gam base)+  calcLower part = pi_star n (upperClass $ defaultSetPartition part) where+    n = partitionWeight part++--------------------------------------------------------------------------------++-- | We know from the pushforward property of CSM clsses that @(pi_* upperCSM) = sum (chi * openCSM)@.+-- we can use this to recursively compute the CSM classes of the open loci+--+openCSM :: ChernBase base => Partition -> ZMod (Gam base)+openCSM = polyCache2 calcOpenCSM where++  calcOpenCSM :: ChernBase base => Partition -> ZMod (Gam base)+  calcOpenCSM part = ZMod.invScale thisCoeff (pushdown `ZMod.sub` smaller) where+    n = partitionWeight part+    pushdown  = lowerClass part+    smaller   = ZMod.linComb [ (c , openCSM q) | (q,c) <- Map.assocs theClosure ]+    (thisCoeff,theClosure) = preimageView part  -- closureView' part++-- | To compute the CSM of the closed loci, we just some over the open strata+-- in the closure. +closedCSM :: ChernBase base => Partition -> ZMod (Gam base)+closedCSM = polyCache2 calcClosedCSM where++  calcClosedCSM :: ChernBase base => Partition -> ZMod (Gam base)+  calcClosedCSM part =+    ZMod.sum [ openCSM q | q <- Set.toList (closureSet part) ]  +      +--------------------------------------------------------------------------------++{-+equivDualClass :: Partition -> ZMod Gam+equivDualClass part = filterGrade (codim part) (closedCSM part)++equivOpenEuler :: Partition -> ZMod Gam+equivOpenEuler part = filterGrade (partitionWeight part) (openCSM part)++equivClosedEuler :: Partition -> ZMod Gam+equivClosedEuler part = filterGrade (partitionWeight part) (closedCSM part)+-}++--------------------------------------------------------------------------------
+ src/Math/RootLoci/CSM/Equivariant/Umbral.hs view
@@ -0,0 +1,214 @@++-- | The umbral formula for the open CSM classes.+--+-- The formula is the following:+--+-- > A(mu)    = 1 / aut(mu) * prod_i Theta(mu_i)+-- > Theta(p) = ( (1 + beta*s) (alpha+t)^p - (1 + alpha*s) (beta+t)^p ) / ( alpha - beta )+--+-- and the umbral subtitution resulting in the CSM class (at least for @length(mu)>=3@) is:+--+-- > t^j  ->  P_j(m)+-- > s^k  ->  (n-3)(n-3-1)(...n-3-k+1) * Q(n-3-k)+--+-- Note that Theta(p) is actually a (symmetric) polynomial in @alpha@ and @beta@; furthermore+-- it's linear in s and degree p in t. ++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances, ScopedTypeVariables #-}+module Math.RootLoci.CSM.Equivariant.Umbral where++--------------------------------------------------------------------------------++import Math.Combinat.Classes+import Math.Combinat.Numbers+import Math.Combinat.Partitions.Integer++import Data.Array.IArray++import qualified Data.Set as Set++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.CSM.Equivariant.PushForward ( tau , piStarTableAff , piStarTableProj )+import Math.RootLoci.CSM.Equivariant.Ordered     ( formulaQPoly )++import qualified Math.RootLoci.CSM.Equivariant.Direct as Direct++--------------------------------------------------------------------------------+-- * The umbral variables++-- | A monomial @s^k * t^j@+data ST +  = ST !Int !Int+  deriving (Eq,Ord,Show)++instance Monoid ST where+  mempty = ST 0 0 +  (ST s1 t1) `mappend` (ST s2 t2) = ST (s1+s2) (t1+t2)++instance Pretty ST where+  pretty st = case st of+    ST 0 0 -> "" +    ST e 0 -> showVarPower "s" e+    ST 0 f -> showVarPower "t" f+    ST e f -> showVarPower "s" e ++ "*" ++ showVarPower "t" f++prettyMixedST :: forall b c. (Pretty b, Num c, Eq c, IsSigned c, Show c) => FreeMod (FreeMod c b) ST -> String+prettyMixedST = prettyFreeMod'' prettyInner pretty where++  prettyInner :: FreeMod c b -> String+  prettyInner = paren . pretty++--------------------------------------------------------------------------------+-- * The umbral formula++-- | @Theta(p)@ is defined by the formula+--+-- > Theta(p) = ( (1 + beta*s) (alpha+t)^p - (1 + alpha*s) (beta+t)^p ) / ( alpha - beta )+--+-- This is actually a polynomial in @alpha@,@beta@,@s@,@t@, also symmetric in @alpha@ and @beta@+--+theta :: ChernBase base => Int -> FreeMod (ZMod base) ST+theta p +  | p >= 1    = ZMod.fromList (term0 ++ term1) +  | otherwise = error "theta: non-positive input"+  where+ +    term0 =  [ (ST 0 i , ZMod.scale (binomial p i) (                         tau (p-i-1)) ) | i<-[0..p-1] ]+    term1 =  [ (ST 1 i , ZMod.scale (binomial p i) (ZMod.mulMonom c2_monom $ tau (p-i-2)) ) | i<-[0..p-2] ] +          ++ [ (ST 1 p , ZMod.konst (-1) ) ]++    c2_monom = select0 (alphaBeta,c2)++-- | Same as 'theta' but with rational coefficients+thetaQ :: ChernBase b => Int -> FreeMod (QMod b) ST+thetaQ p = ZMod.mapCoeff (ZMod.mapCoeff fromIntegral) (theta p)++-------------------------------------------------------------------------------- ++-- | This is just @prod_i Theta_{mu_i}@+integralUmbralFormula :: ChernBase base => Partition -> FreeMod (ZMod base) ST +integralUmbralFormula (Partition ps) = ZMod.product [ theta p | p <- ps ]++-- | This is @1/aut(mu) * prod_i Theta_{mu_i}@+umbralFormula :: ChernBase base => Partition -> FreeMod (QMod base) ST +umbralFormula mu@(Partition ps) = result where+ +  result = ZMod.mapCoeff (ZMod.scale (1 / autmu))+         $ ZMod.product [ thetaQ p | p <- ps ]++  autmu :: Rational+  autmu = fromIntegral (aut mu)++--------------------------------------------------------------------------------+-- * The affine CSM++-- | The polynomial to be substituted in the place of @s^k*t^j@:+--+-- > s^k*t^j  ->  P_j(m) * Q_k(n-3-k) * (n-3)_k+--+-- where @n = length(mu)@ and @m = weight(mu)@.+--+umbralSubstPolyAff :: ChernBase base => Partition -> ST -> ZMod base+umbralSubstPolyAff part = fun where++  n = numberOfParts part+  m = weight part+  tablePPoly = piStarTableAff m++  fun (ST k j) +    | k >= -3 && k <= n-3 && j >= 0 && j <= m  = ZMod.scale falling (qpoly `ZMod.mul` ppoly)+    | otherwise                                = ZMod.zero+    where+      falling :: Integer+      falling = product [ fromIntegral (n-3-i) | i<-[0..k-1] ]++      qpoly   = formulaQPoly (n-3-k)+      ppoly   = tablePPoly ! j++-- | The (affine) umbral substitution+umbralSubstitutionAff :: (ChernBase base) => Partition -> FreeMod (ZMod base) ST -> ZMod base+umbralSubstitutionAff part input = output where++  output   = ZMod.sum [ ab `ZMod.mul` (substfun st) | (st,ab) <- ZMod.toList input ]+  substfun = umbralSubstPolyAff part++-- | CSM of the open stratums from the umbral the formula+umbralAffOpenCSM :: ChernBase base => Partition -> ZMod base   +umbralAffOpenCSM = polyCache1 calc where++  -- the current umbral formula only works for @n >= 3@ ??+  calc mu +    | n < 3     = forgetGamma (Direct.directOpenCSM mu)+    | otherwise = ZMod.invScale (aut mu)+                $ umbralSubstitutionAff mu+                $ integralUmbralFormula mu+    where+      n = numberOfParts mu++-- | Sum over the strata in the closure+umbralAffClosedCSM :: ChernBase base => Partition -> ZMod base   +umbralAffClosedCSM = polyCache1 calc where+  +  calc :: ChernBase base => Partition -> ZMod base+  calc part = ZMod.sum [ umbralAffOpenCSM q | q <- Set.toList (closureSet part) ] ++--------------------------------------------------------------------------------+-- * The projective CSM++-- | The polynomial to be substituted in the place of @s^k*t^j@:+--+-- > s^k*t^j  ->  P_j(m) * Q_k(n-3-k) * (n-3)_k+--+-- where @n = length(mu)@ and @m = weight(mu)@.+--+umbralSubstPolyProj :: forall base. ChernBase base => Partition -> ST -> ZMod (Gam base)+umbralSubstPolyProj part = fun where++  n = numberOfParts part+  m = weight part+  tablePPoly = piStarTableProj m++  fun (ST k j) +    | k >= -3 && k <= n-3 && j >= 0 && j <= m  = ZMod.scale falling (qpoly `ZMod.mul` ppoly)+    | otherwise                                = ZMod.zero+    where+      falling :: Integer+      falling = product [ fromIntegral (n-3-i) | i<-[0..k-1] ]++      qpoly   = injectZMod (formulaQPoly (n-3-k)) :: ZMod (Gam base)+      ppoly   = tablePPoly ! j                    :: ZMod (Gam base)+++-- | The (projective) umbral substitution+umbralSubstitutionProj :: (ChernBase base) => Partition -> FreeMod (ZMod base) ST -> ZMod (Gam base)+umbralSubstitutionProj part input = output where++  output   = ZMod.sum [ injectZMod ab `ZMod.mul` (substfun st) | (st,ab) <- ZMod.toList input ]+  substfun = umbralSubstPolyProj part++-- | CSM of the open stratums from the umbral the formula (for @length(mu) >= 3@)+umbralOpenCSM :: ChernBase base => Partition -> ZMod (Gam base)+umbralOpenCSM = polyCache2 calc where++  -- the current umbral formula only works for @n >= 3@ ??+  calc mu +    | n < 3     = Direct.directOpenCSM mu     +    | otherwise = ZMod.invScale (aut mu)+                $ umbralSubstitutionProj mu+                $ integralUmbralFormula mu+    where+      n = numberOfParts mu++-- | Sum over the strata in the closure+umbralClosedCSM :: ChernBase base => Partition -> ZMod (Gam base)+umbralClosedCSM = polyCache2 calc where+  +  calc :: ChernBase base => Partition -> ZMod (Gam base)+  calc part = ZMod.sum [ umbralOpenCSM q | q <- Set.toList (closureSet part) ] ++--------------------------------------------------------------------------------
+ src/Math/RootLoci/CSM/Projective.hs view
@@ -0,0 +1,218 @@++-- | Compute the non-equivariant CSM in @P^n@ recursively++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}+module Math.RootLoci.CSM.Projective +  ( -- * Pushforwards+    delta_star+  , pi_star+    -- * Easy things+  , tangentChernClass+  , smallestOrbitCSM+    -- * CSM calculation+  , upperCSM , lowerCSM+  , openCSM  , closedCSM+    -- * extracting coefficients+  , highestCoeff_ , lowestCoeff_+  , highestCoeff  , lowestCoeff +  ) +  where++--------------------------------------------------------------------------------++import Data.List+import Data.Maybe++import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set+import Math.Combinat.Sets++import qualified Data.Map as Map ; import Data.Map (Map)+import qualified Data.Set as Set ; import Data.Set (Set)++import Data.Array.IArray+import Data.Array (Array)++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.Algebra.FreeMod as ZMod+ +--------------------------------------------------------------------------------++{-+  +we have maps+* Delta_nu : Q^d -> Q^n+* pi : Q^n -> P^n+  +-}++--------------------------------------------------------------------------------+-- * The order-forgetting map @pi : Q^n -> P^n@++pi_star_1 :: Int -> HS -> (G,Integer)+pi_star_1 n (HS hs) = (gk,c) where+  c  = factorial (n - length hs) +  gk = G (length hs)++-- | The pushforward map @pi_*@ along @pi@.+--+-- A (cohomology) group generator above is a subset (=product) of H-s, which we map to+-- a group generator below. This defines the map on the cohomology ring by additive extension.+--+pi_star +  :: Int           -- ^ the number of points @m@ (with multiplicity)+  -> ZMod HS       -- ^ the cohomoly class \"up\"+  -> ZMod G+pi_star n = ZMod.flatMap (sing . pi_star_1 n) where +  sing (b,c) = ZMod.singleton b c++--------------------------------------------------------------------------------+-- * The diagonal maps @Delta_{\nu} : Q^d -> Q^n@+  +delta_star_1 :: Partition -> US -> ZMod HS+delta_star_1 part@(Partition ps) (US us) = ZMod.histogram almost where++  n = sum    ps+  d = length ps+  +  idxtable = linearIndices part+      +  -- inner lists = monoms+  -- outer lists = linear combination of monoms+  -- now we want to multiply those together+  stuff :: [[[H]]]+  stuff = (map . map . map) H (go 1 idxtable)+  +  almost :: [HS]+  almost = map (HS . concat) $ listTensor stuff     -- this does the multiplication of terms+  +  uis = [ i | U i <- us ]+    +  go :: Int -> [[Int]] -> [[[Int]]]+  go _ []       = []+  go k (is:iss) = this : go (k+1) iss where+    this = if k `elem` uis+      then [is]                     -- "sigma_k"+      else chooseN1 is              -- "sigma_(k-1)"+  +-- | A group generator on the left is a subset (=product) of U-s, which+-- we map to a linear combinaton of H-s. This is then extended additively+-- to the cohomology ring.+--+delta_star :: Partition -> ZMod US -> ZMod HS+delta_star part = ZMod.flatMap (delta_star_1 part)++--------------------------------------------------------------------------------+-- * Easy things++-- | The total Chern class of the tangent bundle of @Q^d = P^1 x P^1 x ... x P^1@+--+-- This is just the product of @(1+2u_i)@-s for @i=[1..d]@+--+tangentChernClass :: Int -> ZMod US+tangentChernClass d = ZMod.fromList $ concatMap worker [0..d] where+  worker k = map (\xs -> (US (map U xs) , 2^k)) (choose_ k d)++-- | The CSM of the smallest orbit: 1 point with multiplicity @n@,+-- which is just the rational normal curve in @P^n@.+--+smallestOrbitCSM :: Int -> ZMod G+smallestOrbitCSM n = ZMod.fromList +  [ (G (n-1) ,     fromIntegral n) +  , (G  n    , 2 * fromIntegral n) +  ] ++--------------------------------------------------------------------------------+-- * CSM calculation++-- | We know that:+-- +-- > csm(im(Delta) = Delta_* c(TQ^d)+-- > c(TQ^d) = (1+2*u1) (1+2*u2) ... (1+2*ud)+--+-- From these, we can compute @csm(im(Delta_nu))@ recursively+--+upperCSM :: Partition -> ZMod HS+upperCSM = pcache calc where++  calc part@(Partition ps) = (delta_star part) (tangentChernClass d) where+    d = length ps++-- | A formula for @pi_*(csm(im(delta)))@. This should satisfy+--+-- > lowerCSM part = pi_star n (upperCSM part)+--+lowerCSM :: Partition -> ZMod G+lowerCSM = pcache calc where++  calc part@(Partition ps) = zmod where+    d = length ps+    n = sum ps+    zmod = ZMod.fromList+      [ ( G (n-d+r) , coeff )+      | r<-[0..d]+      , let coeff = factorial (d-r) * 2^r * symPolyNum (d-r) (map fi ps)+      ]+  +    fi :: Int -> Integer+    fi = fromIntegral++check_lower_upper :: Int -> Bool+check_lower_upper n = and [ pi_star n (upperCSM p) == lowerCSM p | p <- partitions n ]++-- | Cached CSM computation of the open strata+openCSM :: Partition -> ZMod G+openCSM = pcache calcOpenCSM where++  -- | we know that (pi_* upperCSM) = sum (chi * openCSM)+  calcOpenCSM :: Partition -> ZMod G+  calcOpenCSM part = ZMod.invScale thisCoeff (pushdown - smaller) where+    n = partitionWeight part+    pushdown  = lowerCSM part -- pi_star n (upperCSM part) +    smaller   = ZMod.linComb [ (c , openCSM q) | (q,c) <- Map.assocs theClosure ]+    (thisCoeff,theClosure) = preimageView part++-- | To get the CSM of the closed strata, we just sum over the open strata contained+-- in the closure.++closedCSM :: Partition -> ZMod G +closedCSM = pcache calcClosedCSM where  ++  calcClosedCSM :: Partition -> ZMod G+  calcClosedCSM part = ZMod.sum [ openCSM q | q <- Set.toList (closureSet part) ]++--------------------------------------------------------------------------------++lowestCoeff_ :: ZMod G -> Integer+lowestCoeff_ = snd . lowestCoeff++highestCoeff_ :: ZMod G -> Integer+highestCoeff_ = snd . highestCoeff++lowestCoeff :: ZMod G -> (G,Integer)+lowestCoeff = fromJust . ZMod.findMinTerm +-- lowestCoeff = head . ZMod.toList ++highestCoeff :: ZMod G -> (G,Integer)+highestCoeff = fromJust . ZMod.findMaxTerm+-- highestCoeff = last . ZMod.toList ++--------------------------------------------------------------------------------++{-+check_degree :: Partition -> Bool+check_degree p = hilbert p == lowestCoeff_ (closedCSM p)++check_euler_degree :: Partition -> Bool+check_euler_degree p@(Partition ps) = hilbert p == ((csmToEuler n $ closedCSM p) !! d) where+  d = length ps+  n = sum ps+-}++--------------------------------------------------------------------------------+
+ src/Math/RootLoci/Classic.hs view
@@ -0,0 +1,91 @@++-- | Classical results: +--+-- * Hilbert's degree formula+--+-- * some enumarative geometry computations by Schubert+--++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}++module Math.RootLoci.Classic where++--------------------------------------------------------------------------------++import Data.List++import Control.Monad++import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer+import Math.Combinat.Sets++--------------------------------------------------------------------------------++-- | Codimension of a strata. This is simply @(sum mu_i) - length(mu)@.+codim :: Partition -> Int+codim (Partition ps) = sum ps - length ps++-- | Dimension of the strata. @dim = length(mu)@.+dimension :: Partition -> Int+dimension (Partition ps) = length ps++--------------------------------------------------------------------------------+-- * Hilbert formula++-- | Hilbert's formula for the degree of a stratum+hilbert :: Partition -> Integer+hilbert part@(Partition ps) = div numer denom where++  n = sum    ps+  d = length ps++  numer = factorial d * product (map fi ps)          -- d! * prod (nu_i)+  denom = product (map (factorial . snd) ies)        -- prod (e_r!)+ +  ies = toExponentialForm part      -- (r,e_r) pairs+   +  fi :: Int -> Integer+  fi = fromIntegral++-- | Hilbert's degree formula, another version (as a sanity test).+hilbert2 :: Partition -> Integer+hilbert2 part@(Partition ps) = div numer denom where++  -- this is from FNR, opposite notation (d and n are swapped!)+  -- just to be really sure about the formula :)++  n = sum es+  d = sum [ i*ei | (i,ei) <- toExponentialForm part ]+  es =    [ ei   | (i,ei) <- toExponentialForm part ]++  numer = factorial n * product [ (fi i)^ei | (i,ei) <- toExponentialForm part ]+  denom = product [ factorial ei | (i,ei) <- toExponentialForm part ]++  fi :: Int -> Integer+  fi = fromIntegral+   +-- check_hilbert2 :: Bool   +-- check_hilbert2 = and [ hilbert p == hilbert2 p | n<-[0..20] , p<-partitions n ]++--------------------------------------------------------------------------------+-- * Schubert++-- | Number of 4-tangent lines to a generic degree @d@ surface +quadTangentLines :: Int -> Integer+quadTangentLines d0+  | d < 8     = 0+  | otherwise = d * (d - 4) * (d - 5) * (d - 6) * (d - 7) * (d^3 + 6*d^2 + 7*d - 30)+  where+    d = fromIntegral d0 :: Integer++-- | Number of lines meeting a generic degree @d@ surface at point with 5x multiplicity+quintFlexLines :: Int -> Integer+quintFlexLines d0+  | d < 5     = 0+  | otherwise = error "quintFlexLines"+  where+    d = fromIntegral d0 :: Integer++--------------------------------------------------------------------------------
+ src/Math/RootLoci/Dual/Localization.hs view
@@ -0,0 +1,115 @@++-- | Localization formula for the dual class from:+--+-- L. M. Feher, A. Nemethi, R. Rimanyi: Coincident root loci of binary forms;+-- Michigan Math. J. Volume 54, Issue 2 (2006), 375--392.+--+-- Note: This formula is in the form of /rational function/ (as opposed to +-- a polynomial). Since we don't have polynomial factorization implemented here,+-- instead we /evaluate/ it substituting different rational numbers+-- into @alpha@ and @beta@, and then use Lagrange interpolation to figure+-- out the result (we know a priori that it is a homogenenous polynomial+-- in @alpha@ and @beta@).++module Math.RootLoci.Dual.Localization where++--------------------------------------------------------------------------------++import Control.Monad++import Data.List+import Data.Ratio++import Math.Combinat.Numbers+import Math.Combinat.Sets+import Math.Combinat.Sign+import Math.Combinat.Partitions++import qualified Data.Map as Map++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra+import Math.RootLoci.Classic++--------------------------------------------------------------------------------++-- | The localization formula as a string which Mathematica can parse+localizeMathematica :: Partition -> String+localizeMathematica (Partition xs) = formula where+  n   = sum xs+  d   = length xs+  ies = [ (head ys, length ys) | ys <- group (sort xs) ]+  es  = map snd ies++  paren str = '(' : str ++ ")"+  wt j = paren $ show j ++ "a+" ++ show (n-j) ++ "b"++  sumOver = listTensor [ [0..e] | e<-es ] +  formula = global ++ " * " ++ paren (intercalate " + " (map term sumOver)) ++  global = intercalate "*" [ wt j | j<-[0..n] ] ++ " / (b-a)^" ++ show d++  rkonst ss = product [ factorial s * factorial (e-s) | (s,e) <- zip ss es ]+  konst  ss = show (paritySignValue (sum ss)) +            ++ "/" ++ show (rkonst ss)              +  denom  ss = show n ++ "*a - " ++ show (sum [ i*s | ((i,e),s) <- zip ies ss ]) ++ "*(a-b)"+  term   ss = konst ss ++ " / " ++ paren (denom ss)++--------------------------------------------------------------------------------++-- | The localization formula evaluated at given values of @a@ and @b@+localizeEval :: Fractional q => Partition -> q -> q -> q+localizeEval (Partition xs) a b = formula where+  n   = sum xs+  d   = length xs+  ies = [ (head ys, length ys) | ys <- group (sort xs) ]+  es  = map snd ies++  wt j = fromIntegral j * a + fromIntegral (n-j) * b++  sumOver = listTensor [ [0..e] | e<-es ] +  formula = global * sum (map term sumOver)++  global = product [ wt j | j<-[0..n] ] / (b-a)^d++  rkonst ss = product [ factorial s * factorial (e-s) | (s,e) <- zip ss es ]+  konst  ss = fromIntegral (paritySignValue (sum ss)) +            / fromIntegral (rkonst ss)              +  denom  ss = fromIntegral n * a - fromIntegral (sum [ i*s | ((i,e),s) <- zip ies ss ]) * (a-b)+  term   ss = konst ss / denom ss++--------------------------------------------------------------------------------++-- | The dual class, recovered via from the localization formula via Lagrange+-- interpolation+localizeDual :: Partition -> ZMod AB+localizeDual part = ZMod.mapBase worker $ localizeInterpolatedZ part where+  c = codim part+  worker (X i) = AB (c-i) i ++-- | We can use Lagrange interpolation to express the dual class from the+-- localization formula (since we know a priori that the result is a homogeneous+-- polynomial in @a@ and @b@)+--+localizeInterpolatedQ :: Partition -> QMod X+localizeInterpolatedQ part@(Partition xs) = final where+  codim = sum xs - length xs+  bs = map fromIntegral [ 2..codim+2 ]    :: [Rational]+  qs = [ localizeEval part 1 b | b<-bs ] :: [Rational]+  final = lagrangeInterp' (zip bs qs)++localizeInterpolatedZ :: Partition -> ZMod X+localizeInterpolatedZ = (ZMod.mapCoeff f . localizeInterpolatedQ) where+  f :: Rational -> Integer+  f q = case denominator q of+          1 -> numerator q+          _ -> error "non-integral coefficient in the result"++--------------------------------------------------------------------------------++{-+main = do+  forM_ (partitions 9) $ \part@(Partition xs) -> do+    putStrLn $ "X" ++ concatMap show xs ++ " = Factor[ " ++ tp_local_mathematica part ++ " ]"+-}
+ src/Math/RootLoci/Dual/Restriction.hs view
@@ -0,0 +1,223 @@+++-- | Formula for the dual cohomology class of the /cones/ over the strata (sometimes called Thom polynomial) +-- in terms of the Chern classes @c1@ and @c2@, from the author's MSc thesis.+--+-- Note that the dual class agress with the lowest degree part of the CSM class.+--+-- See: Balazs Komuves: Thom Polynomials via Restriction Equations; MSc thesis, ELTE, 2003+--++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances, ScopedTypeVariables #-}+module Math.RootLoci.Dual.Restriction where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio++import Control.Monad++import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set+import Math.Combinat.Sets+import Math.Combinat.Tuples++import qualified Data.Set as Set ; import Data.Set (Set)++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra+import Math.RootLoci.Classic+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++--------------------------------------------------------------------------------+-- * The dual class++-- | The affine Thom polynomial formula from my MSc thesis+affineDualMSc :: Partition -> ZMod Chern+affineDualMSc part@(Partition ps) = ++  case ps of+    []            -> error "affine_tp_msc: empty partition"+    [n]           -> ZMod.fromList [ ( Chern (n-d-2*j) j , rat2int $ single j ) | j<-[ 0 .. div (n-d) 2] ] +    [a,b] | a==b  -> ZMod.fromList [ ( Chern (n-d-2*j) j , rat2int $ double j ) | j<-[ 0 .. div (n-d) 2] ] +    otherwise     -> ZMod.fromList [ ( Chern (n-d-2*j) j , rat2int $ lambda j ) | j<-[ 0 .. div (n-d) 2] ] ++  where++    n = sum ps+    d = length ps++    p = div  n    2+    q = div (n-1) 2++    rat2int r = case denominator r of+      1 -> numerator r+      _ -> error "lambda_j: not integer"++    lambda j = (fi n / 2)^(n-2*q) * fi (doubleFactorial (n-2))^2 * s where+      s = sum+        [ negateIfOdd (n + p + j + lpsi) $ bigTheta j nphi * (fi (2*nphi-n) / fi n)^(d-2) / (fi $ aut phi * aut psi)+        | (phi,psi) <- Set.toList (divideIntoTwoNonEmpty part)+        , let nphi = sum $ fromPartition phi+        , let npsi = sum $ fromPartition psi+        , let lphi = length $ fromPartition phi+        , let lpsi = length $ fromPartition psi+        ] ++    gamma :: Int -> Rational+    gamma k +      | 2*k == n   = 0 +      | otherwise  = fi (k*(k-n)) / fi ((2*k-n)*(2*k-n))++    bigTheta :: Int -> Int -> Rational+    bigTheta j k +      | 2*k == n   = 0 +      | otherwise  = gamma k * smallTheta j k++    smallTheta :: Int -> Int -> Rational+    smallTheta j k = sympoly (q-1-j) [ gamma i | i<-[1..q] , i/=k, i/=n-k ]+   +    fi :: Integral a => a -> Rational+    fi = fromIntegral++    sqj :: Int -> Rational+    sqj j = sympoly (q-j) [ gamma i | i<-[1..q] ]++    sympoly :: Int -> [Rational] -> Rational+    sympoly k xs = sum [ product ys | ys <- choose k xs ]++    -- S(n)+    single j = fi (factorial n) / (product [ gamma i | i<-[1..q] ])+             * negateIfOdd j (sqj j) ++    -- S(p,p)+    double j = fi (doubleFactorial n)^2 / 4 +             * negateIfOdd (q+j) (sqj j) +++--------------------------------------------------------------------------------+-- * Degree++-- | Compute the projective degree from the affine equivariant dual +-- (which can be checked against Hilbert's formula)+-- +-- This is just a simple substition:+--+-- > alpha  ->  1/n+-- > beta   ->  1/n+--+-- or in terms of Chern classes:+--+-- > c1     ->  2/n+-- > c2     ->  1/n^2+--+projDegreeFromDual+  :: Int             -- ^ number of points = dimension of the projective space @P^n@+  -> ZMod Chern      -- ^ dual class+  -> Integer         -- ^ degree+projDegreeFromDual n zm = fromRat s where ++  s :: Rational+  s = sum [ fromIntegral c * c1^e * c2^f  | (Chern e f, c) <- ZMod.toList zm ]++  c1 = 2 / fromIntegral  n    :: Rational+  c2 = 1 / fromIntegral (n*n) :: Rational++-- | Compute the degree of the strata via the formula for the dual class+degreeMSc :: Partition -> Integer+degreeMSc part = projDegreeFromDual (partitionWeight part) (affineDualMSc part)++{-++check_msc_degree :: Bool+check_msc_degree = and+  [ msc_degree part == hilbert part | n<-[1..12] , part <- partitions n ]+-}++--------------------------------------------------------------------------------+-- * extract the dual class from the CSM class ++-- | The dual class of the closure agress with the lowest degree part of the CSM class.+dualClassFromProjCSM :: forall base. ChernBase base => ZMod (Gam base) -> ZMod base+dualClassFromProjCSM csm = dualClassFromAffCSM (ZMod.filterBase nogamma csm) where+  nogamma :: Gam base -> Maybe base+  nogamma (Gam k ab) = if k==0 then Just ab else Nothing++dualClassFromAffCSM :: ChernBase base => ZMod base -> ZMod base+dualClassFromAffCSM csm = filterGrade min_degree csm where+  min_degree = minimum $ map grade $ map fst $ ZMod.toList csm++--------------------------------------------------------------------------------+-- * Lemma 9.1.3++{-+test_lemma_913 = and+  [ lemma913 p h +  | n<-[1..10] +  , p@(Partition ps)<-partitions n+  , let d=length ps+  , h<-[0..d]+  ]++test_lemma_913' =  +  [ (lemma913' p h,(p,h),(d,n))+  | n<-[1..10] +  , p@(Partition ps)<-partitions n+  , let d=length ps+  , h<-[0..d]+  ]+-}++-- | Checks if Lemma 9.1.3 from the thesis is true for the given inputs+lemma913 :: Partition -> Int -> Bool+lemma913 part h = (a==b) where +  (a,b) = lemma913' part h ++  lemma913' :: Partition -> Int -> (Rational, Rational)+  lemma913' part@(Partition ps) h = ( lhs , rhs ) where++    n = sum ps+    d = length ps++    rhs | h == d  = tr (factorial d) * product (map fi ps)+        | h <  d  = 0+        | h >  d  = -666++    lhs = sum+      [ negateIfOdd (length rs) $  (fi (2 * sum qs - n) / 2)^h * (tr $ aut part) / (tr $ aut phi * aut psi)+      | ( phi@(Partition qs) , psi@(Partition rs) ) <- Set.toList (divideIntoTwo part)+      ]++    fi :: Int -> Rational+    fi = fromIntegral++    tr :: Integer -> Rational+    tr = fromIntegral  +++--------------------------------------------------------------------------------+-- * helper functions++-- | Different ways to divide a partition into two +divideIntoTwo :: Partition -> Set (Partition,Partition)+divideIntoTwo (Partition ps) = Set.fromList $ map f (binaryTuples d) where++  d    = length ps+  f ts = ( g ts , g (map not ts) )+  g ts = Partition [ k | (b,k) <- zip ts ps , b ]++  -- nonempty (p,q) = not (isEmptyPartition p) && not (isEmptyPartition q)++-- | Different ways to divide a partition into two /nonempty/ partitions+divideIntoTwoNonEmpty :: Partition -> Set (Partition,Partition)+divideIntoTwoNonEmpty p = Set.delete x $ Set.delete y $ divideIntoTwo p where+  x = (emptyPartition,p)+  y = (p,emptyPartition)++--------------------------------------------------------------------------------+
+ src/Math/RootLoci/Geometry.hs view
@@ -0,0 +1,13 @@++-- | Re-exporting the Geometry.* modules++module Math.RootLoci.Geometry+  ( module Math.RootLoci.Geometry.Cohomology+  , module Math.RootLoci.Geometry.Forget+  , module Math.RootLoci.Geometry.Mobius+  )+  where++import Math.RootLoci.Geometry.Cohomology+import Math.RootLoci.Geometry.Forget+import Math.RootLoci.Geometry.Mobius
+ src/Math/RootLoci/Geometry/Cohomology.hs view
@@ -0,0 +1,279 @@++-- | Bases in the cohomology of the spaces appearing in the computations.+--+-- We have three different spaces: +--+-- * @Q^n = P^1 x P^1 x ... x P^1@ (@n@ times; @m = length lambda@)+--+-- * @Q^m = P^1 x P^1 x ... x P^1 x P^1@ (@m@ times, @m = sum lambda >= n@)+-- +-- * @P^m = P(Sym^m C^2)@+--+-- Furthermore, we have @GL2@ acting naturally on these spaces.+--++{-# LANGUAGE +      BangPatterns, TypeSynonymInstances, FlexibleInstances, DeriveFunctor, +      ScopedTypeVariables, Rank2Types +  #-}+module Math.RootLoci.Geometry.Cohomology where++--------------------------------------------------------------------------------++import Data.List+import Data.Monoid++import Math.Combinat.Numbers++import qualified Data.Map as Map+import qualified Data.Set as Set++import qualified Math.RootLoci.Algebra.FreeMod as ZMod+import Math.RootLoci.Algebra.FreeMod ( ZMod , FreeMod(..) , unFreeMod )++import Math.RootLoci.Algebra.SymmPoly +import Math.RootLoci.Misc.Pretty++--------------------------------------------------------------------------------+-- * The non-equivariant case++-- | a (ring) generator of @H^*(Q^n)@ (note that @u_i^2 = 0@)+newtype U = U Int deriving (Eq,Ord,Show)++-- | (a ring) generator of @H^*(Q^m)@ (note that @h_i^2 = 0@)+newtype H = H Int deriving (Eq,Ord,Show)++-- | the generator of @H^*(P^n)@ (with @g^(n+1) = 0@)+newtype G = G Int deriving (Eq,Ord,Show)++-- | A monomial in @u_i@ (encoded as a subset of @[1..n]@, as @u_i^2=0@)+newtype US = US [U] deriving (Eq,Ord,Show)++-- | A monomial in @h_i@ (encoded as a subset of @[1..m]@, as @h_i^2=0@)+newtype HS = HS [H] deriving (Eq,Ord,Show)++--------------------------------------------------------------------------------++instance Monoid US where+  mempty = US []+  (US us1) `mappend` (US us2) = +    if nub us3 == us3+      then US us3+      else error "[U]/monoid: duplicate indices"+    where+      us3 = sort (us1 ++ us2)++instance Monoid HS where+  mempty = HS []+  (HS hs1) `mappend` (HS hs2) = +    if nub hs3 == hs3+      then HS hs3+      else error "[H]/monoid: duplicate indices"+    where+      hs3 = sort (hs1 ++ hs2)++instance Monoid G where+  mempty = G 0+  (G e) `mappend` (G f) = G (e+f)+ +--------------------------------------------------------------------------------++instance Pretty G where+  pretty (G e) = "g^" ++ show e++instance Pretty H where+  pretty (H i) = "h" ++ show i++instance Pretty U where+  pretty (U i) = "u" ++ show i++instance Pretty HS where+  pretty (HS []) = ""+  pretty (HS hs) = intercalate "*" (map pretty hs)++instance Pretty US where+  pretty (US []) = ""+  pretty (US us) = intercalate "*" (map pretty us)++--------------------------------------------------------------------------------++instance Graded U where grade _ = 1+instance Graded H where grade _ = 1+instance Graded G where grade (G g) = g+instance Graded HS where grade (HS js) = length js+instance Graded US where grade (US js) = length js++instance Graded ab => Graded (Omega ab) where grade (Omega us ab) = length us + grade ab+instance Graded ab => Graded (Eta   ab) where grade (Eta   hs ab) = length hs + grade ab+instance Graded ab => Graded (Gam   ab) where grade (Gam   g  ab) = g + grade ab++--------------------------------------------------------------------------------+-- * The equivariant case++-- | A monomial generator of @Z[alpha,beta;u1,u2,...,u_nd]/(...)@, +-- the cohomology ring of @Q^n@. +--+-- The encoding is that the list is the list of indices of @u@ which appear.+data Omega ab = Omega ![Int] !ab deriving (Eq,Ord,Show,Functor)++-- | A monomial generator of @Z[alpha,beta;eta1,eta2...eta_m]/(...)@,+-- he cohomology ring of @Q^m@. +--+-- The encoding is that the list is the list of indices of @eta@ which appear.+data Eta ab = Eta ![Int] !ab deriving (Eq,Ord,Show,Functor)++-- | A monomial generator of @Z[alpha,beta;gamma]/(...)@,+-- the cohomology ring of @P^m@. +data Gam ab = Gam !Int !ab deriving (Eq,Ord,Show,Functor)++--------------------------------------------------------------------------------++-- | Class of monomial bases which form modules over the @H^*(BGL2)@+class Functor f => Equivariant f where +  injectMonom  :: x -> f x+  projectMonom :: f x -> x++instance Equivariant Omega where +  injectMonom = Omega [] +  projectMonom (Omega _ ab) = ab++instance Equivariant Eta where +  injectMonom = Eta [] +  projectMonom (Eta _ ab) = ab++instance Equivariant Gam where  +  injectMonom = Gam 0  +  projectMonom (Gam _ ab) = ab++injectZMod :: (Equivariant f, ChernBase base, Ord (f base)) => ZMod base -> ZMod (f base)+injectZMod = ZMod.mapBase injectMonom++forgetGamma :: Ord base => ZMod (Gam base) -> ZMod base +forgetGamma = ZMod.filterBase f where+  f (Gam k ab) = case k of+    0 -> Just ab+    _ -> Nothing++forgetEquiv :: ChernBase base => ZMod (Gam base) -> ZMod G+forgetEquiv = ZMod.filterBase f where+  f (Gam k ab) = if (ab == mempty) +    then Just (G k)+    else Nothing++--------------------------------------------------------------------------------+-- * Conversion between different bases++convertOmega   +  :: (Ord ab, Ord cd)+  => (ZMod ab -> ZMod cd) +  -> ZMod (Omega ab) -> ZMod (Omega cd)+convertOmega = convertEach f g Omega where+  f (Omega xs _ ) = xs+  g (Omega _  ab) = ab++convertEta+  :: (Ord ab, Ord cd)+  => (ZMod ab -> ZMod cd) +  -> ZMod (Eta ab) -> ZMod (Eta cd)+convertEta = convertEach f g Eta where+  f (Eta xs _ ) = xs+  g (Eta _  ab) = ab++convertGam+  :: (Ord ab, Ord cd)+  => (ZMod ab -> ZMod cd) +  -> ZMod (Gam ab) -> ZMod (Gam cd)+convertGam = convertEach f g Gam where+  f (Gam k _ ) = k+  g (Gam _ ab) = ab++-- | A generic function which can convert the @GL2@ representations+convertEach +  :: forall f x y ab cd. (Functor f, Ord ab, Ord cd, Ord (f ab), Ord (f cd), Ord x) +  => (forall y. f y -> x)+  -> (forall y. f y -> y)+  -> (forall y. x -> y -> f y)+  -> (ZMod    ab  -> ZMod    cd )+  ->  ZMod (f ab) -> ZMod (f cd)+convertEach selx sely build convert src = tgt where+  tgt    = ZMod.sum [ worker layer | layer <- layers ]+  layers = Set.toList $ Set.map selx $ Map.keysSet $ unFreeMod src :: [x]+  worker layer +    = FreeMod+    $ Map.mapKeys (build layer)+    $ unFreeMod+    $ convert+    $ FreeMod+    $ Map.mapKeys sely +    $ Map.filterWithKey (\k _ -> selx k == layer) +    $ unFreeMod src++--------------------------------------------------------------------------------++-- | This is a hack to reuse the same pushforward code+unsafeEtaToOmega :: Ord ab => FreeMod coeff (Eta ab) -> FreeMod coeff (Omega ab)+unsafeEtaToOmega = ZMod.mapBase f where+  f (Eta js ab) = Omega js ab++unsafeOmegaToEta :: Ord ab => FreeMod coeff (Omega ab) -> FreeMod coeff (Eta ab)+unsafeOmegaToEta = ZMod.mapBase f where+  f (Omega js ab) = Eta js ab++--------------------------------------------------------------------------------++instance Monoid ab => Monoid (Omega ab) where+  mempty = Omega [] mempty+  (Omega as ab1) `mappend` (Omega bs ab2) = +    if nub cs == cs+      then Omega cs (ab1 <> ab2)+      else error "Omega/monoid: duplicate indices"+    where+      cs = sort (as ++ bs)++instance Monoid ab => Monoid (Eta ab) where+  mempty = Eta [] mempty+  (Eta fs ab1) `mappend` (Eta gs ab2) = +    if nub hs == hs+      then Eta hs (ab1 <> ab2)+      else error "Eta/monoid: duplicate indices"+    where+      hs = sort (fs ++ gs)++instance Monoid ab => Monoid (Gam ab) where+  mempty = Gam 0 mempty+  (Gam e ab1) `mappend` (Gam f ab2) = Gam (e+f) (ab1 <> ab2)++--------------------------------------------------------------------------------++instance (Pretty ab, Monoid ab, Eq ab) => Pretty (Gam ab) where+  pretty (Gam 0 ab) = pretty ab+  pretty (Gam g ab)+    | ab == mempty  = "g^" ++ show g+    | otherwise     = "g^" ++ show g ++ "*" ++ pretty ab++instance (Pretty ab, Monoid ab, Eq ab) => Pretty (Eta ab) where+  pretty eta = +    case eta of+      (Eta [] ab)       -> pretty ab +      (Eta is ab)   +        | ab == mempty  -> hs is+        | otherwise     -> hs is ++ "*" ++ pretty ab +    where+      hs is = case is of+        [] -> ""+        _  -> intercalate "*" [ "h" ++ show i | i<-is ]++instance (Pretty ab, Monoid ab, Eq ab) => Pretty (Omega ab) where+  pretty omega = +    case omega of+      (Omega [] ab)       -> pretty ab +      (Omega is ab)    +        | ab == mempty    -> us is+        | otherwise       -> us is ++ "*" ++ pretty ab +    where+      us is = case is of+        [] -> ""+        _  -> intercalate "*" [ "u" ++ show i | i<-is ]++--------------------------------------------------------------------------------+
+ src/Math/RootLoci/Geometry/Forget.hs view
@@ -0,0 +1,111 @@++-- | Geometry of the degree @n!@ finite map @pi@, which just forgets the order points:+--+-- > pi : Q^n = P^1 x P^1 x ... x P^1  ->  P^n = P(Sym^n C^2)+--+-- It's clear that the degree of @pi@ restricted to an open stratum corresponding to+-- a partition @mu@ is the multinomial coefficient corresponding to @n `choose` mu@.+--+-- It is also not hard to see that the degree of @pi@ restricted to the intersection+-- of the open stratum corresponding to @mu@ with the image of the diagonal map +-- corresponding to @nu@ equals the number of \"automorphisms\" @aut(mu) = prod (e_i!)@+-- where @mu = (1^e1 2^e2 ... k^ek)@ and the number of ways @nu@ is refinement of @mu@.+--+-- Note that for @nu=(1,1...1)@ the multinomial agrees with the number of refinements.+--+-- This module contains functions to compute these numbers.+--++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}+module Math.RootLoci.Geometry.Forget where++--------------------------------------------------------------------------------++import Data.List++import Control.Monad++import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set+import Math.Combinat.Sets++import qualified Data.Map.Strict as Map+import Data.Map.Strict (Map)++import qualified Data.Map.Lazy as LMap++import qualified Data.Set as Set+import Data.Set (Set)++import Data.Array.IArray+import Data.Array.Unboxed+import Data.Array (Array)++import Math.RootLoci.Misc.Common+import Math.RootLoci.Misc.PTable+-- import Math.RootLoci.Geometry.Mobius++--------------------------------------------------------------------------------++-- | Given a partition, we list all coarser partitions together+-- with the number of ways the input is a refinement of the+-- coarser partition.+--+-- TODO: at the moment this is just a synonym for 'countCoarseningsNaive' ...+--+countCoarsenings :: Partition -> Map Partition Integer+countCoarsenings = countCoarseningsNaive++-- | Count coarsenings (with multiplicities) which are shorter by just 1.+countDirectCoarsenings :: Partition -> Map Partition Integer+countDirectCoarsenings part = Map.fromListWith (+) list where+  list =  +    [ ( fromExponentialFrom ((i1+i2,1):(i1,e1-1):(i2,e2-1):rest) , fromIntegral (e1*e2) )+    | ( (i1,e1):(i2,e2):[] , rest ) <- choose' 2 ies+    ] +++    [ ( fromExponentialFrom ((2*i,1):(i,e-2):rest) , binomial e 2 )+    | ( (i,e):[] , rest ) <- choose' 1 ies+    , e >= 2+    ]+  ies = toExponentialForm part +  +--------------------------------------------------------------------------------++-- | Naive (very slow) implementation of 'countCoarsenings'.+countCoarseningsNaive :: Partition -> Map Partition Integer+countCoarseningsNaive = pcache count where++  count (Partition ps) = histogram (map f setps) where+    d     = length ps+    setps = map fromSetPartition $ setPartitions d :: [[[Int]]]+    arr   = listArray (1,d) ps :: UArray Int Int+    f iss = mkPartition [ sum [ arr ! k | k <- is ] | is <- iss ]++-- | Given a partition @nu@, we stratify the image of the +-- corresponding diagonal @Delta_nu@ as usual, and list+-- the degree of @pi@ restricted to these strata+--+-- This is just counting the coarsenings, multiplied by+-- the number of \"automorphisms\" of the partition.+--+countPreimage :: Partition -> Map Partition Integer+countPreimage = pcache compute where+  compute part = Map.mapWithKey f (countCoarsenings part) +  f q c = c * aut q++-- | The preimage counts, but the partition itself is separated out.+preimageView :: Partition -> (Integer, Map Partition Integer)+preimageView part = unsafeDeleteLookup part (countPreimage part) ++--------------------------------------------------------------------------------++-- | The preimage @pi^-1(x)@ of a point under the map +-- @pi : Q^n -> P^n@ is just a multinomial coefficient+countFullPreimage :: Partition -> Integer+countFullPreimage part@(Partition ps) = multinomial ps +  +--------------------------------------------------------------------------------++
+ src/Math/RootLoci/Geometry/Mobius.hs view
@@ -0,0 +1,224 @@++-- | Mobius inversion for the coarsening poset of partitions++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-}+module Math.RootLoci.Geometry.Mobius +  ( Partition(..) +  -- * The refinement poset of partitions+  , coarserThan , finerThan+  , (.==.) , (./=.) , (.<=.) , (.>=.) , (.<.) , (.>.) +  -- * closures+  , fastClosure , fastAntiClosure+  , closureSet , closureSet'+  -- * Mobius function+  , zetaOf , mobiusOf+  -- * helpers+  , firstLevelDown , firstLevelUp  +  -- * set partitions+  , closureSetOfSetPartition+  , firstLevelDownSetP+  )+  where++--------------------------------------------------------------------------------++import Data.List++import qualified Data.Map.Strict as Map ; import Data.Map.Strict (Map)+import qualified Data.Set        as Set ; import Data.Set        (Set)++import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set+import Math.Combinat.Sets++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra+import Math.RootLoci.Misc++--------------------------------------------------------------------------------++{-+indicator :: Bool -> Integer+indicator b = if b then 1 else 0++kronecker' :: Partition -> ZMod Partition+kronecker' p = ZMod.singleton p 1++kronecker :: Partition -> Partition -> Integer+kronecker p q = indicator (p .==. q)++zeta :: Partition -> Partition -> Integer+zeta p q = indicator (p .<=. q)+-}++--------------------------------------------------------------------------------+-- * Mobius function++-- | Zeta function of the refinement poset+zetaOf :: Partition -> ZMod Partition+zetaOf = pcache calc where+  calc p = ZMod.fromList $ map (\p -> (p,1)) $ Set.toList $ closureSet p++-- | Mobius function of the refinement poset+mobiusOf :: Partition -> ZMod Partition+mobiusOf = pcache calc where+  calc    p = ZMod.sub (ZMod.singleton p 1) (smaller p)+  smaller p = ZMod.sum [ mobiusOf q | q <- Set.toList (closureSet' p) ]++--------------------------------------------------------------------------------+-- * The refinement poset of partitions++coarserThan :: Partition -> Partition -> Bool+coarserThan p q = Set.member p (closureSet q)++finerThan :: Partition -> Partition -> Bool+finerThan q p = coarserThan p q++(.<=.) :: Partition -> Partition -> Bool+(.<=.) = coarserThan++(.>=.) :: Partition -> Partition -> Bool+(.>=.) = finerThan++(.==.) :: Partition -> Partition -> Bool+(.==.) = (==)++(./=.) :: Partition -> Partition -> Bool+(./=.) = (/=)++(.<.) :: Partition -> Partition -> Bool+(.<.) p q = (p .<=. q) && (p /= q) ++(.>.) :: Partition -> Partition -> Bool+(.>.) p q = (p .>=. q) && (p /= q) ++--------------------------------------------------------------------------------+-- | Efficient first level merge/split++insertRevSorted :: Int -> [Int] -> [Int]+insertRevSorted x = go where+  go yys@(y:ys) = if x >= y then x : yys else y : go ys+  go []         = [x]++insertRevSorted2 :: Int -> Int -> [Int] -> [Int]+insertRevSorted2 x y = insertRevSorted x . insertRevSorted y++-- | Example: +-- +-- > insertGroup [3,3] [[5,5,5],[4],[1,1,1,1]] == [5,5,5,4,3,3,1,1,1,1]+--+insertGroup_ :: [Int] -> [[Int]] -> [Int]+insertGroup_ zs@(z:_) = go where+  go (xs@(x:_):rest) = if z >= x then zs ++ xs ++ concat rest +                                 else xs ++ go rest+  go ([]      :rest) = go rest+  go []              = zs+insertGroup_ [] = concat++-- | These satisfy:+--+-- > concat . insertGroup what == insertGroup_ what+--+insertGroup :: [Int] -> [[Int]] -> [[Int]]+insertGroup zs@(z:_) = go where+  go (xs@(x:_):rest) = if z >= x then zs : xs : rest +                                 else xs : go rest+  go ([]      :rest) = go rest+  go []              = [zs]+insertGroup [] = id++insertGroup2_ :: [Int] -> [Int] -> [[Int]] -> [Int]+insertGroup2_ xs ys = insertGroup_ xs . insertGroup ys++insertGroup2 :: [Int] -> [Int] -> [[Int]] -> [[Int]]+insertGroup2 xs ys = insertGroup xs . insertGroup ys++choose1 :: [a] -> [(a,[a])]+choose1 (x:xs) = (x,xs) : [ (y,x:ys) | (y,ys) <- choose1 xs ]+choose1 []     = []++choose2 :: [a] -> [(a,a,[a])]+choose2 (x:xs) =  [ (x,y,ys  ) |   (y,ys) <- choose1 xs ]+               ++ [ (y,z,x:zs) | (y,z,zs) <- choose2 xs ]+choose2 []     =  []++-- | Merging two parts+firstLevelDown :: Partition -> [Partition]+firstLevelDown (Partition ps) = (one ++ two) where+  gs  = group ps+  one = [ Partition $ insertRevSorted (x+y) (insertGroup_  zs    rest) | ((x:y:zs)     ,rest) <- choose1 gs ]+  two = [ Partition $ insertRevSorted (x+y) (insertGroup2_ xs ys rest) | ((x:xs),(y:ys),rest) <- choose2 gs ]++-- | Splitting one part into two+firstLevelUp :: Partition -> [Partition]+firstLevelUp (Partition ps) = one where+  gs  = group ps+  one = [ Partition $ insertRevSorted2 x (z-x) (insertGroup_  zs rest) | ((z:zs),rest) <- choose1 gs , x<-[1..div z 2] ]++-- | Sanity check+firstLevelDownNaive :: Partition -> [Partition]+firstLevelDownNaive (Partition ps) = unique [ mkPartition ( x+y : zs ) | ([x,y],zs) <- choose' 2 ps ]++firstLevelUpNaive :: Partition -> [Partition]+firstLevelUpNaive (Partition ps) = unique [ mkPartition ( x : z-x : zs ) | ([z],zs) <- choose' 1 ps , x<-[1..z-1] ]++checkDown :: Partition -> Bool+checkDown p = (sort (firstLevelDown p) == firstLevelDownNaive p)++checkUp :: Partition -> Bool+checkUp p = (sort (firstLevelUp p) == firstLevelUpNaive p)++--------------------------------------------------------------------------------++-- | Fast computation of a single closure+fastClosure :: Partition -> Set Partition+fastClosure p = go Set.empty [p] where+  go !acc (p:ps) = case Set.member p acc of+    True  -> go acc ps+    False -> go (Set.insert p acc) (firstLevelDown p ++ ps)+  go !acc []     = acc++-- | Fast computation of a single \"anticlosure\" (opposite poset)+fastAntiClosure :: Partition -> Set Partition+fastAntiClosure p = go Set.empty [p] where+  go !acc (p:ps) = case Set.member p acc of+    True  -> go acc ps+    False -> go (Set.insert p acc) (firstLevelUp p ++ ps)+  go !acc []     = acc++--------------------------------------------------------------------------------++-- | Caches and reuses all closures (lazily), this is the fastest version+closureSet :: Partition -> Set Partition +closureSet = cached where+  cached = monoCache calc +  calc p = go (Set.singleton p) (firstLevelDown p) where+    go !acc (p:ps) = case Set.member p acc of+      True  -> go acc ps+      False -> go (Set.union acc (cached p)) ps+    go !acc []     = acc++-- | The closure without the stratum itself+closureSet' :: Partition -> Set Partition+closureSet' p = Set.delete p (closureSet p)++--------------------------------------------------------------------------------+-- * set partitions++firstLevelDownSetP :: SetPartition -> [SetPartition]+firstLevelDownSetP (SetPartition ps) =+  [ toSetPartition ( (x++y) : zs ) | ([x,y],zs) <- choose' 2 ps ]+  +closureSetOfSetPartition :: SetPartition -> Set SetPartition  +closureSetOfSetPartition = cached where+  cached = monoCache calc+  calc p = go (Set.singleton p) (firstLevelDownSetP p) where+    go !acc (p:ps) = case Set.member p acc of+      True  -> go acc ps+      False -> go (Set.union acc (cached p)) ps+    go !acc []     = acc+ +--------------------------------------------------------------------------------+ +
+ src/Math/RootLoci/Misc.hs view
@@ -0,0 +1,13 @@++-- | Re-exporting the Misc.* modules++module Math.RootLoci.Misc +  ( module Math.RootLoci.Misc.Common +  , module Math.RootLoci.Misc.PTable+  , module Math.RootLoci.Misc.Pretty +  )+  where++import Math.RootLoci.Misc.Common +import Math.RootLoci.Misc.PTable+import Math.RootLoci.Misc.Pretty
+ src/Math/RootLoci/Misc/Common.hs view
@@ -0,0 +1,132 @@++-- | Some auxilary functions++{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances, DeriveFunctor #-}+module Math.RootLoci.Misc.Common where++--------------------------------------------------------------------------------++import Data.List+import Data.Monoid+import Data.Ratio++import Control.Monad++import Math.Combinat.Numbers+import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer +import Math.Combinat.Partitions.Set+import Math.Combinat.Sets++import qualified Data.Map.Strict as Map+import Data.Map (Map)++-- import qualified Math.RootLoci.Algebra.FreeMod as ZMod+-- import Math.RootLoci.Algebra.SymmPoly+-- import Math.RootLoci.Geometry.Cohomology++--------------------------------------------------------------------------------+-- * Pairs++data Pair a +  = Pair a a +  deriving (Eq,Ord,Show,Functor)++--------------------------------------------------------------------------------+-- * Lists++{-# SPECIALIZE unique :: [Partition] -> [Partition] #-}+unique :: Ord a => [a] -> [a]+unique = map head . group . sort++-- | Synonym for histogram+count :: Ord b => [b] -> Map b Integer+count = histogram++{-# SPECIALIZE histogram :: [Partition] -> Map Partition Integer #-}+histogram :: Ord b => [b] -> Map b Integer+histogram xs = foldl' f Map.empty xs where+  f old x = Map.insertWith (+) x 1 old++--------------------------------------------------------------------------------+-- * Maps+  +deleteLookup :: Ord a => a -> Map a b -> (Maybe b, Map a b)+deleteLookup k table = (Map.lookup k table, Map.delete k table)  ++unsafeDeleteLookup :: Ord a => a -> Map a b -> (b, Map a b)+unsafeDeleteLookup k table = (fromJust (Map.lookup k table), Map.delete k table) where+  fromJust mb = case mb of+    Just y  -> y+    Nothing -> error "unsafeDeleteLookup: key not found"++--------------------------------------------------------------------------------+-- * Partitions++-- | @aut(mu)@ is the number of symmetries of the partition mu:+--+-- > aut(mu) = prod_r (e_r)!+--+-- where @mu = (1^e1 2^e2 .. k^ek)@+aut :: Partition -> Integer+aut part = product $ map factorial es where+  es = map snd $ toExponentialForm part   ++--------------------------------------------------------------------------------+-- * Set partitions+ +-- | Makes set partition from a partition (simply filling up from left to right)+-- with the shape giving back the input partition+defaultSetPartition :: Partition -> SetPartition+defaultSetPartition = SetPartition . linearIndices++-- | Produce linear indices from a partition @nu@ (to encode the diagonal map @Delta_nu@).+linearIndices :: Partition -> [[Int]]+linearIndices (Partition ps) = go 0 ps where+  go _  []     = []+  go !k (a:as) = [k+1..k+a] : go (k+a) as++--------------------------------------------------------------------------------+-- * Signs++class IsSigned a where+  signOf :: a -> Maybe Sign++signOfNum :: (Ord a, Num a) => a -> Maybe Sign +signOfNum x = case compare x 0 of+  LT -> Just Minus+  GT -> Just Plus+  EQ -> Nothing++instance IsSigned Int      where signOf = signOfNum+instance IsSigned Integer  where signOf = signOfNum+instance IsSigned Rational where signOf = signOfNum++--------------------------------------------------------------------------------+-- * Numbers++fromRat :: Rational -> Integer+fromRat r = case denominator r of+  1 -> numerator r+  _ -> error "fromRat: not an integer"    ++safeDiv :: Integer -> Integer -> Integer+safeDiv a b = case divMod a b of+  (q,0) -> q+  (q,r) -> error $ "saveDiv: " ++ show a ++ " = " ++ show b ++ " * " ++ show q ++ " + " ++ show r++--------------------------------------------------------------------------------+-- * Combinatorics++-- | Chooses (n-1) elements out of n+chooseN1 :: [a] -> [[a]]+chooseN1 = go where+  go (x:xs) = xs : map (x:) (go xs)+  go []     = []+  +symPolyNum :: Num a => Int -> [a] -> a+symPolyNum k xs = sum' (map prod' $ choose k xs) where+  sum'  = foldl' (+) 0+  prod' = foldl' (*) 1++--------------------------------------------------------------------------------
+ src/Math/RootLoci/Misc/PTable.hs view
@@ -0,0 +1,130 @@++-- | Infinite lazy partition tables (used for caching).+--+-- We cache almost all computations (which would be otherwise typically +-- executed many times); this really helps performance.+--++{-# LANGUAGE Rank2Types #-} +module Math.RootLoci.Misc.PTable where++--------------------------------------------------------------------------------++import Data.List++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set++import qualified Data.Map.Lazy as LMap++import Math.RootLoci.Algebra.SymmPoly++--------------------------------------------------------------------------------+-- * Finite lazy partition tables++newtype PTable a = PTable (LMap.Map Partition a)++createPTable :: (Partition -> a) -> Int -> PTable a+createPTable f n = PTable $ LMap.fromList [ (p, f p) | p <- partitions n ]++lookupPTable :: Partition -> PTable a -> a+lookupPTable p (PTable lmap) = case LMap.lookup p lmap of+  Just y  -> y+  Nothing -> error "lookupPTable"++--------------------------------------------------------------------------------+-- * Infinite lazy partition tables++newtype PSeries a = PSeries [PTable a]+  +createPSeries :: (Partition -> a) -> PSeries a+createPSeries f = PSeries [ createPTable f n | n<-[0..] ]++lookupPSeries :: Partition -> PSeries a -> a+lookupPSeries part (PSeries series) = lookupPTable part (series !! partitionWeight part)++--------------------------------------------------------------------------------+-- * Finite lazy set-partition tables++newtype SetPTable a = SetPTable (LMap.Map SetPartition a)++createSetPTable :: (SetPartition -> a) -> Int -> SetPTable a+createSetPTable f n = SetPTable $ LMap.fromList [ (p, f p) | p <- setPartitions n ]++lookupSetPTable :: SetPartition -> SetPTable a -> a+lookupSetPTable p (SetPTable lmap) = case LMap.lookup p lmap of+  Just y  -> y+  Nothing -> error "lookupSetPTable"++--------------------------------------------------------------------------------+-- * Infinite lazy set-partition tables++newtype SetPSeries a = SetPSeries [SetPTable a]+  +createSetPSeries :: (SetPartition -> a) -> SetPSeries a+createSetPSeries f = SetPSeries [ createSetPTable f n | n<-[0..] ]++lookupSetPSeries :: SetPartition -> SetPSeries a -> a+lookupSetPSeries setp (SetPSeries series) = lookupSetPTable setp (series !! setpWeight setp) where+  setpWeight (SetPartition ps) = foldl' (+) 0 (map length ps)++--------------------------------------------------------------------------------+-- * polymorphic caching ++polyCache1 +  :: (CacheKey key) +  => (forall base. ChernBase base => key -> f base)     -- ^ polymorphic function to be cached+  -> (forall base. ChernBase base => key -> f base)+polyCache1 calc = \key -> select1 (cacheAB key, cacheChern key)  where+  cacheAB    = monoCache $ \k -> spec1' ChernRoot  (calc k)+  cacheChern = monoCache $ \k -> spec1' ChernClass (calc k)++polyCache2 +  :: (CacheKey key) +  => (forall base. ChernBase base => key -> f (g base))     -- ^ polymorphic function to be cached+  -> (forall base. ChernBase base => key -> f (g base))+polyCache2 calc = \key -> select2 (cacheAB key, cacheChern key)  where+  cacheAB    = monoCache $ \k -> spec2' ChernRoot  (calc k)+  cacheChern = monoCache $ \k -> spec2' ChernClass (calc k)++polyCache3 +  :: (CacheKey key) +  => (forall base. ChernBase base => key -> f (g (h base)))     -- ^ polymorphic function to be cached+  -> (forall base. ChernBase base => key -> f (g (h base)))+polyCache3 calc = \key -> select3 (cacheAB key, cacheChern key)  where+  cacheAB    = monoCache $ \k -> spec3' ChernRoot  (calc k)+  cacheChern = monoCache $ \k -> spec3' ChernClass (calc k)++--------------------------------------------------------------------------------+-- * monomorphic caching ++class CacheKey key where+  monoCache :: (key -> a) -> (key -> a)++instance CacheKey Int          where  monoCache = icache+instance CacheKey Partition    where  monoCache = pcache+instance CacheKey SetPartition where  monoCache = setpcache++--------------------------------------------------------------------------------+-- * individual caching functions++pcache :: (Partition -> a) -> (Partition -> a)+pcache calc = lkp where+  lkp p = lookupPSeries p table +  table = createPSeries calc++setpcache :: (SetPartition -> a) -> (SetPartition -> a)+setpcache calc = lkp where+  lkp setp = lookupSetPSeries setp table +  table    = createSetPSeries calc++icache :: (Int -> a) -> (Int -> a)+icache calc = \n -> (table !! n) where+  table = [ calc i | i <- [0..]  ]++icache' :: a -> Int -> (Int -> a) -> (Int -> a)+icache' dflt fstidx calc = \n -> if n < fstidx then dflt else (table !! (n-fstidx)) where+  table = [ calc i | i <- [fstidx..]  ]++--------------------------------------------------------------------------------
+ src/Math/RootLoci/Misc/Pretty.hs view
@@ -0,0 +1,137 @@++{-# LANGUAGE FlexibleInstances #-}++-- | Pretty-printing+ +module Math.RootLoci.Misc.Pretty where++--------------------------------------------------------------------------------++import Data.List++import Math.Combinat.Sign+import Math.Combinat.Partitions.Integer++import qualified Data.Map.Strict as Map+import Data.Map.Strict (Map)++import Math.RootLoci.Algebra.FreeMod ( FreeMod, ZMod, QMod )+import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Misc.Common++--------------------------------------------------------------------------------++class Pretty a where+  pretty :: a -> String++-- instance Pretty a => Pretty (ZMod a) where+--   pretty = prettyZMod pretty++instance (Num c, Eq c, Show c, IsSigned c, Pretty b) => Pretty (FreeMod c b) where+  pretty = prettyFreeMod' True pretty++--------------------------------------------------------------------------------+-- * Pretty printing elements of free modules++-- | Example: @showVarPower "x" 5 == "x^5"@+showVarPower :: String -> Int -> String+showVarPower name expo = case expo of+  0 -> "1"+  1 -> name+  _ -> name ++ "^" ++ show expo++--------------------------------------------------------------------------------++-- | no multiplication sign (ok for mathematica and humans)+prettyZMod_ :: (b -> String) -> ZMod b -> String+prettyZMod_ = prettyFreeMod' False+  +-- | multiplication sign (ok for maple etc)+prettyZMod :: (b -> String) -> ZMod b -> String+prettyZMod = prettyFreeMod' True++--------------------------------------------------------------------------------++prettyFreeMod' +  :: (Num c, Eq c, Show c, IsSigned c) +  => Bool                -- ^ use star for multiplication (@False@ means just concatenation)+  -> (b -> String)       -- ^ show base+  -> FreeMod c b +  -> String+prettyFreeMod' star showBase what = final where+  final = if take 3 stuff == " + " then drop 3 stuff else drop 1 stuff+  stuff = concatMap f (ZMod.toList what) +  f (g,  1) = plus  ++ showBase' g+  f (g, -1) = minus ++ showBase' g+  f (g, c)  = sgn c ++ {- extendStringR 3 -} (show $ abs c) ++ starStr ++ showBase' g+  -- cond (_,c) = (c/=0)+  starStr = if star then "*" else " "+  showBase' g = case showBase g of+    "" -> "1"  -- "(1)"+    s  -> s+  sgn c = case signOf c of+    Just Minus -> minus+    _          -> plus+  plus  = " + "+  minus = " - "++prettyFreeMod'' +  :: (c -> String)    -- ^ show coefficient+  -> (b -> String)    -- ^ show base+  -> FreeMod c b +  -> String+prettyFreeMod'' showCoeff showBase what = result where+  result = intercalate " + " (map f $ ZMod.toList what) +  f (g, c) = showCoeff c ++ starStr ++ showBase' g+  starStr = "*" -- if star then "*" else " "+  showBase' g = case showBase g of+    "" -> "1"  -- "(1)"+    s  -> s++--------------------------------------------------------------------------------+-- * Utility++-- | Put into parentheses+paren :: String -> String+paren s = '(' : s ++ ")"++--------------------------------------------------------------------------------++-- | Exponential form of a partition+expFormString :: Partition -> String+expFormString p = "(" ++ intercalate "," (map f ies) ++ ")" where+  ies = toExponentialForm p+  f (i,e) = show i ++ "^" ++ show e++extendStringL :: Int -> String -> String+extendStringL k s = s ++ replicate (k - length s) ' '++extendStringR :: Int -> String -> String+extendStringR k s = replicate (k - length s) ' ' ++ s++--------------------------------------------------------------------------------+-- * Mathematica-formatted output++class Mathematica a where+  mathematica :: a -> String++instance Mathematica Int where+  mathematica = show++instance Mathematica Integer where+  mathematica = show++instance Mathematica String where+  mathematica = show++instance Mathematica Partition where+  mathematica (Partition ps) = "{" ++ intercalate "," (map show ps) ++ "}"++data Indexed a = Indexed String a++instance Mathematica a => Mathematica (Indexed a) where+  mathematica (Indexed x sub) = "Subscript[" ++ x ++ "," ++ mathematica sub ++ "]"++--------------------------------------------------------------------------------+
+ test/Tests/CSM/Equivariant.hs view
@@ -0,0 +1,91 @@++-- | Tests for the equivariant CSM class+++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.CSM.Equivariant where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.CSM.Equivariant.Ordered   as Ordered+import qualified Math.RootLoci.CSM.Equivariant.Recursive as Recur+import qualified Math.RootLoci.CSM.Equivariant.Direct    as Direct+import qualified Math.RootLoci.CSM.Equivariant.Umbral    as Umbral++import Math.RootLoci.Classic+import Math.RootLoci.CSM.Aluffi++import Tests.Common++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "equivariant (projective) CSM classes"+  [ testGroup "ordered"+      [ testCase "structure lemma"                           (forList [-3..4] "failed"   prop_structure_lemma                 )+      , testCase "rec. ord. distinct = U(n) formula /AB"     (forAllInt   7 "failed" $ prop_Un_equals_recursive ChernRoot   )+      , testCase "rec. ord. distinct = U(n) formula /Chern"  (forAllInt   7 "failed" $ prop_Un_equals_recursive ChernClass  )+      , testCase "recur. upper = sum of direct opens /AB"    (forAllSetp  6 "failed" $ prop_ord_rec_upper_vs_sum_direct ChernRoot   )+      , testCase "recur. upper = sum of direct opens /Chern" (forAllSetp  6 "failed" $ prop_ord_rec_upper_vs_sum_direct ChernClass  )+      , testCase "tangent class = sum of all strata"         (forAllInt   7 "failed" $ prop_tangent_class_equals_sum ChernClass)+      ]+  , testGroup "unordered"+      [ testGroup "direct"+          [ testCase "recursive CSM = direct CSM (open) /AB"      (forAllPart 8 "failed" $ prop_recur_direct_open   ChernRoot  ) +          , testCase "recursive CSM = direct CSM (open) /Chern"   (forAllPart 8 "failed" $ prop_recur_direct_open   ChernClass ) +          , testCase "recursive CSM = direct CSM (closed) /AB"    (forAllPart 8 "failed" $ prop_recur_direct_closed ChernRoot  ) +          , testCase "recursive CSM = direct CSM (closed) /Chern" (forAllPart 8 "failed" $ prop_recur_direct_closed ChernClass ) +          ]+      , testGroup "umbral"+          [ testCase "umbral CSM = direct CSM (open) /AB"        (forAllPart 10 "failed" $ prop_umbral_vs_direct_open   ChernRoot  ) +          , testCase "umbral CSM = direct CSM (open) /Chern"     (forAllPart 10 "failed" $ prop_umbral_vs_direct_open   ChernClass ) +          , testCase "umbral CSM = direct CSM (closed) /AB"      (forAllPart 10 "failed" $ prop_umbral_vs_direct_closed ChernRoot  ) +          , testCase "umbral CSM = direct CSM (closed) /Chern"   (forAllPart 10 "failed" $ prop_umbral_vs_direct_closed ChernClass ) +          ]+      ]+  ]++--------------------------------------------------------------------------------++prop_umbral_vs_direct_open sing part +  = ( spec2' sing (Direct.directOpenCSM part) == spec2' sing (Umbral.umbralOpenCSM part) )++prop_umbral_vs_direct_closed sing part +  = ( spec2' sing (Direct.directClosedCSM part) == spec2' sing (Umbral.umbralClosedCSM part) )++prop_recur_direct_open   sing part = (spec2' sing (Direct.directOpenCSM   part) == spec2' sing (Recur.openCSM   part))+prop_recur_direct_closed sing part = (spec2' sing (Direct.directClosedCSM part) == spec2' sing (Recur.closedCSM part))++--------------------------------------------------------------------------------++-- very slow for n>=5 !! (because we have +3, so 5 -> 8)+prop_structure_lemma n = (Ordered.computeQPolys n == Ordered.formulaQPoly n)++prop_Un_equals_recursive sing n = +  (spec2' sing $ Ordered.computeOpenStratumCSM  n) == +  (spec2' sing $ Ordered.formulaDistinctCSM     n)++prop_ord_rec_upper_vs_sum_direct sing setp =+  (spec2' sing $ unsafeEtaToOmega (Recur.upperClass setp)) ==+  (spec2' sing $ Ordered.computeClosureOfAnyStratumCSM     setp ) ++prop_tangent_class_equals_sum sing n = +  (spec2' sing $ Ordered.tangentChernClass n) == +  (spec2' sing $ ZMod.sum [ Ordered.formulaAnyStratumCSM setp | setp <- setPartitions n ]) +  +--------------------------------------------------------------------------------+
+ test/Tests/CSM/Projective.hs view
@@ -0,0 +1,54 @@++-- | Tests for the non-equivarant CSM classes++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.CSM.Projective where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.CSM.Equivariant.Ordered   as Ordered+import qualified Math.RootLoci.CSM.Equivariant.Direct    as Direct++import Math.RootLoci.Classic++import qualified Math.RootLoci.CSM.Aluffi     as Aluffi+import qualified Math.RootLoci.CSM.Projective as Proj++import Tests.Common++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "non-equivariant CSM classes"+  [ testGroup "unordered"+      [ testCase "aluffi == recursive (open)"       (forAllPart 10 "failed"  prop_aluffi_vs_recursive_open   )+      , testCase "aluffi == recursive (closed)"     (forAllPart 10 "failed"  prop_aluffi_vs_recursive_closed )+      , testCase "aluffi == from equiv  (open)"     (forAllPart 10 "failed"  prop_aluffi_vs_equiv_open       )+      , testCase "aluffi == from equiv. (closed)"   (forAllPart 10 "failed"  prop_aluffi_vs_equiv_closed     )+      , testCase "degree == lowest coeff of CSM"    (forAllPart 10 "failed"  prop_csm_degree                 )+      ]+  ]++prop_aluffi_vs_recursive_open   part = (Proj.openCSM   part == Aluffi.aluffiOpenCSM   part)+prop_aluffi_vs_recursive_closed part = (Proj.closedCSM part == Aluffi.aluffiClosedCSM part)++prop_aluffi_vs_equiv_open   part = (Aluffi.aluffiOpenCSM part   == forgetEquiv (spec2' ChernClass $ Direct.directOpenCSM   part))+prop_aluffi_vs_equiv_closed part = (Aluffi.aluffiClosedCSM part == forgetEquiv (spec2' ChernClass $ Direct.directClosedCSM part))++prop_csm_degree part = (hilbert part == Proj.lowestCoeff_ (Proj.closedCSM part))++--------------------------------------------------------------------------------
+ test/Tests/Common.hs view
@@ -0,0 +1,35 @@++-- | Shared utilities for testing++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.Common where++--------------------------------------------------------------------------------++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Partitions.Set++import Math.RootLoci.Algebra+import Math.RootLoci.Misc++import Test.Tasty.HUnit++--------------------------------------------------------------------------------++forList  :: [a] -> String -> (a -> Bool) -> Assertion+forList xs msg check = assertBool msg $ and [ check x | x <- xs ]++forAllInt  :: Int -> String -> (Int -> Bool) -> Assertion+forAllInt maxn msg check = assertBool msg $ and [ check i | i<-[0..maxn] ]++forAllPart :: Int -> String -> (Partition -> Bool) -> Assertion+forAllPart maxn msg check = assertBool msg $ and [ check p | p <- allPartitions maxn ]++forAllPartPos :: Int -> String -> (Partition -> Bool) -> Assertion+forAllPartPos maxn msg check = assertBool msg $ and [ check p | p <- allPartitions maxn , not (isEmpty p) ]++forAllSetp :: Int -> String -> (SetPartition -> Bool) -> Assertion+forAllSetp maxn msg check = assertBool msg $ and [ check p | k<-[0..maxn] , p <- setPartitions k ]++--------------------------------------------------------------------------------
+ test/Tests/Dual.hs view
@@ -0,0 +1,54 @@++-- | Tests involving the cohomological dual+++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.Dual where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Classes+import Math.Combinat.Partitions++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import qualified Math.RootLoci.Dual.Restriction  as MSc+import qualified Math.RootLoci.Dual.Localization as FNR++import qualified Math.RootLoci.CSM.Equivariant.Direct as Direct++import Math.RootLoci.Classic+import Math.RootLoci.CSM.Aluffi++import Tests.Common++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "dual classes"+  [ testCase "lemma 9.1.3"                         (forAllPartPos 15 "failed" prop_lemma913       ) +  , testCase "proj degree matches Hilbert"         (forAllPartPos 15 "failed" prop_degree_Hilbert )+  , testCase "proj degree matches Aluffi"          (forAllPartPos 14 "failed" prop_degree_Aluffi  )+  , testCase "dual class agrees with localization" (forAllPartPos 14 "failed" prop_msc_equals_local)+  +  , testCase "dual class = lowest deg part of open CSM"   (forAllPartPos 10 "failed" prop_msc_equals_lowest_csm_open  )+  , testCase "dual class = lowest deg part of closed CSM" (forAllPartPos 10 "failed" prop_msc_equals_lowest_csm_closed)+  ]++prop_lemma913 part = and [ MSc.lemma913 part k | k<-[0..m] ] where m = numberOfParts part+ +prop_degree_Hilbert part = (MSc.degreeMSc part == hilbert      part)+prop_degree_Aluffi  part = (MSc.degreeMSc part == aluffiDegree part)++prop_msc_equals_local part = FNR.localizeDual part == chernToAB (MSc.affineDualMSc part)++prop_msc_equals_lowest_csm_open   part = MSc.dualClassFromProjCSM (Direct.directOpenCSM   part) == MSc.affineDualMSc part+prop_msc_equals_lowest_csm_closed part = MSc.dualClassFromProjCSM (Direct.directClosedCSM part) == MSc.affineDualMSc part++--------------------------------------------------------------------------------
+ test/Tests/Pushforward.hs view
@@ -0,0 +1,49 @@++-- | Tests for the push-forward+++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.Pushforward where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Classes+import Math.Combinat.Partitions++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import Math.RootLoci.CSM.Equivariant.PushForward ++import Math.RootLoci.Classic++import Tests.Common++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "pushforward"+  [ testCase "tau definition"                                    (forList [-1..20] "failed" prop_tau_defin                 ) +  , testCase "symm breaking pi_* == recursive formula for P_j"   (forAllInt 20 "failed" prop_symmbreaking_vs_ppolys        ) +  , testCase "affine pi_* == proj pi_* [ gamma -> 0 ] /AB"       (forAllInt 20 "failed" (prop_ppoly_aff_vs_proj ChernRoot ))+  , testCase "affine pi_* == proj pi_* [ gamma -> 0 ] /Chern"    (forAllInt 20 "failed" (prop_ppoly_aff_vs_proj ChernClass))+  ]++prop_symmbreaking_vs_ppolys n = spec3' ChernRoot (piStarTableProj n) == pi_star_table n++prop_ppoly_aff_vs_proj sing n = spec2' sing (piStarTableAff n) == fmap forgetGamma (spec3' sing (piStarTableProj n))++prop_tau_defin n = (tau n * (a - b)) == (apow - bpow) where+  a    = ZMod.generator $ AB  1     0   +  b    = ZMod.generator $ AB  0     1   +  apow = ZMod.generator $ AB (n+1)  0   +  bpow = ZMod.generator $ AB  0    (n+1)++--------------------------------------------------------------------------------
+ test/Tests/RootVsClass/Check.hs view
@@ -0,0 +1,93 @@++-- | Checking polymorphic functions++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies, ScopedTypeVariables #-}+module Tests.RootVsClass.Check where++--------------------------------------------------------------------------------++-- import Data.Proxy+-- import Math.Combinat.Partitions++import Data.Array++import Math.RootLoci.Algebra+import Math.RootLoci.Misc+import Math.RootLoci.Geometry.Cohomology ++import qualified Math.RootLoci.Algebra.FreeMod as ZMod++import Math.RootLoci.CSM.Equivariant.Umbral ( ST )++--------------------------------------------------------------------------------++checkZMod :: (forall b. ChernBase b => ZMod b) -> Bool+checkZMod polymorph+  =  ( abToChern (spec1' ChernRoot polymorph) ==            spec1' ChernClass polymorph )+  && (            spec1' ChernRoot polymorph  == chernToAB (spec1' ChernClass polymorph) ) ++{-+checkZModExt :: forall f. Equivariant f =>  (forall b. ChernBase b => ZMod (f b)) -> Bool+checkZModExt polymorph+  =  ( convertEquiv abToChern (spec2' ChernRoot polymorph) ==                         spec2' ChernClass polymorph  )+  && (                         spec2' ChernRoot polymorph  == convertEquiv chernToAB (spec2' ChernClass polymorph) ) +-}++--------------------------------------------------------------------------------++{-+checkGam   :: (forall b. ChernBase b => ZMod (Gam   b)) -> Bool+checkOmega :: (forall b. ChernBase b => ZMod (Omega b)) -> Bool+checkEta   :: (forall b. ChernBase b => ZMod (Eta   b)) -> Bool++checkGam   = checkZModExt+checkOmega = checkZModExt+checkEta   = checkZModExt+-}++checkOmega :: (forall b. ChernBase b => ZMod (Omega b)) -> Bool+checkOmega polymorph+  =  ( convertOmega abToChern (spec2' ChernRoot polymorph) ==                         spec2' ChernClass polymorph  )+  && (                         spec2' ChernRoot polymorph  == convertOmega chernToAB (spec2' ChernClass polymorph) ) ++checkEta :: (forall b. ChernBase b => ZMod (Eta b)) -> Bool+checkEta polymorph+  =  ( convertEta abToChern (spec2' ChernRoot polymorph) ==                       spec2' ChernClass polymorph  )+  && (                       spec2' ChernRoot polymorph  == convertEta chernToAB (spec2' ChernClass polymorph) ) ++checkGam :: (forall b. ChernBase b => ZMod (Gam b)) -> Bool+checkGam polymorph+  =  ( convertGam abToChern (spec2' ChernRoot polymorph) ==                       spec2' ChernClass polymorph  )+  && (                       spec2' ChernRoot polymorph  == convertGam chernToAB (spec2' ChernClass polymorph) ) ++--------------------------------------------------------------------------------++checkArrZMod :: (forall b. ChernBase b => Array Int (ZMod b)) -> Bool+checkArrZMod polymorph+  =  ( fmap abToChern (spec2' ChernRoot polymorph) ==                 spec2' ChernClass polymorph )+  && (                 spec2' ChernRoot polymorph  == fmap chernToAB (spec2' ChernClass polymorph) ) ++checkArrGam :: (forall b. ChernBase b => Array Int (ZMod (Gam b))) -> Bool+checkArrGam polymorph+  =  ( fmap fwd (spec3' ChernRoot polymorph) ==           spec3' ChernClass polymorph )+  && (           spec3' ChernRoot polymorph  == fmap bwd (spec3' ChernClass polymorph) ) +  where+    fwd = convertGam abToChern+    bwd = convertGam chernToAB+    +--------------------------------------------------------------------------------++{-+checkMixedST :: forall c. (Eq c, Num c) => (forall b. ChernBase b => FreeMod (FreeMod c b) ST) -> Bool+checkMixedST polymorph+  =  ( fwd (spec2' ChernRoot polymorph) ==       spec2' ChernClass polymorph  )+  && (      spec2' ChernRoot polymorph  ==  bwd (spec2' ChernClass polymorph) ) +  where+    fwd :: FreeMod (FreeMod c AB   ) ST -> FreeMod (FreeMod c Chern) ST  +    bwd :: FreeMod (FreeMod c Chern) ST -> FreeMod (FreeMod c AB   ) ST +    fwd = ZMod.mapCoeff abToChern+    bwd = ZMod.mapCoeff chernToAB+-}++--------------------------------------------------------------------------------+
+ test/Tests/RootVsClass/Direct.hs view
@@ -0,0 +1,36 @@++-- | Comparing the Chern root vs. the Chern class versions for stuff in+-- "Math.RootLoci.CSM.Equivariant.Direct"++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.RootVsClass.Direct where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Partitions++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import Math.RootLoci.CSM.Equivariant.Direct ++import Tests.Common+import Tests.RootVsClass.Check++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "direct"+  [ testCase "open CSM"      (forAllPart 9 "failed" prop_directOpenCSM  ) +  , testCase "closed CSM"    (forAllPart 9 "failed" prop_directClosedCSM) +  ]+  +prop_directOpenCSM   part = checkGam (directOpenCSM   part)+prop_directClosedCSM part = checkGam (directClosedCSM part)++--------------------------------------------------------------------------------
+ test/Tests/RootVsClass/Ordered.hs view
@@ -0,0 +1,46 @@++-- | Comparing the Chern root vs. the Chern class versions for stuff in+-- "Math.RootLoci.CSM.Equivariant.Ordered"++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.RootVsClass.Ordered where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Partitions++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import Math.RootLoci.CSM.Equivariant.Ordered ++import Tests.Common+import Tests.RootVsClass.Check++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "ordered"+  [ testCase "tangent Chern class" (forAllInt   9 "failed" prop_tangentChernClass    )+  , testCase "small diagonals"     (forAllInt   9 "failed" prop_smallDiagonal        )+  , testCase "open stratum"        (forAllInt   7 "failed" prop_openStratumCSM       )+  , testCase "any stratum"         (forAllSetp  6 "failed" prop_anyStratumCSM        )+  , testCase "formula for Q-poly"  (forList [-3.. 20] "failed" prop_formulaQPoly         )+  , testCase "formula U(n)"        (forAllInt  10 "failed" prop_formulaDistinctCSM   )+  ]+  +--------------------------------------------------------------------------------++prop_tangentChernClass  n = checkOmega (tangentChernClass  n)+prop_smallDiagonal      n = checkOmega (smallDiagonal      n)+prop_openStratumCSM     n = checkOmega (computeOpenStratumCSM    n)+prop_anyStratumCSM   setp = checkOmega (computeAnyStratumCSM  setp)+prop_formulaDistinctCSM n = checkOmega (formulaDistinctCSM n)+prop_formulaQPoly       n = checkZMod  (formulaQPoly       n)++--------------------------------------------------------------------------------
+ test/Tests/RootVsClass/PushForward.hs view
@@ -0,0 +1,41 @@++-- | Comparing the Chern root vs. the Chern class versions for stuff in+-- "Math.RootLoci.CSM.Equivariant.PushForward"++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.RootVsClass.PushForward where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Partitions++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import Math.RootLoci.CSM.Equivariant.PushForward++import Tests.Common+import Tests.RootVsClass.Check++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "pushforward"+  [ testCase "tau"             (forAllInt 20 "failed" prop_tau    ) +  , testCase "tauEta"          (forAllInt 20 "failed" prop_tauEta ) +  , testCase "piStarTableAff"  (forAllInt 15 "failed" prop_piStarTableAff  ) +  , testCase "piStarTableProj" (forAllInt 15 "failed" prop_piStarTableProj ) +  ]+  +prop_tau    n = checkZMod (tau    n)+prop_tauEta n = checkEta  (tauEta n)++prop_piStarTableAff  n = checkArrZMod (piStarTableAff  n)+prop_piStarTableProj n = checkArrGam  (piStarTableProj n)++--------------------------------------------------------------------------------
+ test/Tests/RootVsClass/Recursive.hs view
@@ -0,0 +1,44 @@++-- | Comparing the Chern root vs. the Chern class versions for stuff in+-- "Math.RootLoci.CSM.Equivariant.Recursive"++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.RootVsClass.Recursive where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Partitions++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import Math.RootLoci.CSM.Equivariant.Recursive ++import Tests.Common+import Tests.RootVsClass.Check++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "recursive"+  [ testCase "upper classes /setp" (forAllSetp 6 "failed" prop_upper_setp    )+  , testCase "upper classes /part" (forAllPart 7 "failed" prop_upper_part    )+  , testCase "lower classes" (forAllPart 7 "failed" prop_lower    )+  , testCase "open CSM"      (forAllPart 8 "failed" prop_openCSM  ) +  , testCase "closed CSM"    (forAllPart 8 "failed" prop_closedCSM) +  ]+  +prop_openCSM   part = checkGam (openCSM   part)+prop_closedCSM part = checkGam (closedCSM part)++prop_upper_setp setp = checkEta (upperClass $ setp)+prop_upper_part part = checkEta (upperClass $ defaultSetPartition part)++prop_lower part = checkGam (lowerClass part)++--------------------------------------------------------------------------------
+ test/Tests/RootVsClass/Umbral.hs view
@@ -0,0 +1,42 @@++-- | Comparing the Chern root vs. the Chern class versions for stuff in+-- "Math.RootLoci.CSM.Equivariant.Umbral"++{-# LANGUAGE Rank2Types, GADTs, TypeFamilies #-}+module Tests.RootVsClass.Umbral where++--------------------------------------------------------------------------------++import Data.Proxy++import Math.Combinat.Partitions++import Math.RootLoci.Algebra+import Math.RootLoci.Geometry+import Math.RootLoci.Misc++import Math.RootLoci.CSM.Equivariant.Umbral++import Tests.Common+import Tests.RootVsClass.Check++import Test.Tasty+import Test.Tasty.HUnit++--------------------------------------------------------------------------------++all_tests = testGroup "umbral"+  [ +    testCase "open affine CSM"      (forAllPart 11 "failed" prop_umbralAffOpenCSM  ) +  , testCase "closed affine CSM"    (forAllPart 11 "failed" prop_umbralAffClosedCSM)     +  -- , testCase "theta"     (forAllPosInt 15 "failed" prop_theta ) +  -- , testCase "thetaQ"    (forAllPosInt 15 "failed" prop_thetaQ) +  ]+  +-- prop_theta  n = checkMixedST (theta  n)+-- prop_thetaQ n = checkMixedST (thetaQ n)++prop_umbralAffOpenCSM   part = checkZMod (umbralAffOpenCSM   part)+prop_umbralAffClosedCSM part = checkZMod (umbralAffClosedCSM part)++--------------------------------------------------------------------------------
+ test/testSuite.hs view
@@ -0,0 +1,49 @@++-- | The test-suite++module Main where++--------------------------------------------------------------------------------++import Test.Tasty++-- import Test.Tasty.HUnit+-- import Test.Tasty.SmallCheck as SC+-- import Test.Tasty.QuickCheck as QC++import qualified Tests.RootVsClass.Ordered +import qualified Tests.RootVsClass.Recursive+import qualified Tests.RootVsClass.Direct+import qualified Tests.RootVsClass.PushForward+import qualified Tests.RootVsClass.Umbral++import qualified Tests.Dual+import qualified Tests.Pushforward++import qualified Tests.CSM.Equivariant+import qualified Tests.CSM.Projective++--------------------------------------------------------------------------------++main = defaultMain tests++tests :: TestTree+tests = testGroup "Tests"  +  [ Tests.Pushforward.all_tests+  , Tests.CSM.Projective.all_tests+  , Tests.CSM.Equivariant.all_tests+  , rootVsClass+  , Tests.Dual.all_tests+  ]++rootVsClass :: TestTree  +rootVsClass = testGroup "chern root vs. chern class tests"+  [ Tests.RootVsClass.Ordered.all_tests+  , Tests.RootVsClass.Recursive.all_tests+  , Tests.RootVsClass.Direct.all_tests+  , Tests.RootVsClass.PushForward.all_tests+  , Tests.RootVsClass.Umbral.all_tests+  ]++--------------------------------------------------------------------------------+