msakai / toysolver / 797

Committed 09 May 2026 11:18AM UTC coverage: 71.463% (-0.4%) from 71.864%

Build # 797

Build Type

push

github

Committed by

web-flow

Commit Message

Merge 2564c0a1a into 36153c78c

Coverage Stats

11046 of 15457 relevant lines covered (71.46%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

83.69

/src/ToySolver/Arith/OmegaTest/Base.hs

{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.OmegaTest.Base
-- Copyright   :  (c) Masahiro Sakai 2011
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- (incomplete) implementation of Omega Test
--
-- References:
--
-- * William Pugh. The Omega test: a fast and practical integer
--   programming algorithm for dependence analysis. In Proceedings of
--   the 1991 ACM/IEEE conference on Supercomputing (1991), pp. 4-13.
--
-- * <http://users.cecs.anu.edu.au/~michaeln/pubs/arithmetic-dps.pdf>
--
-- See also:
--
-- * <http://hackage.haskell.org/package/Omega>
--
-----------------------------------------------------------------------------
module ToySolver.Arith.OmegaTest.Base
    ( Model
    , solve
    , solveQFLIRAConj
    , Options (..)
    , checkRealNoCheck
    , checkRealByFM

    -- * Exported just for testing
    , zmod
    ) where

import Control.Exception (assert)
import Control.Monad
import Data.Default.Class
import Data.List (foldl', foldl1', minimumBy)
import Data.Maybe
import Data.Ord
import Data.Ratio
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.VectorSpace

import ToySolver.Data.OrdRel
import ToySolver.Data.Boolean
import ToySolver.Data.DNF
import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar
import ToySolver.Internal.Util (combineMaybe)
import qualified ToySolver.Arith.FourierMotzkin as FM

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

-- | Options for solving.
--
-- The default option can be obtained by 'def'.
data Options
  = Options
  { optCheckReal :: VarSet -> [LA.Atom Rational] -> Bool
  -- ^ optCheckReal is used for checking whether real shadow is satisfiable.
  --
  -- * If it returns @True@, the real shadow may or may not be satisfiable.
  --
  -- * If it returns @False@, the real shadow must be unsatisfiable
  }

instance Default Options where
  def =
    Options
    { optCheckReal =
        -- checkRealNoCheck
        checkRealByFM
    }

checkRealNoCheck :: VarSet -> [LA.Atom Rational] -> Bool
checkRealNoCheck _ _ = True

checkRealByFM :: VarSet -> [LA.Atom Rational] -> Bool
checkRealByFM vs as = isJust $ FM.solve vs as

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

type ExprZ = LA.Expr Integer

-- | Atomic constraint
data Constr
  = IsNonneg ExprZ
  -- ^ e ≥ 0
  | IsZero ExprZ
  -- ^ e = 0
  deriving (Show, Eq, Ord)

instance Variables Constr where
  vars (IsNonneg t) = vars t
  vars (IsZero t) = vars t

-- 制約集合の単純化
-- It returns Nothing when a inconsistency is detected.
simplify :: [Constr] -> Maybe [Constr]
simplify = fmap concat . mapM f
  where
    f :: Constr -> Maybe [Constr]
    f constr@(IsNonneg e) =
      case LA.asConst e of
        Just x -> guard (x >= 0) >> return []
        Nothing -> return [constr]
    f constr@(IsZero e) =
      case LA.asConst e of
        Just x -> guard (x == 0) >> return []
        Nothing -> return [constr]

leZ, ltZ, geZ, gtZ, eqZ :: ExprZ -> ExprZ -> Constr
-- Note that constants may be floored by division
leZ e1 e2 = IsNonneg (LA.mapCoeff (`div` d) e)
  where
    e = e2 ^-^ e1
    d = abs $ gcd' [c | (c,v) <- LA.terms e, v /= LA.unitVar]
ltZ e1 e2 = (e1 ^+^ LA.constant 1) `leZ` e2
geZ = flip leZ
gtZ = flip ltZ
eqZ e1 e2 =
  case isZero (e1 ^-^ e2) of
    Just c -> c
    Nothing -> IsZero (LA.constant 1) -- unsatisfiable

isZero :: ExprZ -> Maybe Constr
isZero e
  = if LA.coeff LA.unitVar e `mod` d == 0
    then Just $ IsZero $ LA.mapCoeff (`div` d) e
    else Nothing
  where
    d = abs $ gcd' [c | (c,v) <- LA.terms e, v /= LA.unitVar]

applySubst1Constr :: Var -> ExprZ -> Constr -> Constr
applySubst1Constr v e (IsNonneg e2) = LA.applySubst1 v e e2 `geZ` LA.constant 0
applySubst1Constr v e (IsZero e2)   = LA.applySubst1 v e e2 `eqZ` LA.constant 0

instance Eval (Model Integer) Constr Bool where
  eval m (IsNonneg t) = LA.eval m t >= 0
  eval m (IsZero t)   = LA.eval m t == 0

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

-- | (t,c) represents t/c, and c must be >0.
type Rat = (ExprZ, Integer)

{-
(ls,us) represents
{ x | ∀(M,c)∈ls. M/c≤x, ∀(M,c)∈us. x≤M/c }
-}
type BoundsZ = ([Rat],[Rat])

evalBoundsZ :: Model Integer -> BoundsZ -> IntervalZ
evalBoundsZ model (ls,us) =
  foldl' intersectZ univZ $
    [ (Just (ceiling (LA.eval model x % c)), Nothing) | (x,c) <- ls ] ++
    [ (Nothing, Just (floor (LA.eval model x % c))) | (x,c) <- us ]

collectBoundsZ :: Var -> [Constr] -> (BoundsZ, [Constr])
collectBoundsZ v = foldr phi (([],[]),[])
  where
    phi :: Constr -> (BoundsZ,[Constr]) -> (BoundsZ,[Constr])
    phi constr@(IsNonneg t) ((ls,us),xs) =
      case LA.extract v t of
        (c,t') ->
          case c `compare` 0 of
            EQ -> ((ls, us), constr : xs)
            GT -> (((negateV t', c) : ls, us), xs) -- 0 ≤ cx + M ⇔ -M/c ≤ x
            LT -> ((ls, (t', negate c) : us), xs)   -- 0 ≤ cx + M ⇔ x ≤ M/-c
    phi constr@(IsZero _) (bnd,xs) = (bnd, constr : xs) -- we assume v does not appear in constr

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

solve' :: Options -> VarSet -> [Constr] -> Maybe (Model Integer)
solve' opt vs2 xs = simplify xs >>= go vs2
  where
    go :: VarSet -> [Constr] -> Maybe (Model Integer)
    go vs cs | IS.null vs = do
      let m :: Model Integer
          m = IM.empty
      guard $ all (eval m) cs
      return m
    go vs cs | Just (e,cs2) <- extractEq cs = do
      (vs',cs3, mt) <- eliminateEq e vs cs2
      m <- go vs' =<< simplify cs3
      return (mt m)
    go vs cs = do
      guard $ optCheckReal opt vs (map toLAAtom cs)
      if isExact bnd then
        case1
      else
        case1 `mplus` case2

      where
        (v,vs',bnd@(ls,us),rest) = chooseVariable vs cs

        case1 = do
          let zs = [ LA.constant ((a-1)*(b-1)) `leZ` (a *^ d ^-^ b *^ c)
                   | (c,a)<-ls , (d,b)<-us ]
          model <- go vs' =<< simplify (zs ++ rest)
          case pickupZ (evalBoundsZ model bnd) of
            Nothing  -> error "ToySolver.OmegaTest.solve': should not happen"
            Just val -> return $ IM.insert v val model

        case2 = msum
          [ do c <- isZero $ a' *^ LA.var v ^-^ (c' ^+^ LA.constant k)
               model <- go vs (c : cs)
               return model
          | let m = maximum [b | (_,b)<-us]
          , (c',a') <- ls
          , k <- [0 .. (m*a'-a'-m) `div` m]
          ]

isExact :: BoundsZ -> Bool
isExact (ls,us) = and [a==1 || b==1 | (_,a)<-ls , (_,b)<-us]

toLAAtom :: Constr -> LA.Atom Rational
toLAAtom (IsNonneg e) = LA.mapCoeff fromInteger e .>=. LA.constant 0
toLAAtom (IsZero e)   = LA.mapCoeff fromInteger e .==. LA.constant 0

chooseVariable :: VarSet -> [Constr] -> (Var, VarSet, BoundsZ, [Constr])
chooseVariable vs xs = head $ [e | e@(_,_,bnd,_) <- table, isExact bnd] ++ table
  where
    table = [ (v, IS.fromAscList vs', bnd, rest)
            | (v,vs') <- pickup (IS.toAscList vs), let (bnd, rest) = collectBoundsZ v xs
            ]

-- Find an equality constraint (e==0) and extract e with rest of constraints.
extractEq :: [Constr] -> Maybe (ExprZ, [Constr])
extractEq = msum . map f . pickup
  where
    f :: (Constr, [Constr]) -> Maybe (ExprZ, [Constr])
    f (IsZero e, ls) = Just (e, ls)
    f _ = Nothing

-- Eliminate an equality equality constraint (e==0).
eliminateEq :: ExprZ -> VarSet -> [Constr] -> Maybe (VarSet, [Constr], Model Integer -> Model Integer)
eliminateEq e vs cs | Just c <- LA.asConst e = guard (c==0) >> return (vs, cs, id)
eliminateEq e vs cs = do
  let (ak,xk) = minimumBy (comparing (abs . fst)) [(a,x) | (a,x) <- LA.terms e, x /= LA.unitVar]
  if abs ak == 1 then
    case LA.extract xk e of
      (_, e') -> do
        let xk_def = signum ak *^ negateV e'
        return $
           ( IS.delete xk vs
           , [applySubst1Constr xk xk_def c | c <- cs]
           , \model -> IM.insert xk (LA.eval model xk_def) model
           )
  else do
    let m = abs ak + 1
    assert (ak `zmod` m == - signum ak) $ return ()
    let -- sigma is a fresh variable
        sigma = if IS.null vs then 0 else IS.findMax vs + 1
        -- m *^ LA.var sigma == LA.fromTerms [(a `zmod` m, x) | (a,x) <- LA.terms e]
        -- m *^ LA.var sigma == LA.fromTerms [(a `zmod` m, x) | (a,x) <- LA.terms e, x /= xk] ^+^ (ak `zmod` m) *^ LA.var xk
        -- m *^ LA.var sigma == LA.fromTerms [(a `zmod` m, x) | (a,x) <- LA.terms e, x /= xk] ^+^ (- signum ak) *^ LA.var xk
        -- LA.var xk == (- signum ak * m) *^ LA.var sigma ^+^ LA.fromTerms [(signum ak * (a `zmod` m), x) | (a,x) <- LA.terms e, x /= xk]
        xk_def = (- signum ak * m) *^ LA.var sigma ^+^
                   LA.fromTerms [(signum ak * (a `zmod` m), x) | (a,x) <- LA.terms e, x /= xk]
        -- e2 is normalized version of (LA.applySubst1 xk xk_def e).
        e2 = (- abs ak) *^ LA.var sigma ^+^
                LA.fromTerms [((floor (a%m + 1/2) + (a `zmod` m)), x) | (a,x) <- LA.terms e, x /= xk]
    assert (m *^ e2 == LA.applySubst1 xk xk_def e) $ return ()
    let mt :: Model Integer -> Model Integer
        mt model = IM.delete sigma $ IM.insert xk (LA.eval model xk_def) model
    c <- isZero e2
    return (IS.insert sigma (IS.delete xk vs), c : [applySubst1Constr xk xk_def c | c <- cs], mt)

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

type IntervalZ = (Maybe Integer, Maybe Integer)

univZ :: IntervalZ
univZ = (Nothing, Nothing)

intersectZ :: IntervalZ -> IntervalZ -> IntervalZ
intersectZ (l1,u1) (l2,u2) = (combineMaybe max l1 l2, combineMaybe min u1 u2)

pickupZ :: IntervalZ -> Maybe Integer
pickupZ (Nothing,Nothing) = return 0
pickupZ (Just x, Nothing) = return x
pickupZ (Nothing, Just x) = return x
pickupZ (Just x, Just y) = if x <= y then return x else mzero

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

-- | @'solve' opt {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LIA φ@ when
-- such @M@ exists, returns @Nothing@ otherwise.
--
-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
--
solve :: Options -> VarSet -> [LA.Atom Rational] -> Maybe (Model Integer)
solve opt vs cs = msum [solve' opt vs cs | cs <- unDNF dnf]
  where
    dnf = andB (map f cs)
    f (OrdRel lhs op rhs) =
      case op of
        Lt  -> DNF [[a `ltZ` b]]
        Le  -> DNF [[a `leZ` b]]
        Gt  -> DNF [[a `gtZ` b]]
        Ge  -> DNF [[a `geZ` b]]
        Eql -> DNF [[a `eqZ` b]]
        NEq -> DNF [[a `ltZ` b], [a `gtZ` b]]
      where
        (e1,c1) = g lhs
        (e2,c2) = g rhs
        a = c2 *^ e1
        b = c1 *^ e2
    g :: LA.Expr Rational -> (ExprZ, Integer)
    g a = (LA.mapCoeff (\c -> floor (c * fromInteger d)) a, d)
      where
        d = foldl' lcm 1 [denominator c | (c,_) <- LA.terms a]

-- | @'solveQFLIRAConj' {x1,…,xn} φ I@ returns @Just M@ that @M ⊧_LIRA φ@ when
-- such @M@ exists, returns @Nothing@ otherwise.
--
-- * @FV(φ)@ must be a subset of @{x1,…,xn}@.
--
-- * @I@ is a set of integer variables and must be a subset of @{x1,…,xn}@.
--
solveQFLIRAConj :: Options -> VarSet -> [LA.Atom Rational] -> VarSet -> Maybe (Model Rational)
solveQFLIRAConj opt vs cs ivs = listToMaybe $ do
  (cs2, mt) <- FM.projectN rvs cs
  m <- maybeToList $ solve opt ivs cs2
  return $ mt $ IM.map fromInteger m
  where
    rvs = vs `IS.difference` ivs

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

zmod :: Integer -> Integer -> Integer
a `zmod` b = a - b * floor (a % b + 1 / 2)

gcd' :: [Integer] -> Integer
gcd' [] = 1
gcd' xs = foldl1' gcd xs

pickup :: [a] -> [(a,[a])]
pickup [] = []
pickup (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- pickup xs]

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

1	{-# LANGUAGE FlexibleInstances #-}
2	{-# LANGUAGE MultiParamTypeClasses #-}
3	{-# LANGUAGE TypeSynonymInstances #-}
4	{-# OPTIONS_GHC -Wall #-}
5	{-# OPTIONS_HADDOCK show-extensions #-}
6	-----------------------------------------------------------------------------
7	-- \|
8	-- Module : ToySolver.Arith.OmegaTest.Base
9	-- Copyright : (c) Masahiro Sakai 2011
10	-- License : BSD-style
11	--
12	-- Maintainer : masahiro.sakai@gmail.com
13	-- Stability : provisional
14	-- Portability : non-portable
15	--
16	-- (incomplete) implementation of Omega Test
17	--
18	-- References:
19	--
20	-- * William Pugh. The Omega test: a fast and practical integer
21	-- programming algorithm for dependence analysis. In Proceedings of
22	-- the 1991 ACM/IEEE conference on Supercomputing (1991), pp. 4-13.
23	--
24	-- * <http://users.cecs.anu.edu.au/~michaeln/pubs/arithmetic-dps.pdf>
25	--
26	-- See also:
27	--
28	-- * <http://hackage.haskell.org/package/Omega>
29	--
30	-----------------------------------------------------------------------------
31	module ToySolver.Arith.OmegaTest.Base
32	( Model
33	, solve
34	, solveQFLIRAConj
35	, Options (..)
36	, checkRealNoCheck
37	, checkRealByFM
38
39	-- * Exported just for testing
40	, zmod
41	) where
42
43	import Control.Exception (assert)
44	import Control.Monad
45	import Data.Default.Class
46	import Data.List (foldl', foldl1', minimumBy)
47	import Data.Maybe
48	import Data.Ord
49	import Data.Ratio
50	import qualified Data.IntMap as IM
51	import qualified Data.IntSet as IS
52	import Data.VectorSpace
53
54	import ToySolver.Data.OrdRel
55	import ToySolver.Data.Boolean
56	import ToySolver.Data.DNF
57	import qualified ToySolver.Data.LA as LA
58	import ToySolver.Data.IntVar
59	import ToySolver.Internal.Util (combineMaybe)
60	import qualified ToySolver.Arith.FourierMotzkin as FM
61
62	-- ---------------------------------------------------------------------------
63
64	-- \| Options for solving.
65	--
66	-- The default option can be obtained by 'def'.
67	data Options
68	= Options
69	{ optCheckReal :: VarSet -> [LA.Atom Rational] -> Bool	1✔
70	-- ^ optCheckReal is used for checking whether real shadow is satisfiable.
71	--
72	-- * If it returns @True@, the real shadow may or may not be satisfiable.
73	--
74	-- * If it returns @False@, the real shadow must be unsatisfiable
75	}
76
77	instance Default Options where
78	def =	1✔
79	Options	1✔
80	{ optCheckReal =
81	-- checkRealNoCheck
82	checkRealByFM	1✔
83	}
84
85	checkRealNoCheck :: VarSet -> [LA.Atom Rational] -> Bool
86	checkRealNoCheck _ _ = True	×
87
88	checkRealByFM :: VarSet -> [LA.Atom Rational] -> Bool
89	checkRealByFM vs as = isJust $ FM.solve vs as	1✔
90
91	-- ---------------------------------------------------------------------------
92
93	type ExprZ = LA.Expr Integer
94
95	-- \| Atomic constraint
96	data Constr
97	= IsNonneg ExprZ
98	-- ^ e ≥ 0
99	\| IsZero ExprZ
100	-- ^ e = 0
101	deriving (Show, Eq, Ord)	×
102
103	instance Variables Constr where
104	vars (IsNonneg t) = vars t	×
105	vars (IsZero t) = vars t	×
106
107	-- 制約集合の単純化
108	-- It returns Nothing when a inconsistency is detected.
109	simplify :: [Constr] -> Maybe [Constr]
110	simplify = fmap concat . mapM f	1✔
111	where
112	f :: Constr -> Maybe [Constr]
113	f constr@(IsNonneg e) =	1✔
114	case LA.asConst e of	1✔
115	Just x -> guard (x >= 0) >> return []	1✔
116	Nothing -> return [constr]	1✔
117	f constr@(IsZero e) =
118	case LA.asConst e of	1✔
119	Just x -> guard (x == 0) >> return []	×
120	Nothing -> return [constr]	1✔
121
122	leZ, ltZ, geZ, gtZ, eqZ :: ExprZ -> ExprZ -> Constr
123	-- Note that constants may be floored by division
124	leZ e1 e2 = IsNonneg (LA.mapCoeff (`div` d) e)	1✔
125	where
126	e = e2 ^-^ e1	1✔
127	d = abs $ gcd' [c \| (c,v) <- LA.terms e, v /= LA.unitVar]	1✔
128	ltZ e1 e2 = (e1 ^+^ LA.constant 1) `leZ` e2	1✔
129	geZ = flip leZ	1✔
130	gtZ = flip ltZ	1✔
131	eqZ e1 e2 =	1✔
132	case isZero (e1 ^-^ e2) of	1✔
133	Just c -> c	1✔
134	Nothing -> IsZero (LA.constant 1) -- unsatisfiable	1✔
135
136	isZero :: ExprZ -> Maybe Constr
137	isZero e	1✔
138	= if LA.coeff LA.unitVar e `mod` d == 0	1✔
139	then Just $ IsZero $ LA.mapCoeff (`div` d) e	1✔
140	else Nothing	1✔
141	where
142	d = abs $ gcd' [c \| (c,v) <- LA.terms e, v /= LA.unitVar]	1✔
143
144	applySubst1Constr :: Var -> ExprZ -> Constr -> Constr
145	applySubst1Constr v e (IsNonneg e2) = LA.applySubst1 v e e2 `geZ` LA.constant 0	1✔
146	applySubst1Constr v e (IsZero e2) = LA.applySubst1 v e e2 `eqZ` LA.constant 0	1✔
147
148	instance Eval (Model Integer) Constr Bool where
149	eval m (IsNonneg t) = LA.eval m t >= 0	×
150	eval m (IsZero t) = LA.eval m t == 0	×
151
152	-- ---------------------------------------------------------------------------
153
154	-- \| (t,c) represents t/c, and c must be >0.
155	type Rat = (ExprZ, Integer)
156
157	{-
158	(ls,us) represents
159	{ x \| ∀(M,c)∈ls. M/c≤x, ∀(M,c)∈us. x≤M/c }
160	-}
161	type BoundsZ = ([Rat],[Rat])
162
163	evalBoundsZ :: Model Integer -> BoundsZ -> IntervalZ
164	evalBoundsZ model (ls,us) =	1✔
165	foldl' intersectZ univZ $	1✔
166	[ (Just (ceiling (LA.eval model x % c)), Nothing) \| (x,c) <- ls ] ++	1✔
167	[ (Nothing, Just (floor (LA.eval model x % c))) \| (x,c) <- us ]	1✔
168
169	collectBoundsZ :: Var -> [Constr] -> (BoundsZ, [Constr])
170	collectBoundsZ v = foldr phi (([],[]),[])	1✔
171	where
172	phi :: Constr -> (BoundsZ,[Constr]) -> (BoundsZ,[Constr])
173	phi constr@(IsNonneg t) ((ls,us),xs) =	1✔
174	case LA.extract v t of	1✔
175	(c,t') ->
176	case c `compare` 0 of	1✔
177	EQ -> ((ls, us), constr : xs)	1✔
178	GT -> (((negateV t', c) : ls, us), xs) -- 0 ≤ cx + M ⇔ -M/c ≤ x	1✔
179	LT -> ((ls, (t', negate c) : us), xs) -- 0 ≤ cx + M ⇔ x ≤ M/-c	1✔
180	phi constr@(IsZero _) (bnd,xs) = (bnd, constr : xs) -- we assume v does not appear in constr	×
181
182	-- ---------------------------------------------------------------------------
183
184	solve' :: Options -> VarSet -> [Constr] -> Maybe (Model Integer)
185	solve' opt vs2 xs = simplify xs >>= go vs2	1✔
186	where
187	go :: VarSet -> [Constr] -> Maybe (Model Integer)
188	go vs cs \| IS.null vs = do	1✔
189	let m :: Model Integer
190	m = IM.empty	1✔
191	guard $ all (eval m) cs	×
192	return m	1✔
193	go vs cs \| Just (e,cs2) <- extractEq cs = do	1✔
194	(vs',cs3, mt) <- eliminateEq e vs cs2	1✔
195	m <- go vs' =<< simplify cs3	1✔
196	return (mt m)	1✔
197	go vs cs = do	1✔
198	guard $ optCheckReal opt vs (map toLAAtom cs)	1✔
199	if isExact bnd then	1✔
200	case1	1✔
201	else
202	case1 `mplus` case2	1✔
203
204	where
205	(v,vs',bnd@(ls,us),rest) = chooseVariable vs cs	1✔
206
207	case1 = do	1✔
208	let zs = [ LA.constant ((a-1)(b-1)) `leZ` (a ^ d ^-^ b *^ c)	1✔
209	\| (c,a)<-ls , (d,b)<-us ]	1✔
210	model <- go vs' =<< simplify (zs ++ rest)	1✔
211	case pickupZ (evalBoundsZ model bnd) of	1✔
212	Nothing -> error "ToySolver.OmegaTest.solve': should not happen"	×
213	Just val -> return $ IM.insert v val model	1✔
214
215	case2 = msum	1✔
216	[ do c <- isZero $ a' *^ LA.var v ^-^ (c' ^+^ LA.constant k)	1✔
217	model <- go vs (c : cs)	1✔
218	return model	×
219	\| let m = maximum [b \| (_,b)<-us]	1✔
220	, (c',a') <- ls	1✔
221	, k <- [0 .. (m*a'-a'-m) `div` m]	1✔
222	]
223
224	isExact :: BoundsZ -> Bool
225	isExact (ls,us) = and [a==1 \|\| b==1 \| (_,a)<-ls , (_,b)<-us]	1✔
226
227	toLAAtom :: Constr -> LA.Atom Rational
228	toLAAtom (IsNonneg e) = LA.mapCoeff fromInteger e .>=. LA.constant 0	1✔
229	toLAAtom (IsZero e) = LA.mapCoeff fromInteger e .==. LA.constant 0	×
230
231	chooseVariable :: VarSet -> [Constr] -> (Var, VarSet, BoundsZ, [Constr])
232	chooseVariable vs xs = head $ [e \| e@(_,_,bnd,_) <- table, isExact bnd] ++ table	1✔
233	where
234	table = [ (v, IS.fromAscList vs', bnd, rest)	1✔
235	\| (v,vs') <- pickup (IS.toAscList vs), let (bnd, rest) = collectBoundsZ v xs	1✔
236	]
237
238	-- Find an equality constraint (e==0) and extract e with rest of constraints.
239	extractEq :: [Constr] -> Maybe (ExprZ, [Constr])
240	extractEq = msum . map f . pickup	1✔
241	where
242	f :: (Constr, [Constr]) -> Maybe (ExprZ, [Constr])
243	f (IsZero e, ls) = Just (e, ls)	1✔
244	f _ = Nothing	1✔
245
246	-- Eliminate an equality equality constraint (e==0).
247	eliminateEq :: ExprZ -> VarSet -> [Constr] -> Maybe (VarSet, [Constr], Model Integer -> Model Integer)
248	eliminateEq e vs cs \| Just c <- LA.asConst e = guard (c==0) >> return (vs, cs, id)	×
249	eliminateEq e vs cs = do	1✔
250	let (ak,xk) = minimumBy (comparing (abs . fst)) [(a,x) \| (a,x) <- LA.terms e, x /= LA.unitVar]	1✔
251	if abs ak == 1 then	1✔
252	case LA.extract xk e of	1✔
253	(_, e') -> do	1✔
254	let xk_def = signum ak *^ negateV e'	1✔
255	return $	1✔
256	( IS.delete xk vs	1✔
257	, [applySubst1Constr xk xk_def c \| c <- cs]	1✔
258	, \model -> IM.insert xk (LA.eval model xk_def) model	1✔
259	)
260	else do	1✔
261	let m = abs ak + 1	1✔
262	assert (ak `zmod` m == - signum ak) $ return ()	×
263	let -- sigma is a fresh variable
264	sigma = if IS.null vs then 0 else IS.findMax vs + 1	×
265	-- m *^ LA.var sigma == LA.fromTerms [(a `zmod` m, x) \| (a,x) <- LA.terms e]
266	-- m ^ LA.var sigma == LA.fromTerms [(a `zmod` m, x) \| (a,x) <- LA.terms e, x /= xk] ^+^ (ak `zmod` m) ^ LA.var xk
267	-- m ^ LA.var sigma == LA.fromTerms [(a `zmod` m, x) \| (a,x) <- LA.terms e, x /= xk] ^+^ (- signum ak) ^ LA.var xk
268	-- LA.var xk == (- signum ak * m) ^ LA.var sigma ^+^ LA.fromTerms [(signum ak (a `zmod` m), x) \| (a,x) <- LA.terms e, x /= xk]
269	xk_def = (- signum ak * m) *^ LA.var sigma ^+^	1✔
270	LA.fromTerms [(signum ak * (a `zmod` m), x) \| (a,x) <- LA.terms e, x /= xk]	1✔
271	-- e2 is normalized version of (LA.applySubst1 xk xk_def e).
272	e2 = (- abs ak) *^ LA.var sigma ^+^	1✔
273	LA.fromTerms [((floor (a%m + 1/2) + (a `zmod` m)), x) \| (a,x) <- LA.terms e, x /= xk]	1✔
274	assert (m *^ e2 == LA.applySubst1 xk xk_def e) $ return ()	×
275	let mt :: Model Integer -> Model Integer
276	mt model = IM.delete sigma $ IM.insert xk (LA.eval model xk_def) model	1✔
277	c <- isZero e2	1✔
278	return (IS.insert sigma (IS.delete xk vs), c : [applySubst1Constr xk xk_def c \| c <- cs], mt)	1✔
279
280	-- ---------------------------------------------------------------------------
281
282	type IntervalZ = (Maybe Integer, Maybe Integer)
283
284	univZ :: IntervalZ
285	univZ = (Nothing, Nothing)	1✔
286
287	intersectZ :: IntervalZ -> IntervalZ -> IntervalZ
288	intersectZ (l1,u1) (l2,u2) = (combineMaybe max l1 l2, combineMaybe min u1 u2)	1✔
289
290	pickupZ :: IntervalZ -> Maybe Integer
291	pickupZ (Nothing,Nothing) = return 0	1✔
292	pickupZ (Just x, Nothing) = return x	1✔
293	pickupZ (Nothing, Just x) = return x	1✔
294	pickupZ (Just x, Just y) = if x <= y then return x else mzero	×
295
296	-- ---------------------------------------------------------------------------
297
298	-- \| @'solve' opt {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LIA φ@ when
299	-- such @M@ exists, returns @Nothing@ otherwise.
300	--
301	-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
302	--
303	solve :: Options -> VarSet -> [LA.Atom Rational] -> Maybe (Model Integer)
304	solve opt vs cs = msum [solve' opt vs cs \| cs <- unDNF dnf]	1✔
305	where
306	dnf = andB (map f cs)	1✔
307	f (OrdRel lhs op rhs) =	1✔
308	case op of	1✔
309	Lt -> DNF [[a `ltZ` b]]	1✔
310	Le -> DNF [[a `leZ` b]]	1✔
311	Gt -> DNF [[a `gtZ` b]]	1✔
312	Ge -> DNF [[a `geZ` b]]	1✔
313	Eql -> DNF [[a `eqZ` b]]	1✔
314	NEq -> DNF [[a `ltZ` b], [a `gtZ` b]]	×
315	where
316	(e1,c1) = g lhs	1✔
317	(e2,c2) = g rhs	1✔
318	a = c2 *^ e1	1✔
319	b = c1 *^ e2	1✔
320	g :: LA.Expr Rational -> (ExprZ, Integer)
321	g a = (LA.mapCoeff (\c -> floor (c * fromInteger d)) a, d)	1✔
322	where
323	d = foldl' lcm 1 [denominator c \| (c,_) <- LA.terms a]	1✔
324
325	-- \| @'solveQFLIRAConj' {x1,…,xn} φ I@ returns @Just M@ that @M ⊧_LIRA φ@ when
326	-- such @M@ exists, returns @Nothing@ otherwise.
327	--
328	-- * @FV(φ)@ must be a subset of @{x1,…,xn}@.
329	--
330	-- * @I@ is a set of integer variables and must be a subset of @{x1,…,xn}@.
331	--
332	solveQFLIRAConj :: Options -> VarSet -> [LA.Atom Rational] -> VarSet -> Maybe (Model Rational)
333	solveQFLIRAConj opt vs cs ivs = listToMaybe $ do	×
334	(cs2, mt) <- FM.projectN rvs cs	×
335	m <- maybeToList $ solve opt ivs cs2	×
336	return $ mt $ IM.map fromInteger m	×
337	where
338	rvs = vs `IS.difference` ivs	×
339
340	-- ---------------------------------------------------------------------------
341
342	zmod :: Integer -> Integer -> Integer
343	a `zmod` b = a - b * floor (a % b + 1 / 2)	1✔
344
345	gcd' :: [Integer] -> Integer
346	gcd' [] = 1	1✔
347	gcd' xs = foldl1' gcd xs	1✔
348
349	pickup :: [a] -> [(a,[a])]
350	pickup [] = []	1✔
351	pickup (x:xs) = (x,xs) : [(y,x:ys) \| (y,ys) <- pickup xs]	1✔
352
353	-- ---------------------------------------------------------------------------

msakai / toysolver / 797

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous