diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. + diff --git a/README.txt b/README.txt new file mode 100644 --- /dev/null +++ b/README.txt @@ -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: . + + +USAGE: +====== + +sigma-ij -h help +sigma-ij -i3 -j1 -n7 compute $Tp(Sigma^{3,1}(7))$ +sigma-oj -i3 -j1 -n7 -r compute with coefficients in the given ring +sigma-oj -i3 -j1 -n7 -B -b 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 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 diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/sigma-ij.cabal b/sigma-ij.cabal new file mode 100644 --- /dev/null +++ b/sigma-ij.cabal @@ -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 . +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 + diff --git a/src/Math/Algebra/Determinant.hs b/src/Math/Algebra/Determinant.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Determinant.hs @@ -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 + + diff --git a/src/Math/Algebra/ModP.hs b/src/Math/Algebra/ModP.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/ModP.hs @@ -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' + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Algebra/Schur.hs b/src/Math/Algebra/Schur.hs new file mode 100644 --- /dev/null +++ b/src/Math/Algebra/Schur.hs @@ -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 +-} +-------------------------------------------------------------------------------- diff --git a/src/Math/FreeModule/Class.hs b/src/Math/FreeModule/Class.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/Class.hs @@ -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) + +-------------------------------------------------------------------------------- + + diff --git a/src/Math/FreeModule/Helper.hs b/src/Math/FreeModule/Helper.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/Helper.hs @@ -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)) + +-------------------------------------------------------------------------------- diff --git a/src/Math/FreeModule/PP.hs b/src/Math/FreeModule/PP.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/PP.hs @@ -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) ++ ")" + +-------------------------------------------------------------------------------- diff --git a/src/Math/FreeModule/Parser.hs b/src/Math/FreeModule/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/Parser.hs @@ -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 + +-------------------------------------------------------------------------------- + diff --git a/src/Math/FreeModule/PrettyPrint.hs b/src/Math/FreeModule/PrettyPrint.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/PrettyPrint.hs @@ -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 + +-------------------------------------------------------------------------------- + diff --git a/src/Math/FreeModule/SortedList.hs b/src/Math/FreeModule/SortedList.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/SortedList.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/FreeModule/Symbol.hs b/src/Math/FreeModule/Symbol.hs new file mode 100644 --- /dev/null +++ b/src/Math/FreeModule/Symbol.hs @@ -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" + ] + +-------------------------------------------------------------------------------- diff --git a/src/Math/ThomPoly/Formulae.hs b/src/Math/ThomPoly/Formulae.hs new file mode 100644 --- /dev/null +++ b/src/Math/ThomPoly/Formulae.hs @@ -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) + +-------------------------------------------------------------------------------- diff --git a/src/Math/ThomPoly/Shared.hs b/src/Math/ThomPoly/Shared.hs new file mode 100644 --- /dev/null +++ b/src/Math/ThomPoly/Shared.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/Math/ThomPoly/SigmaI.hs b/src/Math/ThomPoly/SigmaI.hs new file mode 100644 --- /dev/null +++ b/src/Math/ThomPoly/SigmaI.hs @@ -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] + +-------------------------------------------------------------------------------- diff --git a/src/Math/ThomPoly/SigmaIJ.hs b/src/Math/ThomPoly/SigmaIJ.hs new file mode 100644 --- /dev/null +++ b/src/Math/ThomPoly/SigmaIJ.hs @@ -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 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 +-} + +-------------------------------------------------------------------------------- + diff --git a/src/Math/ThomPoly/Subs.hs b/src/Math/ThomPoly/Subs.hs new file mode 100644 --- /dev/null +++ b/src/Math/ThomPoly/Subs.hs @@ -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 + +-------------------------------------------------------------------------------- diff --git a/src/cbits/c_det.c b/src/cbits/c_det.c new file mode 100644 --- /dev/null +++ b/src/cbits/c_det.c @@ -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= p) || (mat[i]<0)) { mat[i] = mat[i] % p; } } + + int negative = 0; + + for (int i=0;i= n) || (row[j] == 0) ) { return 0; } + + if (j > i) + { // exchange columns + int64_t *q = row; + for (int k=i;k + +//------------------------------------------------------------------------------ + +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 + diff --git a/src/sigmaij.hs b/src/sigmaij.hs new file mode 100644 --- /dev/null +++ b/src/sigmaij.hs @@ -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 + +-------------------------------------------------------------------------------- +