-- | 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.Algebra.Polynomial.FreeModule 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 :: Partition -> ZMod Chern
affineDualMSc part :: Partition
part@(Partition [Int]
ps) = 

  case [Int]
ps of
    []            -> [Char] -> ZMod Chern
forall a. HasCallStack => [Char] -> a
error [Char]
"affine_tp_msc: empty partition"
    [Int
n]           -> [(Chern, Integer)] -> ZMod Chern
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList [ ( Int -> Int -> Chern
Chern (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
j) Int
j , Ratio Integer -> Integer
forall p. (Eq p, Num p) => Ratio p -> p
rat2int (Ratio Integer -> Integer) -> Ratio Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Ratio Integer
single Int
j ) | Int
j<-[ Int
0 .. Int -> Int -> Int
forall a. Integral a => a -> a -> a
div (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
d) Int
2] ] 
    [Int
a,Int
b] | Int
aInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
b  -> [(Chern, Integer)] -> ZMod Chern
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList [ ( Int -> Int -> Chern
Chern (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
j) Int
j , Ratio Integer -> Integer
forall p. (Eq p, Num p) => Ratio p -> p
rat2int (Ratio Integer -> Integer) -> Ratio Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Ratio Integer
double Int
j ) | Int
j<-[ Int
0 .. Int -> Int -> Int
forall a. Integral a => a -> a -> a
div (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
d) Int
2] ] 
    [Int]
otherwise     -> [(Chern, Integer)] -> ZMod Chern
forall c b. (Eq c, Num c, Ord b) => [(b, c)] -> FreeMod c b
ZMod.fromList [ ( Int -> Int -> Chern
Chern (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
j) Int
j , Ratio Integer -> Integer
forall p. (Eq p, Num p) => Ratio p -> p
rat2int (Ratio Integer -> Integer) -> Ratio Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Ratio Integer
lambda Int
j ) | Int
j<-[ Int
0 .. Int -> Int -> Int
forall a. Integral a => a -> a -> a
div (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
d) Int
2] ] 

  where

    n :: Int
n = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
ps
    d :: Int
d = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
ps

    p :: Int
p = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div  Int
n    Int
2
    q :: Int
q = Int -> Int -> Int
forall a. Integral a => a -> a -> a
div (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
2

    rat2int :: Ratio p -> p
rat2int Ratio p
r = case Ratio p -> p
forall a. Ratio a -> a
denominator Ratio p
r of
      p
1 -> Ratio p -> p
forall a. Ratio a -> a
numerator Ratio p
r
      p
_ -> [Char] -> p
forall a. HasCallStack => [Char] -> a
error [Char]
"lambda_j: not integer"

    lambda :: Int -> Ratio Integer
lambda Int
j = (Int -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi Int
n Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Ratio Integer
2)Ratio Integer -> Int -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
q) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Integer -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi (Int -> Integer
forall a. Integral a => a -> Integer
doubleFactorial (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2))Ratio Integer -> Integer -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Ratio Integer
s where
      s :: Ratio Integer
s = [Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
        [ Int -> Ratio Integer -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
lpsi) (Ratio Integer -> Ratio Integer) -> Ratio Integer -> Ratio Integer
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Ratio Integer
bigTheta Int
j Int
nphi Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* (Int -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nphiInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Int -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi Int
n)Ratio Integer -> Int -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ (Integer -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi (Integer -> Ratio Integer) -> Integer -> Ratio Integer
forall a b. (a -> b) -> a -> b
$ Partition -> Integer
aut Partition
phi Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Partition -> Integer
aut Partition
psi)
        | (Partition
phi,Partition
psi) <- Set (Partition, Partition) -> [(Partition, Partition)]
forall a. Set a -> [a]
Set.toList (Partition -> Set (Partition, Partition)
divideIntoTwoNonEmpty Partition
part)
        , let nphi :: Int
nphi = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Partition -> [Int]
fromPartition Partition
phi
        , let npsi :: Int
npsi = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Partition -> [Int]
fromPartition Partition
psi
        , let lphi :: Int
lphi = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Partition -> [Int]
fromPartition Partition
phi
        , let lpsi :: Int
lpsi = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Partition -> [Int]
fromPartition Partition
psi
        ] 

    gamma :: Int -> Rational
    gamma :: Int -> Ratio Integer
gamma Int
k 
      | Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n   = Ratio Integer
0 
      | Bool
otherwise  = Int -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n)) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Int -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi ((Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n)Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n))

    bigTheta :: Int -> Int -> Rational
    bigTheta :: Int -> Int -> Ratio Integer
bigTheta Int
j Int
k 
      | Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n   = Ratio Integer
0 
      | Bool
otherwise  = Int -> Ratio Integer
gamma Int
k Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Int -> Int -> Ratio Integer
smallTheta Int
j Int
k

    smallTheta :: Int -> Int -> Rational
    smallTheta :: Int -> Int -> Ratio Integer
smallTheta Int
j Int
k = Int -> [Ratio Integer] -> Ratio Integer
sympoly (Int
qInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j) [ Int -> Ratio Integer
gamma Int
i | Int
i<-[Int
1..Int
q] , Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/=Int
k, Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/=Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k ]
   
    fi :: Integral a => a -> Rational
    fi :: a -> Ratio Integer
