-- | 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 ]

--------------------------------------------------------------------------------
-- * Enumerative geometry

-- | The degree of the dual curve is @d(d-1)@
degreeOfDualCurve :: Int -> Integer
degreeOfDualCurve d0 
  | d < 2     = 0
  | otherwise = d*(d-1) 
  where
    d = fromIntegral d0 :: Integer

-- | Number of flex lines to a generic degree @d@ plane curve
numberOfCurveFlexes :: Int -> Integer
numberOfCurveFlexes d0
  | d < 3     = 0
  | otherwise = 3*d*(d-2)
  where
    d = fromIntegral d0 :: Integer

-- | Number of bitangent lines to a generic degree @d@ plane curve
numberOfCurveBiTangents :: Int -> Integer
numberOfCurveBiTangents d0
  | d < 3     = 0
  | otherwise = div ((-3 + d)* (-2 + d)* d* (3 + d)) 2 
  where
    d = fromIntegral d0 :: Integer

-- | Number of 4-tangent lines to a generic degree @d@ surface (Schubert)
numberOfSurface4xTangents :: Int -> Integer
numberOfSurface4xTangents 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
numberOfSurface5xHyperflexes :: Int -> Integer
numberOfSurface5xHyperflexes d0
  | d < 5     = 0
  | otherwise = (35*d^3 - 200*d^2 + 240*d)
  where
    d = fromIntegral d0 :: Integer

-- | Bidegree of bitangent locus of a generic hypersurface
-- 
-- (See: Kathlen Kohn, Bernt Ivar Utstol Nodland, Paolo Tripoli: Secants, bitangents, and their congruences)
--
bidegreeOfSurfaceBiTangents :: Int -> (Integer,Integer)
bidegreeOfSurfaceBiTangents d0 
  | d < 4     = ( 0 , 0 )
  | otherwise = ( div (d*(d-1)*(d-2)*(d-3)) 2 , div (d*(d-2)*(d-3)*(d+3)) 2 )
  where
    d = fromIntegral d0 :: Integer

-- | Bidegree of the flex locus of a generic hypersurface
--
-- (See: Kathlen Kohn, Bernt Ivar Utstol Nodland, Paolo Tripoli: Secants, bitangents, and their congruences)
--
bidegreeOfSurfaceFlexes :: Int -> (Integer,Integer)
bidegreeOfSurfaceFlexes d0
  | d < 4     = ( 0 , 0 ) 
  | otherwise = ( d*(d-1)*(d-3) , 3*d*(d-2) )
  where
    d = fromIntegral d0 :: Integer

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