msakai / toysolver / 604

Committed 06 Feb 2025 03:40AM UTC coverage: 71.288% (-0.07%) from 71.354%

Build # 604

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #145 from msakai/update-stack-lts-202502

Update stack resolvers (2025-02)

Run Details

10696 of 15004 relevant lines covered (71.29%)

0.71 hits per line

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

59.77

/src/ToySolver/Arith/Cooper/Base.hs

{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.Cooper.Base
-- Copyright   :  (c) Masahiro Sakai 2011-2014
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Naive implementation of Cooper's algorithm
--
-- Reference:
--
-- * <http://hagi.is.s.u-tokyo.ac.jp/pub/staff/hagiya/kougiroku/ronri/5.txt>
--
-- * <http://www.cs.cmu.edu/~emc/spring06/home1_files/Presburger%20Arithmetic.ppt>
--
-- See also:
--
-- * <http://hackage.haskell.org/package/presburger>
--
-----------------------------------------------------------------------------
module ToySolver.Arith.Cooper.Base
    (
    -- * Language of presburger arithmetics
      ExprZ
    , Lit (..)
    , QFFormula
    , fromLAAtom
    , (.|.)
    , Model
    , Eval (..)

    -- * Projection
    , project
    , projectN
    , projectCases
    , projectCasesN

    -- * Constraint solving
    , solve
    , solveQFFormula
    , solveQFLIRAConj
    ) where

import Control.Monad
import qualified Data.Foldable as Foldable
import Data.List
import Data.Maybe
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.Ratio
import qualified Data.Semigroup as Semigroup
import Data.Set (Set)
import qualified Data.Set as Set
import Data.VectorSpace hiding (project)

import ToySolver.Data.OrdRel
import ToySolver.Data.Boolean
import ToySolver.Data.BoolExpr (BoolExpr (..))
import qualified ToySolver.Data.BoolExpr as BoolExpr
import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar
import qualified ToySolver.Arith.FourierMotzkin as FM

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

-- | Linear arithmetic expression over integers.
type ExprZ = LA.Expr Integer

fromLAAtom :: LA.Atom Rational -> QFFormula
fromLAAtom (OrdRel a op b) = ordRel op a' b'
  where
    (e1,c1) = toRat a
    (e2,c2) = toRat b
    a' = c2 *^ e1
    b' = c1 *^ e2

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

toRat :: LA.Expr Rational -> Rat
toRat e = seq m $ (LA.mapCoeff f e, m)
  where
    f x = numerator (fromInteger m * x)
    m = foldl' lcm 1 [denominator c | (c,_) <- LA.terms e]

leZ, ltZ, geZ, gtZ :: ExprZ -> ExprZ -> Lit
leZ e1 e2 = e1 `ltZ` (e2 ^+^ LA.constant 1)
ltZ e1 e2 = IsPos $ (e2 ^-^ e1)
geZ = flip leZ
gtZ = flip ltZ

eqZ :: ExprZ -> ExprZ -> QFFormula
eqZ e1 e2 = Atom (e1 `leZ` e2) .&&. Atom (e1 `geZ` e2)

-- | Literals of Presburger arithmetic.
data Lit
    = IsPos ExprZ
    -- ^ @IsPos e@ means @e > 0@
    | Divisible Bool Integer ExprZ
    -- ^
    -- * @Divisible True d e@ means @e@ can be divided by @d@ (i.e. @d | e@)
    -- * @Divisible False d e@ means @e@ can not be divided by @d@ (i.e. @¬(d | e)@)
    deriving (Show, Eq, Ord)

instance Variables Lit where
  vars (IsPos t) = vars t
  vars (Divisible _ _ t) = vars t

instance Complement Lit where
  notB (IsPos e) = e `leZ` LA.constant 0
  notB (Divisible b c e) = Divisible (not b) c e

-- | Quantifier-free formula of Presburger arithmetic.
type QFFormula = BoolExpr Lit

instance IsEqRel (LA.Expr Integer) QFFormula where
  a .==. b = eqZ a b
  a ./=. b = notB $ eqZ a b

instance IsOrdRel (LA.Expr Integer) QFFormula where
  ordRel op lhs rhs =
    case op of
      Le  -> Atom $ leZ lhs rhs
      Ge  -> Atom $ geZ lhs rhs
      Lt  -> Atom $ ltZ lhs rhs
      Gt  -> Atom $ gtZ lhs rhs
      Eql -> lhs .==. rhs
      NEq -> lhs ./=. rhs

-- | @d | e@ means @e@ can be divided by @d@.
(.|.) :: Integer -> ExprZ -> QFFormula
n .|. e = Atom $ Divisible True n e

subst1 :: Var -> ExprZ -> QFFormula -> QFFormula
subst1 x e = fmap f
  where
    f (Divisible b c e1) = Divisible b c $ LA.applySubst1 x e e1
    f (IsPos e1) = IsPos $ LA.applySubst1 x e e1

simplify :: QFFormula -> QFFormula
simplify = BoolExpr.simplify . BoolExpr.fold simplifyLit

simplifyLit :: Lit -> QFFormula
simplifyLit (IsPos e) =
  case LA.asConst e of
    Just c -> if c > 0 then true else false
    Nothing ->
      -- e > 0  <=>  e - 1 >= 0
      -- <=>  LA.mapCoeff (`div` d) (e - 1) >= 0
      -- <=>  LA.mapCoeff (`div` d) (e - 1) + 1 > 0
      Atom $ IsPos $ LA.mapCoeff (`div` d) e1 ^+^ LA.constant 1
  where
    e1 = e ^-^ LA.constant 1
    d  = if null cs then 1 else abs $ foldl1' gcd cs
    cs = [c | (c,x) <- LA.terms e1, x /= LA.unitVar]
simplifyLit lit@(Divisible b c e)
  | LA.coeff LA.unitVar e2 `mod` d /= 0 = if b then false else true
  | c' == 1   = if b then true else false
  | d  == 1   = Atom lit
  | otherwise = Atom $ Divisible b c' e'
  where
    e2 = LA.mapCoeff (`mod` c) e
    d  = abs $ foldl' gcd c [c2 | (c2,x) <- LA.terms e2, x /= LA.unitVar]
    c' = c `checkedDiv` d
    e' = LA.mapCoeff (`checkedDiv` d) e2

instance Eval (Model Integer) Lit Bool where
  eval m (IsPos e) = LA.eval m e > 0
  eval m (Divisible True n e)  = LA.eval m e `mod` n == 0
  eval m (Divisible False n e) = LA.eval m e `mod` n /= 0

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

data Witness = WCase1 Integer ExprZ | WCase2 Integer Integer Integer (Set ExprZ)
  deriving (Show)

instance Eval (Model Integer) Witness Integer where
  eval model (WCase1 c e) = LA.eval model e `checkedDiv` c
  eval model (WCase2 c j delta us)
    | Set.null us' = j `checkedDiv` c
    | otherwise = (j + (((u - delta - 1) `div` delta) * delta)) `checkedDiv` c
    where
      us' = Set.map (LA.eval model) us
      u = Set.findMin us'

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

{-| @'project' x φ@ returns @(ψ, lift)@ such that:

* @⊢_LIA ∀y1, …, yn. (∃x. φ) ↔ ψ@ where @{y1, …, yn} = FV(φ) \\ {x}@, and

* if @M ⊧_LIA ψ@ then @lift M ⊧_LIA φ@.
-}
project :: Var -> QFFormula -> (QFFormula, Model Integer -> Model Integer)
project x formula = (formula', mt)
  where
    xs = projectCases x formula
    formula' = orB' [phi | (phi,_) <- xs]
    mt m = head $ do
      (phi, mt') <- xs
      guard $ eval m phi
      return $ mt' m
    orB' = orB . concatMap f
      where
        f (Or xs) = concatMap f xs
        f x = [x]

{-| @'projectN' {x1,…,xm} φ@ returns @(ψ, lift)@ such that:

* @⊢_LIA ∀y1, …, yn. (∃x1, …, xm. φ) ↔ ψ@ where @{y1, …, yn} = FV(φ) \\ {x1,…,xm}@, and

* if @M ⊧_LIA ψ@ then @lift M ⊧_LIA φ@.
-}
projectN :: VarSet -> QFFormula -> (QFFormula, Model Integer -> Model Integer)
projectN vs2 = f (IS.toList vs2)
  where
    f :: [Var] -> QFFormula -> (QFFormula, Model Integer -> Model Integer)
    f [] formula     = (formula, id)
    f (v:vs) formula = (formula3, mt1 . mt2)
      where
        (formula2, mt1) = project v formula
        (formula3, mt2) = f vs formula2

{-| @'projectCases' x φ@ returns @[(ψ_1, lift_1), …, (ψ_m, lift_m)]@ such that:

* @⊢_LIA ∀y1, …, yn. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_m)@ where @{y1, …, yn} = FV(φ) \\ {x}@, and

* if @M ⊧_LIA ψ_i@ then @lift_i M ⊧_LIA φ@.
-}
projectCases :: Var -> QFFormula -> [(QFFormula, Model Integer -> Model Integer)]
projectCases x formula = do
  (phi, wit) <- projectCases' x formula
  return (phi, \m -> IM.insert x (eval m wit) m)

projectCases' :: Var -> QFFormula -> [(QFFormula, Witness)]
projectCases' x formula = [(phi', w) | (phi,w) <- case1 ++ case2, let phi' = simplify phi, phi' /= false]
  where
    -- eliminate Not, Imply and Equiv.
    formula0 :: QFFormula
    formula0 = pos formula
      where
        pos (Atom a) = Atom a
        pos (And xs) = And (map pos xs)
        pos (Or xs) = Or (map pos xs)
        pos (Not x) = neg x
        pos (Imply x y) = neg x .||. pos y
        pos (Equiv x y) = pos ((x .=>. y) .&&. (y .=>. x))
        pos (ITE c t e) = pos ((c .=>. t) .&&. (Not c .=>. e))

        neg (Atom a) = Atom (notB a)
        neg (And xs) = Or (map neg xs)
        neg (Or xs) = And (map neg xs)
        neg (Not x) = pos x
        neg (Imply x y) = pos x .&&. neg y
        neg (Equiv x y) = neg ((x .=>. y) .&&. (y .=>. x))
        neg (ITE c t e) = neg ((c .=>. t) .&&. (Not c .=>. e))

    -- xの係数の最小公倍数
    c :: Integer
    c = getLCM $ Foldable.foldMap f formula0
      where
         f (IsPos e) = LCM $ fromMaybe 1 (LA.lookupCoeff x e)
         f (Divisible _ _ e) = LCM $ fromMaybe 1 (LA.lookupCoeff x e)

    -- 式をスケールしてxの係数の絶対値をcへと変換し、その後cxをxで置き換え、
    -- xがcで割り切れるという制約を追加した論理式
    formula1 :: QFFormula
    formula1 = simplify $ fmap f formula0 .&&. (c .|. LA.var x)
      where
        f lit@(IsPos e) =
          case LA.lookupCoeff x e of
            Nothing -> lit
            Just a ->
              let s = abs (c `checkedDiv` a)
              in IsPos $ g s e
        f lit@(Divisible b d e) =
          case LA.lookupCoeff x e of
            Nothing -> lit
            Just a ->
              let s = abs (c `checkedDiv` a)
              in Divisible b (s*d) $ g s e

        g :: Integer -> ExprZ -> ExprZ
        g s = LA.mapCoeffWithVar (\c' x' -> if x==x' then signum c' else s*c')

    -- d|x+t という形の論理式の d の最小公倍数
    delta :: Integer
    delta = getLCM $ Foldable.foldMap f formula1
      where
        f (Divisible _ d _) = LCM d
        f (IsPos _) = LCM 1

    -- ts = {t | t < x は formula1 に現れる原子論理式}
    ts :: Set ExprZ
    ts = Foldable.foldMap f formula1
      where
        f (Divisible _ _ _) = Set.empty
        f (IsPos e) =
          case LA.extractMaybe x e of
            Nothing -> Set.empty
            Just (1, e') -> Set.singleton (negateV e') -- IsPos e <=> (x + e' > 0) <=> (-e' < x)
            Just (-1, _) -> Set.empty -- IsPos e <=> (-x + e' > 0) <=> (x < e')
            _ -> error "should not happen"

    -- formula1を真にする最小のxが存在する場合
    case1 :: [(QFFormula, Witness)]
    case1 = [ (subst1 x e formula1, WCase1 c e)
            | j <- [1..delta], t <- Set.toList ts, let e = t ^+^ LA.constant j ]

    -- formula1のなかの x < t を真に t < x を偽に置き換えた論理式
    formula2 :: QFFormula
    formula2 = simplify $ BoolExpr.fold f formula1
      where
        f lit@(IsPos e) =
          case LA.lookupCoeff x e of
            Nothing -> Atom lit
            Just 1    -> false -- IsPos e <=> ( x + e' > 0) <=> -e' < x
            Just (-1) -> true  -- IsPos e <=> (-x + e' > 0) <=>  x  < e'
            _ -> error "should not happen"
        f lit@(Divisible _ _ _) = Atom lit

    -- us = {u | x < u は formula1 に現れる原子論理式}
    us :: Set ExprZ
    us = Foldable.foldMap f formula1
      where
        f (IsPos e) =
          case LA.extractMaybe x e of
            Nothing -> Set.empty
            Just (1, _)   -> Set.empty -- IsPos e <=> (x + e' > 0) <=> -e' < x
            Just (-1, e') -> Set.singleton e' -- IsPos e <=> (-x + e' > 0) <=>  x  < e'
            _ -> error "should not happen"
        f (Divisible _ _ _) = Set.empty

    -- formula1を真にする最小のxが存在しない場合
    case2 :: [(QFFormula, Witness)]
    case2 = [(subst1 x (LA.constant j) formula2, WCase2 c j delta us) | j <- [1..delta]]

{-| @'projectCasesN' {x1,…,xm} φ@ returns @[(ψ_1, lift_1), …, (ψ_n, lift_n)]@ such that:

* @⊢_LIA ∀y1, …, yp. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_n)@ where @{y1, …, yp} = FV(φ) \\ {x1,…,xm}@, and

* if @M ⊧_LIA ψ_i@ then @lift_i M ⊧_LIA φ@.
-}
projectCasesN :: VarSet -> QFFormula -> [(QFFormula, Model Integer -> Model Integer)]
projectCasesN vs2 = f (IS.toList vs2)
  where
    f :: [Var] -> QFFormula -> [(QFFormula, Model Integer -> Model Integer)]
    f [] formula = return (formula, id)
    f (v:vs) formula = do
      (formula2, mt1) <- projectCases v formula
      (formula3, mt2) <- f vs formula2
      return (formula3, mt1 . mt2)

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

newtype LCM a = LCM{ getLCM :: a }

instance Integral a => Semigroup.Semigroup (LCM a) where
  LCM a <> LCM b = LCM $ lcm a b
  stimes = Semigroup.stimesIdempotent

instance Integral a => Monoid (LCM a) where
  mempty = LCM 1

checkedDiv :: Integer -> Integer -> Integer
checkedDiv a b =
  case a `divMod` b of
    (q,0) -> q
    _ -> error "ToySolver.Cooper.checkedDiv: should not happen"

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

-- | @'solveQFFormula' {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LIA φ@ when
-- such @M@ exists, returns @Nothing@ otherwise.
--
-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
--
solveQFFormula :: VarSet -> QFFormula -> Maybe (Model Integer)
solveQFFormula vs formula = listToMaybe $ do
  (formula2, mt) <- projectCasesN vs formula
  let m :: Model Integer
      m = IM.empty
  guard $ eval m formula2
  return $ mt m

-- | @'solve' {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LIA φ@ when
-- such @M@ exists, returns @Nothing@ otherwise.
--
-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
--
solve :: VarSet -> [LA.Atom Rational] -> Maybe (Model Integer)
solve vs cs = solveQFFormula vs $ andB $ map fromLAAtom cs

-- | @'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 :: VarSet -> [LA.Atom Rational] -> VarSet -> Maybe (Model Rational)
solveQFLIRAConj vs cs ivs = listToMaybe $ do
  (cs2, mt) <- FM.projectN rvs cs
  m <- maybeToList $ solve ivs cs2
  return $ mt $ IM.map fromInteger m
  where
    rvs = vs `IS.difference` ivs

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

_testHagiya :: QFFormula
_testHagiya = fst $ project 1 $ andB [c1, c2, c3]
  where
    [x,y,z] = map LA.var [1..3]
    c1 = x .<. (y ^+^ y)
    c2 = z .<. x
    c3 = 3 .|. x

{-
∃ x. x < y+y ∧ z<x ∧ 3|x
⇔
(2y>z+1 ∧ 3|z+1) ∨ (2y>z+2 ∧ 3|z+2) ∨ (2y>z+3 ∧ 3|z+3)
-}

_test3 :: QFFormula
_test3 = fst $ project 1 $ andB [p1,p2,p3,p4]
  where
    x = LA.var 0
    y = LA.var 1
    p1 = LA.constant 0 .<. y
    p2 = 2 *^ x .<. y
    p3 = y .<. x ^+^ LA.constant 2
    p4 = 2 .|. y

{-
∃ y. 0 < y ∧ 2x<y ∧ y < x+2 ∧ 2|y
⇔
(2x < 2 ∧ 0 < x) ∨ (0 < 2x+2 ∧ x < 0)
-}

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

1	{-# OPTIONS_GHC -Wall #-}
2	{-# OPTIONS_HADDOCK show-extensions #-}
3	{-# LANGUAGE FlexibleInstances #-}
4	{-# LANGUAGE MultiParamTypeClasses #-}
5	-----------------------------------------------------------------------------
6	-- \|
7	-- Module : ToySolver.Arith.Cooper.Base
8	-- Copyright : (c) Masahiro Sakai 2011-2014
9	-- License : BSD-style
10	--
11	-- Maintainer : masahiro.sakai@gmail.com
12	-- Stability : provisional
13	-- Portability : non-portable
14	--
15	-- Naive implementation of Cooper's algorithm
16	--
17	-- Reference:
18	--
19	-- * <http://hagi.is.s.u-tokyo.ac.jp/pub/staff/hagiya/kougiroku/ronri/5.txt>
20	--
21	-- * <http://www.cs.cmu.edu/~emc/spring06/home1_files/Presburger%20Arithmetic.ppt>
22	--
23	-- See also:
24	--
25	-- * <http://hackage.haskell.org/package/presburger>
26	--
27	-----------------------------------------------------------------------------
28	module ToySolver.Arith.Cooper.Base
29	(
30	-- * Language of presburger arithmetics
31	ExprZ
32	, Lit (..)
33	, QFFormula
34	, fromLAAtom
35	, (.\|.)
36	, Model
37	, Eval (..)
38
39	-- * Projection
40	, project
41	, projectN
42	, projectCases
43	, projectCasesN
44
45	-- * Constraint solving
46	, solve
47	, solveQFFormula
48	, solveQFLIRAConj
49	) where
50
51	import Control.Monad
52	import qualified Data.Foldable as Foldable
53	import Data.List
54	import Data.Maybe
55	import qualified Data.IntMap as IM
56	import qualified Data.IntSet as IS
57	import Data.Ratio
58	import qualified Data.Semigroup as Semigroup
59	import Data.Set (Set)
60	import qualified Data.Set as Set
61	import Data.VectorSpace hiding (project)
62
63	import ToySolver.Data.OrdRel
64	import ToySolver.Data.Boolean
65	import ToySolver.Data.BoolExpr (BoolExpr (..))
66	import qualified ToySolver.Data.BoolExpr as BoolExpr
67	import qualified ToySolver.Data.LA as LA
68	import ToySolver.Data.IntVar
69	import qualified ToySolver.Arith.FourierMotzkin as FM
70
71	-- ---------------------------------------------------------------------------
72
73	-- \| Linear arithmetic expression over integers.
74	type ExprZ = LA.Expr Integer
75
76	fromLAAtom :: LA.Atom Rational -> QFFormula
77	fromLAAtom (OrdRel a op b) = ordRel op a' b'	1✔
78	where
79	(e1,c1) = toRat a	1✔
80	(e2,c2) = toRat b	1✔
81	a' = c2 *^ e1	1✔
82	b' = c1 *^ e2	1✔
83
84	-- \| (t,c) represents t/c, and c must be >0.
85	type Rat = (ExprZ, Integer)
86
87	toRat :: LA.Expr Rational -> Rat
88	toRat e = seq m $ (LA.mapCoeff f e, m)	1✔
89	where
90	f x = numerator (fromInteger m * x)	1✔
91	m = foldl' lcm 1 [denominator c \| (c,_) <- LA.terms e]	1✔
92
93	leZ, ltZ, geZ, gtZ :: ExprZ -> ExprZ -> Lit
94	leZ e1 e2 = e1 `ltZ` (e2 ^+^ LA.constant 1)	1✔
95	ltZ e1 e2 = IsPos $ (e2 ^-^ e1)	1✔
96	geZ = flip leZ	1✔
97	gtZ = flip ltZ	1✔
98
99	eqZ :: ExprZ -> ExprZ -> QFFormula
100	eqZ e1 e2 = Atom (e1 `leZ` e2) .&&. Atom (e1 `geZ` e2)	1✔
101
102	-- \| Literals of Presburger arithmetic.
103	data Lit
104	= IsPos ExprZ
105	-- ^ @IsPos e@ means @e > 0@
106	\| Divisible Bool Integer ExprZ
107	-- ^
108	-- * @Divisible True d e@ means @e@ can be divided by @d@ (i.e. @d \| e@)
109	-- * @Divisible False d e@ means @e@ can not be divided by @d@ (i.e. @¬(d \| e)@)
110	deriving (Show, Eq, Ord)	×
111
112	instance Variables Lit where
113	vars (IsPos t) = vars t	×
114	vars (Divisible _ _ t) = vars t	×
115
116	instance Complement Lit where
117	notB (IsPos e) = e `leZ` LA.constant 0	×
118	notB (Divisible b c e) = Divisible (not b) c e	×
119
120	-- \| Quantifier-free formula of Presburger arithmetic.
121	type QFFormula = BoolExpr Lit
122
123	instance IsEqRel (LA.Expr Integer) QFFormula where
124	a .==. b = eqZ a b	1✔
125	a ./=. b = notB $ eqZ a b	×
126
127	instance IsOrdRel (LA.Expr Integer) QFFormula where	×
128	ordRel op lhs rhs =	1✔
129	case op of	1✔
130	Le -> Atom $ leZ lhs rhs	1✔
131	Ge -> Atom $ geZ lhs rhs	1✔
132	Lt -> Atom $ ltZ lhs rhs	1✔
133	Gt -> Atom $ gtZ lhs rhs	1✔
134	Eql -> lhs .==. rhs	1✔
135	NEq -> lhs ./=. rhs	×
136
137	-- \| @d \| e@ means @e@ can be divided by @d@.
138	(.\|.) :: Integer -> ExprZ -> QFFormula
139	n .\|. e = Atom $ Divisible True n e	1✔
140
141	subst1 :: Var -> ExprZ -> QFFormula -> QFFormula
142	subst1 x e = fmap f	1✔
143	where
144	f (Divisible b c e1) = Divisible b c $ LA.applySubst1 x e e1	1✔
145	f (IsPos e1) = IsPos $ LA.applySubst1 x e e1	1✔
146
147	simplify :: QFFormula -> QFFormula
148	simplify = BoolExpr.simplify . BoolExpr.fold simplifyLit	1✔
149
150	simplifyLit :: Lit -> QFFormula
151	simplifyLit (IsPos e) =	1✔
152	case LA.asConst e of	1✔
153	Just c -> if c > 0 then true else false	1✔
154	Nothing ->
155	-- e > 0 <=> e - 1 >= 0
156	-- <=> LA.mapCoeff (`div` d) (e - 1) >= 0
157	-- <=> LA.mapCoeff (`div` d) (e - 1) + 1 > 0
158	Atom $ IsPos $ LA.mapCoeff (`div` d) e1 ^+^ LA.constant 1	1✔
159	where
160	e1 = e ^-^ LA.constant 1	1✔
161	d = if null cs then 1 else abs $ foldl1' gcd cs	×
162	cs = [c \| (c,x) <- LA.terms e1, x /= LA.unitVar]	1✔
163	simplifyLit lit@(Divisible b c e)
164	\| LA.coeff LA.unitVar e2 `mod` d /= 0 = if b then false else true	×
165	\| c' == 1 = if b then true else false	×
166	\| d == 1 = Atom lit	1✔
167	\| otherwise = Atom $ Divisible b c' e'	×
168	where
169	e2 = LA.mapCoeff (`mod` c) e	1✔
170	d = abs $ foldl' gcd c [c2 \| (c2,x) <- LA.terms e2, x /= LA.unitVar]	1✔
171	c' = c `checkedDiv` d	1✔
172	e' = LA.mapCoeff (`checkedDiv` d) e2	1✔
173
174	instance Eval (Model Integer) Lit Bool where
175	eval m (IsPos e) = LA.eval m e > 0	1✔
176	eval m (Divisible True n e) = LA.eval m e `mod` n == 0	×
177	eval m (Divisible False n e) = LA.eval m e `mod` n /= 0	×
178
179	-- ---------------------------------------------------------------------------
180
181	data Witness = WCase1 Integer ExprZ \| WCase2 Integer Integer Integer (Set ExprZ)
182	deriving (Show)	×
183
184	instance Eval (Model Integer) Witness Integer where
185	eval model (WCase1 c e) = LA.eval model e `checkedDiv` c	1✔
186	eval model (WCase2 c j delta us)
187	\| Set.null us' = j `checkedDiv` c	1✔
188	\| otherwise = (j + (((u - delta - 1) `div` delta) * delta)) `checkedDiv` c	×
189	where
190	us' = Set.map (LA.eval model) us	1✔
191	u = Set.findMin us'	1✔
192
193	-- ---------------------------------------------------------------------------
194
195	{-\| @'project' x φ@ returns @(ψ, lift)@ such that:
196
197	* @⊢_LIA ∀y1, …, yn. (∃x. φ) ↔ ψ@ where @{y1, …, yn} = FV(φ) \\ {x}@, and
198
199	* if @M ⊧_LIA ψ@ then @lift M ⊧_LIA φ@.
200	-}
201	project :: Var -> QFFormula -> (QFFormula, Model Integer -> Model Integer)
202	project x formula = (formula', mt)	×
203	where
204	xs = projectCases x formula	×
205	formula' = orB' [phi \| (phi,_) <- xs]	×
206	mt m = head $ do	×
207	(phi, mt') <- xs	×
208	guard $ eval m phi	×
209	return $ mt' m	×
210	orB' = orB . concatMap f	×
211	where
212	f (Or xs) = concatMap f xs	×
213	f x = [x]	×
214
215	{-\| @'projectN' {x1,…,xm} φ@ returns @(ψ, lift)@ such that:
216
217	* @⊢_LIA ∀y1, …, yn. (∃x1, …, xm. φ) ↔ ψ@ where @{y1, …, yn} = FV(φ) \\ {x1,…,xm}@, and
218
219	* if @M ⊧_LIA ψ@ then @lift M ⊧_LIA φ@.
220	-}
221	projectN :: VarSet -> QFFormula -> (QFFormula, Model Integer -> Model Integer)
222	projectN vs2 = f (IS.toList vs2)	×
223	where
224	f :: [Var] -> QFFormula -> (QFFormula, Model Integer -> Model Integer)
225	f [] formula = (formula, id)	×
226	f (v:vs) formula = (formula3, mt1 . mt2)	×
227	where
228	(formula2, mt1) = project v formula	×
229	(formula3, mt2) = f vs formula2	×
230
231	{-\| @'projectCases' x φ@ returns @[(ψ_1, lift_1), …, (ψ_m, lift_m)]@ such that:
232
233	* @⊢_LIA ∀y1, …, yn. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_m)@ where @{y1, …, yn} = FV(φ) \\ {x}@, and
234
235	* if @M ⊧_LIA ψ_i@ then @lift_i M ⊧_LIA φ@.
236	-}
237	projectCases :: Var -> QFFormula -> [(QFFormula, Model Integer -> Model Integer)]
238	projectCases x formula = do	1✔
239	(phi, wit) <- projectCases' x formula	1✔
240	return (phi, \m -> IM.insert x (eval m wit) m)	1✔
241
242	projectCases' :: Var -> QFFormula -> [(QFFormula, Witness)]
243	projectCases' x formula = [(phi', w) \| (phi,w) <- case1 ++ case2, let phi' = simplify phi, phi' /= false]	1✔
244	where
245	-- eliminate Not, Imply and Equiv.
246	formula0 :: QFFormula
247	formula0 = pos formula	1✔
248	where
249	pos (Atom a) = Atom a	1✔
250	pos (And xs) = And (map pos xs)	1✔
251	pos (Or xs) = Or (map pos xs)	×
252	pos (Not x) = neg x	×
253	pos (Imply x y) = neg x .\|\|. pos y	×
254	pos (Equiv x y) = pos ((x .=>. y) .&&. (y .=>. x))	×
255	pos (ITE c t e) = pos ((c .=>. t) .&&. (Not c .=>. e))	×
256
257	neg (Atom a) = Atom (notB a)	×
258	neg (And xs) = Or (map neg xs)	×
259	neg (Or xs) = And (map neg xs)	×
260	neg (Not x) = pos x	×
261	neg (Imply x y) = pos x .&&. neg y	×
262	neg (Equiv x y) = neg ((x .=>. y) .&&. (y .=>. x))	×
263	neg (ITE c t e) = neg ((c .=>. t) .&&. (Not c .=>. e))	×
264
265	-- xの係数の最小公倍数
266	c :: Integer
267	c = getLCM $ Foldable.foldMap f formula0	1✔
268	where
269	f (IsPos e) = LCM $ fromMaybe 1 (LA.lookupCoeff x e)	1✔
270	f (Divisible _ _ e) = LCM $ fromMaybe 1 (LA.lookupCoeff x e)	×
271
272	-- 式をスケールしてxの係数の絶対値をcへと変換し、その後cxをxで置き換え、
273	-- xがcで割り切れるという制約を追加した論理式
274	formula1 :: QFFormula
275	formula1 = simplify $ fmap f formula0 .&&. (c .\|. LA.var x)	1✔
276	where
277	f lit@(IsPos e) =	1✔
278	case LA.lookupCoeff x e of	1✔
279	Nothing -> lit	1✔
280	Just a ->
281	let s = abs (c `checkedDiv` a)	1✔
282	in IsPos $ g s e	1✔
283	f lit@(Divisible b d e) =
284	case LA.lookupCoeff x e of	1✔
285	Nothing -> lit	×
286	Just a ->
287	let s = abs (c `checkedDiv` a)	1✔
288	in Divisible b (s*d) $ g s e	1✔
289
290	g :: Integer -> ExprZ -> ExprZ
291	g s = LA.mapCoeffWithVar (\c' x' -> if x==x' then signum c' else s*c')	1✔
292
293	-- d\|x+t という形の論理式の d の最小公倍数
294	delta :: Integer
295	delta = getLCM $ Foldable.foldMap f formula1	1✔
296	where
297	f (Divisible _ d _) = LCM d	1✔
298	f (IsPos _) = LCM 1	1✔
299
300	-- ts = {t \| t < x は formula1 に現れる原子論理式}
301	ts :: Set ExprZ
302	ts = Foldable.foldMap f formula1	1✔
303	where
304	f (Divisible _ _ _) = Set.empty	1✔
305	f (IsPos e) =
306	case LA.extractMaybe x e of	1✔
307	Nothing -> Set.empty	1✔
308	Just (1, e') -> Set.singleton (negateV e') -- IsPos e <=> (x + e' > 0) <=> (-e' < x)	1✔
309	Just (-1, _) -> Set.empty -- IsPos e <=> (-x + e' > 0) <=> (x < e')	1✔
310	_ -> error "should not happen"	×
311
312	-- formula1を真にする最小のxが存在する場合
313	case1 :: [(QFFormula, Witness)]
314	case1 = [ (subst1 x e formula1, WCase1 c e)	1✔
315	\| j <- [1..delta], t <- Set.toList ts, let e = t ^+^ LA.constant j ]	1✔
316
317	-- formula1のなかの x < t を真に t < x を偽に置き換えた論理式
318	formula2 :: QFFormula
319	formula2 = simplify $ BoolExpr.fold f formula1	1✔
320	where
321	f lit@(IsPos e) =	1✔
322	case LA.lookupCoeff x e of	1✔
323	Nothing -> Atom lit	1✔
324	Just 1 -> false -- IsPos e <=> ( x + e' > 0) <=> -e' < x	1✔
325	Just (-1) -> true -- IsPos e <=> (-x + e' > 0) <=> x < e'	1✔
326	_ -> error "should not happen"	×
327	f lit@(Divisible _ _ _) = Atom lit	1✔
328
329	-- us = {u \| x < u は formula1 に現れる原子論理式}
330	us :: Set ExprZ
331	us = Foldable.foldMap f formula1	1✔
332	where
333	f (IsPos e) =	1✔
334	case LA.extractMaybe x e of	1✔
335	Nothing -> Set.empty	1✔
336	Just (1, _) -> Set.empty -- IsPos e <=> (x + e' > 0) <=> -e' < x	×
337	Just (-1, e') -> Set.singleton e' -- IsPos e <=> (-x + e' > 0) <=> x < e'	1✔
338	_ -> error "should not happen"	×
339	f (Divisible _ _ _) = Set.empty	1✔
340
341	-- formula1を真にする最小のxが存在しない場合
342	case2 :: [(QFFormula, Witness)]
343	case2 = [(subst1 x (LA.constant j) formula2, WCase2 c j delta us) \| j <- [1..delta]]	1✔
344
345	{-\| @'projectCasesN' {x1,…,xm} φ@ returns @[(ψ_1, lift_1), …, (ψ_n, lift_n)]@ such that:
346
347	* @⊢_LIA ∀y1, …, yp. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_n)@ where @{y1, …, yp} = FV(φ) \\ {x1,…,xm}@, and
348
349	* if @M ⊧_LIA ψ_i@ then @lift_i M ⊧_LIA φ@.
350	-}
351	projectCasesN :: VarSet -> QFFormula -> [(QFFormula, Model Integer -> Model Integer)]
352	projectCasesN vs2 = f (IS.toList vs2)	1✔
353	where
354	f :: [Var] -> QFFormula -> [(QFFormula, Model Integer -> Model Integer)]
355	f [] formula = return (formula, id)	1✔
356	f (v:vs) formula = do	1✔
357	(formula2, mt1) <- projectCases v formula	1✔
358	(formula3, mt2) <- f vs formula2	1✔
359	return (formula3, mt1 . mt2)	1✔
360
361	-- ---------------------------------------------------------------------------
362
363	newtype LCM a = LCM{ getLCM :: a }	1✔
364
365	instance Integral a => Semigroup.Semigroup (LCM a) where	×
366	LCM a <> LCM b = LCM $ lcm a b	1✔
367	stimes = Semigroup.stimesIdempotent	×
368
369	instance Integral a => Monoid (LCM a) where	×
370	mempty = LCM 1	1✔
371
372	checkedDiv :: Integer -> Integer -> Integer
373	checkedDiv a b =	1✔
374	case a `divMod` b of	1✔
375	(q,0) -> q	1✔
376	_ -> error "ToySolver.Cooper.checkedDiv: should not happen"	×
377
378	-- ---------------------------------------------------------------------------
379
380	-- \| @'solveQFFormula' {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LIA φ@ when
381	-- such @M@ exists, returns @Nothing@ otherwise.
382	--
383	-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
384	--
385	solveQFFormula :: VarSet -> QFFormula -> Maybe (Model Integer)
386	solveQFFormula vs formula = listToMaybe $ do	1✔
387	(formula2, mt) <- projectCasesN vs formula	1✔
388	let m :: Model Integer
389	m = IM.empty	1✔
390	guard $ eval m formula2	1✔
391	return $ mt m	1✔
392
393	-- \| @'solve' {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LIA φ@ when
394	-- such @M@ exists, returns @Nothing@ otherwise.
395	--
396	-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
397	--
398	solve :: VarSet -> [LA.Atom Rational] -> Maybe (Model Integer)
399	solve vs cs = solveQFFormula vs $ andB $ map fromLAAtom cs	1✔
400
401	-- \| @'solveQFLIRAConj' {x1,…,xn} φ I@ returns @Just M@ that @M ⊧_LIRA φ@ when
402	-- such @M@ exists, returns @Nothing@ otherwise.
403	--
404	-- * @FV(φ)@ must be a subset of @{x1,…,xn}@.
405	--
406	-- * @I@ is a set of integer variables and must be a subset of @{x1,…,xn}@.
407	--
408	solveQFLIRAConj :: VarSet -> [LA.Atom Rational] -> VarSet -> Maybe (Model Rational)
409	solveQFLIRAConj vs cs ivs = listToMaybe $ do	×
410	(cs2, mt) <- FM.projectN rvs cs	×
411	m <- maybeToList $ solve ivs cs2	×
412	return $ mt $ IM.map fromInteger m	×
413	where
414	rvs = vs `IS.difference` ivs	×
415
416	-- ---------------------------------------------------------------------------
417
418	_testHagiya :: QFFormula
419	_testHagiya = fst $ project 1 $ andB [c1, c2, c3]	×
420	where
421	[x,y,z] = map LA.var [1..3]	×
422	c1 = x .<. (y ^+^ y)	×
423	c2 = z .<. x	×
424	c3 = 3 .\|. x	×
425
426	{-
427	∃ x. x < y+y ∧ z<x ∧ 3\|x
428	⇔
429	(2y>z+1 ∧ 3\|z+1) ∨ (2y>z+2 ∧ 3\|z+2) ∨ (2y>z+3 ∧ 3\|z+3)
430	-}
431
432	_test3 :: QFFormula
433	_test3 = fst $ project 1 $ andB [p1,p2,p3,p4]	×
434	where
435	x = LA.var 0	×
436	y = LA.var 1	×
437	p1 = LA.constant 0 .<. y	×
438	p2 = 2 *^ x .<. y	×
439	p3 = y .<. x ^+^ LA.constant 2	×
440	p4 = 2 .\|. y	×
441
442	{-
443	∃ y. 0 < y ∧ 2x<y ∧ y < x+2 ∧ 2\|y
444	⇔
445	(2x < 2 ∧ 0 < x) ∨ (0 < 2x+2 ∧ x < 0)
446	-}
447
448	-- ---------------------------------------------------------------------------

msakai / toysolver / 604

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