fi = a -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral

    sqj :: Int -> Rational
    sqj :: Int -> Ratio Integer
sqj Int
j = Int -> [Ratio Integer] -> Ratio Integer
sympoly (Int
qInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
j) [ Int -> Ratio Integer
gamma Int
i | Int
i<-[Int
1..Int
q] ]

    sympoly :: Int -> [Rational] -> Rational
    sympoly :: Int -> [Ratio Integer] -> Ratio Integer
sympoly Int
k [Ratio Integer]
xs = [Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ [Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Ratio Integer]
ys | [Ratio Integer]
ys <- Int -> [Ratio Integer] -> [[Ratio Integer]]
forall a. Int -> [a] -> [[a]]
choose Int
k [Ratio Integer]
xs ]

    -- S(n)
    single :: Int -> Ratio Integer
single Int
j = Integer -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi (Int -> Integer
forall a. Integral a => a -> Integer
factorial Int
n) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ ([Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [ Int -> Ratio Integer
gamma Int
i | Int
i<-[Int
1..Int
q] ])
             Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Int -> Ratio Integer -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd Int
j (Int -> Ratio Integer
sqj Int
j) 

    -- S(p,p)
    double :: Int -> Ratio Integer
double Int
j = Integer -> Ratio Integer
forall a. Integral a => a -> Ratio Integer
fi (Int -> Integer
forall a. Integral a => a -> Integer
doubleFactorial Int
n)Ratio Integer -> Integer -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2 Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Ratio Integer
4 
             Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Int -> Ratio Integer -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd (Int
qInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j) (Int -> Ratio Integer
sqj Int
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 :: Int -> ZMod Chern -> Integer
projDegreeFromDual Int
n ZMod Chern
zm = Ratio Integer -> Integer
fromRat Ratio Integer
s where 

  s :: Rational
  s :: Ratio Integer
s = [Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Integer -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
c Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Ratio Integer
c1Ratio Integer -> Int -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
e Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* Ratio Integer
c2Ratio Integer -> Int -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
f  | (Chern Int
e Int
f, Integer
c) <- ZMod Chern -> [(Chern, Integer)]
forall c b. FreeMod c b -> [(b, c)]
ZMod.toList ZMod Chern
zm ]

  c1 :: Ratio Integer
c1 = Ratio Integer
2 Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Int -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral  Int
n    :: Rational
  c2 :: Ratio Integer
c2 = Ratio Integer
1 Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Int -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) :: Rational

-- | Compute the degree of the strata via the formula for the dual class
degreeMSc :: Partition -> Integer
degreeMSc :: Partition -> Integer
degreeMSc Partition
part = Int -> ZMod Chern -> Integer
projDegreeFromDual (Partition -> Int
partitionWeight Partition
part) (Partition -> ZMod Chern
affineDualMSc Partition
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 :: ZMod (Gam base) -> ZMod base
dualClassFromProjCSM ZMod (Gam base)
csm = ZMod base -> ZMod base
forall base. ChernBase base => ZMod base -> ZMod base
dualClassFromAffCSM ((Gam base -> Maybe base) -> ZMod (Gam base) -> ZMod base
forall a b c.
(Ord a, Ord b, Eq c, Num c) =>
(a -> Maybe b) -> FreeMod c a -> FreeMod c b
ZMod.mapMaybeBase Gam base -> Maybe base
nogamma ZMod (Gam base)
csm) where
  nogamma :: Gam base -> Maybe base
  nogamma :: Gam base -> Maybe base
nogamma (Gam Int
k base
ab) = if Int
kInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 then base -> Maybe base
forall a. a -> Maybe a
Just base
ab else Maybe base
forall a. Maybe a
Nothing

dualClassFromAffCSM :: ChernBase base => ZMod base -> ZMod base
dualClassFromAffCSM :: ZMod base -> ZMod base
dualClassFromAffCSM ZMod base
csm = Int -> ZMod base -> ZMod base
forall b. (Ord b, Graded b) => Int -> ZMod b -> ZMod b
filterGrade Int
min_degree ZMod base
csm where
  min_degree :: Int
min_degree = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ (base -> Int) -> [base] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map base -> Int
forall a. Graded a => a -> Int
grade ([base] -> [Int]) -> [base] -> [Int]
forall a b. (a -> b) -> a -> b
$ ((base, Integer) -> base) -> [(base, Integer)] -> [base]
forall a b. (a -> b) -> [a] -> [b]
map (base, Integer) -> base
forall a b. (a, b) -> a
fst ([(base, Integer)] -> [base]) -> [(base, Integer)] -> [base]
forall a b. (a -> b) -> a -> b
$ ZMod base -> [(base, Integer)]
forall c b. FreeMod c b -> [(b, c)]
ZMod.toList ZMod base
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 :: Partition -> Int -> Bool
lemma913 Partition
part Int
h = (Ratio Integer
aRatio Integer -> Ratio Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Ratio Integer
b) where 
  (Ratio Integer
a,Ratio Integer
b) = Partition -> Int -> (Ratio Integer, Ratio Integer)
lemma913' Partition
part Int
h 

  lemma913' :: Partition -> Int -> (Rational, Rational)
  lemma913' :: Partition -> Int -> (Ratio Integer, Ratio Integer)
lemma913' part :: Partition
part@(Partition [Int]
ps) Int
h = ( Ratio Integer
lhs , Ratio Integer
rhs ) where

    n :: Int
n = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
ps
    d :: Int
d = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
ps

    rhs :: Ratio Integer
rhs | Int
h Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d  = Integer -> Ratio Integer
tr (Int -> Integer
forall a. Integral a => a -> Integer
factorial Int
d) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* [Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Int -> Ratio Integer) -> [Int] -> [Ratio Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Ratio Integer
fi [Int]
ps)
        | Int
h Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
d  = Ratio Integer
0
        | Int
h Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>  Int
d  = -Ratio Integer
666

    lhs :: Ratio Integer
lhs = [Ratio Integer] -> Ratio Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
      [ Int -> Ratio Integer -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b -> b
negateIfOdd ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
rs) (Ratio Integer -> Ratio Integer) -> Ratio Integer -> Ratio Integer
forall a b. (a -> b) -> a -> b
$  (Int -> Ratio Integer
fi (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
qs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ Ratio Integer
2)Ratio Integer -> Int -> Ratio Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
h Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Num a => a -> a -> a
* (Integer -> Ratio Integer
tr (Integer -> Ratio Integer) -> Integer -> Ratio Integer
forall a b. (a -> b) -> a -> b
$ Partition -> Integer
aut Partition
part) Ratio Integer -> Ratio Integer -> Ratio Integer
forall a. Fractional a => a -> a -> a
/ (Integer -> Ratio Integer
tr (Integer -> Ratio Integer) -> Integer -> Ratio Integer
forall a b. (a -> b) -> a -> b
$ Partition -> Integer
aut Partition
phi Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Partition -> Integer
aut Partition
psi)
      | ( phi :: Partition
phi@(Partition [Int]
qs) , psi :: Partition
psi@(Partition [Int]
rs) ) <- Set (Partition, Partition) -> [(Partition, Partition)]
forall a. Set a -> [a]
Set.toList (Partition -> Set (Partition, Partition)
divideIntoTwo Partition
part)
      ]

    fi :: Int -> Rational
    fi :: Int -> Ratio Integer
fi = Int -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral

    tr :: Integer -> Rational
    tr :: Integer -> Ratio Integer
tr = Integer -> Ratio Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral  


--------------------------------------------------------------------------------
-- * helper functions

-- | Different ways to divide a partition into two 
divideIntoTwo :: Partition -> Set (Partition,Partition)
divideIntoTwo :: Partition -> Set (Partition, Partition)
divideIntoTwo (Partition [Int]
ps) = [(Partition, Partition)] -> Set (Partition, Partition)
forall a. Ord a => [a] -> Set a
Set.fromList ([(Partition, Partition)] -> Set (Partition, Partition))
-> [(Partition, Partition)] -> Set (Partition, Partition)
forall a b. (a -> b) -> a -> b
$ ([Bool] -> (Partition, Partition))
-> [[Bool]] -> [(Partition, Partition)]
forall a b. (a -> b) -> [a] -> [b]
map [Bool] -> (Partition, Partition)
f (Int -> [[Bool]]
binaryTuples Int
d) where

  d :: Int
d    = [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
ps
  f :: [Bool] -> (Partition, Partition)
f [Bool]
ts = ( [Bool] -> Partition
g [Bool]
ts , [Bool] -> Partition
g ((Bool -> Bool) -> [Bool] -> [Bool]
forall a b. (a -> b) -> [a] -> [b]
map Bool -> Bool
not [Bool]
ts) )
  g :: [Bool] -> Partition
g [Bool]
ts = [Int] -> Partition
Partition [ Int
k | (Bool
b,Int
k) <- [Bool] -> [Int] -> [(Bool, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Bool]
ts [Int]
ps , Bool
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 :: Partition -> Set (Partition, Partition)
divideIntoTwoNonEmpty Partition
p = (Partition, Partition)
-> Set (Partition, Partition) -> Set (Partition, Partition)
forall a. Ord a => a -> Set a -> Set a
Set.delete (Partition, Partition)
x (Set (Partition, Partition) -> Set (Partition, Partition))
-> Set (Partition, Partition) -> Set (Partition, Partition)
forall a b. (a -> b) -> a -> b
$ (Partition, Partition)
-> Set (Partition, Partition) -> Set (Partition, Partition)
forall a. Ord a => a -> Set a -> Set a
Set.delete (Partition, Partition)
y (Set (Partition, Partition) -> Set (Partition, Partition))
-> Set (Partition, Partition) -> Set (Partition, Partition)
forall a b. (a -> b) -> a -> b
$ Partition -> Set (Partition, Partition)
divideIntoTwo Partition
p where
  x :: (Partition, Partition)
x = (Partition
emptyPartition,Partition
p)
  y :: (Partition, Partition)
y = (Partition
p,Partition
emptyPartition)

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