test: new linear solver benchmark by Marc

This commit is contained in:
Sebastian Ullrich 2021-12-02 17:03:35 +01:00
parent b6d5cd8155
commit 5869283694
4 changed files with 442 additions and 0 deletions

View file

@ -0,0 +1,51 @@
50 50
27 -144 2 -281 4 172 7 -284 8 16 9 -158 12 -220 13 25 15 -125 16 -123 19 -53 21 -260 22 0 23 -184 28 -251 29 89 31 -24 32 -27 37 -290 40 -83 44 -258 45 253 46 -242 47 -198 48 -2 49 244 50 204409 0
27 -139 2 134 5 -7 6 3 7 256 8 110 9 101 10 -125 12 -275 14 -75 16 142 20 186 21 -295 25 -296 27 -154 28 -107 33 -153 34 -17 35 -137 36 0 39 24 41 -163 42 289 43 -128 44 -251 46 156 47 -88912 0
31 -84 4 270 5 227 6 -24 7 -78 8 -155 9 -135 10 -196 11 -174 14 280 16 -2 17 166 18 143 21 278 24 223 25 -237 26 110 28 -175 29 24 31 -220 32 86 34 -94 36 10 37 73 38 82 39 -189 42 258 43 27 45 -221 48 -30 50 80123 0
20 -77 1 -221 3 94 6 217 16 -9 20 139 21 -63 26 -25 28 -253 29 -26 33 42 35 90 38 -61 40 -40 41 156 44 -248 46 121 48 114 49 -18 50 -97596 0
28 55 1 255 2 261 3 135 4 172 10 79 13 -83 14 216 15 -74 17 65 18 141 19 -1 20 192 21 -32 23 -294 24 266 25 118 27 50 29 287 32 180 33 14 34 281 37 222 38 217 40 208 46 40 47 49 49 -31347 0
33 -11 1 -300 3 243 4 -213 5 247 6 172 7 -100 8 -179 9 19 10 197 11 65 13 206 15 11 18 292 19 -50 22 22 23 212 24 -86 25 -219 27 -28 28 274 31 -14 32 214 33 46 37 193 38 -251 40 -174 43 82 44 112 45 -224 48 -33 49 -255 50 -72266 0
24 -122 3 -49 5 -123 6 31 10 135 11 -181 15 -185 18 -107 19 291 23 -121 24 -220 25 17 27 171 30 -55 31 93 35 -217 36 -72 37 -48 38 -130 39 154 42 -257 44 -279 45 157 47 -24256 0
25 10 2 -275 4 -197 6 198 8 -288 9 50 13 155 16 264 18 -292 21 -142 23 -128 24 52 25 23 28 85 29 13 30 -74 33 -5 34 -168 35 80 37 -182 38 -284 40 -225 42 76 43 -203 45 2432 0
30 -123 1 -140 3 -87 4 170 5 277 6 52 7 100 8 -66 9 -234 10 -168 11 -84 12 -209 13 169 14 -155 15 -78 16 -2 17 253 18 -209 19 201 20 163 21 264 23 17 24 125 25 247 26 162 27 -201 35 -36 37 -5 47 133 50 178291 0
27 -19 1 140 3 -120 6 -42 7 177 8 44 9 257 12 -215 15 8 16 -258 17 -214 18 -155 19 87 20 -233 21 178 25 35 27 -74 29 -141 30 236 32 203 33 -200 35 239 37 176 39 43 41 -276 48 -103 49 20229 0
28 36 1 136 6 -201 7 155 8 -300 9 95 11 89 12 94 13 -21 15 172 16 20 18 29 19 152 20 -53 23 -165 24 24 25 222 26 135 27 53 29 -5 32 193 34 196 36 38 39 94 40 -246 41 -34 43 4 44 -10827 0
27 238 2 278 3 139 6 -140 7 126 8 263 9 243 11 -150 12 -87 13 180 16 -151 17 207 20 264 21 133 22 -204 23 121 29 -299 34 259 35 -52 38 95 39 -251 40 35 41 -127 43 263 45 -32 46 -226 50 -298973 0
26 -101 4 80 5 12 7 -84 9 55 12 259 13 116 15 165 18 -43 19 95 21 121 22 -141 25 -261 27 45 30 151 31 211 35 -285 37 -82 38 -14 39 255 40 -56 41 67 42 -210 45 -296 46 -267 49 -58757 0
28 -267 1 122 3 142 4 -213 5 39 6 274 8 -116 9 -72 13 -250 14 -176 16 -162 17 -110 19 -123 21 -123 22 -102 24 -205 28 128 29 -184 30 198 33 115 35 -213 36 -245 37 18 40 -54 42 -73 45 7 46 -144 50 66742 0
18 5 7 142 9 -23 11 289 12 230 14 203 16 0 18 180 19 32 20 287 31 -158 32 123 36 259 40 9 42 193 47 202 48 -175 49 49417 0
31 54 1 -160 3 30 4 -5 6 -246 8 -102 9 -229 10 -66 11 -238 12 30 14 -88 15 59 17 106 18 195 22 -270 25 28 26 -262 27 258 29 62 31 188 34 178 35 -38 36 170 37 -50 38 140 39 176 42 216 44 55 45 -122 47 84 50 84718 0
26 -134 3 279 5 193 6 -75 10 195 11 56 12 5 13 72 14 -46 15 -42 16 -167 17 -258 18 119 19 143 20 242 29 -144 30 -281 32 -104 33 -15 35 211 37 205 39 141 40 -19 43 274 47 -268 49 90385 0
27 -298 5 -9 10 -38 11 -47 13 -4 14 -274 15 104 16 118 17 166 19 253 20 -121 21 -163 23 12 24 296 29 -137 30 59 31 143 33 206 34 -191 35 -88 37 -180 39 174 40 192 44 96 46 -283 47 250 50 123002 0
29 218 1 -76 6 45 8 -275 9 58 10 59 13 -236 14 185 15 241 16 282 17 234 18 -29 19 64 22 271 24 250 25 172 29 -244 30 268 33 -209 37 -212 38 -102 39 217 40 -162 43 -261 44 153 45 41 46 -234 48 -214 50 -139007 0
26 295 2 71 3 -251 4 -235 6 69 8 299 10 -198 15 -179 19 260 21 -106 23 -230 24 258 26 156 27 -213 29 -113 35 267 36 -214 38 15 39 -111 40 -276 42 -36 46 -8 47 82 48 -230 49 -45 50 -84778 0
30 -58 2 120 3 -258 4 65 5 -230 6 -154 7 -267 8 -216 11 -297 13 96 16 156 17 -121 19 -64 20 -264 21 292 26 -243 28 -48 29 -152 31 141 34 39 35 -16 36 269 37 -191 41 -160 42 -250 44 110 45 -16 46 167 48 203 49 -17836 0
27 -231 1 50 3 66 4 6 5 -111 6 30 8 -150 12 286 14 -261 16 21 17 97 18 -70 19 247 23 25 24 -1 27 -127 32 -107 34 -162 35 185 37 183 39 160 44 -241 45 147 47 -157 48 256 49 230 50 8609 0
29 278 1 -172 2 44 7 -285 8 284 11 -243 13 22 14 -132 19 50 22 147 23 198 24 -135 25 -157 27 -82 32 190 34 -150 35 144 36 -206 39 33 40 -89 41 -299 42 95 43 -111 44 174 45 110 46 -150 48 294 49 175 50 182491 0
28 -269 1 188 4 -45 6 15 7 120 9 146 10 -274 11 267 15 73 17 -221 20 -69 21 178 22 -7 23 -161 26 -147 27 -250 28 -210 29 -259 33 -71 34 287 36 -72 37 -187 41 295 43 135 44 257 46 -236 49 -221 50 186243 0
23 -291 1 284 2 -105 4 -237 5 -171 9 131 11 232 13 247 16 27 17 22 23 67 24 -155 26 19 27 55 28 257 29 -49 33 244 34 203 38 157 40 -85 44 188 45 72 49 -49835 0
26 179 1 83 3 -8 4 42 8 -154 10 -127 12 193 13 133 14 -51 15 74 18 175 23 -109 24 -230 28 -81 29 274 31 187 32 103 34 -164 37 -105 38 -21 39 22 43 290 45 -164 46 -259 49 -254 50 -1705 0
24 111 1 122 2 142 5 -16 11 132 12 -8 14 243 17 -217 19 -121 23 115 25 -298 26 -186 27 -121 28 -186 29 197 31 242 35 2 41 276 42 -161 45 -119 46 239 47 285 49 -41 50 -289044 0
24 -256 1 -4 4 -277 5 -234 6 5 10 -8 13 -63 15 160 16 -178 21 -125 23 291 24 -295 25 238 28 -262 29 -222 30 6 31 -275 33 286 34 -209 38 -36 40 47 42 186 45 -53 46 277639 0
25 -150 1 282 4 102 6 -61 8 -92 9 -47 10 -82 11 34 12 -257 16 77 20 226 21 78 22 -24 25 150 27 -12 28 261 29 -38 31 203 39 -134 40 -255 42 -244 43 104 45 -110 46 -221 49 50648 0
26 -69 2 -187 4 -4 8 -181 10 -222 11 -295 13 -297 18 -55 20 113 21 -157 22 -108 24 136 25 -102 27 -41 28 -13 30 -272 31 -75 32 40 33 -175 34 131 38 102 39 137 43 -141 45 287 47 147 49 -131506 0
25 -235 1 -139 2 -65 3 193 5 -34 6 -124 8 28 12 204 14 -136 15 160 16 221 21 291 24 10 26 -108 30 -222 33 187 34 54 38 109 39 -169 41 -73 42 191 44 113 46 294 48 -57 50 49727 0
26 -48 4 152 5 75 7 -92 9 237 11 -4 15 -27 16 -262 18 -3 19 -247 21 212 22 -184 24 285 25 193 26 -195 30 218 31 86 32 198 35 208 39 -81 43 252 44 -7 45 -149 48 -64 49 89 50 41586 0
22 -169 3 -294 6 23 8 -246 15 -261 17 -198 20 -228 22 -106 23 -57 24 -155 26 250 28 -229 30 -228 32 235 33 143 34 238 35 -284 39 96 40 128 45 -172 46 -257 47 -150791 0
26 -290 1 253 3 4 5 192 6 -211 8 -78 9 -248 11 -13 12 -230 15 140 16 6 19 -136 20 -287 22 -281 24 -57 27 -272 28 -111 29 205 36 -179 37 -112 39 33 42 30 45 -260 47 -68 48 68 49 -39862 0
29 153 2 72 3 -18 4 128 5 -145 6 122 7 79 8 258 9 205 12 259 15 -242 16 -284 17 -36 18 132 19 194 20 228 21 -190 24 178 26 -46 31 227 32 293 33 156 34 -198 36 -130 38 186 41 -46 47 114 48 -47 49 -89775 0
31 180 1 -22 2 85 4 -218 6 -209 9 -244 10 54 11 282 12 -158 14 285 15 -202 16 -158 18 -61 19 73 21 -72 22 279 23 38 24 107 25 -112 26 -297 28 -15 29 -227 30 -276 35 100 37 -233 40 -248 43 236 44 -76 47 52 48 203 49 -103843 0
30 -100 2 -68 3 7 4 -104 6 -87 7 -276 8 53 9 19 11 -233 12 -294 14 110 15 -246 16 45 17 78 18 25 21 95 24 27 25 -150 28 247 33 215 35 200 36 86 37 16 40 -233 43 117 46 -157 47 185 48 22 49 -230 50 13768 0
21 137 4 -294 6 108 8 65 9 125 11 -203 14 -210 15 -20 19 25 21 -152 22 -238 24 -9 28 96 29 -276 30 -176 32 271 34 -16 40 244 41 -115 43 -226 50 -29238 0
32 268 1 -87 2 -191 4 -71 5 -147 9 14 12 159 13 23 16 251 18 -6 19 217 20 -277 21 75 23 29 24 13 25 -110 26 24 27 28 28 -132 29 -231 30 218 31 251 32 92 33 122 38 204 39 74 40 3 43 161 45 -250 47 7 48 -36 49 39470 0
30 -288 3 190 4 -128 5 -108 7 -239 9 13 10 -64 11 102 12 167 16 -128 17 -138 18 177 19 263 20 201 23 161 26 181 27 95 28 -178 30 -125 31 192 33 -156 34 -205 36 -168 38 -128 40 -274 41 -200 42 -190 44 -99 47 99 49 -100595 0
28 -92 2 -250 3 -13 4 214 5 -179 7 -42 10 249 11 -231 12 -149 14 -28 15 273 16 165 22 -107 25 -48 28 -128 29 -293 32 -35 34 103 35 -168 36 -98 38 -26 40 217 41 62 42 -240 44 53 46 -155 47 -187 48 20275 0
26 117 2 51 5 -237 6 93 8 58 9 -273 10 143 15 -178 17 58 18 -197 21 224 27 124 28 -243 30 -85 31 -114 32 196 33 185 34 -165 35 -287 36 -131 39 30 41 -236 43 108 44 142 45 -13 47 110196 0
26 -148 1 106 2 221 6 -114 10 209 11 184 14 17 16 212 17 198 20 158 22 139 24 -99 26 -198 27 226 28 136 31 -159 32 -273 35 262 37 39 39 107 40 -210 42 -55 43 67 44 -269 47 69 50 98508 0
33 -206 1 247 2 7 3 129 6 -249 7 164 8 184 9 159 10 29 13 211 14 -261 16 -162 17 189 20 49 21 123 22 111 23 -136 24 -156 26 -148 27 234 28 -272 29 284 31 249 35 222 36 -115 37 -190 39 -130 40 -263 41 99 45 30 46 73 48 227 49 -190958 0
26 119 1 -116 2 -125 3 37 4 -86 5 181 7 130 11 -141 13 205 14 -202 15 283 16 -203 17 -7 18 -71 20 -293 22 -6 27 52 29 -215 31 21 37 99 38 -191 39 231 41 -282 43 262 44 -99 46 -98035 0
29 -297 1 -28 4 -92 6 -88 8 -38 10 173 12 115 13 -24 14 67 16 237 18 -121 19 101 21 229 23 -31 25 19 27 185 29 76 33 -155 34 -113 37 -32 38 146 39 178 42 -47 44 15 46 256 47 188 48 -91 49 -228 50 -85530 0
24 -22 1 -51 3 -43 8 26 13 -20 14 -25 15 -293 19 237 20 -9 21 -49 23 65 26 95 27 204 28 -113 29 -119 31 164 33 -72 36 300 39 121 40 27 41 -277 42 -230 47 -131 48 11135 0
21 190 1 144 2 182 6 175 8 -1 9 63 11 -86 13 152 14 -48 15 -274 16 -298 20 -109 30 79 32 77 35 -185 37 116 38 -181 41 275 42 -89 45 -211 46 -111362 0
27 -225 3 -218 5 -122 7 -8 12 -32 14 169 20 -266 24 42 26 -48 27 -233 28 -159 31 282 33 -187 34 -10 35 103 37 103 39 221 40 -21 41 101 42 130 43 -24 44 113 45 206 46 -148 47 -268 49 -259 50 -54796 0
23 -222 1 109 2 202 8 152 10 146 12 -285 14 183 15 -5 21 -8 22 66 25 -139 27 93 28 -265 32 -4 37 -261 38 267 40 -47 41 -171 42 -189 43 9 44 202 45 174 46 38730 0

388
tests/bench/liasolver.lean Normal file
View file

@ -0,0 +1,388 @@
/-
Linear Diophantine equation solver
Author: Marc Huisinga
-/
import Std.Data.HashMap
open Std
namespace Int
def roundedDiv (a b : Int) : Int := do
let mut div := a / b
let sgn := if a ≥ 0 ∧ b ≥ 0 a < 0 ∧ b < 0 then 1 else -1
let rest := (a % b).natAbs
if 2*rest ≥ b.natAbs then
div := div + sgn
return div
def mod' (a b : Int) : Int :=
a - b*(a.roundedDiv b)
end Int
namespace Std.AssocList
def map (f : α → β → δ) : AssocList α β → AssocList α δ
| AssocList.nil => AssocList.nil
| AssocList.cons k v t => AssocList.cons k (f k v) (map f t)
def filter (p : α → β → Bool) : AssocList α β → AssocList α β
| AssocList.nil => AssocList.nil
| AssocList.cons k v t =>
if p k v then
AssocList.cons k v (filter p t)
else
filter p t
end Std.AssocList
namespace Std.HashMap
variable [BEq α] [Hashable α]
@[inline] protected def forIn {δ : Type w} {m : Type w → Type w'} [Monad m]
(as : HashMap α β) (init : δ) (f : (α × β) → δ → m (ForInStep δ)) : m δ := do
as.val.buckets.val.forIn init fun bucket acc => do
let (done, v) ← bucket.forIn (false, acc) fun v (_, acc) => do
let r ← f v acc
match r with
| ForInStep.done r' =>
return ForInStep.done (true, r')
| ForInStep.yield r' =>
return ForInStep.yield (false, r')
if done then
return ForInStep.done v
else
return ForInStep.yield v
instance : ForIn m (HashMap α β) (α × β) where
forIn := HashMap.forIn
def modify! [Inhabited β] (xs : HashMap α β) (k : α) (f : β → β) : HashMap α β :=
let v := xs.find! k
xs.erase k |>.insert k (f v)
def any (xs : HashMap α β) (p : α → β → Bool) : Bool := do
for (k, v) in xs do
if p k v then
return true
return false
def mapValsM [Monad m] (f : β → m γ) (xs : HashMap α β) : m (HashMap α γ) :=
mkHashMap (capacity := xs.size) |> xs.foldM fun acc k v => do acc.insert k (←f v)
def mapVals (f : β → γ) (xs : HashMap α β) : HashMap α γ :=
mkHashMap (capacity := xs.size) |> xs.fold fun acc k v => acc.insert k (f v)
def fastMapVals (f : α → β → β) (xs : HashMap α β) : HashMap α β :=
let size := xs.val.size
let buckets := xs.val.buckets.val.map (·.map f)
⟨⟨size, ⟨buckets, sorry⟩⟩, sorry⟩
def filter (p : α → β → Bool) (xs : HashMap α β) : HashMap α β :=
let buckets := xs.val.buckets.val.map (·.filter p)
let size := buckets.foldl (fun acc bucket => bucket.foldl (fun acc _ _ => acc + 1) acc) 0
⟨⟨size, ⟨buckets, sorry⟩⟩, sorry⟩
def getAny? (x : HashMap α β) : Option (α × β) := do
for bucket in x.val.buckets.val do
for (k, v) in bucket do
return some (k, v)
return none
end Std.HashMap
structure Equation where
id : Nat
coeffs : HashMap Nat Int
const : Int
deriving Inhabited
def gcd (coeffs : HashMap Nat Int) : Nat :=
let coeffs := coeffs.mapVals (·.natAbs)
let coeffsContent := coeffs.toArray
match coeffsContent with
| #[] => panic! "Cannot calculate GCD of empty list of coefficients"
| #[(i, x)] => x
| coeffsContent =>
coeffsContent[0].2.gcd coeffsContent[1].2
|> coeffs.fold fun acc k v => acc.gcd v
namespace Equation
def preprocess? (e : Equation) : Option Equation := do
let gcd : Int := gcd e.coeffs
if e.const % gcd ≠ 0 then
return none
return some { e with
coeffs := e.coeffs.fastMapVals fun _ coeff => coeff / gcd
const := e.const / gcd }
def subst (fromEq toEq : Equation) (varIdx : Nat) : Equation := do
-- varIdx ≡ k
-- fromEq ≡ sₖxₖ + ∑ i ∈ V_fromEq\{k}. aᵢxᵢ = A
-- ⇔ xₖ = sₖA - ∑ i ∈ V_fromEq\{k}. sₖaᵢxᵢ
-- toEq ≡ bₖxₖ + ∑ i ∈ V_toEq\{k}. bᵢxᵢ = B
-- ⇝ B = bₖ(sₖA - ∑ i ∈ V_fromEq\{k}. sₖaᵢxᵢ) + ∑ i ∈ V_toEq\{k}. bᵢxᵢ
-- = bₖsₖA - ∑ i ∈ V_fromEq\{k}. sₖbₖaᵢxᵢ + ∑ i ∈ V_toEq\{k}. bᵢxᵢ
-- = bₖsₖA + ∑ i ∈ V_fromEq\V_toEq. -sₖbₖaᵢxᵢ
-- + ∑ i ∈ V_toEq\{k} ∩ V_fromEq. (bᵢ - sₖbₖaᵢ)xᵢ
-- + ∑ i ∈ V_toEq\V_fromEq. bᵢxᵢ
-- ⇔ B - bₖsₖA = + ∑ i ∈ X. -sₖbₖaᵢxᵢ
-- + ∑ i ∈ Y. (bᵢ - sₖbₖaᵢ)xᵢ
-- + ∑ i ∈ Z. bᵢxᵢ
-- with X, Y, Z defined as above, X Y Z = (V_fromEq V_toEq)\{k}
-- and X, Y, Z pairwise disjoint
let A := fromEq.const
let B := toEq.const
let V_fromEq := fromEq.coeffs
let V_toEq := toEq.coeffs
let k := varIdx
let sₖ := V_fromEq.find! k
let bₖ := V_toEq.find! k
let mut V_toEq := V_toEq.fastMapVals fun i bᵢ =>
match V_fromEq.find? i with
| none =>
bᵢ
| some aᵢ =>
bᵢ - sₖ*bₖ*aᵢ
for (i, aᵢ) in V_fromEq do
if ¬V_toEq.contains i then
V_toEq := V_toEq.insert i (-sₖ*bₖ*aᵢ)
V_toEq := V_toEq.filter fun i bᵢ => i ≠ k ∧ bᵢ ≠ 0
let B' := B - bₖ*sₖ*A
{ toEq with coeffs := V_toEq, const := B' }
def normalize (e : Equation) : Equation := do
if e.coeffs.size ≠ 1 then
return e
let (i, c) := e.coeffs.getAny?.get!
return { e with
coeffs := e.coeffs.insert i 1
const := e.const / c }
def invert (e : Equation) : Equation :=
{ e with
coeffs := e.coeffs.fastMapVals fun _ coeff => (-1)*coeff
const := (-1)*e.const }
def reorganizeFor (e : Equation) (varIdx : Nat) : Equation := do
let singletonCoeff := e.coeffs.find! varIdx
let mut e := { e with coeffs := e.coeffs.fastMapVals fun _ coeff => (-1)*coeff }
if singletonCoeff = -1 then
e := e.invert
{ e with coeffs := e.coeffs.erase varIdx }
def findSingleton? (e : Equation) : Option (Nat × Int) := do
for (i, coeff) in e.coeffs do
if coeff = 1 coeff = -1 then
return some (i, coeff)
return none
def findAbsMinimumCoeff? (e : Equation) : Option (Nat × Int) := do
let mut r? : Option (Nat × Int) := none
for (i, coeff) in e.coeffs do
match r? with
| none =>
r? := some (i, coeff)
| some (i', coeff') =>
if coeff.natAbs < coeff'.natAbs then
r? := some (i, coeff)
return r?
end Equation
structure Problem where
equations : HashMap Nat Equation
solvedEquations : HashMap Nat Equation
nEquations : Nat
nVars : Nat
deriving Inhabited
def preprocess? (eqs : HashMap Nat Equation) : Option (HashMap Nat Equation) :=
OptionM.run <| eqs.mapValsM (·.preprocess?)
def eliminateSingleton (p : Problem) (singletonEq : Equation) (varIdx : Nat) : Problem := do
let mut eqsWithVarIdx : Array Nat := #[]
for (id, eq) in p.equations do
if eq.coeffs.contains varIdx then
eqsWithVarIdx := eqsWithVarIdx.push id
let mut equations := p.equations
for id in eqsWithVarIdx do
if id == singletonEq.id then
continue
equations := equations.modify! id fun eq => singletonEq.subst eq varIdx |>.normalize
equations := equations.erase singletonEq.id
let solvedEquations := p.solvedEquations.insert varIdx <| singletonEq.reorganizeFor varIdx
return { p with
equations := equations
solvedEquations := solvedEquations }
partial def eliminateSingletons (p : Problem) : Problem := do
let mut r? : Option (Equation × Nat) := none
for (id, eq) in p.equations do
match eq.findSingleton? with
| none =>
continue
| some (varIdx, coeff) =>
r? := some (eq, varIdx)
match r? with
| none =>
return p
| some (eq, varIdx) =>
let p := eliminateSingleton p eq varIdx
return eliminateSingletons p
def addAuxEquation (p : Problem) : Problem := do
let mut E? : Option Equation := none
let mut k? : Option Nat := none
let mut aₖ? : Option Int := none
for (_, eq) in p.equations do
match eq.findAbsMinimumCoeff?, aₖ? with
| none, _ => continue
| some (k', aₖ'), none =>
E? := some eq
k? := some k'
aₖ? := some aₖ'
| some (k', aₖ'), some aₖ =>
if aₖ'.natAbs < aₖ.natAbs then
E? := some eq
k? := some k'
aₖ? := some aₖ'
let mut E := E?.get!
let k := k?.get!
let mut aₖ := aₖ?.get!
if aₖ < 0 then
aₖ := -aₖ
E := E.invert
let m := aₖ + 1
let σIdx := p.nVars
let newEqCoeffs := E.coeffs.fastMapVals (fun _ coeff => coeff.mod' m)
|>.insert σIdx (-m)
|>.filter (fun _ coeff => coeff ≠ 0)
let newEqConst := E.const.mod' m
let newEq : Equation := ⟨p.nEquations, newEqCoeffs, newEqConst⟩
let E'coeffs := E.coeffs.filter (fun i _ => i ≠ k)
|>.fastMapVals (fun _ aᵢ => aᵢ.roundedDiv m + aᵢ.mod' m)
|>.insert σIdx (-aₖ)
|>.filter (fun _ coeff => coeff ≠ 0)
let c := E.const
let E'const := c.roundedDiv m + c.mod' m
let E' := { E with coeffs := E'coeffs, const := E'const }.normalize
let equations' := p.equations.insert E'.id E' |>.insert newEq.id newEq
let p' : Problem := { p with
equations := equations'
nVars := p.nVars + 1
nEquations := p.nEquations + 1 }
return eliminateSingleton p' newEq k
inductive Solution
| unsat
| sat (assignment : Array Int)
deriving Inhabited
partial def readSolution? (p : Problem) : Option Solution := do
if p.equations.any (fun _ eq => eq.coeffs.size ≠ 0) then
return none
if p.equations.any (fun _ eq => eq.const ≠ 0) then
return some Solution.unsat
let mut assignment : Array (Option Int) := mkArray p.nVars none
for i in [0:p.nVars] do
assignment := readSolution i assignment
return Solution.sat <| assignment.map (·.get!)
where
readSolution (varIdx : Nat) (assignment : Array (Option Int)) : Array (Option Int) := do
match p.solvedEquations.find? varIdx with
| none =>
return assignment.set! varIdx (some 0)
| some eq =>
let mut assignment := assignment
let mut r := eq.const
for (i, coeff) in eq.coeffs do
if assignment[i].isNone then
assignment := readSolution i assignment
r := r + coeff*assignment[i].get!
return assignment.set! varIdx (some r)
partial def solveProblem' (p : Problem) : Solution := do
match readSolution? p with
| some solution => return solution
| none =>
let p := eliminateSingletons p
match readSolution? p with
| some solution => return solution
| none =>
let p := addAuxEquation p
return solveProblem' p
def isSatAssignment (p : Problem) (assignment : Array Int) : Bool :=
¬ p.equations.any fun _ (eq : Equation) => do
let mut r := 0
for (i, coeff) in eq.coeffs do
r := r + coeff*assignment[i]
return r ≠ eq.const
def solveProblem (p : Problem) : Solution :=
let nVars := p.nVars
match solveProblem' p with
| Solution.unsat =>
Solution.unsat
| Solution.sat assignment =>
let assignment' := assignment.extract 0 nVars
if isSatAssignment p assignment' then
Solution.sat assignment'
else
Solution.unsat
def error (msg : String) : IO α :=
throw <| IO.userError s!"Error: {msg}."
def Array.ithVal (xs : Array String) (i : Nat) (name : String) : IO Int := do
let some unparsed ← xs.get? i
| error s!"Missing {name}"
let some parsed ← String.toInt? unparsed
| error s!"Invalid {name}: `{unparsed}`"
return parsed
def main (args : List String) : IO UInt32 := do
let some path ← args.head?
| error "Usage: liasolver <input file>"
let lines ← IO.FS.lines path <&> Array.filter (¬·.isEmpty)
let some headerLine ← lines.get? 0
| error "No header line"
let header := headerLine.splitOn.toArray
let nEquations ← header.ithVal 0 "amount of equations"
let nVars ← header.ithVal 1 "amount of variables"
let mut equations : HashMap Nat Equation := mkHashMap
for line in lines[1:] do
let elems := line.splitOn.toArray
let nTerms ← elems.ithVal 0 "amount of equation terms"
let 0 ← elems.ithVal (elems.size - 1) "end of line symbol"
| error "Non-zero end of line symbol"
let const ← elems.ithVal (elems.size - 2) "constant value"
let mut coeffs := mkHashMap
for i in [1:elems.size-2:2] do
let coeff ← elems.ithVal i "coefficient"
let varIdx ← elems.ithVal (i + 1) "variable index"
if varIdx < 1 then
error "Invalid variable index"
let varIdx := varIdx.toNat - 1
if coeff ≠ 0 then
coeffs := coeffs.insert varIdx coeff
if coeffs.size ≠ 0 then
equations := equations.insert equations.size ⟨equations.size, coeffs, const⟩
match preprocess? equations with
| none =>
IO.println "UNSAT"
| some equations' =>
let problem : Problem := ⟨equations', mkHashMap, equations'.size, nVars.natAbs⟩
match solveProblem problem with
| Solution.unsat =>
IO.println "UNSAT"
| Solution.sat assignment =>
IO.println "SAT"
IO.println <| String.intercalate " " <| assignment.data.map toString
return 0

View file

@ -0,0 +1 @@
ex-50-50-1.leq

View file

@ -0,0 +1,2 @@
SAT
-104 -253 65 75 -271 269 234 -125 59 -291 23 -150 -127 -232 11 -66 -199 133 -51 -120 -141 276 24 6 -85 28 240 54 -182 -160 128 14 -212 269 -154 28 134 -125 -49 192 14 -130 69 -53 -157 264 103 48 -271 184