module Algebra.Algorithms.ZeroDim
(
solveM, solve', solveViaCompanion, solveLinear,
radical, isRadical,
fglm, fglmMap,
solveWith, univPoly,
reduction,
matrixRep, subspMatrix,
vectorRep
) where
import Algebra.Algorithms.FGLM
import Algebra.Algorithms.Groebner
import Algebra.Instances ()
import Algebra.Internal hiding (OLt)
import qualified Algebra.Matrix as AM
import Algebra.Prelude.Core
import Algebra.Ring.Polynomial.Quotient
import Control.Lens hiding ((:<))
import Control.Monad.Loops (whileM_)
import Control.Monad.Random hiding (next)
import Control.Monad.Reader (runReaderT)
import Control.Monad.ST (runST)
import Data.Complex (Complex (..), magnitude)
import Data.Convertible (Convertible, convert)
import qualified Data.Matrix as M
import Data.Maybe (fromJust)
import Data.Ord (comparing)
import Data.Reflection (Reifies)
import qualified Data.Sized.Builtin as SV
import Data.STRef.Strict (newSTRef)
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import qualified Numeric.Algebra as NA
import qualified Numeric.LinearAlgebra as LA
import qualified Prelude as P
import Proof.Propositional (IsTrue (Witness))
solveM :: forall m r ord n.
(Normed r, Ord r, MonadRandom m, Field r, CoeffRing r, KnownNat n,
IsMonomialOrder n ord, Convertible r Double,
(0 :< n) ~ 'True)
=> Ideal (OrderedPolynomial r ord n)
-> m [Sized n (Complex Double)]
solveM ideal =
withKnownNat (sSucc (sing :: SNat n)) $
reifyQuotient (radical ideal) $ \pxy ->
case standardMonomials' pxy of
Just bs -> step 10 (length bs)
Nothing -> error "Given ideal is not zero-dimensional"
where
step bd len = do
coeffs <-
replicateM (sNatToInt (sSucc (sing :: SNat n))) $ getRandomR (bd, bd)
let vars = one : SV.toList allVars
f = sum $ zipWith (.*.) (map (NA.fromInteger :: Integer -> r) coeffs) vars
case solveWith f ideal of
Nothing -> step (bd*2) len
Just sols -> return sols
solveWith :: forall r n ord. (DecidableZero r, Normed r, Ord r, Field r, CoeffRing r,
(0 :< n) ~ 'True, IsMonomialOrder n ord,
KnownNat n, Convertible r Double)
=> OrderedPolynomial r ord n
-> Ideal (OrderedPolynomial r ord n)
-> Maybe [Sized n (Complex Double)]
solveWith f0 i0 =
withKnownNat (sSucc (sing :: SNat n)) $
reifyQuotient (radical i0) $ \pxy ->
let ideal = gBasis' pxy
Just base = map (leadingMonomial . quotRepr) <$> standardMonomials' pxy
in case elemIndex one base of
Nothing -> Just []
Just cind ->
let f = modIdeal' pxy f0
vars = sortBy (flip $ comparing snd) $
map (\on -> (on, leadingMonomial $ var on `asTypeOf` f0)) $
enumOrdinal $ sArity' f0
inds = flip map vars $ second $ \b ->
case findIndex (==b) base of
Just ind -> Right ind
Nothing ->
let Just g = find ((==b) . leadingMonomial) ideal
r = leadingCoeff g
answer = mapCoeff toComplex $ injectCoeff (recip r) * (toPolynomial (leadingTerm g) g)
in Left answer
mf = AM.fromLists $ map (map toComplex) $ matrixRep f
(_, evecs) = LA.eig $ LA.tr mf
calc vec =
let c = vec LA.! cind
phi (idx, Right nth) acc = acc & ix idx .~ (vec LA.! nth) P./ c
phi (idx, Left g) acc = acc & ix idx .~ substWith (*) acc g
in if c == 0
then Nothing
else Just $ foldr ( phi) (SV.replicate (sArity' f0) (error "indec!")) inds
in sequence $ map calc $ LA.toColumns evecs
solve' :: forall r n ord.
(Field r, CoeffRing r, KnownNat n, (0 :< n) ~ 'True,
IsMonomialOrder n ord, Convertible r Double)
=> Double
-> Ideal (OrderedPolynomial r ord n)
-> [Sized n (Complex Double)]
solve' err ideal =
reifyQuotient ideal $ \ii ->
if gBasis' ii == [one]
then []
else
let vs = map (nub . LA.toList . LA.eigenvalues . AM.fromLists . map (map toComplex). matrixRep . modIdeal' ii) $
SV.toList allVars
mul p q = toComplex p * q
in [ xs
| xs0 <- sequence vs
, let xs = SV.unsafeFromList' xs0
, all ((< err) . magnitude . substWith mul xs) $ generators ideal
]
where
_ = Witness :: IsTrue (0 :< n)
subspMatrix :: (Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord)
=> Ordinal n -> Ideal (OrderedPolynomial r ord n) -> M.Matrix r
subspMatrix on ideal =
let poly = univPoly on ideal
v = var on `asTypeOf` head (generators ideal)
dim = fromIntegral $ totalDegree' poly
cfs = [negate $ coeff (leadingMonomial $ pow v (j :: Natural)) poly | j <- [0..fromIntegral (dim 1)]]
in (M.fromLists [replicate (dim 1) zero]
M.<->
fmap unwrapAlgebra (M.identity (dim 1))) M.<|> M.colVector (V.fromList cfs)
solveViaCompanion :: forall r ord n.
(Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord, Convertible r Double)
=> Double
-> Ideal (OrderedPolynomial r ord n)
-> [Sized n (Complex Double)]
solveViaCompanion err ideal =
if calcGroebnerBasis ideal == [one]
then []
else
let vs = map (nub . LA.toList . LA.eigenvalues . LA.fromLists . matToLists . fmap toComplex . flip subspMatrix ideal) $
enumOrdinal (sing :: SNat n)
mul p q = toComplex p * q
in [ xs
| xs0 <- sequence vs
, let xs = SV.unsafeFromList' xs0
, all ((< err) . magnitude . substWith mul xs) $ generators ideal
]
matToLists :: M.Matrix a -> [[a]]
matToLists mat = [ V.toList $ M.getRow i mat | i <- [1.. M.nrows mat] ]
matrixRep :: (DecidableZero t, Eq t, Field t, KnownNat n, IsMonomialOrder n order,
Reifies ideal (QIdeal (OrderedPolynomial t order n)))
=> Quotient (OrderedPolynomial t order n) ideal -> [[t]]
matrixRep f =
case standardMonomials of
Just [] -> []
Just bases ->
let anss = map (quotRepr . (f *)) bases
in transpose $ map (\a -> map (flip coeff a . leadingMonomial . quotRepr) bases) anss
Nothing -> error "Not finite dimension"
toComplex :: Convertible a Double => a -> Complex Double
toComplex a = convert a :+ 0
reduction :: (CoeffRing r, KnownNat n, IsMonomialOrder n ord, Field r)
=> Ordinal n -> OrderedPolynomial r ord n -> OrderedPolynomial r ord n
reduction on f =
let df = diff on f
in snd $ head $ f `divPolynomial` calcGroebnerBasis (toIdeal [f, df])
univPoly :: forall r ord n. (Ord r, Field r, CoeffRing r, KnownNat n, IsMonomialOrder n ord)
=> Ordinal n
-> Ideal (OrderedPolynomial r ord n)
-> OrderedPolynomial r ord n
univPoly nth ideal =
reifyQuotient ideal $ \pxy ->
if gBasis' pxy == [one]
then one
else let x = var nth
p0 : pows = [fmap WrapAlgebra $ vectorRep $ modIdeal' pxy (pow x i)
| i <- [0:: Natural ..]
| _ <- zero : fromJust (standardMonomials' pxy) ]
step m ~(p : ps) =
case solveLinear m p of
Nothing -> step (m M.<|> M.colVector p) ps
Just ans ->
let cur = fromIntegral $ V.length ans :: Natural
in
pow x cur sum (zipWith (.*.) (fmap unwrapAlgebra $ V.toList ans)
[pow x i | i <- [0 :: Natural .. cur P.- 1]])
in step (M.colVector p0) pows
solveLinear :: (Ord r, P.Fractional r)
=> M.Matrix r
-> V.Vector r
-> Maybe (V.Vector r)
solveLinear mat vec =
if ( uRank u) < uRank u' || M.diagProd u == 0 || uRank u < M.ncols mat
then Nothing
else let ans = M.getCol 1 $ p P.* M.colVector vec
lsol = solveL ans
cfs = M.getCol 1 $ q P.* M.colVector ( solveU lsol)
in Just cfs
where
Just (u, l, p, q, _, _) = M.luDecomp' mat
Just (u', _,_, _, _, _) = M.luDecomp' (mat M.<|> M.colVector vec)
uRank = V.foldr (\a acc -> if a /= 0 then acc + 1 else acc) (0 :: Int) . M.getDiag
solveL v = V.create $ do
let stop = min (M.ncols l) (M.nrows l)
mv <- MV.replicate (M.ncols l) 0
forM_ [0..stop 1] $ \i -> do
MV.write mv i $ v V.! i
forM_ [0,1..min (i1) (M.ncols l 1)] $ \j -> do
a <- MV.read mv i
b <- MV.read mv j
MV.write mv i $ a P.- (l M.! (i + 1, j + 1)) P.* b
return mv
solveU v = V.create $ do
let stop = min (M.ncols u) (M.nrows u)
mv <- MV.replicate (M.ncols u) 0
forM_ [stop 1, stop 2 .. 0] $ \ i -> do
MV.write mv i $ v V.! i
forM_ [i+1,i+2..M.ncols u1] $ \j -> do
a <- MV.read mv i
b <- MV.read mv j
MV.write mv i $ a P.- (u M.! (i+1, j+1)) P.* b
a0 <- MV.read mv i
MV.write mv i $ a0 P./ (u M.! (i+1, i+1))
return mv
radical :: forall r ord n . (Ord r, CoeffRing r,
KnownNat n, Field r, IsMonomialOrder n ord)
=> Ideal (OrderedPolynomial r ord n) -> Ideal (OrderedPolynomial r ord n)
radical ideal =
let gens = map (\on -> reduction on $ univPoly on ideal) $ enumOrdinal (sing :: SNat n)
in toIdeal $ calcGroebnerBasis $ toIdeal $ generators ideal ++ gens
isRadical :: forall r ord n. (Ord r, CoeffRing r, KnownNat n,
(0 :< n) ~ 'True,
Field r, IsMonomialOrder n ord)
=> Ideal (OrderedPolynomial r ord n) -> Bool
isRadical ideal =
let gens = map (\on -> reduction on $ univPoly on ideal) $
enumOrdinal (sing :: SNat n)
in all (`isIdealMember` ideal) gens
where
_ = Witness :: IsTrue (0 :< n)
fglm :: (Ord r, KnownNat n, Field r,
IsMonomialOrder n ord, (0 :< n) ~ 'True)
=> Ideal (OrderedPolynomial r ord n)
-> ([OrderedPolynomial r Lex n], [OrderedPolynomial r Lex n])
fglm ideal = reifyQuotient ideal $ \pxy ->
fglmMap (\f -> vectorRep $ modIdeal' pxy f)
fglmMap :: forall k ord n. (Ord k, Field k, (0 :< n) ~ 'True,
IsMonomialOrder n ord, CoeffRing k, KnownNat n)
=> (OrderedPolynomial k ord n -> V.Vector k)
-> ( [OrderedPolynomial k Lex n]
, [OrderedPolynomial k Lex n]
)
fglmMap l = runST $ do
env <- FGLMEnv l <$> newSTRef [] <*> newSTRef [] <*> newSTRef Nothing <*> newSTRef one
flip runReaderT env $ do
mainLoop
whileM_ toContinue $ nextMonomial >> mainLoop
(,) <$> look gLex <*> (map (changeOrder Lex) <$> look bLex)
mainLoop :: (DecidableZero r, Ord r, KnownNat n, Field r, IsMonomialOrder n o)
=> Machine s r o n ()
mainLoop = do
m <- look monomial
let f = toPolynomial (one, changeMonomialOrderProxy Proxy m)
lx <- image f
bs <- mapM image =<< look bLex
let mat = foldr (M.<|>) (M.fromList 0 0 []) $ map (M.colVector . fmap WrapAlgebra) bs
cond | null bs = if V.all (== zero) lx
then Just $ V.replicate (length bs) 0
else Nothing
| otherwise = solveLinear mat (fmap WrapAlgebra lx)
case cond of
Nothing -> do
proced .== Nothing
bLex %== (f : )
Just cs -> do
bps <- look bLex
let g = changeOrder Lex $ f sum (zipWith (.*.) (V.toList $ fmap unwrapAlgebra cs) bps)
proced .== Just (changeOrder Lex f)
gLex %== (g :)
toContinue :: forall s r o n.
((0 :< n) ~ 'True, Ord r,
KnownNat n, Field r)
=> Machine s r o n Bool
toContinue = do
mans <- look proced
case mans of
Nothing -> return True
Just g -> do
let xLast = P.maximum allVars `asTypeOf` g
return $ not $ leadingMonomial g `isPowerOf` leadingMonomial xLast
where
_ = Witness :: IsTrue (0 :< n)
nextMonomial :: forall s r ord n.
(CoeffRing r, KnownNat n) => Machine s r ord n ()
nextMonomial = do
m <- look monomial
gs <- map leadingMonomial <$> look gLex
let next = fst $ maximumBy (comparing snd) $
[ (OrderedMonomial monom, fromEnum od)
| od <- enumOrdinal (sing :: SNat n)
, let monom = beta (getMonomial m) od
, all (not . (`divs` OrderedMonomial monom)) gs
]
monomial .== next
beta :: Monomial n -> Ordinal n -> Monomial n
beta xs o@(OLt k) =
let n = sizedLength xs
in withRefl (lneqSuccLeq k n) $
withRefl (plusComm (n %:- sSucc k) (sSucc k)) $
withRefl (minusPlus n (sSucc k) Witness) $
(SV.take (sSucc k) $ xs & ix o +~ 1) SV.++ SV.replicate (n %:- sSucc k) 0
beta _ _ = error "beta: Bug in ghc!"