packages feed

sigma-ij (empty) → 0.2

raw patch · 22 files changed

+2510/−0 lines, 22 filesdep +arraydep +basedep +combinatsetup-changed

Dependencies added: array, base, combinat, containers, optparse-applicative, parsec2, random, time

Files

+ LICENSE view
@@ -0,0 +1,29 @@+Copyright (c) 2010, 2016, 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.+
+ README.txt view
@@ -0,0 +1,34 @@++This is a program to compute Thom polynomials of second-order +Thom-Boardman singularities $Sigma^{i,j}(n)$.++The computation is based on the localization method described in +the author's PhD thesis: <http://renyi.hu/~komuves/phdthesis.pdf>.+++USAGE:+======++sigma-ij -h                        help+sigma-ij -i3 -j1 -n7               compute $Tp(Sigma^{3,1}(7))$+sigma-oj -i3 -j1 -n7 -r<RING>      compute with coefficients in the given ring+sigma-oj -i3 -j1 -n7 -B<N> -b<n>   compute the n-th (out of N) part+sigma-oj -i3 -j1 -n7 -rZp          compute in the (baked-in) prime field Zp+sigma-oj -i3 -j1 -n7 -o<FILE>      change the output file++Supported rings:+ * rationals + * integers (remark: the division-free determinant algorithm often fails)+ * Zp, a baked-in prime field ++The -B and -b options are useful to parallelize the computation over +many computers.+ ++TODO:+=====++ - better (and faster) prime field implementation(s)+ - allow arbitrary prime fields instead of just a baked-in one+ - pivoting for the Bareiss (division-free) determinant algorithm+ - implement explicit formula for j=1
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ sigma-ij.cabal view
@@ -0,0 +1,78 @@++Name:                sigma-ij+Version:             0.2+Synopsis:            Thom polynomials of second order Thom-Boardman singularities+Description:         A program to compute Thom polynomials of second order Thom-Boardman +                     singularities, using the localization method described in the+                     author's PhD thesis <http://renyi.hu/~komuves/phdthesis.pdf>.+License:             BSD3+License-file:        LICENSE+Author:              Balazs Komuves+Copyright:           (c) 2010, 2016 Balazs Komuves+Maintainer:          bkomuves (plus) hackage (at) gmail (dot) com+Homepage:            http://code.haskell.org/~bkomuves/+Stability:           Experimental+Category:            Math+Tested-With:         GHC == 7.10.3+Cabal-Version:       >= 1.18+Build-Type:          Simple++--------------------------------------------------------------------------------++extra-source-files:  src/cbits/c_det.c+                     src/cbits/c_det.h+                     README.txt++--------------------------------------------------------------------------------++Executable sigma-ij++  hs-source-dirs:      src+  main-is:             sigmaij.hs++  Build-Depends:       base >= 4 && < 5, array >= 0.5, containers >= 0.5, random,+                       time, parsec2, optparse-applicative, +                       combinat >= 0.2.8++  -- cabal gets confused if the executable is in the same source tree...+  c-sources:           src/cbits/c_det.c    +  cc-options:          -std=c99 ++  Default-Language:    Haskell2010++--------------------------------------------------------------------------------++Library  ++  hs-source-dirs:      src+  +  c-sources:           src/cbits/c_det.c+  cc-options:          -std=c99 ++  exposed-modules:     Math.ThomPoly.SigmaI+                       Math.ThomPoly.SigmaIJ+                       Math.ThomPoly.Formulae+                       Math.ThomPoly.Shared+                       Math.ThomPoly.Subs+                       Math.Algebra.Schur+                       Math.Algebra.Determinant+                       Math.Algebra.ModP+                       Math.FreeModule.Class+                       Math.FreeModule.Helper+                       Math.FreeModule.Parser+                       Math.FreeModule.PP+                       Math.FreeModule.PrettyPrint+                       Math.FreeModule.Symbol+                       Math.FreeModule.SortedList++  Build-Depends:       base >= 4 && < 5, array >= 0.5, containers >= 0.5, random,+                       time, parsec2, optparse-applicative, +                       combinat >= 0.2.8++  Default-Extensions:  CPP, BangPatterns, ScopedTypeVariables+  Other-Extensions:    TypeFamilies, ForeignFunctionInterface++  Default-Language:    Haskell2010++  ghc-options:         -fwarn-tabs -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-unused-imports+
+ src/Math/Algebra/Determinant.hs view
@@ -0,0 +1,291 @@++-- | Determinants.+--+-- TODO: specialized prime fields; fast C implementation; pivoting for Bareiss+--++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, +             FlexibleInstances, TypeSynonymInstances,+             ForeignFunctionInterface+  #-}+module Math.Algebra.Determinant where++--------------------------------------------------------------------------------++import Control.Monad+import Control.Monad.ST++import Data.Array.Base+import Data.Array.IArray+import Data.Array.MArray+import Data.Array.Unsafe+import Data.Array.ST++import Data.List+import Data.Ratio+import Data.STRef++import Data.Bits+import Data.Word+import Data.Int++import Foreign.C+import Foreign.Ptr+import Foreign.Marshal+import System.IO.Unsafe as Unsafe++import System.Random++import Debug.Trace+import GHC.IO ( unsafeIOToST )++import Math.Algebra.ModP++--------------------------------------------------------------------------------+-- * matrices++type Matrix a = Array (Int,Int) a++printMatrix :: Show a => Matrix a -> IO ()+printMatrix = putStrLn . showMatrix++showMatrix :: Show a => Matrix a -> String+showMatrix = unlines . showMatrix'++showMatrix' :: Show a => Matrix a -> [String]+showMatrix' mat = map mkRow (transpose cols) where+  ((1,1),(n,m)) = bounds mat+  cols = map extend [ [ show (mat!(i,j)) | i<-[1..n] ] | j<-[1..m] ]++  mkRow strs = "[ " ++ intercalate " " strs ++ " ]"++  extend :: [String] -> [String]+  extend xs = map f xs where+    n = maximum (map length xs)+    f s = replicate (n - length s) ' ' ++ s++--------------------------------------------------------------------------------+-- * a type class for determinants++class (Eq a, Num a, Show a) => Determinant a where +  determinant :: Matrix a -> a++instance Determinant Integer  where determinant = bareissDeterminantFullRank+instance Determinant Int      where determinant = bareissDeterminantFullRank+instance Determinant Rational where determinant = gaussElimDeterminant+instance Determinant Zp       where determinant = gaussElimDeterminantInt64++--------------------------------------------------------------------------------+-- * C implementation of determinant in a prime field (gaussian elimination, fitting into 64 bit)++foreign import ccall "c_det.h inv_modp" c_inv_modp :: Int64 -> Int64 -> Int64+foreign import ccall "c_det.h det_modp" c_det_modp :: Int64 -> CInt  -> Ptr Int64 -> IO Int64++fastDetModP :: Int64 -> Matrix Int64 -> Int64+fastDetModP p mat = Unsafe.unsafePerformIO $ ioFastDetModP p mat++ioFastDetModP :: Int64 -> Matrix Int64 -> IO Int64+ioFastDetModP p mat = do+  let ((1,1),(n,_)) = bounds mat+  withArray (elems mat) $ \ptr -> c_det_modp p (fromIntegral n :: CInt) ptr++gaussElimDeterminantInt64 :: Matrix Zp -> Zp+gaussElimDeterminantInt64 mat = +  Unsafe.unsafePerformIO $ do+    let pp = fromIntegral p :: Int64+    let ((1,1),(n,_)) = bounds mat+        xs = map (fromIntegral . fromZp) (elems mat) :: [Int64]+    d <- withArray xs $ \ptr -> c_det_modp pp (fromIntegral n :: CInt) ptr+    return $ Zp $ fromIntegral d++--------------------------------------------------------------------------------+-- * Bareiss determinant algorithm++type STMatrix s a = STArray s (Int,Int) a++-- | Works only if the top-left minors all have nonzero determinants+{-# SPECIALIZE bareissDeterminantFullRank :: Matrix Integer -> Integer #-}+{-# SPECIALIZE bareissDeterminantFullRank :: Matrix Int     -> Int     #-}+bareissDeterminantFullRank :: forall a . Integral a => Matrix a -> a+bareissDeterminantFullRank mat = ++  if n>0 +    then runST $ do+      ar1   <- thaw mat       :: ST s (STMatrix s a)  +      ar2   <- newArray_ siz  :: ST s (STMatrix s a)+      last  <- newSTRef 1     :: ST s (STRef s a)+      (ar,_) <- foldM (worker last) (ar1,ar2) [1..n-1] +      readArray ar (n,n) +    else 1  -- determinant of the empty matrix is 1++  where ++    siz@((1,1),(n,_)) = bounds mat++    unsafeReadArray :: STMatrix s a -> (Int,Int) -> ST s a+    unsafeReadArray ar ij = unsafeRead ar (index siz ij)++    unsafeWriteArray :: STMatrix s a -> (Int,Int) -> a -> ST s ()+    unsafeWriteArray ar ij x = unsafeWrite ar (index siz ij) x++    worker :: STRef s a -> (STMatrix s a, STMatrix s a)  -> Int -> ST s (STMatrix s a, STMatrix s a)+    worker last (ar1,ar2) !k = do+      q <- readSTRef last             ++      when (q==0) $ unsafeIOToST $ do+        putStrLn "divison by zero while computing the determinant..."++      forM_ [k+1..n] $ \(!i) -> +        forM_ [k+1..n] $ \(!j) -> do+          a <- unsafeReadArray ar1 (k,k)+          b <- unsafeReadArray ar1 (i,k)+          c <- unsafeReadArray ar1 (k,j)+          d <- unsafeReadArray ar1 (i,j)+          unsafeWriteArray ar2 (i,j) $ (a*d - b*c) `div` q      +      unsafeReadArray ar1 (k,k) >>= writeSTRef last +      return (ar2,ar1)++--------------------------------------------------------------------------------+-- * Gaussian elimination++{-# SPECIALIZE gaussElimDeterminant :: Matrix Rational -> Rational #-}+{-# SPECIALIZE gaussElimDeterminant :: Matrix Zp       -> Zp       #-}+gaussElimDeterminant :: forall a. (Eq a, Show a, Fractional a) => Matrix a -> a+gaussElimDeterminant mat =  ++  if n <= 0 +    then 1             -- determinant of the empty matrix is 1+    else runST $ do+      -- unsafeIOToST (printMatrix mat >> putStrLn "")+      neg <- newSTRef False +      arr <- thaw mat :: ST s (STMatrix s a)  +      worker neg arr 1++  where ++    siz@((1,1),(n,_)) = bounds mat++    unsafeReadArray :: STMatrix s a -> (Int,Int) -> ST s a+    unsafeReadArray !ar !ij = unsafeRead ar (index siz ij)++    unsafeWriteArray :: STMatrix s a -> (Int,Int) -> a -> ST s ()+    unsafeWriteArray !ar !ij !x = unsafeWrite ar (index siz ij) x++    finish :: STRef s Bool -> STMatrix s a -> ST s a+    finish !neg !arr = do+      diag <- sequence [ unsafeReadArray arr (i,i) | i<-[1..n] ]        +      b    <- readSTRef neg+      return $ if b +        then negate $ product diag+        else          product diag++    worker :: STRef s Bool -> STMatrix s a -> Int -> ST s a+    worker !neg !arr !i = if i >= n +      then finish neg arr+      else do+        ps <- sequence [ unsafeReadArray arr (i,j) | j<-[i..n] ]+        case findIndex (/=0) ps of+          Nothing    -> return 0                    -- no pivot -> line is full zero -> determinant is zero+          Just pivot -> cont neg arr i (i+pivot)++    cont :: STRef s Bool -> STMatrix s a -> Int -> Int -> ST s a+    cont !neg !arr !i !pivot = do+--      printST (i,pivot)+      when (pivot > i) $ xchg neg arr i pivot+      p <- unsafeReadArray arr (i,i)+      forM_ [i+1..n] $ \k -> do+        q <- unsafeReadArray arr (k,i)+        unsafeWriteArray arr (k,i) 0+        let z = q / p+        forM_ [i+1..n] $ \j -> do+          a <- unsafeReadArray arr (i,j)+          b <- unsafeReadArray arr (k,j)+          unsafeWriteArray arr (k,j) (b - a*z)              +      worker neg arr (i+1)  ++    xchg :: STRef s Bool -> STMatrix s a -> Int -> Int -> ST s ()+    xchg !neg !arr !i !j = do+      modifySTRef neg not             -- exchanging two rows flip the sign of the determinant+      forM_ [i..n] $ \k -> do+        a <- unsafeReadArray arr (k,i)+        b <- unsafeReadArray arr (k,j)        +        unsafeWriteArray arr (k,j) a+        unsafeWriteArray arr (k,i) b++--------------------------------------------------------------------------------+-- * naive determinant algorithm (for testing purposes)++naiveDeterminant :: forall a. (Num a) => Matrix a -> a+naiveDeterminant mat+  | n <= 0    = 1+  | n == 1    = mat!(1,1)+  | n == 2    = mat!(1,1) * mat!(2,2) - mat!(1,2) * mat!(2,1)+  | otherwise = worker [1..n] [1..n]+  where++    siz@((1,1),(n,_)) = bounds mat++    signs = cycle [True,False]++    worker []     []     = 1+    worker [a]    [b]    = mat!(a,b)+    worker [a,b]  [p,q]  = mat!(a,p) * mat!(b,q) -  mat!(a,q) * mat!(b,p)+    worker (i:is) js     = foldl' (+) 0 (zipWith f signs js) where+      f b j = if b +        then          mat!(i,j) * worker is (js\\[j])+        else negate $ mat!(i,j) * worker is (js\\[j])+++--------------------------------------------------------------------------------+-- * random matrices++mkSquareMatrix :: (Int -> Int -> a) -> Int -> Matrix a+mkSquareMatrix f n = array ((1,1),(n,n)) [ ((i,j) , f i j ) | i<-[1..n] , j<-[1..n] ]++testMatrix :: Num a => Int -> Matrix a+testMatrix n = mkSquareMatrix f n where+  f i j = fromIntegral +        $ 3 + i*i*i - j*j + (4*i*j + 3*i + 5*j + 7) + xor (13+i) (17+j) where++randomMatrix :: (Random a, Num a) => Int -> IO (Matrix a)+randomMatrix = randomMatrix' 10++randomMatrix' :: (Random a, Num a) => a -> Int -> IO (Matrix a)+randomMatrix' bnd n = do+  xs <- replicateM (n*n) (randomRIO (-bnd,bnd))+  return $ listArray ((1,1),(n,n)) xs++printST :: Show a => a -> ST s ()+printST x = unsafeIOToST (print x)++--------------------------------------------------------------------------------+-- * testing++test = do+  forM_ [1..10] $ \n -> do+    putStrLn $ "testing matrices of size " ++ show n ++ " x " ++ show n ++ "..."+    replicateM_ 100 $ do+      imat <- randomMatrix n :: IO (Matrix Integer)+      let mat = fmap fromInteger imat :: Matrix Rational++      let a = naiveDeterminant     mat+          b = gaussElimDeterminant mat++      let ia = naiveDeterminant  imat :: Integer+          amodp = mkZp $ fromIntegral (mod ia (fromIntegral p))++      let c  = gaussElimDeterminant (fmap mkZp imat)+          d0 = fastDetModP (fromIntegral p) (fmap (\a -> fromIntegral (mod a (fromIntegral p))) imat)+          d  = fromIntegral d0 :: Zp++      when (a/=b) $ do+        putStrLn "\nERROR!"+        print (a,b)+        print imat++      when (c/=d || d/=amodp) $ do+        putStrLn "\nC ERROR!"+        print (c,d,amodp)+        print imat++
+ src/Math/Algebra/ModP.hs view
@@ -0,0 +1,117 @@++-- | Prime fields. +--+-- TODO: do it properly; and fast implementation for specialized prime fields+--++{-# LANGUAGE BangPatterns #-}+module Math.Algebra.ModP where++--------------------------------------------------------------------------------++import Data.Bits+import Data.Ratio+import Data.Int++--------------------------------------------------------------------------------++-- | @2^31-1@ is a prime (in practice this seems to be significantly faster than @2^63-25@)+p :: Int64+p = 2^31 - 1      ++-- p = 20551         -- max coefficient in 3/1/8 is 20460+-- p = 2^31 - 1      -- @2^31-1@ +-- p = 2^33 - 9      -- @2^33-9@+-- p = 2^62 - 57     -- @2^62-57@+-- p = 2^63 - 25     -- @2^63-25@+++--------------------------------------------------------------------------------++newtype Zp = Zp Int64 deriving (Eq, Show)++fromZp :: Zp -> Int+fromZp (Zp k) = fromIntegral k++mkZp :: Integral a => a -> Zp+mkZp n = Zp (mod (fromIntegral n) p)++--------------------------------------------------------------------------------++instance Num Zp where+  (+)          = addZp +  (-)          = subZp +  (*)          = mulZp +  fromInteger  = mkZp . fromInteger+  abs          = id+  signum _     = Zp 1++instance Fractional Zp where+  recip (Zp a)   = mkZp $ invZp_euclid a+  a / b          = a * recip b+  fromRational r = fromInteger (numerator r) / fromInteger (denominator r)++--------------------------------------------------------------------------------++addZp :: Zp -> Zp -> Zp+addZp (Zp a) (Zp b) +  | c <  0    = Zp (c - p)               -- overflow+  | c >= p    = Zp (c - p)+  | otherwise = Zp  c+  where+    c = a + b++subZp :: Zp -> Zp -> Zp+subZp (Zp a) (Zp b) = Zp (if b<=a then a-b else a+p-b)++mulZp :: Zp -> Zp -> Zp+mulZp (Zp a0) (Zp b0) = Zp (fromInteger c) where+  a = fromIntegral a0 :: Integer                    -- because Int can overflow :(+  b = fromIntegral b0 :: Integer+  c = mod (a * b) (fromIntegral p)++-- | Inverse using the binary Euclidean algorithm +invZp_euclid :: Int64 -> Int64+invZp_euclid a +  | a == 0     = 0+  | otherwise  = go 1 0 a p+  where+  +    modp :: Int64 -> Int64+    modp n = mod n p++    halfp1 = shiftR (p+1) 1++    go :: Int64 -> Int64 -> Int64 -> Int64 -> Int64+    go !x1 !x2 !u !v +      | u==1       = x1+      | v==1       = x2+      | otherwise  = stepU x1 x2 u v++    stepU :: Int64 -> Int64 -> Int64 -> Int64 -> Int64+    stepU !x1 !x2 !u !v = if even u +      then let u'  = shiftR u 1+               x1' = if even x1 then shiftR x1 1 else shiftR x1 1 + halfp1+           in  stepU x1' x2 u' v+      else     stepV x1  x2 u  v++    stepV :: Int64 -> Int64 -> Int64 -> Int64 -> Int64+    stepV !x1 !x2 !u !v = if even v+      then let v'  = shiftR v 1+               x2' = if even x2 then shiftR x2 1 else shiftR x2 1 + halfp1+           in  stepV x1 x2' u v' +      else     final x1 x2  u v++    final :: Int64 -> Int64 -> Int64 -> Int64 -> Int64+    final !x1 !x2 !u !v = if u>=v++      then let u'  = u-v+               x1' = if x1 >= x2 then modp (x1-x2) else modp (x1+p-x2)               +           in  go x1' x2  u' v ++      else let v'  = v-u+               x2' = if x2 >= x1 then modp (x2-x1) else modp (x2+p-x1)+           in  go x1  x2' u  v'++--------------------------------------------------------------------------------+
+ src/Math/Algebra/Schur.hs view
@@ -0,0 +1,164 @@++-- | Schur polynomials++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns #-}+module Math.Algebra.Schur where++--------------------------------------------------------------------------------++import Control.Monad+import Control.Monad.ST++import Data.Array.Base+import Data.Array.IArray+import Data.Array.MArray+import Data.Array.Unsafe+import Data.Array.ST++import Data.List+import Data.Ratio+import Data.STRef++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Sets++import qualified Data.Map as Map++import Debug.Trace+import GHC.IO ( unsafeIOToST )++import Math.Algebra.Determinant+import Math.Algebra.ModP++--------------------------------------------------------------------------------++-- segre :: Num a => Int -> [a] -> a+-- segre k xs = sum $ map product $ combine k xs++--------------------------------------------------------------------------------+-- * Elementary and complete symmetric polynomials++-- | Precalc chern classes+elemSymmArray :: forall a . Num a => [a] -> Array Int a+elemSymmArray xs = +  runST $ do+    ar <- newArray (1,n) 0 :: ST s (STArray s Int a)+    mapM_ (worker ar) (zip [1..n] xs)+    unsafeFreeze ar+  where+    n = length xs+    worker ar (i,x) = +      forM_ [i,i-1..1] $ \j -> do+        a  <- lkp ar    j+        b  <- lkp ar (  j - 1 )+        writeArray ar j (a + x*b)+    lkp ar j = if j>=1 +      then  readArray ar j +      else  return 1++-- | Precalc segre classes+completeSymmArray :: forall a . Num a => Int -> [a] -> Array Int a+completeSymmArray m xs = +  runST $ do+    ar <- newArray ((1,1),(n,m)) 0 :: ST s (STArray s (Int,Int) a)+    mapM_ (worker ar) (zip [1..n] xs)+    ys <- forM [1..m] $ \j -> readArray ar (n,j)+    return $ listArray (1,m) ys+  where+    n = length xs++    worker :: (STArray s (Int,Int) a) -> (Int,a) -> ST s ()+    worker ar (i,x) = +      forM_ [1..m] $ \j -> do+        a  <- lkp ar (i-1) (j  )+        b  <- lkp ar (i  ) (j-1)+        writeArray ar (i,j) (a + x*b)++    lkp ar i j +      | j>=1 && i>=1  = readArray ar (i,j)+      | j==0          = return 1+      | i==0          = return 0++--------------------------------------------------------------------------------+-- * Schur polynomials+          +schurMatrixChern :: Num a => (Int -> a) -> Partition -> Matrix a+schurMatrixChern c shape = schurMatrixSegre c (dualPartition shape)++schurMatrixSegre :: Num a => (Int -> a) -> Partition -> Matrix a+schurMatrixSegre s shape = matrix where+  matrix = array ((1,1),(n,n)) entries+  n = height (dualPartition shape)+  f k  | k  <  0  =  0+       | k  == 0  =  1+       | k  >  0  =  s k +  entries = [ ( (i,j) , f (k + j - i) ) | (i,k) <- zip [1..n] shape' , j<-[1..n] ]+  shape' = fromPartition shape ++ repeat 0++--------------------------------------------------------------------------------++{-++{-# SPECIALIZE schurDeterminantChern :: (Int -> Integer) -> Partition -> Integer #-}+{-# SPECIALIZE schurDeterminantSegre :: (Int -> Integer) -> Partition -> Integer #-}  ++-- | Jacobi-Trudi formula+schurDeterminantChern :: Integral a => (Int -> a) -> Partition -> a+schurDeterminantChern chern = bareissDeterminantFullRank . schurMatrixChern chern++schurDeterminantSegre :: Integral a => (Int -> a) -> Partition -> a+schurDeterminantSegre segre = bareissDeterminantFullRank . schurMatrixSegre segre++schurFromChernArray :: Integral a => Array Int a -> Partition -> a+schurFromChernArray ar part = schurDeterminantChern f part where+  (1,n) = bounds ar+  f k | k<=n  =  ar!k+      | k> n  =  0++schurFromSegreArray :: Integral a => Array Int a -> Partition -> a+schurFromSegreArray ar part = schurDeterminantSegre f part where+  (1,n) = bounds ar+  f k | k<=n = ar!k+      | k>n  = error $ "schur-segre " ++ show k ++ " " ++ show n ++ " " ++ show part+-}++--------------------------------------------------------------------------------++schurDeterminantChern :: (Determinant a) => (Int -> a) -> Partition -> a+schurDeterminantChern chern = determinant . schurMatrixChern chern++schurDeterminantSegre :: (Determinant a)  => (Int -> a) -> Partition -> a+schurDeterminantSegre segre = determinant . schurMatrixSegre segre++schurFromChernArray :: (Determinant a)  => Array Int a -> Partition -> a+schurFromChernArray ar part = schurDeterminantChern f part where+  (1,n) = bounds ar+  f k | k<=n  =  ar!k+      | k> n  =  0++schurFromSegreArray :: (Determinant a)  => Array Int a -> Partition -> a+schurFromSegreArray ar part = schurDeterminantSegre f part where+  (1,n) = bounds ar+  f k | k<=n = ar!k+      | k>n  = error $ "schur-segre " ++ show k ++ " " ++ show n ++ " " ++ show part+       +--------------------------------------------------------------------------------+{-++-- * caching++makeSegreSchurCache :: forall s. Array Int Integer -> ST s (Partition -> ST s Integer)+makeSegreSchurCache ar = do+  cacheRef <- newSTRef Map.empty :: ST s (STRef s (Map.Map Partition Integer))+  let fun !part = do+        table <- readSTRef cacheRef+        case Map.lookup part table of+          Just y   -> return y+          Nothing  -> do+            let y = schurFromSegreArray ar part+            writeSTRef cacheRef $! Map.insert part y table+            return y+  return fun+-}+--------------------------------------------------------------------------------
+ src/Math/FreeModule/Class.hs view
@@ -0,0 +1,94 @@++-- | Class interface to different free module implementations.+--+-- Free modules are like maps from a base type to a numeric type,+-- with the additional invariant that the values are never zero.++{-# LANGUAGE TypeFamilies, FlexibleContexts, CPP #-}+module Math.FreeModule.Class where++--------------------------------------------------------------------------------  ++-- | generic baseMap implementation, converts to list and back.+baseMap  :: (FreeModule x, FreeModule y, Coeff x ~ Coeff y) => (Base x -> Base y) -> x -> y +baseMap f = fromList . map h . toList where h (b,c) = (f b, c)++-- | generic coeffMap implementation, converts to list and back.+coeffMap :: (FreeModule x, FreeModule y, Base x ~ Base y) => (Coeff x -> Coeff y) -> x -> y +coeffMap g = fromList . map h . toList  where h (b,c) = (b, g c)++--------------------------------------------------------------------------------  ++class (Ord (Base a), Eq (Coeff a), Num (Coeff a)) => FreeModule a where+  type Base  a :: *+  type Coeff a :: *++  isZero    :: a -> Bool+  zero      :: a+  fromBase  :: Base a -> a+  fromTerm  :: Base a -> Coeff a -> a+  (^+^)     :: a -> a -> a+  (^-^)     :: a -> a -> a+  neg       :: a -> a+  scalarMul :: Coeff a -> a -> a++  -- | We should call the function even when the given base is present  +  -- only in one of the arguments! So that @unionWith (-)@ works correctly.+  unionWith  :: (Coeff a -> Coeff a -> Coeff a) -> a -> a -> a+  +  coeff      :: Base a -> a -> Coeff a+  +  size       :: a -> Int +  minTerm    :: a -> (Base a, Coeff a)+  maxTerm    :: a -> (Base a, Coeff a)+  +  -- | split into two approximately equal parts @x@ and @y@, such that+  -- @maxTerm x < minTerm y@+  split      :: a -> (a, a)+  -- | we assume that @maxTerm x < minTerm y@+  unsafeJoin :: a -> a -> a+  +  toList     :: a -> [(Base a, Coeff a)]+  fromList   :: [(Base a, Coeff a)] -> a+  +  fromAscendingList :: [(Base a, Coeff a)] -> a+    +  isZero x = (size x == 0)+  neg x    = scalarMul (-1) x+  x ^+^ y  = unionWith (+) x y+  x ^-^ y  = unionWith (-) x y -- x ^+^ (neg y)+  fromAscendingList = fromList+  fromBase b = fromTerm b 1+  fromTerm b c = scalarMul c (fromBase b)+  +--------------------------------------------------------------------------------    ++(*^) :: FreeModule a => Coeff a -> a -> a +(*^) = scalarMul++(^*) :: FreeModule a => a -> Coeff a -> a+(^*) = flip scalarMul++infixl 6 ^+^+infixl 6 ^-^++infixl 7 *^+infixl 7 ^*++--------------------------------------------------------------------------------  ++lookupTerm :: FreeModule a => Base a -> a -> Maybe (Base a, Coeff a)+lookupTerm b x = +  case coeff b x of+    0 -> Nothing+    c -> Just (b,c)++minTermMaybe :: FreeModule a => a -> Maybe (Base a, Coeff a)+minTermMaybe x = if isZero x then Nothing else Just (minTerm x)++maxTermMaybe :: FreeModule a  => a -> Maybe (Base a, Coeff a)+maxTermMaybe x = if isZero x then Nothing else Just (maxTerm x)+     +--------------------------------------------------------------------------------  ++
+ src/Math/FreeModule/Helper.hs view
@@ -0,0 +1,26 @@+++-- | misc helper functions++module Math.FreeModule.Helper where++--------------------------------------------------------------------------------++import Data.Ord+import Data.List++--------------------------------------------------------------------------------++(<#>) :: (a -> b) -> (c -> d) -> (a,c) -> (b,d)+(f<#>g) (x,y) = (f x, g y)+ +equating :: Eq b => (a -> b) -> a -> a -> Bool+equating f x y = (f x == f y)+ +sortByFst :: Ord b => [(b,c)] -> [(b,c)]+sortByFst = sortBy (comparing fst)++filterNotZero :: (Eq c, Num c) => [(b,c)] -> [(b,c)]+filterNotZero = filter (\(b,c) -> (c/=0))++--------------------------------------------------------------------------------
+ src/Math/FreeModule/PP.hs view
@@ -0,0 +1,52 @@++-- | More concrete prettyprinting.+-- this should be a separated package.+ +{-# LANGUAGE TypeFamilies, FlexibleContexts, FlexibleInstances, CPP #-}+module Math.FreeModule.PP where++--------------------------------------------------------------------------------++import Data.Ratio++import Math.FreeModule.Class+import Math.FreeModule.Symbol+import Math.FreeModule.PrettyPrint++import qualified Math.FreeModule.SortedList as SL++import Math.Algebra.ModP++--------------------------------------------------------------------------------++class Pretty a where+  pretty :: a -> String+  +pp :: Pretty a => a -> IO ()+pp = putStrLn . pretty++--------------------------------------------------------------------------------++instance Pretty Symbol where pretty = showSymbol++--instance (Pretty b, Ord b, Real c, Show c) => Pretty (SL.FreeMod b c) where+--  pretty = bracket (prettyPrintRealWith pretty)++instance (Pretty b, Ord b) => Pretty (SL.FreeMod b Zp) where+  pretty = bracket (prettyPrintArbWith pretty showZp)++showZp :: Zp -> String+showZp (Zp n) = show n++instance (Pretty b, Ord b) => Pretty (SL.FreeMod b Integer) where+  pretty = bracket (prettyPrintRealWith' show pretty)++instance (Pretty b, Ord b) => Pretty (SL.FreeMod b Rational) where+  pretty = bracket (prettyPrintRealWith' showRational pretty)++showRational :: Rational -> String+showRational r = if denominator r == 1 +  then show (numerator r)+  else "(" ++ show (numerator r) ++ "/" ++ show (denominator r) ++ ")"+  +--------------------------------------------------------------------------------
+ src/Math/FreeModule/Parser.hs view
@@ -0,0 +1,131 @@++{-# LANGUAGE TypeFamilies, CPP #-}+module Math.FreeModule.Parser where++--------------------------------------------------------------------------------++import Control.Monad+import Text.ParserCombinators.Parsec++import Math.FreeModule.Class+import Math.FreeModule.Symbol++--------------------------------------------------------------------------------++type Par s a = GenParser Char s a++--------------------------------------------------------------------------------++-- | Parses @\"alpha[5]\"@ style symbols+symbolP :: Par s Symbol+symbolP = do+  n <- many1 alphaNum+  i <- option Nothing $ do+    char '['+    xs <- many1 digit+    char ']'+    return $ Just (read xs :: Int)+  return (Symbol n i)+  +-- | Parses @\"e2\"@ style symbols+symbolP' :: Par s Symbol+symbolP' = do+  n <- many1 letter+  i <- option Nothing $ do+    xs <- many1 digit+    return $ Just (read xs :: Int)+  return (Symbol n i)++--------------------------------------------------------------------------------++integerP :: Par s Integer+integerP = do+  s <- option 1 signP+  xs <- many1 digit+  return $ s * (read xs)+  +--------------------------------------------------------------------------------++signP :: Num a => Par s a +signP = do+  c <- oneOf "+-"+  return $ case c of { '+' -> 1 ; '-' -> (-1) }++betweenSpaces :: Par s a -> Par s a+betweenSpaces p = do+  spaces+  x <- p+  spaces+  return x  ++--------------------------------------------------------------------------------++notEmpty :: GenParser tok st a -> GenParser tok st a+notEmpty parser = do+  pos1 <- getPosition+  x <- parser+  pos2 <- getPosition+  if (pos1 == pos2)+    then fail "empty"+    else return x++-- this is useful for exterior algebras, for example. +freeModuleP' :: FreeModule a => Par s (Base a,Coeff a) -> Par s (Coeff a) -> Par s a+freeModuleP' baseP coeffP = try p <|> q where+  p = betweenSpaces (string "0") >> eof >> return zero+  q = liftM fromList $ do+    xs <- liftM helper $ many1 (termP baseP coeffP) +    spaces+    eof+    return xs+  helper = map $ \((b,c1),c2) -> (b,c1*c2)+    +freeModuleP :: FreeModule a => Par s (Base a) -> Par s (Coeff a) -> Par s a+freeModuleP baseP coeffP = try p <|> q where+  p = betweenSpaces (string "0") >> eof >> return zero+  q = liftM fromList $ do+    xs <- many1 (termP baseP coeffP) +    spaces+    eof+    return xs+    +termP :: Num c => Par s b -> Par s c -> Par s (b,c)+termP baseP coeffP = +  do+    s <- option 1 (betweenSpaces signP)  +    (b,c) <- try q <|> p+    return (b,s*c)+  where +    p = do+      b <- notEmpty baseP  +      return (b,1)+    q = do+      c <- coeffP+      optional (betweenSpaces (char '*'))+      b <- baseP  +      return (b,c)+{-+  s <- option 1 (betweenSpaces signP)  +  c <- option 1 $ do+    c <- coeffP+    optional (betweenSpaces (char '*'))+    return c+  b <- baseP  +  return (b,s*c)+-}  ++--------------------------------------------------------------------------------+  +parseLinearExpr :: (FreeModule a, Base a ~ Symbol, Coeff a ~ Integer) => String -> a+parseLinearExpr = parseFreeModule symbolP integerP++parseFreeModule :: FreeModule a => Parser (Base a) -> Parser (Coeff a) -> String -> a+parseFreeModule baseP coeffP s =+  case runParser p () "input" s of+    Left err -> error (show err)+    Right x  -> x+  where +    p = freeModuleP baseP coeffP+  +--------------------------------------------------------------------------------+
+ src/Math/FreeModule/PrettyPrint.hs view
@@ -0,0 +1,56 @@++{-# LANGUAGE TypeFamilies, FlexibleContexts, CPP #-}+module Math.FreeModule.PrettyPrint where++--------------------------------------------------------------------------------++import Math.FreeModule.Class++--------------------------------------------------------------------------------++bracket :: (a -> String) -> a -> String+bracket f x = "(" ++ f x ++ ")"++-- | Print stuff with real (eg integral or rational) coefficients+prettyPrintRealWith +  :: (FreeModule x, Real (Coeff x), Show (Coeff x)) +  => (Base x -> String) -> x -> String+prettyPrintRealWith showBase x = s where+  y = toList x+  s = if isZero x +    then "0"+    else if take 3 t == " + " then drop 3 t else t+  t = concatMap h y+  h (b,c) = (if c<0 then " - " else " + ") ++ show (abs c) ++ t+    where t = case showBase b of+            "" -> "" +            xs -> "*" ++ xs++prettyPrintRealWith'+  :: (FreeModule x, Real (Coeff x), Show (Coeff x)) +  => (Coeff x -> String) -> (Base x -> String) -> x -> String+prettyPrintRealWith' showCoeff showBase x = s where+  y = toList x+  s = if isZero x +    then "0"+    else if take 3 t == " + " then drop 3 t else t+  t = concatMap h y+  h (b,c) = (if c<0 then " - " else " + ") ++ showCoeff (abs c) ++ t+    where t = case showBase b of+            "" -> "" +            xs -> "*" ++ xs+               +-- | Print stuff with arbitrary coefficients+prettyPrintArbWith +  :: (FreeModule x) +  => (Base x -> String) -> (Coeff x -> String) -> x -> String+prettyPrintArbWith showBase showCoeff x = s where+  y = toList x+  s = if isZero x +    then "0"+    else drop 3 t+  t = concatMap h y+  h (b,c) = " + " ++ showCoeff c ++ "*" ++ showBase b++--------------------------------------------------------------------------------+
+ src/Math/FreeModule/SortedList.hs view
@@ -0,0 +1,105 @@++-- | Free modules implemented as sorted lists of @(base,coeff)@ pairs.+-- The functions 'coeff', 'maxTerm', 'split', 'unsafeJoin' are slow +-- in this implementation.++{-# LANGUAGE TypeFamilies, DeriveFunctor #-}+module Math.FreeModule.SortedList+  ( module Math.FreeModule.Class  +  , baseMap+  , coeffMap+  , FreeMod+  , ZModule+  , QModule+  )+  where++--------------------------------------------------------------------------------++import Data.List+import Data.Ord++import Math.FreeModule.Class hiding (baseMap,coeffMap)+import Math.FreeModule.PrettyPrint+import Math.FreeModule.Helper++--------------------------------------------------------------------------------++newtype FreeMod b c = S [(b,c)] deriving (Eq,Ord,Show,Functor)++type ZModule b = FreeMod b Integer+type QModule b = FreeMod b Rational++--------------------------------------------------------------------------------++-- hackish solution to implementation-specific baseMap/coeffMap:+-- import this module only, which hides the generic implementation+baseMap :: Ord b => (a -> b) -> FreeMod a c -> FreeMod b c+baseMap = sortedlistBaseMap++coeffMap :: (c -> d) -> FreeMod b c -> FreeMod b d+coeffMap = sortedlistCoeffMap++sortedlistBaseMap :: Ord b => (a -> b) -> FreeMod a c -> FreeMod b c+sortedlistBaseMap  f (S xs) = S (sortByFst (map (f<#>id) xs))++sortedlistCoeffMap :: (c -> d) -> FreeMod b c -> FreeMod b d+sortedlistCoeffMap g (S xs) = S (map (id<#>g) xs)++-- does not work?+{- RULES "baseMap/SortedList"  baseMap  = slBaseMap  -}+{- RULES "coeffMap/SortedList" coeffMap = slCoeffMap -}++--------------------------------------------------------------------------------++instance (Ord b, Eq c, Num c) => FreeModule (FreeMod b c) where++  type Base  (FreeMod b c) = b+  type Coeff (FreeMod b c) = c+  +  isZero (S xs) = case xs of { [] -> True ; _ -> False }+  zero = S []+  fromBase b   = S [(b,1)]+  fromTerm b c = S [(b,c)]+  scalarMul c (S xs) = S (map (id<#>(*c)) xs)++  coeff b (S xs) = case lookup b xs of+    Nothing -> 0+    Just c  -> c+  +  unionWith f (S xs) (S ys) = S (unionWorker f xs ys)+  +  size (S xs) = length xs+  +  minTerm (S xs) = case xs of +    [] -> error "minTerm: empty"+    _  -> head xs+  maxTerm (S xs) = case xs of +    [] -> error "maxTerm: empty"+    _  -> last xs+    +  split (S xs) = (S ys, S zs) where (ys,zs) = splitAt (length xs `div` 2) xs+  unsafeJoin (S xs) (S ys) = S (xs++ys)+  +  toList (S xs) = xs+  fromList xs = S $ filterNotZero $ collapse $ sortByFst $ xs where+    collapse = map f . groupBy (equating fst) +    f xs = (fst (head xs), sum (map snd xs))+  fromAscendingList = S++--------------------------------------------------------------------------------+  +unionWorker :: (Ord b, Eq c, Num c) => (c -> c -> c) -> [(b,c)] -> [(b,c)] -> [(b,c)]+unionWorker f xs [] = map (\(b,x) -> (b, f x 0)) xs+unionWorker f [] ys = map (\(b,y) -> (b, f 0 y)) ys+unionWorker f xxs@(x@(b1,c1):xs) yys@(y@(b2,c2):ys) = +  case compare b1 b2 of+    LT -> g b1 c1 0  (unionWorker f xs  yys)+    GT -> g b2 0  c2 (unionWorker f xxs ys ) +    EQ -> g b1 c1 c2 (unionWorker f xs  ys )+  where+    g b c1 c2 rest = case f c1 c2 of+      0 -> rest+      c -> (b,c) : rest+        +--------------------------------------------------------------------------------
+ src/Math/FreeModule/Symbol.hs view
@@ -0,0 +1,90 @@++-- | Possibly indexed symbols.++module Math.FreeModule.Symbol where++--------------------------------------------------------------------------------++import Data.Set (Set)+import qualified Data.Set as Set++--------------------------------------------------------------------------------++data Symbol = Symbol+  { _name  :: String+  , _index :: Maybe Int+  }+  deriving (Eq,Ord,Show)++--------------------------------------------------------------------------------++-- | Shows the symbols in @\"alpha[5]\"@ style+showSymbol :: Symbol -> String+showSymbol (Symbol name idx) = case idx of+  Just j  -> name ++ "[" ++ show j ++ "]"+  Nothing -> name++-- | Shows the symbols in @\"alpha5\"@ style+showSymbol' :: Symbol -> String+showSymbol' (Symbol name idx) = case idx of+  Just j  -> name ++ show j +  Nothing -> name++-- | Shows the symbols in @\"\\alpha_{5}\"@ style+showSymbolLatex :: Symbol -> String+showSymbolLatex (Symbol name idx) = +  case idx of+    Just j  -> name' ++ "_{" ++ show j ++ "}"+    Nothing -> name'+  where +    name' = if Set.member name latexGreek+      then '\\' : name+      else name+      +--------------------------------------------------------------------------------++latexGreek :: Set String+latexGreek = Set.fromList (latexSmallGreek ++ latexCapitalGreek)++latexSmallGreek :: [String]+latexSmallGreek =+  [ "alpha"+  , "beta"+  , "gamma"+  , "delta"+  , "epsilon"+  , "zeta"+  , "eta"+  , "theta"+  , "iota"+  , "kappa"+  , "lambda"+  , "mu"+  , "nu"+  , "xi"+  , "pi"+  , "rho"+  , "sigma"+  , "tau"+  , "upsilon"+  , "phi"+  , "chi"+  , "psi"+  , "omega"+  ]+  +latexCapitalGreek :: [String]+latexCapitalGreek =+  [ "Gamma"+  , "Delta"+  , "Theta"+  , "Lambda"+  , "Xi"+  , "Pi"+  , "Sigma"+  , "Upsilon"+  , "Phi"+  , "Psi"+  ]+  +--------------------------------------------------------------------------------
+ src/Math/ThomPoly/Formulae.hs view
@@ -0,0 +1,17 @@++-- | For some special cases (like @j=1@), we have explicit formulae.+--+-- TODO: implement them! (this is just a placeholder)+++module Math.ThomPoly.Formulae where++--------------------------------------------------------------------------------++data SigmaI1 = SigmaI1+  { _i :: !Int+  , _n :: !Int+  }+  deriving (Eq,Show)++--------------------------------------------------------------------------------
+ src/Math/ThomPoly/Shared.hs view
@@ -0,0 +1,215 @@++-- | Shared code++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, PackageImports,+             TypeSynonymInstances, FlexibleInstances, FlexibleContexts,+             ExistentialQuantification +  #-}+module Math.ThomPoly.Shared where++--------------------------------------------------------------------------------++import Data.List+import Data.Ratio+import Data.Proxy++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Sets++import Math.FreeModule.Symbol+import Math.FreeModule.SortedList+import Math.FreeModule.PrettyPrint+import Math.FreeModule.PP+-- import FreeModule.Parser++import Math.Algebra.ModP+import Math.Algebra.Schur+import Math.Algebra.Determinant++--------------------------------------------------------------------------------+-- * Rings and fields++data AnyRing = forall r. CoeffRing r => AnyRing (Proxy r)++solveAny :: Problem problem => AnyRing -> Batch -> problem -> FreeMod Schur Integer+solveAny anyring batch prob = case anyring of+  AnyRing pxy -> solveAndProject pxy batch prob++ringZZ, ringQQ, ringZp :: AnyRing+ringZZ = AnyRing (Proxy :: Proxy Integer )+ringQQ = AnyRing (Proxy :: Proxy Rational)+ringZp = AnyRing (Proxy :: Proxy Zp      )++--------------------------------------------------------------------------------++class +  ( Eq a , Num a , Show a , Determinant a+  , Eq (FieldOfFractions a) , Show (FieldOfFractions a) , Fractional (FieldOfFractions a)+  , Pretty (Term a)+  ) => CoeffRing a +  where+    type FieldOfFractions a :: *+    embed    :: a -> FieldOfFractions a+    project  :: Proxy a -> FieldOfFractions a -> Maybe a+    toBigInt :: Proxy a -> FieldOfFractions a -> Maybe Integer++ratToInt :: Rational -> Maybe Integer+ratToInt x = case denominator x of { 1 -> Just (numerator x) ; _ -> Nothing }++--------------------------------------------------------------------------------++instance CoeffRing Integer where+  type FieldOfFractions Integer = Rational+  embed      = fromInteger+  project  _ = ratToInt+  toBigInt _ = ratToInt++instance CoeffRing Rational where+  type FieldOfFractions Rational = Rational+  embed      = id+  project  _ = Just+  toBigInt _ = ratToInt++instance CoeffRing Zp where+  type FieldOfFractions Zp = Zp+  embed      = id+  project  _ = Just+  toBigInt _ = Just . fromIntegral . fromZp++unsafeProject :: CoeffRing c => Proxy c -> FieldOfFractions c -> Integer+unsafeProject pxy x = case toBigInt pxy x of+  Just y  -> y+  Nothing -> error "cannot project back result"  ++--------------------------------------------------------------------------------+-- * Thom polynomial problems++class Problem problem where+  baseFName :: problem -> String+  calcStats :: problem -> Stats+  solve     :: CoeffRing coeff => Proxy coeff -> Batch -> problem -> FreeMod Schur (FieldOfFractions coeff)++fullFName :: Problem problem => Batch -> problem -> FilePath+fullFName batch prob = baseFName prob ++ batchSuffix batch ++ ".txt"++solveAndProject +  :: forall problem coeff. (Problem problem, CoeffRing coeff)+  => Proxy coeff -> Batch -> problem -> FreeMod Schur Integer+solveAndProject pxy batch prob = coeffMap (unsafeProject pxy) $ solve pxy batch prob where++--------------------------------------------------------------------------------++-- | \"Statistics\" of a problem+data Stats  = Stats +  { _codim0   :: !Int        -- ^ codimension (for @m=n@)+  , _mu       :: !Int        -- ^ algebraic multiplicity (minus 1)+  , _maxPairs :: !Int        -- ^ maximum number of possible non-zero coefficients+  }+  deriving Show++-------------------------------------------------------------------------------+-- * Batches++data Batch = Batch +  { _whichBatch :: !Int +  , _nBatches   :: !Int +  }+  deriving Show++defaultBatch :: Batch+defaultBatch = Batch 1 1++selectBatch :: Batch -> [a] -> [a]+selectBatch (Batch a b) xs  +  | a < 1     = error "selectBatch: a<1"+  | a > b     = error "selectBatch: a>b"+  | b == 1    = xs+  | otherwise = take bsize $ drop ((a-1)*bsize) $ xs+  where+    n     = length xs+    (q,r) = divMod n b+    bsize = case r of+      0 -> q+      _ -> q+1++batchSuffix :: Batch -> String+batchSuffix (Batch a b)+  | b == 1    = ""+  | otherwise = "_batch" ++ show a ++ "of" ++ show b++{-+-- sanity test+testBatch'' b n = concat [ selectBatch (Batch i b) [1..n] | i<-[1..b] ] == [1..n]+testBatch'  b   = and [ testBatch'' b n | n<-[0..1000] ]+testBatch       = and [ testBatch'  b   | b<-[1..100 ] ] +-}+ +--------------------------------------------------------------------------------+-- * Misc++-- type CoeffRing = Zp  -- Rational -- Integer++type Term coeff = FreeMod Symbol coeff++alpha :: CoeffRing coeff => Int -> Term coeff+alpha i = fromBase $ Symbol "alpha" (Just i) ++newtype Schur = Schur Partition deriving (Eq,Ord,Show)++instance Pretty Schur where +  pretty (Schur part) = 's' : show (fromPartition part)++--------------------------------------------------------------------------------+-- * Evaluate++evaluate :: (Num a, FreeModule x) => (Base x -> Coeff x -> a) -> x -> a+evaluate f = sum . map (uncurry f) . toList+         +--------------------------------------------------------------------------------+-- * Signed partitions++-- | Pairs of partition with weights of fix difference, given by +-- the third parameter, @ofs=|pos|-|neg|@, and complementary length;+-- the first giving the positive deviation compared to the box of (m-n+i)*i,+-- and the second giving the negative one.+-- Picture:+--+-- >         m-n+i                    n-i+-- >    +------------------+----------------+---------++-- >    |                  |          _____/          |+-- >  i |      lambda      |  pos____/                |+-- >    |                  |    /                     |+-- >    |                  |   /                      |+-- > mu +................._|__/                       | mu+-- >    |              __/ |       C lambda ~         |+-- >    |          ___/    |                          |+-- >    |         /  neg   |                          |+-- >    +--------+---------+--------------------------++-- >                       m+--+-- The length ("width" in /combinat-speak/, unfortunately) of the partitions+-- are less than mu; the "height" (first element) of @pos@ is at most @(n-i)@,+-- the height of @neg@ is unlimited (well, it is limited by @(mu-1)*(n-i)@ of course).+--+-- Actually, in the case of sigmaij, length(pos)<=i and length(neg)<=mu-i !+--+partitionPairs :: Int -> Int -> Int -> Int -> [(Partition,Partition)]+partitionPairs mu n i ofs = +  [ (pos,neg) +  | d <- [0..i*(n-i)] +  , pos <- partitions' (n-i,i) d+  , let l = width pos+  , neg <- partitions' (d,mu-i) (d-ofs)+  ]++-- | Given the parameters @(m-n+i,mu) (pos,neg)@, this computes @lambda@+-- in the picture above. +posnegPairToPartition :: (Int,Int) -> (Partition,Partition) -> Partition+posnegPairToPartition (h,w) (pos,neg) = toPartitionUnsafe xs where+  xs = zipWith (+) ys (replicate w h)+  ys = pos' ++ replicate (w - width pos - width neg) 0 ++ map negate (reverse neg')+  pos' = fromPartition pos+  neg' = fromPartition neg++--------------------------------------------------------------------------------
+ src/Math/ThomPoly/SigmaI.hs view
@@ -0,0 +1,155 @@++-- | Calculates the Thom polynomial of @Sigma^{i}@ with localization and the substitution trick+-- (for sanity testing only, as we know the answer anyway)+--++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, PackageImports #-}+module Math.ThomPoly.SigmaI where++--------------------------------------------------------------------------------++import Control.Monad+import Control.Monad.ST+import Data.STRef++import Data.Array.IArray+import Data.Array.Unsafe+import Data.Array.ST++import Data.List+import Data.Ratio+import Data.Proxy++import Debug.Trace+import GHC.IO ( unsafeIOToST )++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Sets++import Math.FreeModule.Symbol+import Math.FreeModule.SortedList+-- import Math.FreeModule.PrettyPrint+-- import Math.FreeModule.PP+-- import Math.FreeModule.Parser++import Math.Algebra.ModP+import Math.Algebra.Schur++import Math.ThomPoly.Subs+import Math.ThomPoly.Shared++--------------------------------------------------------------------------------++instance Problem SigmaI where+  calcStats = statsI+  solve     = sigmai+  baseFName (SigmaI i n)  = "sigmai__i" ++ show i ++ "_n" ++ show n+  +--------------------------------------------------------------------------------+-- * @Sigma^{i}@++data SigmaI = SigmaI+  { _i :: !Int            -- ^ corank of the differential+  , _n :: !Int            -- ^ source dimension+  }+ deriving (Eq,Show)++-- | We need @n >= mu@ with this method+smallestI :: Int -> SigmaI+smallestI i = SigmaI i i++-- | The codimension of @Sigma^{i}(n,m)@ is @codim = i*(m-n+i)@+codim :: SigmaI -> Int -> Int+codim (SigmaI i n) m = i * (m-n+i)  ++-- | There is a sign in the localization formula.+signCorrection :: SigmaI -> Int+signCorrection (SigmaI i n) = i*(n-i)++statsI :: SigmaI -> Stats+statsI prob@(SigmaI i n) = +  Stats +    { _mu       = i +    , _codim0   = codim prob n +    , _maxPairs = length posneg+    } +  where+    posneg = partitionPairs mu n i 0+    mu = i+++--------------------------------------------------------------------------------+-- @Sigma^i@  ++type Fixpoint1 = [Int]   + +sigmai :: CoeffRing coeff => Proxy coeff -> Batch -> SigmaI -> FreeMod Schur (FieldOfFractions coeff)+sigmai pxy batch problem@(SigmaI i n) = sigmai' pxy problem (selectBatch batch posneg) where+  posneg = partitionPairs mu n i 0+  mu     = i++sigmai' +  :: forall coeff. CoeffRing coeff +  => Proxy coeff -> SigmaI -> [(Partition,Partition)] -> FreeMod Schur (FieldOfFractions coeff)+sigmai' _ problem@(SigmaI i n) posneg = result where++  result = runST stuff  ++  stuff :: forall s. ST s (FreeMod Schur (FieldOfFractions coeff))+  stuff = do+    starr <- newArray (1,nparts) 0 :: ST s (STArray s Int (FieldOfFractions coeff))+    forM_ fixpoints (worker starr)+    arr <- unsafeFreeze starr :: ST s (Array Int (FieldOfFractions coeff))+    let g (j,x) = ( Schur (renormLambdaArr!j) , x ) +        bcs = map g (assocs arr)+    return (fromList bcs)+    +  renormLambdaArr = +    listArray (1,nparts) +      [ posnegPairToPartition (   i,mu) (pos,neg) | (pos,neg) <- posneg ]+        :: Array Int Partition+  complLambdaArr = +    listArray (1,nparts) +      [ posnegPairToPartition ( n-i,mu) (neg,pos) | (pos,neg) <- posneg ]+        :: Array Int Partition++{-+  subs :: Term -> Integer+  subs = evaluate f where +    f (Symbol "alpha" (Just i)) coeff = coeff * q^(i-1)+    q = 1 + fromIntegral n :: Integer+-}++  subs :: Term coeff -> coeff+  subs = evaluate f where +    f (Symbol "alpha" (Just i)) coeff = coeff * fromInteger (subsTable!i)++  subsTable = getSubs n++  worker :: STArray s Int (FieldOfFractions coeff) -> Fixpoint1 -> ST s ()+  worker arr fixpoint = do+    let sol = map subs $ solution fixpoint+        tng = map subs $ tangent  fixpoint+        z = product tng+        chern = elemSymmArray sol+    forM_ [1..nparts] $ \j -> do+      let clambda = complLambdaArr ! j+          y = schurFromChernArray chern clambda+      readArray arr j >>= \x -> writeArray arr j (x + correctTheSign (embed y / embed z))+      return ()++  correctTheSign :: FieldOfFractions coeff -> FieldOfFractions coeff+  correctTheSign = if signCorrection problem < 0 then negate else id+  +  solution :: Fixpoint1 -> [Term coeff]+  solution = map alpha +  +  tangent :: Fixpoint1 -> [Term coeff]+  tangent xs = [ alpha j ^-^ alpha i | i<-xs, j<-ys ] where ys = [1..n] \\ xs+  +  mu     = i+  nparts = length posneg+  fixpoints = choose i [1..n] +  +--------------------------------------------------------------------------------
+ src/Math/ThomPoly/SigmaIJ.hs view
@@ -0,0 +1,230 @@++-- | Calculates the Thom polynomial of @Sigma^{ij}@ with localization +-- and the substitution trick++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, PackageImports #-}+module Math.ThomPoly.SigmaIJ where++--------------------------------------------------------------------------------++import Control.Monad+import Control.Monad.ST+import Data.STRef++import Data.Array.IArray+import Data.Array.Unsafe+import Data.Array.ST++import Data.List+import Data.Ratio+import Data.Proxy++import Debug.Trace+import GHC.IO ( unsafeIOToST )++import System.Mem+import System.IO++import Math.Combinat.Classes+import Math.Combinat.Partitions.Integer+import Math.Combinat.Sets++import Math.FreeModule.Symbol+import Math.FreeModule.SortedList+import Math.FreeModule.PrettyPrint+import Math.FreeModule.PP+-- import Math.FreeModule.Parser++import Math.Algebra.ModP+import Math.Algebra.Schur++import Math.ThomPoly.Subs+import Math.ThomPoly.Shared++--------------------------------------------------------------------------------++instance Problem SigmaIJ where+  calcStats = statsIJ+  solve     = sigmaij+  baseFName (SigmaIJ i j n)  = "sigmaij__i" ++ show i ++ "_j" ++ show j ++ "_n" ++ show n+  +--------------------------------------------------------------------------------+-- * @Sigma^{ij}@++data SigmaIJ = SigmaIJ+  { _i :: !Int     -- ^ the index @i@+  , _j :: !Int     -- ^ the index @j@+  , _n :: !Int     -- ^ the source dimension @n@+  }+ deriving (Eq,Show)++-- | We need @n >= mu@ with this method+smallestIJ :: (Int,Int) -> SigmaIJ+smallestIJ ij@(i,j) = SigmaIJ i j (calcMu ij)++-- | The codimension of @Sigma^{i,j}(n,m)@+codim :: SigmaIJ -> Int -> Int+codim (SigmaIJ i j n) m = calcMu (i,j) * (m-n+i)  - (i-j)*j++-- | There is a sign in the localization formula.+signCorrection :: SigmaIJ -> Int+signCorrection (SigmaIJ i j n) = (-1)^p where+  p = n*mu + i*(j-mu)-j*j +  mu = calcMu (i,j)++-- | computes the (shifted) algebraic multiplicity @mu = i + (j `o` i)@+calcMu :: (Int,Int) -> Int +calcMu (i,j) = i + (j `o` i)++-- | Signed pairs of partitions appearing in the Thom polynomial of @Sigma^{ij}@+listPosNeg :: SigmaIJ -> [(Partition,Partition)]+listPosNeg (SigmaIJ i j n) = list where+  list = partitionPairs mu n i (-j*(i-j))+  mu   = i + (j `o` i)++statsIJ :: SigmaIJ -> Stats+statsIJ prob@(SigmaIJ i j n) = Stats +  { _mu       = calcMu (i,j) +  , _codim0   = codim prob n+  , _maxPairs = length $ listPosNeg prob+  }++--------------------------------------------------------------------------------++-- | A fixed point   +data Fixpoint2 = Fix2 +  { _ii  :: [Int] +  , _jj  :: [Int] +  , _ioj :: [(Int,Int)] -- ioj = jj `o` ii+  , _kk  :: [Int]       -- kk = nn\ii+  , _ss  :: [Int]       -- ioj resze+  , _rr  :: [Int]       -- nn\\ii resze+  }+  deriving Show++-- | dimension of a \"half-symmetric tensor product\"+o :: Int -> Int -> Int  +j `o` i = +  if j<=i +    then div (j*(j+1)) 2 + j*(i-j)+    else error "half-symmetric tensor product [dim]: error"+    +-- | \"half-symmetric tensor product\"+--+-- > length (js `oo` is) == (length js) `o` (length is)+--+oo :: [Int] -> [Int] -> [(Int,Int)]+jj `oo` ii = +  if and [ j `elem` ii | j<-jj ] +    then map (\[x,y]->(x,y)) (choose 2 jj) ++ +         [ (j,j) | j<-jj ] +++         [ (j,i) | j<-jj, i<-ii_minus_jj ] +    else error "half-symmetric tensor product [list]: error"+  where+    ii_minus_jj = ii \\ jj ++--------------------------------------------------------------------------------++sigmaij :: CoeffRing coeff => Proxy coeff -> Batch -> SigmaIJ -> FreeMod Schur (FieldOfFractions coeff)+sigmaij pxy batch problem@(SigmaIJ i j n) = sigmaij' pxy problem (selectBatch batch posneg) where+  posneg = partitionPairs mu n i (-j*(i-j))+  mu     = i + (j `o` i)+    +sigmaij' +  :: forall coeff. CoeffRing coeff +  => Proxy coeff -> SigmaIJ -> [(Partition,Partition)] -> FreeMod Schur (FieldOfFractions coeff)+sigmaij' _ problem@(SigmaIJ i j n) posneg = {- if n<mu then error "n<mu" else -} result where++  result = runST stuff  ++  phi (j,i) = alpha j ^+^ alpha i++  stuff :: forall s. ST s (FreeMod Schur (FieldOfFractions coeff))+  stuff = do+    starr <- newArray (1,nparts) 0 :: ST s (STArray s Int (FieldOfFractions coeff))+    +    forM_ (choose i nn) $ \ii -> do+      let ni = nn \\ ii+          tng1' = [ alpha b ^-^ alpha a | a<-ii, b<-ni ] +          sol1' = [ alpha a | a<-ii] +          tng1 = map subs tng1'+          sol1 = map subs sol1'+      forM_ (choose j ii) $ \jj -> do+        let ij = ii \\ jj    :: [Int]+            ioj = jj `oo` ii :: [(Int,Int)]          +            tng2' = [ alpha b ^-^ alpha a | a<-jj, b<-ij ]+            tng2  = map subs tng2'+        forM_ [0..mu'] $ \k -> do+          forM_ (choose k ioj) $ \ss -> do     -- ss is 'coim'+            forM_ (choose k ni) $ \rr -> do    -- rr is 'im'+              let ker   = ioj \\ ss+                  coker = ni \\ rr+                  tng3' =  [ alpha b ^-^ phi   a | a<-ss  , b<-rr    ]+                        ++ [ phi   a ^-^ phi   b | a<-ss  , b<-ker   ] -- itt van az elojel!+                        ++ [ alpha b ^-^ alpha a | a<-rr  , b<-coker ] +                        ++ [ alpha b ^-^ phi   a | a<-ker , b<-coker ] +                  tng3 = map subs tng3'++              let tng123' = tng1' ++ tng2' ++ tng3'+                  tng123  = tng1  ++ tng2  ++ tng3+                  z = product tng123+                  sol2 = map subs +                       $ [ phi a | a<-ker ] ++ [ alpha b | b<-rr]  ++              when (z==0) $ unsafeIOToST $ do+                putStrLn $ "error: zero denominator!"+                putStrLn $ "substitution table: " ++ show (elems subsTable)+                forM_ (zip tng123 tng123') $ \(a,p) -> do+                  when (a==0) $ putStrLn (pretty p ++ " == 0")+                     +              let sol = sol1 ++ sol2+                  -- chern = elemSymmArray sol+                  segre = completeSymmArray (i*(n-i)+j*(i-j)+mu+(n-i)) sol++              -- cachedSchur <- makeSegreSchurCache segre+                  +              forM_ [1..nparts] $ \j -> do+                let clambda = complLambdaArr ! j+                -- let y = (if odd k then negate else id) (schurFromChernArray chern clambda)+                let y = (if odd k then negate else id) (schurFromSegreArray segre clambda)+                x <- readArray starr j +                x `seq` y `seq` z `seq` writeArray starr j (x + correctTheSign (embed y / embed z))+                return ()+    +    arr <- unsafeFreeze starr :: ST s (Array Int (FieldOfFractions coeff))+    let g (j,x) = ( Schur (renormLambdaArr!j) , x ) +        bcs = map g (assocs arr)+    return (fromList bcs)++  correctTheSign :: FieldOfFractions coeff -> FieldOfFractions coeff+  correctTheSign = if signCorrection problem < 0 then negate else id+    +  nn  = [1..n] +  mu' = j `o` i +  mu  = i + mu'+  nparts = length posneg+    +  renormLambdaArr = +    listArray (1,nparts) +      [ posnegPairToPartition (   i,mu) (pos,neg) | (pos,neg) <- posneg ]+        :: Array Int Partition+  complLambdaArr = +    listArray (1,nparts) +      [ posnegPairToPartition ( n-i,mu) (neg,pos) | (pos,neg) <- posneg ]+        :: Array Int Partition++  subs :: Term coeff -> coeff+  subs = evaluate f where +    f (Symbol "alpha" (Just i)) coeff = coeff * fromInteger (subsTable!i)++  subsTable = getSubsNum n++{-+  subs :: Term -> Integer+  subs = evaluate f where +    f (Symbol "alpha" (Just i)) coeff = coeff * q^(i-1)+    q = 1 + fromIntegral n :: Integer+-}++--------------------------------------------------------------------------------+
+ src/Math/ThomPoly/Subs.hs view
@@ -0,0 +1,91 @@++-- | We need to find (small) integer substitution such+-- that the denominator in our formula never vanishes.+--+-- That is, for @n@ we need to find @n@ integers @a[i]@ such that:+--+-- * @0  /=  a[i] -  a[j]@+--   +-- * @0  /=  a[i] - (a[j] + a[k])@+--+-- * @0  /=  (a[i] + a[j]) - (a[k] + a[l])@+--++{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}+module Math.ThomPoly.Subs where++--------------------------------------------------------------------------------++import Data.Array+import Data.List++import Control.Monad+import System.Random+import System.IO.Unsafe as Unsafe++import Math.Combinat.Sets++import Math.Algebra.ModP++--------------------------------------------------------------------------------+-- * Lazily cached tables of substitutions++getSubs :: Int -> Array Int Integer+getSubs n = theSubsTable!!n++getSubsNum :: Num a => Int -> Array Int a+getSubsNum n = fmap fromInteger (theSubsTable!!n)++getSubsZp :: Int -> Array Int Zp+getSubsZp n = fmap mkZp (theSubsTable!!n)++-- | We cache a substitution table+theSubsTable :: [Array Int Integer]+theSubsTable = [ listArray (1,n) (Unsafe.unsafePerformIO (findSubs n)) | n <- [0..] ]++--------------------------------------------------------------------------------+-- * Find substitutions++-- | Select @k@ elements from a list in all possible orders+choosePerm :: Int -> [a] -> [[a]]+choosePerm n xs = concatMap Data.List.permutations (choose n xs)++-- | Checks if a substitution satisfies the constraints+{-# SPECIALIZE checkSubs :: [Int]     -> Bool #-}+{-# SPECIALIZE checkSubs :: [Integer] -> Bool #-}+checkSubs :: forall a. (Eq a, Num a) => [a] -> Bool+checkSubs input = ok2 && ok3 && ok4 where++  n   = length input+  nn  = [1..n] +  arr = listArray (1,n) input :: Array Int a++  ok2 = and [ 0 /= (arr!i - arr!j)         | [i,j] <- choose 2 nn ] ++  ok3 = and [ 0 /= (arr!i - arr!j - arr!k) | [i,j,k] <- choosePerm 3 nn ] &&+        and [ 0 /= (arr!i - 2 * arr!j    ) | [i,j]   <- choosePerm 2 nn ]++  ok4 = and [ 0 /= (arr!i + arr!j - arr!k - arr!l) | [i,j,k,l] <- choosePerm 4 nn ] && +        and [ 0 /= (arr!i + arr!j - 2 * arr!k    ) | [i,j,k]   <- choosePerm 3 nn ] ++--------------------------------------------------------------------------------++findSubsZp :: Int -> IO [Zp]+findSubsZp = liftM (map mkZp) . findSubs++-- | Find random substitution which satisfies the constraints+findSubs :: Int -> IO [Integer]+findSubs n = go 25 where++  go !bound = tryN 100 bound ++  tryN 0    !bound = go (div (bound*3) 2)+  tryN !cnt !bound = do+    subs <- replicateM n $ randomRIO (-bound,bound)+    case checkSubs subs of+      False -> tryN (cnt-1) bound+      True  -> do+        putStrLn $ "good substitution found! " ++ show subs+        return subs++--------------------------------------------------------------------------------
+ src/cbits/c_det.c view
@@ -0,0 +1,135 @@++#include "c_det.h"++typedef __int128  int128_t;++// -----------------------------------------------------------------------------++// we assume a and b are already mod p+inline int64_t sub_modp( int64_t p , int64_t a , int64_t b )+{+  if (b <= a) +    { return (a - b); }+  else+    { return (a + p - b); }+}++inline int64_t mul_modp( int64_t p0 , int64_t a0 , int64_t b0 )+{+  int128_t p = p0;+  int128_t a = a0;+  int128_t b = b0;+  int128_t c = a*b;+  c = c % p;+  return ((int64_t)c);+}++// -----------------------------------------------------------------------------++int64_t euclid( int64_t p , int64_t x1_ , int64_t x2_ , int64_t u_ , int64_t v_ )+{+  int64_t halfp1 = (p + 1) >> 1;++  int64_t x1 = x1_; +  int64_t x2 = x2_;+  int64_t u  = u_;+  int64_t v  = v_;++  while( (u!=1) && (v!=1) )+  {+    while (!(u & 1))+    { // u even+      u = u >> 1;+      if (x1 & 1) { /* x1 odd */  x1 = (x1 >> 1) + halfp1; } else { x1 = x1 >> 1; }+    }+      +    while (!(v & 1))+    { // v even+      v = v >> 1;+      if (x2 & 1) { /* x2 odd */  x2 = (x2 >> 1) + halfp1; } else { x2 = x2 >> 1; }+    }++    if (u >= v)+    {+      u  = u - v;+      if ( x1 >= x2 ) { x1 = (x1 - x2); } else { x1 = (x1 + p - x2); } +    }+    else +    {+      v  = v - u;+      if ( x2 >= x1)  { x2 = (x2 - x1); } else { x2 = (x2 + p - x1); } +    }++  }++  if (u==1) { return x1; }+  if (v==1) { return x2; }+  return 0;                       // shouldn't happen+}++// -----------------------------------------------------------------------------++inline int64_t div_modp( int64_t p , int64_t a , int64_t b )+{+  // return mul_modp( p , a , inv_modp( p , b ) );+  return euclid( p , a , 0 , b , p );}++// mod p inverse using the binary Euclidean algorithm +int64_t inv_modp( int64_t p , int64_t a )+{+  return euclid( p , 1 , 0 , a , p );+}++// -----------------------------------------------------------------------------++// determinant mod p (64 bit), using Gauss elimination+int64_t det_modp(int64_t p, int n, int64_t *mat)+{+  // safety first+  for (int i=0;i<n*n;i++) { if ((mat[i] >= p) || (mat[i]<0)) { mat[i] = mat[i] % p; } }++  int negative = 0;++  for (int i=0;i<n-1;i++)+  {+    int64_t *row = mat + i*n;+   +    // find pivot element+    int j; +    for (j=i;j<n;j++) { if (row[j] != 0) break; }+    if ( (j >= n) || (row[j] == 0) ) { return 0; }+     +    if (j > i)+    { // exchange columns+      int64_t *q = row;      +      for (int k=i;k<n;k++)+      { +        int64_t x;+        x    = q[i];+        q[i] = q[j];+        q[j] = x;+        q += n;     +      }+      negative = negative ^ 1;     // track the sign changes+    }++    // zero out the i-th column+    int64_t *q = row + n;+    for (int k=i+1;k<n;k++)+    { +      int64_t m = div_modp( p , q[i] , row[i] );+      q[i] = 0;  +      for (int l=i+1;l<n;l++) +      { +        q[l] = sub_modp( p , q[l] , mul_modp( p , m , row[l] ) );+      }+      q += n;+    }+  }++  int64_t det = mat[0];+  for (int i=1;i<n;i++) { det = mul_modp( p , det , mat[i*(n+1)] ); }++  if ((negative) && (det!=0)) { return (p-det); } else { return det; }+}+
+ src/cbits/c_det.h view
@@ -0,0 +1,20 @@++//------------------------------------------------------------------------------++#ifndef C_DET_H_INCLUDED+#define C_DET_H_INCLUDED++#include <stdint.h>++//------------------------------------------------------------------------------++int64_t inv_modp( int64_t p , int64_t a );+int64_t div_modp( int64_t p , int64_t a , int64_t b );+int64_t det_modp( int64_t p , int n , int64_t *mat );++int64_t euclid( int64_t p , int64_t x1_ , int64_t x2_ , int64_t u_ , int64_t v_ );++//------------------------------------------------------------------------------++#endif // C_DET_H_INCLUDED+
+ src/sigmaij.hs view
@@ -0,0 +1,377 @@++-- | Calculates Thom polynomial of Sigma^{ij} with localization +-- and the substituion trick+--+-- Some example usages:+--+-- > sigma-ij -h                           # help+-- > sigma-ij -i3 -j2 -n7                  # compute @Tp(Sigma^{3,2}(7))@ in the default ring+-- > sigma-ij -i3 -j2 -n7  -rZp            # compute @Tp(Sigma^{3,2}(7))@ in the hard-coded prime field+-- > sigma-ij -i3 -j2 -n10 -rZp -b3 -B10   # compute the 3rd part of @Tp(Sigma^{3,2}(10))@ divided into 10 pieces+--+-- The task can be parallezied using the @-B@ and @-b@ options+--++{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, PackageImports, PatternGuards #-}+module Main where++--------------------------------------------------------------------------------++import Data.Char+import Data.List+import Data.Ratio+import Data.Monoid++import Control.Monad+import Control.Applicative+import Control.Concurrent+import Control.Concurrent.MVar++import System.Environment+import System.Mem+import System.IO+import System.Exit++import "time" Data.Time.Clock.POSIX++import Math.Combinat.Numbers.Primes++import Math.FreeModule.Symbol+import Math.FreeModule.SortedList+import Math.FreeModule.PrettyPrint+import Math.FreeModule.PP+-- import FreeModule.Parser++import Math.Algebra.ModP+import Math.Algebra.Schur++import Math.ThomPoly.Shared+import Math.ThomPoly.SigmaI   as SigmaI+import Math.ThomPoly.SigmaIJ  as SigmaIJ+import Math.ThomPoly.Formulae as Formulae++import Options.Applicative+    +--------------------------------------------------------------------------------++data Config = Config+  { _problem   :: !AnyProblem  +  , _tgtDim    :: !(Maybe Int)+  , _ring      :: !Ring +  , _outFile   :: !(Maybe FilePath)+  , _batch     :: !(Maybe Batch)+  , _printStat :: !Bool+  , _dry       :: !Bool+  , _timeout   :: !(Maybe Int)+  }+  deriving Show++--------------------------------------------------------------------------------+  +run :: Config -> IO ()+run config = do+  -- print config+  void $ mbTimeout (_timeout config) $ do++    let problem = _problem config+        batch   = maybe defaultBatch id (_batch config)+        ring    = selectRing (_ring config)++    when (_printStat config) $ do+      print $ case problem of+        PI  si  -> calcStats si +        PIJ sij -> calcStats sij++    let fname = case _outFile config of+          Just fname -> fname+          Nothing    -> case problem of+            PI  si  -> fullFName batch si +            PIJ sij -> fullFName batch sij        ++    let answer = case problem of+          PI  si  -> solveAny ring batch si +          PIJ sij -> solveAny ring batch sij++    let text   = pretty answer++    answer `seq` do+      unless (_dry config) $ writeFile fname text++--------------------------------------------------------------------------------+-- * configuration++data AnyProblem +  = PI  !SigmaI+  | PIJ !SigmaIJ+  | PI1 !SigmaI1            -- ^ we have an explicit formula for @Sigma^{i,1}@+  deriving Show++-- | Coefficient ring +data Ring+  = Integers+  | Rationals+  | HardCodedZp             -- ^ temporary+  | PrimeField !Integer+  | SpecPrime  !Int         -- ^ special primes just below @2^k@ for @k=7,15,31,63@+  deriving Show++selectRing :: Ring -> AnyRing+selectRing r = case r of+  Integers    -> ringZZ+  Rationals   -> ringQQ+  HardCodedZp -> ringZp++-- | Primes close to the bounds of (signed) machine words.+specPrimes :: [(Int,Integer)]+specPrimes = +  [ ( 7   , 2^7   -  1 )+  , ( 15  , 2^15  - 19 )+  , ( 31  , 2^31  -  1 )+  , ( 63  , 2^63  - 25 )+--  , ( 127 , 2^127 -  1 )+--  , ( 255 , 2^255 - 19 )+  ]++--------------------------------------------------------------------------------++maybeRead :: Read a => String -> Maybe a+maybeRead s = case reads s of +  [(x,"")] -> Just x+  _        -> Nothing++parseRing :: String -> Either String Ring+parseRing str0 +  | str `elem` ["zz","integer" ,"integers"  ]   =  Right Integers+  | str `elem` ["qq","rational","rationals" ]   =  Right Rationals+  | str `elem` ["zp","primefield" ]             =  Right HardCodedZp+--  | take 2+  where+    str = map toLower str0 ++--------------------------------------------------------------------------------++class Validate a where+  isValid :: a -> Maybe String+  +instance Validate Batch where+  isValid (Batch a b) +    | b < 1           =  Just "the number of batches B should be at least 1"+    | a < 1 || a > b  =  Just "batch index b should be between 1 and B"+    | otherwise       =  Nothing++instance Validate Ring where+  isValid r = case r of+    PrimeField p -> if isProbablyPrime p +      then Nothing +      else Just "order of the finite field should be a prime"+    SpecPrime  q -> case lookup q specPrimes of +      Nothing     -> Just "unimplemented special prime field (BITS should be one of 7, 15, 31 or 63)"+      _           -> Nothing+    _ -> Nothing       ++instance Validate AnyProblem where++  isValid problem = case problem of+    PI (SigmaI i n)+      | i < 1     -> Just "the index I should be at least 1"+      | n < 1     -> Just "the source dimension N should be at least 1"+      | otherwise -> Nothing+    PIJ (SigmaIJ i j n)+      | i < 1          -> Just "the index I should be at least 1"+      | j < 1 || j > i -> Just "the index J should be between 1 and I"+      | n < 1          -> Just "the source dimension N should be at least 1"+      | otherwise      -> Nothing+    PI1 (SigmaI1 i n)+      | i < 1     -> Just "the index I should be at least 1"+      | n < 1     -> Just "the source dimension N should be at least 1"+      | otherwise -> Nothing++--------------------------------------------------------------------------------+-- * option parsing++configOpt :: Parser Config+configOpt = Config +  <$> problemOpt +  <*> mOpt +  <*> ringNameOpt+  <*> outOpt +  <*> batchOpt +  <*> statFlag +  <*> dryFlag +  <*> timeoutOpt++problemOpt :: Parser AnyProblem+problemOpt = f <$> iOpt <*> jOpt <*> nOpt where+  f i mbj n = case mbj of+    Nothing -> PI $ SigmaI i n +    Just j  -> case j of +      0 -> PI  $ SigmaI  i   n+      _ -> PIJ $ SigmaIJ i j n++batchOpt :: Parser (Maybe Batch)+batchOpt = f <$> whichBatchOpt <*> nbatchOpt where+  f a b +    | a >= 1 && a <= b  =  if b > 1 then Just (Batch a b) else Nothing+    | otherwise         =  error "the batch index should be between 1 and B"++ringNameOpt :: Parser Ring+ringNameOpt = option (eitherReader parseRing)+  (  long    "ring"+  <> short   'r'+  <> metavar "R"+  <> value   Rationals+  <> help    "The coefficient ring (or field) R we compute in, for example a prime field" +  <> completeWith [ "Integers" , "Rationals" , "ZZ" , "QQ"+                  , "PrimeField" , "Zp"+--                  , "Zp7bit" , "Zp15bit" , "Zp31bit" , "Zp63bit"+                  ]+  <> showDefault+  )++timeoutOpt :: Parser (Maybe Int)+timeoutOpt = option (Just <$> auto)+  (  long    "timeout"+  <> short   't'+  <> metavar "TIMEOUT"+  <> value   Nothing+  <> help    "Timeout (specified in minutes)" +  )++primeOpt :: Parser Integer+primeOpt = option auto+  (  long    "prime"+  <> short   'p'+  <> metavar "P"+  <> value   1000000007+  <> help    "The order of the prime field" +  )++bitsOpt :: Parser Int+bitsOpt = option auto+  (  long    "bits"+  <> short   'q'+  <> metavar "BITS"+  <> value   63+  <> help    "Number of bits in the order of a special prime fields" +  )++statFlag :: Parser Bool+statFlag = switch+  (  long  "stats"+  <> short 's'+  <> help  "print \"statistics\" (codimension, algebraic multiplicity)" +  )++dryFlag :: Parser Bool+dryFlag = switch+  (  long  "dry"+  <> help  "do not write the result into a file" +  )++outOpt :: Parser (Maybe FilePath)+outOpt = option (Just <$> str)+  (  long    "output"+  <> short   'o'+  <> metavar "FILE"+  <> value   Nothing+  <> help    "Write output to FILE (use --dry to skip)" +  )++nbatchOpt :: Parser Int+nbatchOpt = option auto+  (  long    "nbatches"+  <> short   'B'+  <> metavar "B"+  <> value   1+  <> help    "number of batches" +  )++whichBatchOpt :: Parser Int+whichBatchOpt = option auto+  (  long "batch"+  <> short   'b'+  <> metavar "b"+  <> value   1+  <> help    "which batch to run (from 1 to B)" +  )++iOpt :: Parser Int+iOpt = option auto+  (  short   'i'+  <> metavar "I"+  <> help    "first Thom-Boardman index (I)"+  <> noArgError (ErrorMsg "specifying I is mandatory")+  )++jOpt :: Parser (Maybe Int)+jOpt = option (Just <$> auto)+  (  short   'j'+  <> metavar "J"+  <> value   Nothing+  <> help    "second Thom-Boardman index (J)"+  )++nOpt :: Parser Int+nOpt = option auto+  (  short   'n'+  <> metavar "N"+  <> help    "source dimension (N)"+  )++mOpt :: Parser (Maybe Int)+mOpt = option (Just <$> auto)+  (  short   'm'+  <> metavar "N"+  <> value   Nothing+  <> help    "target dimension (M, optional)"+  <> hidden+  )++--------------------------------------------------------------------------------++main :: IO ()+main = execParser opts >>= run where+  opts = info (helper <*> configOpt)+    (  fullDesc+    <> progDesc shortDesc+    <> header   longDesc+    )+  shortDesc = "Thom polynomials of second order Thom-Boardman singularities"+  longDesc  = "A program computing Thom polynomials of second order Thom-Boardman singularities"++--------------------------------------------------------------------------------+-- * timeout++-- | argument: number of minutes+mbTimeout :: Maybe Int -> IO a -> IO (Maybe a)+mbTimeout mb action = +  case mb of +    Nothing      -> Just <$> action+    Just minutes -> do+      mv <- newEmptyMVar +      t0 <- getPOSIXTime+      threadid <- forkIO $ do+        y <- action +        putMVar mv $! y+      wait mv t0 minutes threadid++  where      +    wait mv t0 minutes threadid = do+      let seconds = minutes * 60+      let loop = do+            threadDelay 1000000     -- wait 1 sec+            mb <- tryTakeMVar mv+            case mb of+              Just y  -> return $ Just y+              Nothing -> do+                t <- getPOSIXTime+                if t - t0 < fromIntegral seconds+                  then loop+                  else do+                    putStrLn $ "timeout after " ++ show minutes ++ " minutes"+                    killThread threadid+                    return Nothing+      loop+      +--------------------------------------------------------------------------------+