toysolver 0.0.5 → 0.0.6

raw patch · 248 files changed

+2603/−2437 lines, 248 filesdep +hashabledep +type-level-numbersdep ~basedep ~containersdep ~finite-field

Dependencies added: hashable, type-level-numbers

Dependency ranges changed: base, containers, finite-field, lattices

Files

• benchmarks/UF250.1065.100/uf250-01.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-010.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-0100.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-011.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-012.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-013.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-014.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-015.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-016.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-017.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-018.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-019.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-02.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-020.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-021.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-022.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-023.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-024.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-025.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-026.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-027.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-028.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-029.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-03.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-030.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-031.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-032.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-033.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-034.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-035.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-036.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-037.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-038.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-039.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-04.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-040.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-041.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-042.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-043.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-044.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-045.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-046.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-047.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-048.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-049.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-05.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-050.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-051.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-052.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-053.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-054.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-055.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-056.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-057.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-058.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-059.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-06.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-060.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-061.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-062.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-063.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-064.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-065.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-066.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-067.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-068.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-069.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-07.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-070.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-071.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-072.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-073.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-074.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-075.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-076.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-077.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-078.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-079.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-08.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-080.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-081.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-082.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-083.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-084.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-085.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-086.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-087.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-088.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-089.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-09.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-090.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-091.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-092.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-093.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-094.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-095.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-096.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-097.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-098.cnf +0/−3
• benchmarks/UF250.1065.100/uf250-099.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-01.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-010.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-0100.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-011.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-012.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-013.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-014.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-015.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-016.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-017.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-018.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-019.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-02.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-020.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-021.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-022.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-023.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-024.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-025.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-026.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-027.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-028.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-029.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-03.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-030.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-031.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-032.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-033.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-034.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-035.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-036.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-037.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-038.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-039.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-04.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-040.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-041.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-042.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-043.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-044.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-045.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-046.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-047.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-048.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-049.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-05.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-050.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-051.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-052.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-053.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-054.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-055.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-056.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-057.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-058.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-059.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-06.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-060.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-061.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-062.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-063.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-064.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-065.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-066.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-067.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-068.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-069.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-07.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-070.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-071.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-072.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-073.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-074.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-075.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-076.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-077.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-078.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-079.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-08.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-080.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-081.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-082.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-083.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-084.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-085.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-086.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-087.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-088.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-089.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-09.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-090.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-091.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-092.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-093.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-094.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-095.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-096.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-097.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-098.cnf +0/−3
• benchmarks/UUF250.1065.100/uuf250-099.cnf +0/−3
• src/Algorithm/CAD.hs +111/−107
• src/Algorithm/CongruenceClosure.hs +25/−40
• src/Algorithm/ContiTraverso.hs +14/−10
• src/Algorithm/FOLModelFinder.hs +18/−16
• src/Algorithm/Simplex2.hs +53/−51
• src/Converter/LP2SMT.hs +2/−1
• src/Converter/MaxSAT2LP.hs +2/−2
• src/Converter/PB2LP.hs +16/−14
• src/Converter/SAT2LP.hs +2/−2
• src/Data/AlgebraicNumber/Real.hs +29/−28
• src/Data/AlgebraicNumber/Root.hs +29/−31
• src/Data/LA.hs +34/−33
• src/Data/Polyhedron.hs +5/−4
• src/Data/Polynomial.hs +59/−680
• src/Data/Polynomial/Base.hs +802/−0
• src/Data/Polynomial/Factorization/FiniteField.hs +44/−31
• src/Data/Polynomial/Factorization/Hensel.hs +147/−0
• src/Data/Polynomial/Factorization/Integer.hs +8/−121
• src/Data/Polynomial/Factorization/Kronecker.hs +132/−0
• src/Data/Polynomial/Factorization/Rational.hs +12/−14
• src/Data/Polynomial/Factorization/SquareFree.hs +29/−14
• src/Data/Polynomial/Factorization/Zassenhaus.hs +171/−0
• src/Data/Polynomial/GBasis.hs +0/−150
• src/Data/Polynomial/GroebnerBasis.hs +144/−0
• src/Data/Polynomial/Interpolation/Lagrange.hs +5/−4
• src/Data/Polynomial/RootSeparation/Graeffe.hs +13/−12
• src/Data/Polynomial/RootSeparation/Sturm.hs +11/−10
• src/Data/Sign.hs +88/−44
• src/Data/Var.hs +6/−5
• src/SAT.hs +7/−2
• src/SAT/PBO/UnsatBased.hs +11/−10
• src/SAT/TseitinEncoder.hs +5/−3
• src/SAT/Types.hs +23/−22
• src/Text/GurobiSol.hs +2/−1
• src/Text/LPFile.hs +10/−8
• src/Text/MPSFile.hs +9/−7
• src/Text/SDPFile.hs +7/−6
• src/Util.hs +2/−1
• src/maxsatverify.hs +42/−0
• test/TestAReal.hs +30/−29
• test/TestAReal2.hs +3/−2
• test/TestCongruenceClosure.hs +34/−0
• test/TestContiTraverso.hs +10/−9
• test/TestPolynomial.hs +259/−208
• test/TestQE.hs +1/−1
• toysat/toysat.hs +64/−56
• toysolver.cabal +18/−5
• toysolver/toysolver.hs +55/−43

@@ -1071,6 +1071,3 @@ -184 203 126 0 -249 81 -231 0 141 231 25 0-%-0-

@@ -1071,6 +1071,3 @@ -175 141 143 0 8 -52 35 0 224 -163 24 0-%-0-

@@ -1071,6 +1071,3 @@ -139 197 50 0 -84 -240 18 0 37 -16 5 0-%-0-

@@ -1071,6 +1071,3 @@ 105 -128 -205 0 -107 1 60 0 -184 70 190 0-%-0-

@@ -1071,6 +1071,3 @@ 202 -84 73 0 180 -214 -52 0 -102 -97 222 0-%-0-

@@ -1071,6 +1071,3 @@ -210 128 -190 0 -226 -56 -127 0 118 -104 -67 0-%-0-

@@ -1071,6 +1071,3 @@ 144 -172 -173 0 -159 104 -66 0 -1 44 -243 0-%-0-

@@ -1071,6 +1071,3 @@ -188 129 -18 0 -104 87 223 0 204 -32 54 0-%-0-

@@ -1071,6 +1071,3 @@ 60 -245 94 0 192 172 -177 0 -162 -180 59 0-%-0-

@@ -1071,6 +1071,3 @@ -173 -84 103 0 98 14 118 0 192 209 -170 0-%-0-

@@ -1071,6 +1071,3 @@ 131 -47 -220 0 137 -177 100 0 -131 45 36 0-%-0-

@@ -1071,6 +1071,3 @@ -161 73 -121 0 124 -134 -172 0 -188 -227 -113 0-%-0-

@@ -1071,6 +1071,3 @@ 79 -33 -202 0 -40 -134 212 0 -239 146 250 0-%-0-

@@ -1071,6 +1071,3 @@ -110 -143 151 0 80 89 -94 0 -116 129 186 0-%-0-

@@ -1071,6 +1071,3 @@ -217 -100 63 0 -83 67 -85 0 17 -159 -223 0-%-0-

@@ -1071,6 +1071,3 @@ -207 66 193 0 -112 174 25 0 -122 27 -131 0-%-0-

@@ -1071,6 +1071,3 @@ 154 131 -208 0 204 210 -82 0 -116 158 137 0-%-0-

@@ -1071,6 +1071,3 @@ 66 208 -187 0 -61 149 -165 0 -136 153 -22 0-%-0-

@@ -1071,6 +1071,3 @@ 192 -234 -203 0 -182 -82 143 0 247 -205 -221 0-%-0-

@@ -1071,6 +1071,3 @@ 37 157 -34 0 176 -134 -81 0 -228 -110 207 0-%-0-

@@ -1071,6 +1071,3 @@ -178 -72 -132 0 232 185 -209 0 -229 -231 -113 0-%-0-

@@ -1071,6 +1071,3 @@ -94 -220 65 0 83 74 84 0 -81 124 197 0-%-0-

@@ -1071,6 +1071,3 @@ -180 237 204 0 91 130 239 0 7 181 25 0-%-0-

@@ -1071,6 +1071,3 @@ 215 101 177 0 -209 140 240 0 160 -191 127 0-%-0-

@@ -1071,6 +1071,3 @@ 93 82 -6 0 76 101 241 0 4 -93 -108 0-%-0-

@@ -1071,6 +1071,3 @@ -140 102 -203 0 59 -132 73 0 52 -14 -111 0-%-0-

@@ -1071,6 +1071,3 @@ 170 -213 63 0 142 -127 55 0 72 -78 -121 0-%-0-

@@ -1071,6 +1071,3 @@ 179 53 214 0 -58 223 -29 0 129 102 35 0-%-0-

@@ -1071,6 +1071,3 @@ -245 -179 -191 0 217 191 -67 0 226 167 -227 0-%-0-

@@ -1071,6 +1071,3 @@ 42 96 -17 0 -134 -238 45 0 223 129 185 0-%-0-

@@ -1071,6 +1071,3 @@ -139 -222 -183 0 -144 106 -107 0 175 154 -11 0-%-0-

@@ -1071,6 +1071,3 @@ -183 77 -186 0 22 -69 -168 0 21 -95 -127 0-%-0-

@@ -1071,6 +1071,3 @@ 157 -172 2 0 214 -230 120 0 -156 84 -40 0-%-0-

@@ -1071,6 +1071,3 @@ 15 -189 -13 0 209 -199 78 0 195 -66 -84 0-%-0-

@@ -1071,6 +1071,3 @@ -57 -223 13 0 181 70 -7 0 -232 69 186 0-%-0-

@@ -1071,6 +1071,3 @@ 48 150 96 0 194 -105 -169 0 -34 -76 246 0-%-0-

@@ -1071,6 +1071,3 @@ 220 -146 232 0 -127 39 -18 0 -107 11 86 0-%-0-

@@ -1071,6 +1071,3 @@ 93 136 205 0 -94 171 -145 0 160 101 83 0-%-0-

@@ -1071,6 +1071,3 @@ 207 6 227 0 217 98 101 0 132 -9 -157 0-%-0-

@@ -1071,6 +1071,3 @@ 154 21 165 0 3 221 -85 0 -179 151 31 0-%-0-

@@ -1071,6 +1071,3 @@ -216 52 172 0 74 224 37 0 62 73 17 0-%-0-

@@ -1071,6 +1071,3 @@ -94 39 142 0 146 -27 -3 0 -71 -25 93 0-%-0-

@@ -1071,6 +1071,3 @@ 46 193 -37 0 142 -15 -138 0 -161 -133 -62 0-%-0-

@@ -1071,6 +1071,3 @@ -145 -160 -138 0 142 -133 100 0 135 61 -3 0-%-0-

@@ -1071,6 +1071,3 @@ -50 -200 209 0 131 170 -3 0 -146 207 -147 0-%-0-

@@ -1071,6 +1071,3 @@ -137 -139 -216 0 -128 107 47 0 184 -231 -169 0-%-0-

@@ -1071,6 +1071,3 @@ -160 -43 -37 0 201 -200 28 0 -152 -185 -41 0-%-0-

@@ -1071,6 +1071,3 @@ 207 -7 2 0 162 156 -24 0 106 -99 -32 0-%-0-

@@ -1071,6 +1071,3 @@ 193 -197 -182 0 -115 -241 189 0 152 -86 206 0-%-0-

@@ -1071,6 +1071,3 @@ -25 122 -80 0 137 57 -50 0 -233 220 -12 0-%-0-

@@ -1071,6 +1071,3 @@ -198 -242 -117 0 140 -225 -227 0 42 -76 45 0-%-0-

@@ -1071,6 +1071,3 @@ 185 -156 36 0 200 -109 -159 0 -35 249 -194 0-%-0-

@@ -1071,6 +1071,3 @@ 187 -193 239 0 203 -143 -5 0 -27 131 -115 0-%-0-

@@ -1071,6 +1071,3 @@ -219 -247 148 0 11 -59 151 0 -41 -21 36 0-%-0-

@@ -1071,6 +1071,3 @@ -190 -228 -26 0 -13 -144 -76 0 -248 191 180 0-%-0-

@@ -1071,6 +1071,3 @@ -66 10 67 0 -211 -136 -150 0 -30 -190 78 0-%-0-

@@ -1071,6 +1071,3 @@ 76 -99 -30 0 212 70 177 0 -55 57 -250 0-%-0-

@@ -1071,6 +1071,3 @@ 160 124 50 0 108 -200 -157 0 -82 -133 -2 0-%-0-

@@ -1071,6 +1071,3 @@ -250 -157 -98 0 102 -143 -179 0 -48 -19 -132 0-%-0-

@@ -1071,6 +1071,3 @@ 154 188 65 0 174 -100 -213 0 -131 77 26 0-%-0-

@@ -1071,6 +1071,3 @@ -121 245 105 0 152 188 76 0 -115 143 250 0-%-0-

@@ -1071,6 +1071,3 @@ -160 95 15 0 -120 -214 -209 0 -107 -195 -93 0-%-0-

@@ -1071,6 +1071,3 @@ 91 66 249 0 144 51 -127 0 214 3 -166 0-%-0-

@@ -1071,6 +1071,3 @@ 48 93 -115 0 -49 4 120 0 176 -126 135 0-%-0-

@@ -1071,6 +1071,3 @@ -108 -37 226 0 2 103 -83 0 208 82 -154 0-%-0-

@@ -1071,6 +1071,3 @@ 227 -238 -94 0 111 211 -180 0 90 92 142 0-%-0-

@@ -1071,6 +1071,3 @@ -44 235 -156 0 -22 -223 14 0 161 -205 -149 0-%-0-

@@ -1071,6 +1071,3 @@ -14 -143 -62 0 169 198 250 0 -234 -157 10 0-%-0-

@@ -1071,6 +1071,3 @@ -162 -83 -217 0 -226 177 120 0 160 200 -134 0-%-0-

@@ -1071,6 +1071,3 @@ 136 -195 -202 0 32 -220 -16 0 -210 180 -12 0-%-0-

@@ -1071,6 +1071,3 @@ 44 109 -55 0 -145 -68 171 0 43 -239 -120 0-%-0-

@@ -1071,6 +1071,3 @@ -232 -84 -171 0 -198 -250 83 0 218 -238 -157 0-%-0-

@@ -1071,6 +1071,3 @@ -141 173 123 0 -68 -229 221 0 -128 225 -219 0-%-0-

@@ -1071,6 +1071,3 @@ 177 -120 -137 0 155 -249 -128 0 -186 154 -165 0-%-0-

@@ -1071,6 +1071,3 @@ 250 123 -151 0 -80 -90 -135 0 215 -83 200 0-%-0-

@@ -1071,6 +1071,3 @@ -13 -34 -25 0 198 55 -8 0 -149 -152 -27 0-%-0-

@@ -1071,6 +1071,3 @@ 118 -82 -168 0 8 93 -154 0 130 -220 -159 0-%-0-

@@ -1071,6 +1071,3 @@ -166 95 73 0 -154 213 -2 0 -203 -137 230 0-%-0-

@@ -1071,6 +1071,3 @@ 67 -10 23 0 -122 134 -167 0 75 185 -193 0-%-0-

@@ -1071,6 +1071,3 @@ 149 -59 118 0 -136 -96 229 0 164 -69 -107 0-%-0-

@@ -1071,6 +1071,3 @@ -241 58 168 0 113 122 178 0 185 -245 85 0-%-0-

@@ -1071,6 +1071,3 @@ 160 -158 -216 0 44 -34 -247 0 -16 36 139 0-%-0-

@@ -1071,6 +1071,3 @@ -133 -119 86 0 -10 -191 -103 0 -154 250 -219 0-%-0-

@@ -1071,6 +1071,3 @@ 39 -248 -209 0 -115 -90 125 0 174 -90 51 0-%-0-

@@ -1071,6 +1071,3 @@ -226 -68 110 0 15 -169 -241 0 142 -124 -215 0-%-0-

@@ -1071,6 +1071,3 @@ 179 147 209 0 124 -132 -241 0 -239 -163 -247 0-%-0-

@@ -1071,6 +1071,3 @@ 113 31 100 0 1 -228 61 0 197 -28 152 0-%-0-

@@ -1071,6 +1071,3 @@ -129 -163 -98 0 -164 -222 49 0 235 -8 195 0-%-0-

@@ -1071,6 +1071,3 @@ -103 84 170 0 -109 -86 93 0 110 76 -102 0-%-0-

@@ -1071,6 +1071,3 @@ -166 -103 -249 0 166 -36 208 0 144 -137 174 0-%-0-

@@ -1071,6 +1071,3 @@ -105 -159 -180 0 -150 244 179 0 192 -16 -154 0-%-0-

@@ -1071,6 +1071,3 @@ -202 183 -195 0 -12 -172 -67 0 -108 57 240 0-%-0-

@@ -1071,6 +1071,3 @@ -99 75 -160 0 -99 108 53 0 238 92 -53 0-%-0-

@@ -1071,6 +1071,3 @@ 2 -6 48 0 198 181 244 0 241 9 146 0-%-0-

@@ -1071,6 +1071,3 @@ -247 -88 176 0 233 -205 -132 0 -71 248 -127 0-%-0-

@@ -1071,6 +1071,3 @@ -211 -233 -30 0 37 4 -125 0 47 109 -151 0-%-0-

@@ -1071,6 +1071,3 @@ -51 -208 156 0 137 -43 126 0 -1 250 22 0-%-0-

@@ -1071,6 +1071,3 @@ -210 -167 -181 0 -52 -173 -60 0 -43 29 235 0-%-0-

@@ -1071,6 +1071,3 @@ -95 109 -150 0 -196 -8 16 0 51 -15 91 0-%-0-

@@ -1071,6 +1071,3 @@ -24 -144 9 0 137 -120 -228 0 -24 -117 158 0-%-0-

@@ -1071,6 +1071,3 @@ -206 107 -223 0 -159 -243 -13 0 118 -9 57 0-%-0-

@@ -1071,6 +1071,3 @@ 207 128 -96 0 -17 92 -195 0 -198 -228 -247 0-%-0-

@@ -1071,6 +1071,3 @@ 119 41 -152 0 -99 -187 -14 0 -207 -121 171 0-%-0-

@@ -1071,6 +1071,3 @@ 208 75 64 0 -171 163 207 0 -76 19 158 0-%-0-

@@ -1071,6 +1071,3 @@ 138 -111 201 0 -53 25 -175 0 115 230 -97 0-%-0-

@@ -1071,6 +1071,3 @@ -45 -64 56 0 9 220 150 0 149 32 -74 0-%-0-

@@ -1071,6 +1071,3 @@ -172 -98 139 0 48 83 -97 0 -8 239 -61 0-%-0-

@@ -1071,6 +1071,3 @@ 51 63 187 0 239 210 243 0 -30 -198 -169 0-%-0-

@@ -1071,6 +1071,3 @@ 49 -213 -32 0 -222 -4 -28 0 -106 -165 56 0-%-0-

@@ -1071,6 +1071,3 @@ 21 61 77 0 219 55 -237 0 61 -248 1 0-%-0-

@@ -1071,6 +1071,3 @@ -240 237 69 0 95 209 157 0 146 -149 224 0-%-0-

@@ -1071,6 +1071,3 @@ -140 -97 99 0 -73 23 -77 0 -96 4 -6 0-%-0-

@@ -1071,6 +1071,3 @@ 51 53 198 0 -94 13 -217 0 -211 188 -43 0-%-0-

@@ -1071,6 +1071,3 @@ -228 74 165 0 195 -7 74 0 129 190 76 0-%-0-

@@ -1071,6 +1071,3 @@ -200 195 -214 0 -230 -184 -113 0 -248 -58 -125 0-%-0-

@@ -1071,6 +1071,3 @@ 126 -144 164 0 -79 114 116 0 -62 -55 84 0-%-0-

@@ -1071,6 +1071,3 @@ -34 -122 129 0 126 21 171 0 4 11 -83 0-%-0-

@@ -1071,6 +1071,3 @@ -127 -38 194 0 -246 -66 128 0 -195 -175 -33 0-%-0-

@@ -1071,6 +1071,3 @@ 230 -161 5 0 -210 123 118 0 156 25 -115 0-%-0-

@@ -1071,6 +1071,3 @@ -93 -133 43 0 85 -152 -57 0 89 7 230 0-%-0-

@@ -1071,6 +1071,3 @@ -2 -236 -116 0 241 156 55 0 -99 59 190 0-%-0-

@@ -1071,6 +1071,3 @@ 48 231 -178 0 138 -58 42 0 -45 -62 -54 0-%-0-

@@ -1071,6 +1071,3 @@ -170 250 -47 0 208 -245 106 0 -225 91 98 0-%-0-

@@ -1071,6 +1071,3 @@ -216 45 -179 0 -30 -236 -59 0 244 -124 -189 0-%-0-

@@ -1071,6 +1071,3 @@ -81 -124 -207 0 131 -27 -39 0 37 113 -245 0-%-0-

@@ -1071,6 +1071,3 @@ 70 -180 -238 0 172 -231 -87 0 -42 -85 -54 0-%-0-

@@ -1071,6 +1071,3 @@ -166 -151 243 0 106 -84 -140 0 -246 144 -24 0-%-0-

@@ -1071,6 +1071,3 @@ 131 -80 244 0 114 -215 153 0 -42 -54 -55 0-%-0-

@@ -1071,6 +1071,3 @@ -47 38 223 0 -215 207 -196 0 -23 56 -167 0-%-0-

@@ -1071,6 +1071,3 @@ -51 22 -17 0 -208 96 -215 0 -10 7 125 0-%-0-

@@ -1071,6 +1071,3 @@ 152 13 -34 0 86 -184 -99 0 142 149 226 0-%-0-

@@ -1071,6 +1071,3 @@ 79 -106 -8 0 -93 -21 232 0 45 13 -30 0-%-0-

@@ -1071,6 +1071,3 @@ -74 -200 -39 0 134 -198 88 0 232 211 214 0-%-0-

@@ -1071,6 +1071,3 @@ -19 -208 164 0 246 -225 241 0 -159 40 -227 0-%-0-

@@ -1071,6 +1071,3 @@ 77 -158 170 0 -157 28 6 0 137 116 118 0-%-0-

@@ -1071,6 +1071,3 @@ -131 46 -87 0 27 -39 143 0 -93 6 -247 0-%-0-

@@ -1071,6 +1071,3 @@ 228 -144 66 0 150 227 66 0 135 61 47 0-%-0-

@@ -1071,6 +1071,3 @@ -8 49 237 0 -118 -103 -39 0 26 23 -231 0-%-0-

@@ -1071,6 +1071,3 @@ 219 85 31 0 -176 139 -28 0 27 74 -72 0-%-0-

@@ -1071,6 +1071,3 @@ -88 136 19 0 -10 -170 58 0 16 113 -146 0-%-0-

@@ -1071,6 +1071,3 @@ -61 234 193 0 -5 -115 99 0 217 193 7 0-%-0-

@@ -1071,6 +1071,3 @@ -176 84 -64 0 -105 -39 -227 0 -149 141 -246 0-%-0-

@@ -1071,6 +1071,3 @@ -25 7 215 0 10 230 -138 0 154 -176 188 0-%-0-

@@ -1071,6 +1071,3 @@ -138 -237 81 0 -27 -199 121 0 -248 -65 156 0-%-0-

@@ -1071,6 +1071,3 @@ 61 -60 -106 0 -38 86 62 0 -145 107 228 0-%-0-

@@ -1071,6 +1071,3 @@ 168 -3 175 0 1 -233 -51 0 -185 196 145 0-%-0-

@@ -1071,6 +1071,3 @@ 69 -171 -109 0 223 247 -196 0 56 145 -167 0-%-0-

@@ -1071,6 +1071,3 @@ 235 -132 133 0 -138 144 -109 0 192 138 179 0-%-0-

@@ -1071,6 +1071,3 @@ -217 239 47 0 -175 68 34 0 -16 195 65 0-%-0-

@@ -1071,6 +1071,3 @@ -46 247 -15 0 213 -82 -189 0 -17 -148 -86 0-%-0-

@@ -1071,6 +1071,3 @@ 66 -170 -70 0 -98 169 141 0 181 190 -103 0-%-0-

@@ -1071,6 +1071,3 @@ 201 -93 -215 0 185 -195 107 0 86 69 -143 0-%-0-

@@ -1071,6 +1071,3 @@ -123 -124 48 0 139 86 171 0 27 177 -108 0-%-0-

@@ -1071,6 +1071,3 @@ -93 161 201 0 -141 -186 1 0 -132 -156 101 0-%-0-

@@ -1071,6 +1071,3 @@ 136 77 167 0 75 -117 24 0 -146 -184 -70 0-%-0-

@@ -1071,6 +1071,3 @@ 73 -194 134 0 89 -75 81 0 -195 -95 232 0-%-0-

@@ -1071,6 +1071,3 @@ -223 -4 152 0 147 173 181 0 -176 150 54 0-%-0-

@@ -1071,6 +1071,3 @@ 10 -219 -211 0 208 216 -91 0 127 99 34 0-%-0-

@@ -1071,6 +1071,3 @@ 147 160 -115 0 141 -76 101 0 -117 30 212 0-%-0-

@@ -1071,6 +1071,3 @@ -124 5 229 0 -86 -136 -131 0 -241 -47 -19 0-%-0-

@@ -1071,6 +1071,3 @@ 138 224 23 0 230 -178 116 0 119 116 223 0-%-0-

@@ -1071,6 +1071,3 @@ -168 -183 -173 0 119 -102 3 0 214 -98 34 0-%-0-

@@ -1071,6 +1071,3 @@ -18 47 -148 0 -13 -91 -54 0 -175 -12 -90 0-%-0-

@@ -1071,6 +1071,3 @@ 80 180 15 0 187 -185 -27 0 31 45 -29 0-%-0-

@@ -1071,6 +1071,3 @@ 249 11 -93 0 -10 203 85 0 -239 -21 -214 0-%-0-

@@ -1071,6 +1071,3 @@ -231 -185 -109 0 17 -80 -236 0 -88 -163 -71 0-%-0-

@@ -1071,6 +1071,3 @@ -151 30 -111 0 35 15 -92 0 134 -228 -186 0-%-0-

@@ -1071,6 +1071,3 @@ -91 10 173 0 161 191 210 0 -141 243 218 0-%-0-

@@ -1071,6 +1071,3 @@ -68 28 -147 0 33 -179 -164 0 -148 -132 -154 0-%-0-

@@ -1071,6 +1071,3 @@ -145 155 233 0 -202 31 -144 0 53 42 -165 0-%-0-

@@ -1071,6 +1071,3 @@ -107 -38 -211 0 129 211 -47 0 -23 120 -196 0-%-0-

@@ -1071,6 +1071,3 @@ 53 -25 -183 0 242 -16 -223 0 -171 241 127 0-%-0-

@@ -1071,6 +1071,3 @@ 46 -24 176 0 -217 130 -154 0 68 -118 -64 0-%-0-

@@ -1071,6 +1071,3 @@ -176 -83 67 0 102 -13 -3 0 92 119 -70 0-%-0-

@@ -1071,6 +1071,3 @@ -15 38 -187 0 -201 -245 -59 0 -173 -42 38 0-%-0-

@@ -1071,6 +1071,3 @@ -20 104 246 0 94 79 136 0 -88 -37 -173 0-%-0-

@@ -1071,6 +1071,3 @@ -193 -199 -220 0 -186 184 -178 0 -105 -72 -247 0-%-0-

@@ -1071,6 +1071,3 @@ 196 132 -91 0 -76 52 -178 0 191 -77 227 0-%-0-

@@ -1071,6 +1071,3 @@ -17 -194 -202 0 -233 -70 -176 0 126 92 -57 0-%-0-

@@ -1071,6 +1071,3 @@ 9 -171 -192 0 -100 -234 163 0 213 162 -118 0-%-0-

@@ -1071,6 +1071,3 @@ -75 -231 -2 0 98 -197 189 0 -32 155 246 0-%-0-

@@ -1071,6 +1071,3 @@ -85 76 -241 0 171 216 71 0 228 -2 -15 0-%-0-

@@ -1071,6 +1071,3 @@ -151 126 -192 0 -171 -49 11 0 -60 -248 195 0-%-0-

@@ -1071,6 +1071,3 @@ -79 -11 172 0 -30 67 37 0 -131 36 110 0-%-0-

@@ -1071,6 +1071,3 @@ 150 -187 -52 0 205 52 -155 0 -250 51 -105 0-%-0-

@@ -1071,6 +1071,3 @@ 161 -89 -229 0 127 232 212 0 -126 -151 161 0-%-0-

@@ -1071,6 +1071,3 @@ -223 -141 -154 0 -8 -54 242 0 -31 161 99 0-%-0-

@@ -1071,6 +1071,3 @@ 13 -221 -72 0 -85 175 -48 0 -77 44 217 0-%-0-

@@ -1071,6 +1071,3 @@ -99 195 -168 0 201 -153 -197 0 -57 -39 -125 0-%-0-

@@ -1071,6 +1071,3 @@ -10 210 238 0 215 124 -100 0 -48 224 141 0-%-0-

@@ -1071,6 +1071,3 @@ -59 -218 231 0 112 8 -102 0 100 235 -126 0-%-0-

@@ -1071,6 +1071,3 @@ 198 -13 203 0 3 -42 -170 0 -124 47 16 0-%-0-

@@ -1071,6 +1071,3 @@ -94 25 182 0 2 -120 -134 0 239 69 -190 0-%-0-

@@ -1071,6 +1071,3 @@ -66 -65 -178 0 198 228 151 0 -233 -162 -189 0-%-0-

@@ -1071,6 +1071,3 @@ 215 199 -237 0 -112 239 155 0 -105 -94 203 0-%-0-

@@ -1071,6 +1071,3 @@ 90 -152 -197 0 -23 -75 -46 0 221 140 -139 0-%-0-

@@ -1071,6 +1071,3 @@ 97 -83 124 0 -16 76 -26 0 36 -54 174 0-%-0-

@@ -1071,6 +1071,3 @@ 156 8 101 0 91 -220 98 0 228 -200 -109 0-%-0-

@@ -1071,6 +1071,3 @@ 103 206 204 0 -33 -147 -15 0 182 -189 -160 0-%-0-

@@ -1071,6 +1071,3 @@ -140 -187 46 0 -166 -6 138 0 239 -111 -223 0-%-0-

@@ -29,7 +29,6 @@   -- * Basic data structures     Point (..)   , Cell (..)-  , module Data.Sign    -- * Projection   , project@@ -50,7 +49,9 @@ import Data.List import Data.Maybe import Data.Ord+import Data.Map (Map) import qualified Data.Map as Map+import Data.Set (Set) import qualified Data.Set as Set import Text.Printf import Text.PrettyPrint.HughesPJClass@@ -58,9 +59,12 @@ import Data.ArithRel import qualified Data.AlgebraicNumber.Real as AReal import Data.DNF-import Data.Polynomial-import Data.Sign+import Data.Polynomial (Polynomial, UPolynomial, X (..), PrettyVar, PrettyCoeff)+import qualified Data.Polynomial as P+import Data.Sign (Sign (..))+import qualified Data.Sign as Sign + import Debug.Trace  -- ---------------------------------------------------------------------------@@ -84,7 +88,7 @@  -- --------------------------------------------------------------------------- -type SignConf c = [(Cell c, Map.Map (UPolynomial c) Sign)]+type SignConf c = [(Cell c, Map (UPolynomial c) Sign)]  emptySignConf :: SignConf c emptySignConf =@@ -99,9 +103,9 @@     f :: SignConf c -> [String]     f = concatMap $ \(cell, m) -> showCell cell : g m -    g :: Map.Map (UPolynomial c) Sign -> [String]+    g :: Map (UPolynomial c) Sign -> [String]     g m =-      [printf "  %s: %s" (prettyShow p) (showSign s) | (p, s) <- Map.toList m]+      [printf "  %s: %s" (prettyShow p) (Sign.symbol s) | (p, s) <- Map.toList m]  -- --------------------------------------------------------------------------- @@ -112,70 +116,70 @@   -> UPolynomial k   -> (k, Integer, UPolynomial k) mr p q-  | n >= m    = assert (constant (bm^(n-m+1)) * p == q * l + r && m > deg r) $ (bm, n-m+1, r)+  | n >= m    = assert (P.constant (bm^(n-m+1)) * p == q * l + r && m > P.deg r) $ (bm, n-m+1, r)   | otherwise = error "mr p q: not (deg p >= deg q)"   where-    x = var X-    n = deg p-    m = deg q-    (bm, _) = leadingTerm grlex q+    x = P.var X+    n = P.deg p+    m = P.deg q+    bm = P.lc P.grlex q     (l,r) = f p n      f :: UPolynomial k -> Integer -> (UPolynomial k, UPolynomial k)     f p n       | n==m =-          let l = constant an-              r = constant bm * p - constant an * q-          in assert (constant (bm^(n-m+1)) * p == q*l + r && m > deg r) $ (l, r)+          let l = P.constant an+              r = P.constant bm * p - P.constant an * q+          in assert (P.constant (bm^(n-m+1)) * p == q*l + r && m > P.deg r) $ (l, r)       | otherwise =-          let p'     = (constant bm * p - constant an * x^(n-m) * q)+          let p'     = (P.constant bm * p - P.constant an * x^(n-m) * q)               (l',r) = f p' (n-1)-              l      = l' + constant (an*bm^(n-m)) * x^(n-m)-          in assert (n > deg p') $-             assert (constant (bm^(n-m+1)) * p == q*l + r && m > deg r) $ (l, r)+              l      = l' + P.constant (an*bm^(n-m)) * x^(n-m)+          in assert (n > P.deg p') $+             assert (P.constant (bm^(n-m+1)) * p == q*l + r && m > P.deg r) $ (l, r)       where-        an = coeff (mmFromList [(X, n)]) p+        an = P.coeff (P.var X `P.mpow` n) p  test_mr_1 :: (Coeff Int, Integer, UPolynomial (Coeff Int))-test_mr_1 = mr (toUPolynomialOf p 3) (toUPolynomialOf q 3)+test_mr_1 = mr (P.toUPolynomialOf p 3) (P.toUPolynomialOf q 3)   where-    a = var 0-    b = var 1-    c = var 2-    x = var 3+    a = P.var 0+    b = P.var 1+    c = P.var 2+    x = P.var 3     p = a*x^(2::Int) + b*x + c     q = 2*a*x + b  test_mr_2 :: (Coeff Int, Integer, UPolynomial (Coeff Int))-test_mr_2 = mr (toUPolynomialOf p 3) (toUPolynomialOf p 3)+test_mr_2 = mr (P.toUPolynomialOf p 3) (P.toUPolynomialOf p 3)   where-    a = var 0-    b = var 1-    c = var 2-    x = var 3+    a = P.var 0+    b = P.var 1+    c = P.var 2+    x = P.var 3     p = a*x^(2::Int) + b*x + c  -- ---------------------------------------------------------------------------  type Coeff v = Polynomial Rational v -type M v = StateT (Map.Map (Polynomial Rational v) (Set.Set Sign)) []+type M v = StateT (Map (Polynomial Rational v) (Set Sign)) [] -runM :: M v a -> [(a, Map.Map (Polynomial Rational v) (Set.Set Sign))]+runM :: M v a -> [(a, Map (Polynomial Rational v) (Set Sign))] runM m = runStateT m Map.empty  assume :: (Ord v, Show v, PrettyVar v) => Polynomial Rational v -> [Sign] -> M v () assume p ss =-  if deg p == 0+  if P.deg p <= 0     then do-      let c = coeff mmOne p-      guard $ signOf c `elem` ss-    else do      -      let (c,_) = leadingTerm grlex p-          p' = mapCoeff (/c) p+      let c = P.coeff P.mone p+      guard $ Sign.signOf c `elem` ss+    else do+      let c  = P.lc P.grlex p+          p' = P.mapCoeff (/c) p       m <- get       let ss1 = Map.findWithDefault (Set.fromList [Neg, Zero, Pos]) p' m-          ss2 = Set.intersection ss1 $ Set.fromList $ [s `signDiv` signOf c | s <- ss]+          ss2 = Set.intersection ss1 $ Set.fromList $ [s `Sign.div` Sign.signOf c | s <- ss]       guard $ not $ Set.null ss2       put $ Map.insert p' ss2 m @@ -185,24 +189,24 @@   -> [([(Polynomial Rational v, [Sign])], [Cell (Polynomial Rational v)])] project cs = [ (guess2cond gs, cells) | (cells, gs) <- result ]   where-    result :: [([Cell (Polynomial Rational v)], Map.Map (Polynomial Rational v) (Set.Set Sign))]+    result :: [([Cell (Polynomial Rational v)], Map (Polynomial Rational v) (Set Sign))]     result = runM $ do       forM_ cs $ \(p,ss) -> do-        when (1 > deg p) $ assume (coeff mmOne p) ss+        when (1 > P.deg p) $ assume (P.coeff P.mone p) ss       conf <- buildSignConf (map fst cs)       let satCells = [cell | (cell, m) <- conf, cell /= Point NegInf, cell /= Point PosInf, ok m]       guard $ not $ null satCells       return satCells -    ok :: Map.Map (UPolynomial (Polynomial Rational v)) Sign -> Bool+    ok :: Map (UPolynomial (Polynomial Rational v)) Sign -> Bool     ok m = and [checkSign m p ss | (p,ss) <- cs]       where         checkSign m p ss =-          if 1 > deg p +          if 1 > P.deg p              then True -- already assumed             else (m Map.! p) `elem` ss -    guess2cond :: Map.Map (Polynomial Rational v) (Set.Set Sign) -> [(Polynomial Rational v, [Sign])]+    guess2cond :: Map (Polynomial Rational v) (Set Sign) -> [(Polynomial Rational v, [Sign])]     guess2cond gs = [(p, Set.toList ss)  | (p, ss) <- Map.toList gs]  buildSignConf@@ -211,20 +215,20 @@   -> M v (SignConf (Polynomial Rational v)) buildSignConf ps = do   ps2 <- collectPolynomials (Set.fromList ps)-  let ts = sortBy (comparing deg) (Set.toList ps2)+  let ts = sortBy (comparing P.deg) (Set.toList ps2)   foldM (flip refineSignConf) emptySignConf ts  collectPolynomials   :: (Ord v, Show v, PrettyVar v)-  => Set.Set (UPolynomial (Polynomial Rational v))-  -> M v (Set.Set (UPolynomial (Polynomial Rational v)))+  => Set (UPolynomial (Polynomial Rational v))+  -> M v (Set (UPolynomial (Polynomial Rational v))) collectPolynomials ps = go Set.empty (f ps)   where-    f = Set.filter (\p -> deg p > 0) +    f = Set.filter (\p -> P.deg p > 0)      go result ps | Set.null ps = return result     go result ps = do-      let rs1 = filter (\p -> deg p > 0) [deriv p X | p <- Set.toList ps]-      rs2 <- liftM (filter (\p -> deg p > 0) . map (\(_,_,r) -> r) . concat) $+      let rs1 = filter (\p -> P.deg p > 0) [P.deriv p X | p <- Set.toList ps]+      rs2 <- liftM (filter (\p -> P.deg p > 0) . map (\(_,_,r) -> r) . concat) $         forM [(p1,p2) | p1 <- Set.toList ps, p2 <- Set.toList ps ++ Set.toList result, p1 /= p2] $ \(p1,p2) -> do           ret1 <- zmod p1 p2           ret2 <- zmod p2 p1@@ -238,7 +242,7 @@   -> M v (Polynomial Rational v, Integer) getHighestNonzeroTerm p = go $ sortBy (flip (comparing snd)) cs   where-    cs = [(c, deg mm) | (c,mm) <- terms p]+    cs = [(c, P.deg mm) | (c,mm) <- P.terms p]      go :: [(Polynomial Rational v, Integer)] -> M v (Polynomial Rational v, Integer)     go [] = return (0, -1)@@ -258,8 +262,8 @@   if not (d >= e) || 0 >= e     then return Nothing     else do-      let p' = fromTerms [(pi, mm) | (pi, mm) <- terms p, deg mm <= d]-          q' = fromTerms [(qi, mm) | (qi, mm) <- terms q, deg mm <= e]+      let p' = P.fromTerms [(pi, mm) | (pi, mm) <- P.terms p, P.deg mm <= d]+          q' = P.fromTerms [(qi, mm) | (qi, mm) <- P.terms q, P.deg mm <= e]       return $ Just $ mr p' q'  refineSignConf@@ -270,8 +274,8 @@ refineSignConf p conf = liftM (extendIntervals 0) $ mapM extendPoint conf   where      extendPoint-      :: (Cell (Polynomial Rational v), Map.Map (UPolynomial (Polynomial Rational v)) Sign)-      -> M v (Cell (Polynomial Rational v), Map.Map (UPolynomial (Polynomial Rational v)) Sign)+      :: (Cell (Polynomial Rational v), Map (UPolynomial (Polynomial Rational v)) Sign)+      -> M v (Cell (Polynomial Rational v), Map (UPolynomial (Polynomial Rational v)) Sign)     extendPoint (Point pt, m) = do       s <- signAt pt m       return (Point pt, Map.insert p s m)@@ -279,8 +283,8 @@       extendIntervals       :: Int-      -> [(Cell (Polynomial Rational v), Map.Map (UPolynomial (Polynomial Rational v)) Sign)]-      -> [(Cell (Polynomial Rational v), Map.Map (UPolynomial (Polynomial Rational v)) Sign)]+      -> [(Cell (Polynomial Rational v), Map (UPolynomial (Polynomial Rational v)) Sign)]+      -> [(Cell (Polynomial Rational v), Map (UPolynomial (Polynomial Rational v)) Sign)]     extendIntervals !n (pt1@(Point _, m1) : (Interval lb ub, m) : pt2@(Point _, m2) : xs) =       pt1 : ys ++ extendIntervals n2 (pt2 : xs)       where@@ -300,7 +304,7 @@                           )     extendIntervals _ xs = xs  -    signAt :: Point (Polynomial Rational v) -> Map.Map (UPolynomial (Polynomial Rational v)) Sign -> M v Sign+    signAt :: Point (Polynomial Rational v) -> Map (UPolynomial (Polynomial Rational v)) Sign -> M v Sign     signAt PosInf _ = do       (c,_) <- getHighestNonzeroTerm p       signCoeff c@@ -308,18 +312,18 @@       (c,d) <- getHighestNonzeroTerm p       if even d         then signCoeff c-        else liftM signNegate $ signCoeff c+        else liftM Sign.negate $ signCoeff c     signAt (RootOf q _) m = do       Just (bm,k,r) <- zmod p q-      s1 <- if deg r > 0+      s1 <- if P.deg r > 0             then return $ m Map.! r-            else signCoeff $ coeff mmOne r+            else signCoeff $ P.coeff P.mone r       -- 場合分けを出来るだけ避ける       if even k         then return s1         else do           s2 <- signCoeff bm-          return $ signDiv s1 (signPow s2 k)+          return $ s1 `Sign.div` Sign.pow s2 k      signCoeff :: Polynomial Rational v -> M v Sign     signCoeff c =@@ -329,7 +333,7 @@  -- --------------------------------------------------------------------------- -type Model v = Map.Map v AReal.AReal+type Model v = Map v AReal.AReal  findSample :: Ord v => Model v -> Cell (Polynomial Rational v) -> Maybe AReal.AReal findSample m cell =@@ -359,13 +363,13 @@ evalPoint _ PosInf = PosInf evalPoint m (RootOf p n) = RootOf (AReal.minimalPolynomial a) (AReal.rootIndex a)   where-    a = AReal.realRootsEx (mapCoeff (eval (m Map.!) . mapCoeff fromRational) p) !! n+    a = AReal.realRootsEx (P.mapCoeff (P.eval (m Map.!) . P.mapCoeff fromRational) p) !! n  -- ---------------------------------------------------------------------------  solve   :: forall v. (Ord v, Show v, PrettyVar v)-  => Set.Set v+  => Set v   -> [(Rel (Polynomial Rational v))]   -> Maybe (Model v) solve vs cs0 = solve' vs (map f cs0)@@ -380,18 +384,18 @@  solve'   :: forall v. (Ord v, Show v, PrettyVar v)-  => Set.Set v+  => Set v   -> [(Polynomial Rational v, [Sign])]   -> Maybe (Model v) solve' vs0 cs0 = go (Set.toList vs0) cs0   where     go :: [v] -> [(Polynomial Rational v, [Sign])] -> Maybe (Model v)     go [] cs =-      if and [signOf v `elem` ss | (p,ss) <- cs, let v = eval (\_ -> undefined) p]+      if and [Sign.signOf v `elem` ss | (p,ss) <- cs, let v = P.eval (\_ -> undefined) p]       then Just Map.empty       else Nothing     go (v:vs) cs = listToMaybe $ do-      (cs2, cell:_) <- project [(toUPolynomialOf p v, ss) | (p,ss) <- cs]+      (cs2, cell:_) <- project [(P.toUPolynomialOf p v, ss) | (p,ss) <- cs]       case go vs cs2 of         Nothing -> mzero         Just m -> do@@ -440,7 +444,7 @@ dumpSignConf   :: forall v.      (Ord v, PrettyVar v, Show v)-  => [(SignConf (Polynomial Rational v), Map.Map (Polynomial Rational v) (Set.Set Sign))]+  => [(SignConf (Polynomial Rational v), Map (Polynomial Rational v) (Set Sign))]   -> IO () dumpSignConf x =    forM_ x $ \(conf, as) -> do@@ -454,7 +458,7 @@ test1a :: IO () test1a = mapM_ putStrLn $ showSignConf conf   where-    x = var X+    x = P.var X     ps :: [UPolynomial (Polynomial Rational Int)]     ps = [x + 1, -2*x + 3, x]     [(conf, _)] = runM $ buildSignConf ps@@ -462,7 +466,7 @@ test1b :: Bool test1b = isJust $ solve vs cs   where-    x = var X+    x = P.var X     vs = Set.singleton X     cs = [x + 1 .>. 0, -2*x + 3 .>. 0, x .>. 0] @@ -471,16 +475,16 @@   m <- solve' (Set.singleton X) cs   guard $ and $ do     (p, ss) <- cs-    let val = eval (m Map.!) (mapCoeff fromRational p)-    return $ signOf val `elem` ss+    let val = P.eval (m Map.!) (P.mapCoeff fromRational p)+    return $ Sign.signOf val `elem` ss   where-    x = var X+    x = P.var X     cs = [(x + 1, [Pos]), (-2*x + 3, [Pos]), (x, [Pos])]  test2a :: IO () test2a = mapM_ putStrLn $ showSignConf conf   where-    x = var X+    x = P.var X     ps :: [UPolynomial (Polynomial Rational Int)]     ps = [x^(2::Int)]     [(conf, _)] = runM $ buildSignConf ps@@ -488,7 +492,7 @@ test2b :: Bool test2b = isNothing $ solve vs cs   where-    x = var X+    x = P.var X     vs = Set.singleton X     cs = [x^(2::Int) .<. 0] @@ -497,53 +501,53 @@ test_project :: DNF (Polynomial Rational Int, [Sign]) test_project = DNF $ map fst $ project [(p', [Zero])]   where-    a = var 0-    b = var 1-    c = var 2-    x = var 3+    a = P.var 0+    b = P.var 1+    c = P.var 2+    x = P.var 3     p :: Polynomial Rational Int     p = a*x^(2::Int) + b*x + c-    p' = toUPolynomialOf p 3+    p' = P.toUPolynomialOf p 3  test_project_print :: IO () test_project_print = putStrLn $ showDNF $ test_project  test_project_2 = project [(p, [Zero]), (x, [Pos])]   where-    x = var X+    x = P.var X     p :: UPolynomial (Polynomial Rational Int)     p = x^(2::Int) + 4*x - 10 -test_project_3_print =  dumpProjection $ project [(toUPolynomialOf p 0, [Neg])]+test_project_3_print =  dumpProjection $ project [(P.toUPolynomialOf p 0, [Neg])]   where-    a = var 0-    b = var 1-    c = var 2+    a = P.var 0+    b = P.var 1+    c = P.var 2     p :: Polynomial Rational Int     p = a^(2::Int) + b^(2::Int) + c^(2::Int) - 1  test_solve = solve vs [p .<. 0]   where-    a = var 0-    b = var 1-    c = var 2+    a = P.var 0+    b = P.var 1+    c = P.var 2     vs = Set.fromList [0,1,2]     p :: Polynomial Rational Int     p = a^(2::Int) + b^(2::Int) + c^(2::Int) - 1  test_collectPolynomials-  :: [( Set.Set (UPolynomial (Polynomial Rational Int))-      , Map.Map (Polynomial Rational Int) (Set.Set Sign)+  :: [( Set (UPolynomial (Polynomial Rational Int))+      , Map (Polynomial Rational Int) (Set Sign)       )] test_collectPolynomials = runM $ collectPolynomials (Set.singleton p')   where-    a = var 0-    b = var 1-    c = var 2-    x = var 3+    a = P.var 0+    b = P.var 1+    c = P.var 2+    x = P.var 3     p :: Polynomial Rational Int     p = a*x^(2::Int) + b*x + c-    p' = toUPolynomialOf p 3+    p' = P.toUPolynomialOf p 3  test_collectPolynomials_print :: IO () test_collectPolynomials_print = do@@ -553,33 +557,33 @@     forM_  (Map.toList s) $ \(p, sign) ->       printf "%s %s\n" (prettyShow p) (show sign) -test_buildSignConf :: [(SignConf (Polynomial Rational Int), Map.Map (Polynomial Rational Int) (Set.Set Sign))]-test_buildSignConf = runM $ buildSignConf [toUPolynomialOf p 3]+test_buildSignConf :: [(SignConf (Polynomial Rational Int), Map (Polynomial Rational Int) (Set Sign))]+test_buildSignConf = runM $ buildSignConf [P.toUPolynomialOf p 3]   where-    a = var 0-    b = var 1-    c = var 2-    x = var 3+    a = P.var 0+    b = P.var 1+    c = P.var 2+    x = P.var 3     p :: Polynomial Rational Int     p = a*x^(2::Int) + b*x + c  test_buildSignConf_print :: IO () test_buildSignConf_print = dumpSignConf test_buildSignConf -test_buildSignConf_2 :: [(SignConf (Polynomial Rational Int), Map.Map (Polynomial Rational Int) (Set.Set Sign))]-test_buildSignConf_2 = runM $ buildSignConf [toUPolynomialOf p 0 | p <- ps]+test_buildSignConf_2 :: [(SignConf (Polynomial Rational Int), Map (Polynomial Rational Int) (Set Sign))]+test_buildSignConf_2 = runM $ buildSignConf [P.toUPolynomialOf p 0 | p <- ps]   where-    x = var 0+    x = P.var 0     ps :: [Polynomial Rational Int]     ps = [x + 1, -2*x + 3, x]  test_buildSignConf_2_print :: IO () test_buildSignConf_2_print = dumpSignConf test_buildSignConf_2 -test_buildSignConf_3 :: [(SignConf (Polynomial Rational Int), Map.Map (Polynomial Rational Int) (Set.Set Sign))]-test_buildSignConf_3 = runM $ buildSignConf [toUPolynomialOf p 0 | p <- ps]+test_buildSignConf_3 :: [(SignConf (Polynomial Rational Int), Map (Polynomial Rational Int) (Set Sign))]+test_buildSignConf_3 = runM $ buildSignConf [P.toUPolynomialOf p 0 | p <- ps]   where-    x = var 0+    x = P.var 0     ps :: [Polynomial Rational Int]     ps = [x, 2*x] 

@@ -31,7 +31,8 @@ import Control.Monad import Data.IORef import Data.Maybe-import qualified Data.IntMap as IM+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap  type Var = Int @@ -47,20 +48,20 @@   = Solver   { svVarCounter           :: IORef Int   , svPending              :: IORef [PendingEqn]-  , svRepresentativeTable  :: IORef (IM.IntMap Var) -- 本当は配列が良い-  , svClassList            :: IORef (IM.IntMap [Var])-  , svUseList              :: IORef (IM.IntMap [Eqn1])-  , svLookupTable          :: IORef (IM.IntMap (IM.IntMap Eqn1))+  , svRepresentativeTable  :: IORef (IntMap Var) -- 本当は配列が良い+  , svClassList            :: IORef (IntMap [Var])+  , svUseList              :: IORef (IntMap [Eqn1])+  , svLookupTable          :: IORef (IntMap (IntMap Eqn1))   }  newSolver :: IO Solver newSolver = do   vcnt     <- newIORef 0   pending  <- newIORef []-  rep      <- newIORef IM.empty-  classes  <- newIORef IM.empty-  useList  <- newIORef IM.empty-  lookup   <- newIORef IM.empty+  rep      <- newIORef IntMap.empty+  classes  <- newIORef IntMap.empty+  useList  <- newIORef IntMap.empty+  lookup   <- newIORef IntMap.empty   return $     Solver     { svVarCounter          = vcnt@@ -75,9 +76,9 @@ newVar solver = do   v <- readIORef (svVarCounter solver)   writeIORef (svVarCounter solver) $! v + 1-  modifyIORef (svRepresentativeTable solver) (IM.insert v v)-  modifyIORef (svClassList solver) (IM.insert v [v])-  modifyIORef (svUseList solver) (IM.insert v [])+  modifyIORef (svRepresentativeTable solver) (IntMap.insert v v)+  modifyIORef (svClassList solver) (IntMap.insert v [v])+  modifyIORef (svUseList solver) (IntMap.insert v [])   return v  merge :: Solver -> (FlatTerm, Var) -> IO ()@@ -97,8 +98,8 @@         Nothing -> do           setLookup solver a1' a2' (FTApp a1 a2, a)           modifyIORef (svUseList solver) $-            IM.alter (Just . ((FTApp a1 a2, a) :) . fromMaybe []) a1' .-            IM.alter (Just . ((FTApp a1 a2, a) :) . fromMaybe []) a2'+            IntMap.alter (Just . ((FTApp a1 a2, a) :) . fromMaybe []) a1' .+            IntMap.alter (Just . ((FTApp a1 a2, a) :) . fromMaybe []) a2'  propagate :: Solver -> IO () propagate solver = go@@ -122,20 +123,20 @@         then return ()         else do           clist <- readIORef (svClassList  solver)-          let classA = clist IM.! a'-              classB = clist IM.! b'+          let classA = clist IntMap.! a'+              classB = clist IntMap.! b'           if length classA < length classB             then update a' b' classA classB             else update b' a' classB classA      update a' b' classA classB = do       modifyIORef (svRepresentativeTable solver) $ -        IM.union (IM.fromList [(c,b') | c <- classA])+        IntMap.union (IntMap.fromList [(c,b') | c <- classA])       modifyIORef (svClassList solver) $-        IM.insert b' (classA ++ classB) . IM.delete a'+        IntMap.insert b' (classA ++ classB) . IntMap.delete a'        useList <- readIORef (svUseList solver)-      forM_ (useList IM.! a') $ \(FTApp c1 c2, c) -> do -- FIXME: not exhaustive+      forM_ (useList IntMap.! a') $ \(FTApp c1 c2, c) -> do -- FIXME: not exhaustive         c1' <- getRepresentative solver c1         c2' <- getRepresentative solver c2         ret <- lookup solver c1' c2'@@ -144,7 +145,7 @@             addToPending solver $ Right ((FTApp c1 c2, c), (FTApp d1 d2, d))           Nothing -> do             return ()-      writeIORef (svUseList solver) $ IM.delete a' useList        +      writeIORef (svUseList solver) $ IntMap.delete a' useList          areCongruent :: Solver -> FlatTerm -> FlatTerm -> IO Bool areCongruent solver t1 t2 = do@@ -170,13 +171,13 @@ lookup solver c1 c2 = do   tbl <- readIORef $ svLookupTable solver   return $ do-     m <- IM.lookup c1 tbl-     IM.lookup c2 m+     m <- IntMap.lookup c1 tbl+     IntMap.lookup c2 m  setLookup :: Solver -> Var -> Var -> Eqn1 -> IO () setLookup solver a1 a2 eqn = do   modifyIORef (svLookupTable solver) $-    IM.insertWith IM.union a1 (IM.singleton a2 eqn)+    IntMap.insertWith IntMap.union a1 (IntMap.singleton a2 eqn)  addToPending :: Solver -> PendingEqn -> IO () addToPending solver eqn = modifyIORef (svPending solver) (eqn :)@@ -184,20 +185,4 @@ getRepresentative :: Solver -> Var -> IO Var getRepresentative solver c = do   m <- readIORef $ svRepresentativeTable solver-  return $ m IM.! c--{---------------------------------------------------------------------  Test---------------------------------------------------------------------}--test = do-  solver <- newSolver-  a <- newVar solver-  b <- newVar solver-  c <- newVar solver-  d <- newVar solver-  merge solver (FTConst a, c)-  print =<< areCongruent solver (FTApp a b) (FTApp c d) -- False-  merge solver (FTConst b, d)-  print =<< areCongruent solver (FTApp a b) (FTApp c d) -- True-+  return $ m IntMap.! c

@@ -32,6 +32,7 @@ import Data.Function import qualified Data.IntMap as IM import qualified Data.IntSet as IS+import qualified Data.Map as Map import Data.List import Data.Monoid import Data.Ratio@@ -40,8 +41,9 @@ import Data.ArithRel import qualified Data.LA as LA import Data.OptDir-import Data.Polynomial-import Data.Polynomial.GBasis as GB+import Data.Polynomial (Polynomial, UPolynomial, Monomial, MonomialOrder)+import qualified Data.Polynomial as P+import Data.Polynomial.GroebnerBasis as GB import Data.Var import qualified Algorithm.LPUtil as LPUtil @@ -89,22 +91,22 @@     cmp2 = elimOrdering (IS.fromList vs2) `mappend` elimOrdering (IS.singleton t) `mappend` costOrdering obj `mappend` cmp      gb :: [Polynomial Rational Var]-    gb = GB.basis' GB.defaultOptions cmp2 (product (map var (t:vs2)) - 1 : phi)+    gb = GB.basis' GB.defaultOptions cmp2 (product (map P.var (t:vs2)) - 1 : phi)       where         phi = do           xj <- vs           let aj = [(yi, aij) | (yi,(ai,_)) <- zip vs2 cs, let aij = LA.coeff xj ai]-          return $  product [var yi ^ aij    | (yi, aij) <- aj, aij > 0]-                  - product [var yi ^ (-aij) | (yi, aij) <- aj, aij < 0] * var xj+          return $  product [P.var yi ^ aij    | (yi, aij) <- aj, aij > 0]+                  - product [P.var yi ^ (-aij) | (yi, aij) <- aj, aij < 0] * P.var xj -    yb = product [var yi ^ bi | ((_,bi),yi) <- zip cs vs2]+    yb = product [P.var yi ^ bi | ((_,bi),yi) <- zip cs vs2] -    [(_,z)] = terms (reduce cmp2 yb gb)+    [(_,z)] = P.terms (P.reduce cmp2 yb gb)      m = mkModel (vs++vs2++[t]) z -mkModel :: [Var] -> MonicMonomial Var -> Model Integer-mkModel vs xs = mmToIntMap xs `IM.union` IM.fromList [(x, 0) | x <- vs] +mkModel :: [Var] -> Monomial Var -> Model Integer+mkModel vs xs = IM.fromDistinctAscList (Map.toAscList (P.mindicesMap xs)) `IM.union` IM.fromList [(x, 0) | x <- vs] -- IM.union is left-biased  costOrdering :: LA.Expr Integer -> MonomialOrder Var@@ -116,4 +118,6 @@ elimOrdering :: IS.IntSet -> MonomialOrder Var elimOrdering xs = compare `on` f   where-    f ys = not $ IS.null $ xs `IS.intersection` IM.keysSet (mmToIntMap ys)+    f ys = not $ IS.null $ xs `IS.intersection` ys'+      where+        ys' = IS.fromDistinctAscList [y | (y,_) <- Map.toAscList $ P.mindicesMap ys]

@@ -50,7 +50,9 @@ import Data.IORef import Data.List import Data.Maybe+import Data.Map (Map) import qualified Data.Map as Map+import Data.Set (Set) import qualified Data.Set as Set import Text.Printf @@ -68,7 +70,7 @@ type PSym = String  class Vars a where-  vars :: a -> Set.Set Var+  vars :: a -> Set Var  instance Vars a => Vars [a] where   vars = Set.unions . map vars@@ -168,7 +170,7 @@ toSkolemNF :: forall m. Monad m => (String -> Int -> m FSym) -> Formula -> m [Clause] toSkolemNF skolem phi = f [] Map.empty (toNNF phi)   where-    f :: [Var] -> Map.Map Var Term -> Formula -> m [Clause]+    f :: [Var] -> Map Var Term -> Formula -> m [Clause]     f _ _ T = return []     f _ _ F = return [[]]     f _ s (Atom a) = return [[Pos (substAtom s a)]]@@ -189,15 +191,15 @@       f uvs (Map.insert v (TmApp fsym [TmVar v | v <- reverse uvs]) s) phi     f _ _ _ = error "toSkolemNF: should not happen" -    gensym :: String -> Set.Set Var -> Var+    gensym :: String -> Set Var -> Var     gensym template vs = head [name | name <- names, Set.notMember name vs]       where         names = template : [template ++ show n | n <-[1..]] -    substAtom :: Map.Map Var Term -> Atom -> Atom+    substAtom :: Map Var Term -> Atom -> Atom     substAtom s (PApp p ts) = PApp p (map (substTerm s) ts) -    substTerm :: Map.Map Var Term -> Term -> Term+    substTerm :: Map Var Term -> Term -> Term     substTerm s (TmVar v)    = fromMaybe (TmVar v) (Map.lookup v s)     substTerm s (TmApp f ts) = TmApp f (map (substTerm s) ts) @@ -255,7 +257,7 @@  -- --------------------------------------------------------------------------- -type M = State (Set.Set Var, Int, [SLit])+type M = State (Set Var, Int, [SLit])  flatten :: Clause -> SClause flatten c =@@ -359,9 +361,9 @@ type GroundLit    = GenLit GroundAtom type GroundClause = [GroundLit] -type Subst = Map.Map Var Entity+type Subst = Map Var Entity -enumSubst :: Set.Set Var -> [Entity] -> [Subst]+enumSubst :: Set Var -> [Entity] -> [Subst] enumSubst vs univ = do   ps <- sequence [[(v,e) | e <- univ] | v <- Set.toList vs]   return $ Map.fromList ps@@ -391,25 +393,25 @@     f (Pos (SEq (STmVar x) y)) = if x==y then Nothing else return []     f lit = return [lit] -collectFSym :: SClause -> Set.Set (FSym, Int)+collectFSym :: SClause -> Set (FSym, Int) collectFSym = Set.unions . map f   where-    f :: SLit -> Set.Set (FSym, Int)+    f :: SLit -> Set (FSym, Int)     f (Pos a) = g a     f (Neg a) = g a -    g :: SAtom -> Set.Set (FSym, Int)+    g :: SAtom -> Set (FSym, Int)     g (SEq (STmApp f xs) _) = Set.singleton (f, length xs)     g _ = Set.empty -collectPSym :: SClause -> Set.Set (PSym, Int)+collectPSym :: SClause -> Set (PSym, Int) collectPSym = Set.unions . map f   where-    f :: SLit -> Set.Set (PSym, Int)+    f :: SLit -> Set (PSym, Int)     f (Pos a) = g a     f (Neg a) = g a -    g :: SAtom -> Set.Set (PSym, Int)+    g :: SAtom -> Set (PSym, Int)     g (SPApp p xs) = Set.singleton (p, length xs)     g _ = Set.empty @@ -418,8 +420,8 @@ data Model   = Model   { mUniverse  :: [Entity]-  , mRelations :: Map.Map PSym [[Entity]]-  , mFunctions :: Map.Map FSym [([Entity], Entity)]+  , mRelations :: Map PSym [[Entity]]+  , mFunctions :: Map FSym [([Entity], Entity)]   }  showModel :: Model -> [String]

@@ -98,8 +98,10 @@ import Data.List import Data.Maybe import Data.Ratio+import Data.Map (Map) import qualified Data.Map as Map-import qualified Data.IntMap as IM+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap import Text.Printf import Data.Time import Data.OptDir@@ -120,15 +122,15 @@  data GenericSolver v   = GenericSolver-  { svTableau :: !(IORef (IM.IntMap (LA.Expr Rational)))-  , svLB      :: !(IORef (IM.IntMap v))-  , svUB      :: !(IORef (IM.IntMap v))-  , svModel   :: !(IORef (IM.IntMap v))+  { svTableau :: !(IORef (IntMap (LA.Expr Rational)))+  , svLB      :: !(IORef (IntMap v))+  , svUB      :: !(IORef (IntMap v))+  , svModel   :: !(IORef (IntMap v))   , svVCnt    :: !(IORef Int)   , svOk      :: !(IORef Bool)   , svOptDir  :: !(IORef OptDir) -  , svDefDB  :: !(IORef (Map.Map (LA.Expr Rational) Var))+  , svDefDB  :: !(IORef (Map (LA.Expr Rational) Var))    , svLogger :: !(IORef (Maybe (String -> IO ())))   , svPivotStrategy :: !(IORef PivotStrategy)@@ -143,10 +145,10 @@  newSolver :: SolverValue v => IO (GenericSolver v) newSolver = do-  t <- newIORef (IM.singleton objVar zeroV)-  l <- newIORef IM.empty-  u <- newIORef IM.empty-  m <- newIORef (IM.singleton objVar zeroV)+  t <- newIORef (IntMap.singleton objVar zeroV)+  l <- newIORef IntMap.empty+  u <- newIORef IntMap.empty+  m <- newIORef (IntMap.singleton objVar zeroV)   v <- newIORef 0   ok <- newIORef True   dir <- newIORef OptMin@@ -223,7 +225,7 @@         delta0 = if null ys then 0.1 else minimum ys         f :: Delta Rational -> Rational         f (Delta r k) = r + k * delta0-    liftM (IM.map f) $ readIORef (svModel solver)+    liftM (IntMap.map f) $ readIORef (svModel solver)  {- Largest coefficient rule: original rule suggested by G. Dantzig.@@ -250,7 +252,7 @@ newVar solver = do   v <- readIORef (svVCnt solver)   writeIORef (svVCnt solver) $! v+1-  modifyIORef (svModel solver) (IM.insert v zeroV)+  modifyIORef (svModel solver) (IntMap.insert v zeroV)   return v  assertAtom :: Solver -> LA.Atom Rational -> IO ()@@ -315,7 +317,7 @@     (Just l0', _) | l <= l0' -> return ()     (_, Just u0') | u0' < l -> markBad solver     _ -> do-      modifyIORef (svLB solver) (IM.insert x l)+      modifyIORef (svLB solver) (IntMap.insert x l)       b <- isNonBasicVariable solver x       v <- getValue solver x       when (b && not (l <= v)) $ update solver x l@@ -329,7 +331,7 @@     (_, Just u0') | u0' <= u -> return ()     (Just l0', _) | u < l0' -> markBad solver     _ -> do-      modifyIORef (svUB solver) (IM.insert x u)+      modifyIORef (svUB solver) (IntMap.insert x u)       b <- isNonBasicVariable solver x       v <- getValue solver x       when (b && not (v <= u)) $ update solver x u@@ -345,9 +347,9 @@ setRow :: SolverValue v => GenericSolver v -> Var -> LA.Expr Rational -> IO () setRow solver v e = do   modifyIORef (svTableau solver) $ \t ->-    IM.insert v (LA.applySubst t e) t+    IntMap.insert v (LA.applySubst t e) t   modifyIORef (svModel solver) $ \m -> -    IM.insert v (LA.evalLinear m (toValue 1) e) m  +    IntMap.insert v (LA.evalLinear m (toValue 1) e) m    setOptDir :: GenericSolver v -> OptDir -> IO () setOptDir solver dir = writeIORef (svOptDir solver) dir@@ -365,7 +367,7 @@ isBasicVariable :: GenericSolver v -> Var -> IO Bool isBasicVariable solver v = do   t <- readIORef (svTableau solver)-  return $ v `IM.member` t+  return $ v `IntMap.member` t  isNonBasicVariable  :: GenericSolver v -> Var -> IO Bool isNonBasicVariable solver x = liftM not (isBasicVariable solver x)@@ -568,7 +570,7 @@    -- Upper bounds of θ   -- NOTE: xj 自体の上限も考慮するのに注意-  ubs <- liftM concat $ forM ((xj,1) : IM.toList col) $ \(xi,aij) -> do+  ubs <- liftM concat $ forM ((xj,1) : IntMap.toList col) $ \(xi,aij) -> do     v1 <- getValue solver xi     li <- getLB solver xi     ui <- getUB solver xi@@ -591,7 +593,7 @@    -- Lower bounds of θ   -- NOTE: xj 自体の下限も考慮するのに注意-  lbs <- liftM concat $ forM ((xj,1) : IM.toList col) $ \(xi,aij) -> do+  lbs <- liftM concat $ forM ((xj,1) : IntMap.toList col) $ \(xi,aij) -> do     v1 <- getValue solver xi     li <- getLB solver xi     ui <- getUB solver xi@@ -690,19 +692,19 @@   Extract results --------------------------------------------------------------------} -type RawModel v = IM.IntMap v+type RawModel v = IntMap v  rawModel :: GenericSolver v -> IO (RawModel v) rawModel solver = do   xs <- variables solver-  liftM IM.fromList $ forM xs $ \x -> do+  liftM IntMap.fromList $ forM xs $ \x -> do     val <- getValue solver x     return (x,val)  getObjValue :: GenericSolver v -> IO v getObjValue solver = getValue solver objVar   -type Model = IM.IntMap Rational+type Model = IntMap Rational    {--------------------------------------------------------------------   major function@@ -718,8 +720,8 @@    aj <- getCol solver xj   modifyIORef (svModel solver) $ \m ->-    let m2 = IM.map (\aij -> aij *^ diff) aj-    in IM.insert xj v $ IM.unionWith (^+^) m2 m+    let m2 = IntMap.map (\aij -> aij *^ diff) aj+    in IntMap.insert xj v $ IntMap.unionWith (^+^) m2 m    -- log solver $ printf "after update x%d (%s)" xj (show v)   -- dump solver@@ -728,9 +730,9 @@ pivot solver xi xj = do   modifyIORef' (svNPivot solver) (+1)   modifyIORef' (svTableau solver) $ \defs ->-    case LA.solveFor (LA.var xi .==. (defs IM.! xi)) xj of+    case LA.solveFor (LA.var xi .==. (defs IntMap.! xi)) xj of       Just (Eql, xj_def) ->-        IM.insert xj xj_def . IM.map (LA.applySubst1 xj xj_def) . IM.delete xi $ defs+        IntMap.insert xj xj_def . IntMap.map (LA.applySubst1 xj xj_def) . IntMap.delete xi $ defs       _ -> error "pivot: should not happen"  pivotAndUpdate :: SolverValue v => GenericSolver v -> Var -> Var -> v -> IO ()@@ -745,13 +747,13 @@   m <- readIORef (svModel solver)    aj <- getCol solver xj-  let aij = aj IM.! xi-  let theta = (v ^-^ (m IM.! xi)) ^/ aij+  let aij = aj IntMap.! xi+  let theta = (v ^-^ (m IntMap.! xi)) ^/ aij -  let m' = IM.fromList $-           [(xi, v), (xj, (m IM.! xj) ^+^ theta)] ++-           [(xk, (m IM.! xk) ^+^ (akj *^ theta)) | (xk, akj) <- IM.toList aj, xk /= xi]-  writeIORef (svModel solver) (IM.union m' m) -- note that 'IM.union' is left biased.+  let m' = IntMap.fromList $+           [(xi, v), (xj, (m IntMap.! xj) ^+^ theta)] +++           [(xk, (m IntMap.! xk) ^+^ (akj *^ theta)) | (xk, akj) <- IntMap.toList aj, xk /= xi]+  writeIORef (svModel solver) (IntMap.union m' m) -- note that 'IntMap.union' is left biased.    pivot solver xi xj @@ -761,34 +763,34 @@ getLB :: GenericSolver v -> Var -> IO (Maybe v) getLB solver x = do   lb <- readIORef (svLB solver)-  return $ IM.lookup x lb+  return $ IntMap.lookup x lb  getUB :: GenericSolver v -> Var -> IO (Maybe v) getUB solver x = do   ub <- readIORef (svUB solver)-  return $ IM.lookup x ub+  return $ IntMap.lookup x ub -getTableau :: GenericSolver v -> IO (IM.IntMap (LA.Expr Rational))+getTableau :: GenericSolver v -> IO (IntMap (LA.Expr Rational)) getTableau solver = do   t <- readIORef (svTableau solver)-  return $ IM.delete objVar t+  return $ IntMap.delete objVar t  getValue :: GenericSolver v -> Var -> IO v getValue solver x = do   m <- readIORef (svModel solver)-  return $ m IM.! x+  return $ m IntMap.! x  getRow :: GenericSolver v -> Var -> IO (LA.Expr Rational) getRow solver x = do   -- x should be basic variable or 'objVar'   t <- readIORef (svTableau solver)-  return $! (t IM.! x)+  return $! (t IntMap.! x)  -- aijが非ゼロの列も全部探しているのは効率が悪い-getCol :: SolverValue v => GenericSolver v -> Var -> IO (IM.IntMap Rational)+getCol :: SolverValue v => GenericSolver v -> Var -> IO (IntMap Rational) getCol solver xj = do   t <- readIORef (svTableau solver)-  return $ IM.mapMaybe (LA.lookupCoeff xj) t+  return $ IntMap.mapMaybe (LA.lookupCoeff xj) t  getCoeff :: GenericSolver v -> Var -> Var -> IO Rational getCoeff solver xi xj = do@@ -826,7 +828,7 @@ basicVariables :: GenericSolver v -> IO [Var] basicVariables solver = do   t <- readIORef (svTableau solver)-  return (IM.keys t)+  return (IntMap.keys t)  #if !MIN_VERSION_base(4,6,0) @@ -900,9 +902,9 @@   x0 <- newVar solver   x1 <- newVar solver -  writeIORef (svTableau solver) (IM.fromList [(x1, LA.var x0)])-  writeIORef (svLB solver) (IM.fromList [(x0, toValue 0), (x1, toValue 0)])-  writeIORef (svUB solver) (IM.fromList [(x0, toValue 2), (x1, toValue 3)])+  writeIORef (svTableau solver) (IntMap.fromList [(x1, LA.var x0)])+  writeIORef (svLB solver) (IntMap.fromList [(x0, toValue 0), (x1, toValue 0)])+  writeIORef (svUB solver) (IntMap.fromList [(x0, toValue 2), (x1, toValue 3)])   setObj solver (LA.fromTerms [(-1, x0)])    ret <- optimize solver defaultOptions@@ -916,9 +918,9 @@   x0 <- newVar solver   x1 <- newVar solver -  writeIORef (svTableau solver) (IM.fromList [(x1, LA.var x0)])-  writeIORef (svLB solver) (IM.fromList [(x0, toValue 0), (x1, toValue 0)])-  writeIORef (svUB solver) (IM.fromList [(x0, toValue 2), (x1, toValue 0)])+  writeIORef (svTableau solver) (IntMap.fromList [(x1, LA.var x0)])+  writeIORef (svLB solver) (IntMap.fromList [(x0, toValue 0), (x1, toValue 0)])+  writeIORef (svUB solver) (IntMap.fromList [(x0, toValue 2), (x1, toValue 0)])   setObj solver (LA.fromTerms [(-1, x0)])    checkFeasibility solver@@ -946,10 +948,10 @@ dumpSize :: SolverValue v => GenericSolver v -> IO () dumpSize solver = do   t <- readIORef (svTableau solver)-  let nrows = IM.size t - 1 -- -1 is objVar+  let nrows = IntMap.size t - 1 -- -1 is objVar   xs <- variables solver   let ncols = length xs - nrows-  let nnz = sum [IM.size $ LA.coeffMap xi_def | (xi,xi_def) <- IM.toList t, xi /= objVar]+  let nnz = sum [IntMap.size $ LA.coeffMap xi_def | (xi,xi_def) <- IntMap.toList t, xi /= objVar]   log solver $ printf "%d rows, %d columns, %d non-zeros" nrows ncols nnz  dump :: SolverValue v => GenericSolver v -> IO ()@@ -958,8 +960,8 @@    log solver "Tableau:"   t <- readIORef (svTableau solver)-  log solver $ printf "obj = %s" (show (t IM.! objVar))-  forM_ (IM.toList t) $ \(xi, e) -> do+  log solver $ printf "obj = %s" (show (t IntMap.! objVar))+  forM_ (IntMap.toList t) $ \(xi, e) -> do     when (xi /= objVar) $ log solver $ printf "x%d = %s" xi (show e)    log solver ""

@@ -22,6 +22,7 @@ import Data.List import Data.Ratio import qualified Data.Set as Set+import Data.Map (Map) import qualified Data.Map as Map import System.IO import Text.Printf@@ -54,7 +55,7 @@ -- ------------------------------------------------------------------------  type Var = String-type Env = Map.Map LP.Var Var+type Env = Map LP.Var Var  concatS :: [ShowS] -> ShowS concatS = foldr (.) id

@@ -14,12 +14,12 @@   ( convert   ) where -import qualified Data.Map as Map+import Data.Map (Map) import qualified Text.LPFile as LPFile import qualified Text.MaxSAT as MaxSAT import SAT.Types import qualified Converter.MaxSAT2WBO as MaxSAT2WBO import qualified Converter.PB2LP as PB2LP -convert :: Bool -> MaxSAT.WCNF -> (LPFile.LP, Map.Map LPFile.Var Rational -> Model)+convert :: Bool -> MaxSAT.WCNF -> (LPFile.LP, Map LPFile.Var Rational -> Model) convert useIndicator wcnf = PB2LP.convertWBO useIndicator (MaxSAT2WBO.convert wcnf)

@@ -18,14 +18,16 @@ import Data.Array.IArray import Data.List import Data.Maybe-import qualified Data.IntSet as IS+import Data.IntSet (IntSet)+import qualified Data.IntSet as IntSet import qualified Data.Set as Set+import Data.Map (Map) import qualified Data.Map as Map import qualified Text.PBFile as PBFile import qualified Text.LPFile as LPFile import qualified SAT.Types as SAT -convert :: PBFile.Formula -> (LPFile.LP, Map.Map LPFile.Var Rational -> SAT.Model)+convert :: PBFile.Formula -> (LPFile.LP, Map LPFile.Var Rational -> SAT.Model) convert formula@(obj, cs) = (lp, mtrans (PBFile.pbNumVars formula))   where     lp = LPFile.LP@@ -47,7 +49,7 @@       }      vs1 = collectVariables formula-    vs2 = (Set.fromList . map convVar . IS.toList) $ vs1+    vs2 = (Set.fromList . map convVar . IntSet.toList) $ vs1      (dir,obj2) =       case obj of@@ -90,16 +92,16 @@ convVar :: PBFile.Var -> LPFile.Var convVar x = ("x" ++ show x) -collectVariables :: PBFile.Formula -> IS.IntSet-collectVariables (obj, cs) = IS.unions $ maybe IS.empty f obj : [f s | (s,_,_) <- cs]+collectVariables :: PBFile.Formula -> IntSet+collectVariables (obj, cs) = IntSet.unions $ maybe IntSet.empty f obj : [f s | (s,_,_) <- cs]   where-    f :: PBFile.Sum -> IS.IntSet-    f xs = IS.fromList $ do+    f :: PBFile.Sum -> IntSet+    f xs = IntSet.fromList $ do       (_,ts) <- xs       lit <- ts       return $ abs lit -convertWBO :: Bool -> PBFile.SoftFormula -> (LPFile.LP, Map.Map LPFile.Var Rational -> SAT.Model)+convertWBO :: Bool -> PBFile.SoftFormula -> (LPFile.LP, Map LPFile.Var Rational -> SAT.Model) convertWBO useIndicator formula@(top, cs) = (lp, mtrans (PBFile.wboNumVars formula))   where     lp = LPFile.LP@@ -121,7 +123,7 @@       }      vs1 = collectVariablesWBO formula-    vs2 = ((Set.fromList . map convVar . IS.toList) $ vs1) `Set.union` vs3+    vs2 = ((Set.fromList . map convVar . IntSet.toList) $ vs1) `Set.union` vs3     vs3 = Set.fromList [v | (ts, _) <- cs2, (_, v) <- ts]      obj2 = [LPFile.Term (fromIntegral w) [v] | (ts, _) <- cs2, (w, v) <- ts]@@ -194,16 +196,16 @@   where     lhs_ub = sum [max c 0 | LPFile.Term c _ <- lhs] -collectVariablesWBO :: PBFile.SoftFormula -> IS.IntSet-collectVariablesWBO (_top, cs) = IS.unions [f s | (_,(s,_,_)) <- cs]+collectVariablesWBO :: PBFile.SoftFormula -> IntSet+collectVariablesWBO (_top, cs) = IntSet.unions [f s | (_,(s,_,_)) <- cs]   where-    f :: PBFile.Sum -> IS.IntSet-    f xs = IS.fromList $ do+    f :: PBFile.Sum -> IntSet+    f xs = IntSet.fromList $ do       (_,ts) <- xs       lit <- ts       return $ abs lit -mtrans :: Int -> Map.Map LPFile.Var Rational -> SAT.Model+mtrans :: Int -> Map LPFile.Var Rational -> SAT.Model mtrans nvar m =   array (1, nvar)     [　(i, val)

@@ -14,12 +14,12 @@   ( convert   ) where -import qualified Data.Map as Map+import Data.Map (Map) import qualified Text.LPFile as LPFile import qualified Language.CNF.Parse.ParseDIMACS as DIMACS import qualified SAT.Types as SAT import qualified Converter.PB2LP as PB2LP import qualified Converter.SAT2PB as SAT2PB -convert :: DIMACS.CNF -> (LPFile.LP, Map.Map LPFile.Var Rational -> SAT.Model)+convert :: DIMACS.CNF -> (LPFile.LP, Map LPFile.Var Rational -> SAT.Model) convert cnf = PB2LP.convert (SAT2PB.convert cnf)

@@ -28,7 +28,6 @@    -- * Properties   , minimalPolynomial-  , deg   , isRational   , isAlgebraicInteger   , height@@ -54,9 +53,8 @@ import qualified Text.PrettyPrint.HughesPJClass as PP import Text.PrettyPrint.HughesPJClass (Doc, PrettyLevel, Pretty (..), prettyParen) -import Data.Polynomial+import Data.Polynomial (Polynomial, UPolynomial, X (..)) import qualified Data.Polynomial as P-import qualified Data.Polynomial.Factorization.Rational as FactorQ import qualified Data.Polynomial.RootSeparation.Sturm as Sturm import Data.Interval (Interval, EndPoint (..), (<=..<), (<..<=), (<..<), (<!), (>!)) import qualified Data.Interval as Interval@@ -73,25 +71,25 @@ -- | Real roots of the polynomial in ascending order. realRoots :: UPolynomial Rational -> [AReal] realRoots p = Set.toAscList $ Set.fromList $ do-  (q,_) <- FactorQ.factor p+  (q,_) <- P.factor p   realRoots' q  -- | Real roots of the polynomial in ascending order. realRootsEx :: UPolynomial AReal -> [AReal] realRootsEx p-  | and [isRational c | (c,_) <- terms p] = realRoots $ mapCoeff toRational p-  | otherwise = [a | a <- realRoots (simpARealPoly p), a `isRootOf` p]+  | and [isRational c | (c,_) <- P.terms p] = realRoots $ P.mapCoeff toRational p+  | otherwise = [a | a <- realRoots (simpARealPoly p), a `P.isRootOf` p]  -- p must already be factored. realRoots' :: UPolynomial Rational -> [AReal] realRoots' p = do-  guard $ deg p > 0+  guard $ P.deg p > 0   i <- Sturm.separate p   return $ realRoot' p i  realRoot :: UPolynomial Rational -> Interval Rational -> AReal realRoot p i = -  case [q | (q,_) <- FactorQ.factor p, deg q > 0, Sturm.numRoots q i == 1] of+  case [q | (q,_) <- P.factor p, P.deg q > 0, Sturm.numRoots q i == 1] of     p2:_ -> realRoot' p2 i     []   -> error "Data.AlgebraicNumber.Real.realRoot: invalid interval" @@ -104,7 +102,7 @@ --------------------------------------------------------------------}  isZero :: AReal -> Bool-isZero a = 0 `Interval.member` (interval a) && 0 `isRootOf` minimalPolynomial a+isZero a = 0 `Interval.member` (interval a) && 0 `P.isRootOf` minimalPolynomial a  scaleAReal :: Rational -> AReal -> AReal scaleAReal r a = realRoot' p2 i2@@ -202,23 +200,30 @@   fromInteger = fromRational . toRational  instance Fractional AReal where-  fromRational r = realRoot' (x - constant r) (Interval.singleton r)+  fromRational r = realRoot' (x - P.constant r) (Interval.singleton r)     where-      x = var X+      x = P.var X    recip a     | isZero a  = error "AReal.recip: zero division"     | otherwise = realRoot' p2 i2       where         p2 = rootRecip (minimalPolynomial a)-        i2 = recip (interval a)+        c1 = sturmChain a+        c2 = Sturm.sturmChain p2+        i2 = go (interval a) (Sturm.separate' c2)+        go i1 is2 =+          case [i2 | i2 <- is2, Interval.member 1 (i1 * i2)] of+            [] -> error "AReal.recip: should not happen"+            [i2] -> i2+            is2'  -> go (Sturm.halve' c1 i1) [Sturm.halve' c2 i2 | i2 <- is2']  instance Real AReal where   toRational x     | isRational x =         let p = minimalPolynomial x             a = P.coeff (P.var X) p-            b = P.coeff P.mmOne p+            b = P.coeff P.mone p         in - b / a     | otherwise  = error "toRational: proper algebraic number" @@ -365,29 +370,25 @@ --  -- If the algebraic number's 'minimalPolynomial' has degree @n@, -- then the algebraic number is said to be degree @n@.-instance Degree AReal where-  deg a = deg $ minimalPolynomial a+instance P.Degree AReal where+  deg a = P.deg $ minimalPolynomial a  -- | Whether the algebraic number is a rational. isRational :: AReal -> Bool-isRational x = deg x == 1+isRational x = P.deg x == 1  -- | Whether the algebraic number is a root of a polynomial with integer -- coefficients with leading coefficient @1@ (a monic polynomial). isAlgebraicInteger :: AReal -> Bool-isAlgebraicInteger x = cn * fromIntegral d == 1-  where-    p = minimalPolynomial x-    d = foldl' lcm 1 [denominator c | (c,_) <- terms p]-    (cn,_) = leadingTerm grlex p+isAlgebraicInteger x = abs (P.lc P.grlex (P.pp (minimalPolynomial x))) == 1  -- | Height of the algebraic number.+--+-- The height of an algebraic number is the greatest absolute value of the+-- coefficients of the irreducible and primitive polynomial with integral+-- rational coefficients. height :: AReal -> Integer-height x = maximum [ assert (denominator c' == 1) (abs (numerator c'))-                   | (c,_) <- terms p, let c' = c * fromInteger d ]-  where-    p = minimalPolynomial x-    d = foldl' lcm 1 [denominator c | (c,_) <- terms p]+height x = maximum [abs (numerator c) | (c,_) <- P.terms $ P.pp $ minimalPolynomial x]  -- | root index, satisfying --@@ -412,7 +413,7 @@       p = minimalPolynomial r       appPrec = 10 -instance PrettyCoeff AReal where+instance P.PrettyCoeff AReal where   pPrintCoeff = pPrintPrec   isNegativeCoeff = (0>) @@ -432,4 +433,4 @@ goldenRatio :: AReal goldenRatio = (1 + root5) / 2   where-    [_, root5] = sort $ realRoots' ((var X)^2 - 5)+    [_, root5] = sort $ realRoots' ((P.var X)^2 - 5)

@@ -20,11 +20,13 @@  import Data.List import Data.Maybe+import Data.Map (Map) import qualified Data.Map as Map import qualified Data.Set as Set -import Data.Polynomial-import qualified Data.Polynomial.GBasis as GB+import Data.Polynomial (Polynomial, UPolynomial, X (..))+import qualified Data.Polynomial as P+import qualified Data.Polynomial.GroebnerBasis as GB  type Var = Int @@ -33,11 +35,7 @@ --------------------------------------------------------------------}  normalizePoly :: UPolynomial Rational -> UPolynomial Rational-normalizePoly p-  | c == 1    = p-  | otherwise = mapCoeff (/ c) p-  where-    (c,_) = leadingTerm grlex p+normalizePoly = P.toMonic P.grlex  rootAdd :: UPolynomial Rational -> UPolynomial Rational -> UPolynomial Rational rootAdd p1 p2 = lift2 (+) p1 p2@@ -47,35 +45,35 @@  rootShift :: Rational -> UPolynomial Rational -> UPolynomial Rational rootShift 0 p = p-rootShift r p = normalizePoly $ subst p (\X -> var X - constant r)+rootShift r p = normalizePoly $ P.subst p (\X -> P.var X - P.constant r)  rootScale :: Rational -> UPolynomial Rational -> UPolynomial Rational-rootScale 0 p = var X-rootScale r p = normalizePoly $ subst p (\X -> constant (recip r) * var X)+rootScale 0 p = P.var X+rootScale r p = normalizePoly $ P.subst p (\X -> P.constant (recip r) * P.var X)  rootRecip :: UPolynomial Rational -> UPolynomial Rational-rootRecip p = normalizePoly $ fromTerms [(c, mmFromList [(X, d - deg xs)]) | (c, xs) <- terms p]+rootRecip p = normalizePoly $ P.fromTerms [(c, P.var X `P.mpow` (d - P.deg xs)) | (c, xs) <- P.terms p]   where-    d = deg p+    d = P.deg p  -- 代数的数を係数とする多項式の根を、有理数係数多項式の根として表す rootSimpPoly :: (a -> UPolynomial Rational) -> UPolynomial a -> UPolynomial Rational-rootSimpPoly f p = findPoly (var 0) ps+rootSimpPoly f p = findPoly (P.var 0) ps   where     ys :: [(UPolynomial Rational, Var)]-    ys = zip (Set.toAscList $ Set.fromList [f c | (c, _) <- terms p]) [1..]+    ys = zip (Set.toAscList $ Set.fromList [f c | (c, _) <- P.terms p]) [1..] -    m :: Map.Map (UPolynomial Rational) Var+    m :: Map (UPolynomial Rational) Var     m = Map.fromDistinctAscList ys      p' :: Polynomial Rational Var-    p' = eval (\X -> var 0) (mapCoeff (\c -> var (m Map.! (f c))) p)+    p' = P.eval (\X -> P.var 0) (P.mapCoeff (\c -> P.var (m Map.! (f c))) p)      ps :: [Polynomial Rational Var]-    ps = p' : [subst q (\X -> var x) | (q, x) <- ys]+    ps = p' : [P.subst q (\X -> P.var x) | (q, x) <- ys]  rootNthRoot :: Integer -> UPolynomial Rational -> UPolynomial Rational-rootNthRoot n p = subst p (\X -> (var X)^n)+rootNthRoot n p = P.subst p (\X -> (P.var X)^n)  lift2 :: (forall a. Num a => a -> a -> a)       -> UPolynomial Rational -> UPolynomial Rational -> UPolynomial Rational@@ -86,37 +84,37 @@     b = 1      f_a_b :: Polynomial Rational Var-    f_a_b = f (var a) (var b)+    f_a_b = f (P.var a) (P.var b)      gbase :: [Polynomial Rational Var]-    gbase = [ subst p1 (\X -> var a), subst p2 (\X -> var b) ]              +    gbase = [ P.subst p1 (\X -> P.var a), P.subst p2 (\X -> P.var b) ]                -- ps のもとで c を根とする多項式を求める findPoly :: Polynomial Rational Var -> [Polynomial Rational Var] -> UPolynomial Rational-findPoly c ps = normalizePoly $ sum [constant coeff * (var X) ^ n | (n,coeff) <- zip [0..] coeffs]+findPoly c ps = normalizePoly $ sum [P.constant coeff * (P.var X) ^ n | (n,coeff) <- zip [0..] coeffs]   where       vn :: Var     vn = if Set.null vs then 0 else Set.findMax vs + 1       where-        vs = Set.fromList [x | p <- (c:ps), (_,xs) <- terms p, (x,_) <- mmToList xs]+        vs = Set.fromList [x | p <- (c:ps), (_,xs) <- P.terms p, (x,_) <- P.mindices xs]      coeffs :: [Rational]     coeffs = head $ catMaybes $ [isLinearlyDependent cs2 | cs2 <- inits cs]       where-        cmp = grlex+        cmp = P.grlex         ps' = GB.basis cmp ps-        cs  = iterate (\p -> reduce cmp (c * p) ps') 1+        cs  = iterate (\p -> P.reduce cmp (c * p) ps') 1      isLinearlyDependent :: [Polynomial Rational Var] -> Maybe [Rational]     isLinearlyDependent cs = if any (0/=) sol then Just sol else Nothing       where         cs2 = zip [vn..] cs-        sol = map (\(l,_) -> eval (\_ -> 1) $ reduce cmp2 (var l) gbase2) cs2-        cmp2   = grlex+        sol = map (\(l,_) -> P.eval (\_ -> 1) $ P.reduce cmp2 (P.var l) gbase2) cs2+        cmp2   = P.grlex         gbase2 = GB.basis cmp2 es         es = Map.elems $ Map.fromListWith (+) $ do-          (n,xs) <- terms $ sum [var ln * cn | (ln,cn) <- cs2]-          let xs' = mmToList xs-              ys = mmFromList [(x,m) | (x,m) <- xs', x < vn]-              zs = mmFromList [(x,m) | (x,m) <- xs', x >= vn]-          return (ys, fromMonomial (n,zs))+          (n,xs) <- P.terms $ sum [P.var ln * cn | (ln,cn) <- cs2]+          let xs' = P.mindicesMap xs+              ys = P.mfromIndicesMap $ Map.filterWithKey (\x _ -> x <  vn) xs'+              zs = P.mfromIndicesMap $ Map.filterWithKey (\x _ -> x >= vn) xs'+          return (ys, P.fromTerm (n,zs))

@@ -58,8 +58,9 @@ import Control.DeepSeq import Data.List import Data.Maybe-import qualified Data.IntMap as IM-import qualified Data.IntSet as IS+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap+import qualified Data.IntSet as IntSet import qualified Data.ArithRel as ArithRel import Data.Interval import Data.Var@@ -73,23 +74,23 @@ newtype Expr r   = Expr   { -- | a mapping from variables to coefficients-    coeffMap :: IM.IntMap r+    coeffMap :: IntMap r   } deriving (Eq, Ord)  -- | Create a @Expr@ from a mapping from variables to coefficients.-fromCoeffMap :: (Num r, Eq r) => IM.IntMap r -> Expr r+fromCoeffMap :: (Num r, Eq r) => IntMap r -> Expr r fromCoeffMap m = normalizeExpr (Expr m)  -- | terms contained in the expression. terms :: Expr r -> [(r,Var)]-terms (Expr m) = [(c,v) | (v,c) <- IM.toList m]+terms (Expr m) = [(c,v) | (v,c) <- IntMap.toList m]  -- | Create a @Expr@ from a list of terms. fromTerms :: (Num r, Eq r) => [(r,Var)] -> Expr r-fromTerms ts = fromCoeffMap $ IM.fromListWith (+) [(x,c) | (c,x) <- ts]+fromTerms ts = fromCoeffMap $ IntMap.fromListWith (+) [(x,c) | (c,x) <- ts]  instance Variables (Expr r) where-  vars (Expr m) = IS.delete unitVar (IM.keysSet m)+  vars (Expr m) = IntSet.delete unitVar (IntMap.keysSet m)  instance Show r => Show (Expr r) where   showsPrec d m  = showParen (d > 10) $@@ -110,41 +111,41 @@  asConst :: Num r => Expr r -> Maybe r asConst (Expr m) =-  case IM.toList m of+  case IntMap.toList m of     [] -> Just 0     [(v,x)] | v==unitVar -> Just x     _ -> Nothing  normalizeExpr :: (Num r, Eq r) => Expr r -> Expr r-normalizeExpr (Expr t) = Expr $ IM.filter (0/=) t+normalizeExpr (Expr t) = Expr $ IntMap.filter (0/=) t  -- | variable var :: Num r => Var -> Expr r-var v = Expr $ IM.singleton v 1+var v = Expr $ IntMap.singleton v 1  -- | constant constant :: (Num r, Eq r) => r -> Expr r-constant c = normalizeExpr $ Expr $ IM.singleton unitVar c+constant c = normalizeExpr $ Expr $ IntMap.singleton unitVar c  -- | map coefficients. mapCoeff :: (Num b, Eq b) => (a -> b) -> Expr a -> Expr b-mapCoeff f (Expr t) = Expr $ IM.mapMaybe g t+mapCoeff f (Expr t) = Expr $ IntMap.mapMaybe g t   where     g c = if c' == 0 then Nothing else Just c'       where c' = f c  -- | map coefficients. mapCoeffWithVar :: (Num b, Eq b) => (a -> Var -> b) -> Expr a -> Expr b-mapCoeffWithVar f (Expr t) = Expr $ IM.mapMaybeWithKey g t+mapCoeffWithVar f (Expr t) = Expr $ IntMap.mapMaybeWithKey g t   where     g v c = if c' == 0 then Nothing else Just c'       where c' = f c v  instance (Num r, Eq r) => AdditiveGroup (Expr r) where-  Expr t ^+^ e2 | IM.null t = e2-  e1 ^+^ Expr t | IM.null t = e1+  Expr t ^+^ e2 | IntMap.null t = e2+  e1 ^+^ Expr t | IntMap.null t = e1   e1 ^+^ e2 = normalizeExpr $ plus e1 e2-  zeroV = Expr $ IM.empty+  zeroV = Expr $ IntMap.empty   negateV = ((-1) *^)  instance (Num r, Eq r) => VectorSpace (Expr r) where@@ -154,30 +155,30 @@   c *^ e = mapCoeff (c*) e  plus :: Num r => Expr r -> Expr r -> Expr r-plus (Expr t1) (Expr t2) = Expr $ IM.unionWith (+) t1 t2+plus (Expr t1) (Expr t2) = Expr $ IntMap.unionWith (+) t1 t2  -- | evaluate the expression under the model. evalExpr :: Num r => Model r -> Expr r -> r-evalExpr m (Expr t) = sum [(m' IM.! v) * c | (v,c) <- IM.toList t]-  where m' = IM.insert unitVar 1 m+evalExpr m (Expr t) = sum [(m' IntMap.! v) * c | (v,c) <- IntMap.toList t]+  where m' = IntMap.insert unitVar 1 m  -- | evaluate the expression under the model. evalLinear :: VectorSpace a => Model a -> a -> Expr (Scalar a) -> a-evalLinear m u (Expr t) = sumV [c *^ (m' IM.! v) | (v,c) <- IM.toList t]-  where m' = IM.insert unitVar u m+evalLinear m u (Expr t) = sumV [c *^ (m' IntMap.! v) | (v,c) <- IntMap.toList t]+  where m' = IntMap.insert unitVar u m  lift1 :: VectorSpace x => x -> (Var -> x) -> Expr (Scalar x) -> x-lift1 unit f (Expr t) = sumV [c *^ (g v) | (v,c) <- IM.toList t]+lift1 unit f (Expr t) = sumV [c *^ (g v) | (v,c) <- IntMap.toList t]   where     g v       | v==unitVar = unit       | otherwise   = f v  applySubst :: (Num r, Eq r) => VarMap (Expr r) -> Expr r -> Expr r-applySubst s (Expr m) = sumV (map f (IM.toList m))+applySubst s (Expr m) = sumV (map f (IntMap.toList m))   where     f (v,c) = c *^ (-      case IM.lookup v s of+      case IntMap.lookup v s of         Just tm -> tm         Nothing -> var v) @@ -193,7 +194,7 @@ --   coeff v e == fst (extract v e) -- @ coeff :: Num r => Var -> Expr r -> r-coeff v (Expr m) = IM.findWithDefault 0 v m+coeff v (Expr m) = IntMap.findWithDefault 0 v m  -- | lookup a coefficient of the variable. -- It returns @Nothing@ if the expression does not contain @v@.@@ -201,15 +202,15 @@ --   lookupCoeff v e == fmap fst (extractMaybe v e) -- @ lookupCoeff :: Num r => Var -> Expr r -> Maybe r-lookupCoeff v (Expr m) = IM.lookup v m  +lookupCoeff v (Expr m) = IntMap.lookup v m    -- | @extract v e@ returns @(c, e')@ such that @e == c *^ v ^+^ e'@ extract :: Num r => Var -> Expr r -> (r, Expr r)-extract v (Expr m) = (IM.findWithDefault 0 v m, Expr (IM.delete v m))+extract v (Expr m) = (IntMap.findWithDefault 0 v m, Expr (IntMap.delete v m)) {- -- Alternative implementation which may be faster but allocte more memory extract v (Expr m) = -  case IM.updateLookupWithKey (\_ _ -> Nothing) v m of+  case IntMap.updateLookupWithKey (\_ _ -> Nothing) v m of     (Nothing, _) -> (0, Expr m)     (Just c, m2) -> (c, Expr m2) -}@@ -218,13 +219,13 @@ -- if @e@ contains v, and returns @Nothing@ otherwise. extractMaybe :: Num r => Var -> Expr r -> Maybe (r, Expr r) extractMaybe v (Expr m) =-  case IM.lookup v m of+  case IntMap.lookup v m of     Nothing -> Nothing-    Just c -> Just (c, Expr (IM.delete v m))+    Just c -> Just (c, Expr (IntMap.delete v m)) {- -- Alternative implementation which may be faster but allocte more memory extractMaybe v (Expr m) =-  case IM.updateLookupWithKey (\_ _ -> Nothing) v m of+  case IntMap.updateLookupWithKey (\_ _ -> Nothing) v m of     (Nothing, _) -> Nothing     (Just c, m2) -> Just (c, Expr m2) -}@@ -241,8 +242,8 @@     ts = [if c==1             then showString (env x)             else showsPrec 8 c . showString "*" . showString (env x)-          | (x,c) <- IM.toList m, x /= unitVar] ++-         [showsPrec 7 c | c <- maybeToList (IM.lookup unitVar m)]+          | (x,c) <- IntMap.toList m, x /= unitVar] +++         [showsPrec 7 c | c <- maybeToList (IntMap.lookup unitVar m)]  ----------------------------------------------------------------------------- 

@@ -22,7 +22,8 @@  import Data.List import Data.Ratio-import qualified Data.IntSet as IS+import qualified Data.IntSet as IntSet+import Data.Map (Map) import qualified Data.Map as Map import Data.VectorSpace import Prelude hiding (null)@@ -43,13 +44,13 @@  -- | Intersection of half-spaces data Polyhedron-  = Polyhedron (Map.Map ExprZ IntervalR)+  = Polyhedron (Map ExprZ IntervalR)   | Empty   deriving (Eq)  instance Variables Polyhedron where-  vars (Polyhedron m) = IS.unions [vars e | e <- Map.keys m]-  vars Empty = IS.empty+  vars (Polyhedron m) = IntSet.unions [vars e | e <- Map.keys m]+  vars Empty = IntSet.empty  instance JoinSemiLattice Polyhedron where   join Empty b = b

@@ -1,4 +1,4 @@-{-# LANGUAGE ScopedTypeVariables, TypeFamilies, BangPatterns, DeriveDataTypeable #-}+{-# OPTIONS_GHC -Wall #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Polynomial@@ -7,7 +7,7 @@ --  -- Maintainer  :  masahiro.sakai@gmail.com -- Stability   :  provisional--- Portability :  non-portable (ScopedTypeVariables, TypeFamilies, BangPatterns, DeriveDataTypeable)+-- Portability :  portable -- -- Polynomials --@@ -24,76 +24,84 @@   (   -- * Polynomial type     Polynomial-  , UPolynomial-  , X (..)    -- * Conversion-  , Variables (..)+  , Var (..)   , constant   , terms   , fromTerms   , coeffMap   , fromCoeffMap-  , fromMonomial+  , fromTerm    -- * Query   , Degree (..)-  , leadingTerm+  , Vars (..)+  , lt+  , lc+  , lm   , coeff   , lookupCoeff   , isPrimitive+  , isRootOf    -- * Operations+  , Factor (..)+  , SQFree (..)   , ContPP (..)   , deriv   , integral   , eval-  , evalA-  , evalM   , subst-  , substA-  , substM-  , isRootOf-  , mapVar   , mapCoeff-  , associatedMonicPolynomial+  , toMonic   , toUPolynomialOf-  , polyDiv-  , polyMod-  , polyDivMod-  , polyGCD-  , polyLCM-  , prem-  , polyGCD'-  , polyMDivMod+  , divModMP   , reduce -  -- * Monomial-  , Monomial-  , monomialDegree-  , monomialProd-  , monomialDivisible-  , monomialDiv-  , monomialDeriv-  , monomialIntegral+  -- * Univariate polynomials+  , UPolynomial+  , X (..)+  , UTerm+  , UMonomial+  , div+  , mod+  , divMod+  , divides+  , gcd+  , lcm+  , exgcd+  , pdivMod+  , pdiv+  , pmod+  , gcd'+  , isSquareFree +  -- * Term+  , Term+  , tdeg+  , tmult+  , tdivides+  , tdiv+  , tderiv+  , tintegral+   -- * Monic monomial-  , MonicMonomial-  , mmOne-  , mmFromList-  , mmFromMap-  , mmFromIntMap-  , mmToList-  , mmToMap-  , mmToIntMap-  , mmProd-  , mmDivisible-  , mmDiv-  , mmDeriv-  , mmIntegral-  , mmLCM-  , mmGCD-  , mmMapVar+  , Monomial+  , mone+  , mfromIndices+  , mfromIndicesMap+  , mindices+  , mindicesMap+  , mmult+  , mpow+  , mdivides+  , mdiv+  , mderiv+  , mintegral+  , mlcm+  , mgcd+  , mcoprime    -- * Monomial order   , MonomialOrder@@ -110,637 +118,8 @@   , PrettyVar (..)   ) where -import Prelude hiding (lex)-import Control.Applicative-import Control.DeepSeq-import Control.Exception (assert)-import Control.Monad-import Data.Data-import qualified Data.FiniteField as FF-import Data.Function-import Data.List-import Data.Monoid-import Data.Ratio-import qualified Data.Map as Map-import qualified Data.Set as Set-import qualified Data.IntMap as IM-import Data.Traversable (for, traverse)-import Data.Typeable-import Data.VectorSpace-import qualified Text.PrettyPrint.HughesPJClass as PP-import Text.PrettyPrint.HughesPJClass (Doc, PrettyLevel, Pretty (..), prettyParen)--infixl 7  `polyDiv`, `polyMod`--{---------------------------------------------------------------------  Classes---------------------------------------------------------------------}--class Variables f where-  var       :: Ord v => v -> f v-  variables :: Ord v => f v -> Set.Set v---- | total degree of a given polynomial-class Degree t where-  deg :: t -> Integer--{---------------------------------------------------------------------  Polynomial type---------------------------------------------------------------------}---- | Polynomial over commutative ring r-newtype Polynomial k v = Polynomial{ coeffMap :: Map.Map (MonicMonomial v) k }-  deriving (Eq, Ord, Typeable)--instance (Eq k, Num k, Ord v) => Num (Polynomial k v) where-  (+)      = plus-  (*)      = prod-  negate   = neg-  abs x    = x -- OK?-  signum x = 1 -- OK?-  fromInteger x = constant (fromInteger x)--instance (Eq k, Num k, Ord v) => AdditiveGroup (Polynomial k v) where-  (^+^)   = plus-  zeroV   = zero-  negateV = neg--instance (Eq k, Num k, Ord v) => VectorSpace (Polynomial k v) where-  type Scalar (Polynomial k v) = k-  k *^ p = scale k p--instance (Show v, Ord v, Show k) => Show (Polynomial k v) where-  showsPrec d p  = showParen (d > 10) $-    showString "fromTerms " . shows (terms p)--instance (NFData k, NFData v) => NFData (Polynomial k v) where-  rnf (Polynomial m) = rnf m--instance (Eq k, Num k) => Variables (Polynomial k) where-  var x       = fromMonomial (1, var x)-  variables p = Set.unions $ [variables mm | (_, mm) <- terms p]--instance Degree (Polynomial k v) where-  deg p-    | isZero p  = -1-    | otherwise = maximum [deg mm | (_,mm) <- terms p]--normalize :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v-normalize (Polynomial m) = Polynomial (Map.filter (0/=) m)--asConstant :: Num k => Polynomial k v -> Maybe k-asConstant p =-  case terms p of-    [] -> Just 0-    [(c,xs)] | Map.null (mmToMap xs) -> Just c-    _ -> Nothing--scale :: (Eq k, Num k, Ord v) => k -> Polynomial k v -> Polynomial k v-scale 0 _ = zero-scale 1 p = p-scale a (Polynomial m) = normalize $ Polynomial (Map.map (a*) m)--zero :: (Eq k, Num k, Ord v) => Polynomial k v-zero = Polynomial $ Map.empty--plus :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v -> Polynomial k v-plus (Polynomial m1) (Polynomial m2) = normalize $ Polynomial $ Map.unionWith (+) m1 m2--neg :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v-neg (Polynomial m) = Polynomial $ Map.map negate m--prod :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v -> Polynomial k v-prod a b-  | Just c <- asConstant a = scale c b-  | Just c <- asConstant b = scale c a-prod (Polynomial m1) (Polynomial m2) = normalize $ Polynomial $ Map.fromListWith (+)-      [ (xs1 `mmProd` xs2, c1*c2)-      | (xs1,c1) <- Map.toList m1, (xs2,c2) <- Map.toList m2-      ]--isZero :: Polynomial k v -> Bool-isZero (Polynomial m) = Map.null m---- | construct a polynomial from a constant-constant :: (Eq k, Num k, Ord v) => k -> Polynomial k v-constant c = fromMonomial (c, mmOne)---- | construct a polynomial from a list of monomials-fromTerms :: (Eq k, Num k, Ord v) => [Monomial k v] -> Polynomial k v-fromTerms = normalize . Polynomial . Map.fromListWith (+) . map (\(c,xs) -> (xs,c))--fromCoeffMap :: (Eq k, Num k, Ord v) => Map.Map (MonicMonomial v) k -> Polynomial k v-fromCoeffMap m = normalize $ Polynomial m---- | construct a polynomial from a monomial-fromMonomial :: (Eq k, Num k, Ord v) => Monomial k v -> Polynomial k v-fromMonomial (c,xs) = normalize $ Polynomial $ Map.singleton xs c---- | list of monomials-terms :: Polynomial k v -> [Monomial k v]-terms (Polynomial m) = [(c,xs) | (xs,c) <- Map.toList m]---- | leading term with respect to a given monomial order-leadingTerm :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Monomial k v-leadingTerm cmp p =-  case terms p of-    [] -> (0, mmOne) -- should be error?-    ms -> maximumBy (cmp `on` snd) ms--coeff :: (Num k, Ord v) => MonicMonomial v -> Polynomial k v -> k-coeff xs (Polynomial m) = Map.findWithDefault 0 xs m--lookupCoeff :: Ord v => MonicMonomial v -> Polynomial k v -> Maybe k-lookupCoeff xs (Polynomial m) = Map.lookup xs m--contI :: (Integral r, Ord v) => Polynomial r v -> r-contI 0 = 1-contI p = foldl1' gcd [abs c | (c,_) <- terms p]--ppI :: (Integral r, Ord v) => Polynomial r v -> Polynomial r v-ppI p = mapCoeff f p-  where-    c = contI p-    f x = assert (x `mod` c == 0) $ x `div` c--class ContPP k where-  -- | Content of a polynomial  -  cont :: (Ord v) => Polynomial k v -> k-  -- constructive-algebra-0.3.0 では cont 0 は error になる--  -- | Primitive part of a polynomial-  pp :: (Ord v) => Polynomial k v -> Polynomial k v--instance ContPP Integer where-  cont = contI-  pp   = ppI--instance Integral r => ContPP (Ratio r) where-  {-# SPECIALIZE instance ContPP (Ratio Integer) #-}--  cont 0 = 1-  cont p = foldl1' gcd ns % foldl' lcm 1 ds-    where-      ns = [abs (numerator c) | (c,_) <- terms p]-      ds = [denominator c     | (c,_) <- terms p]  --  pp p = mapCoeff (/ c) p-    where-      c = cont p--isPrimitive :: (Eq k, Num k, ContPP k, Ord v) => Polynomial k v -> Bool-isPrimitive p = isZero p || cont p == 1---- | Formal derivative of polynomials-deriv :: (Eq k, Num k, Ord v) => Polynomial k v -> v -> Polynomial k v-deriv p x = sumV [fromMonomial (monomialDeriv m x) | m <- terms p]---- | Formal integral of polynomials-integral :: (Eq k, Fractional k, Ord v) => Polynomial k v -> v -> Polynomial k v-integral p x = sumV [fromMonomial (monomialIntegral m x) | m <- terms p]---- | Evaluation-eval :: (Num k, Ord v) => (v -> k) -> Polynomial k v -> k-eval env p = sum [c * product [(env x) ^ e | (x,e) <- mmToList xs] | (c,xs) <- terms p]---- | Evaluation-evalA :: forall k v f. (Num k, Ord v, Applicative f) => (v -> f k) -> Polynomial k v -> f k-evalA env p = sum <$> traverse f (terms p)-  where-    f :: Monomial k v -> f k-    f (c,xs) = ((c*) . product) <$> g xs-    g :: MonicMonomial v -> f [k]-    g xs = traverse (\(x,e) -> liftA (^ e) (env x)) (mmToList xs)---- | Evaluation-evalM :: (Num k, Ord v, Monad m) => (v -> m k) -> Polynomial k v -> m k-evalM env p = do-  liftM sum $ forM (terms p) $ \(c,xs) -> do-    rs <- mapM (\(x,e) -> liftM (^ e) (env x)) (mmToList xs)-    return (c * product rs)---- | Substitution or bind-subst-  :: (Eq k, Num k, Ord v1, Ord v2)-  => Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2-subst p s =-  sumV [constant c * product [(s x)^e | (x,e) <- mmToList xs] | (c, xs) <- terms p]---- | Substitution or bind-substA-  :: forall k v1 v2 f. (Eq k, Num k, Ord v1, Ord v2, Applicative f)-  => Polynomial k v1 -> (v1 -> f (Polynomial k v2)) -> f (Polynomial k v2)-substA p s = sumV <$> traverse f (terms p)-  where-    f :: Monomial k v1 -> f (Polynomial k v2)-    f (c,xs) =  ((constant c *) . product) <$> g xs-    g :: MonicMonomial v1 -> f [Polynomial k v2]-    g xs = traverse (\(x,e) -> liftA (^ e) (s x)) (mmToList xs)---- | Substitution or bind-substM-  :: (Eq k, Num k, Ord v1, Ord v2, Monad m)-  => Polynomial k v1 -> (v1 -> m (Polynomial k v2)) -> m (Polynomial k v2)-substM p s = liftM sum $ forM (terms p) $ \(c,xs) -> do-  xs <- forM (mmToList xs) $ \(x,e) -> liftM (^e) (s x)-  return $ constant c * product xs--isRootOf :: (Eq k, Num k) => k -> UPolynomial k -> Bool-isRootOf x p = eval (\_ -> x) p == 0--mapVar :: (Eq k, Num k, Ord v1, Ord v2) => (v1 -> v2) -> Polynomial k v1 -> Polynomial k v2-mapVar f (Polynomial m) = normalize $ Polynomial $ Map.mapKeysWith (+) (mmMapVar f) m--mapCoeff :: (Eq k1, Num k1, Ord v) => (k -> k1) -> Polynomial k v -> Polynomial k1 v-mapCoeff f (Polynomial m) = Polynomial $ Map.mapMaybe g m-  where-    g x = if y == 0 then Nothing else Just y-      where-        y = f x--associatedMonicPolynomial :: (Eq r, Fractional r, Ord v) => MonomialOrder v -> Polynomial r v -> Polynomial r v-associatedMonicPolynomial cmp p-  | c == 0 = p-  | otherwise = mapCoeff (/c) p-  where-    (c,_) = leadingTerm cmp p--toUPolynomialOf :: (Ord k, Num k, Ord v) => Polynomial k v -> v -> UPolynomial (Polynomial k v)-toUPolynomialOf p v = fromTerms $ do-  (c,mm) <- terms p-  let m = mmToMap mm-  return ( fromTerms [(c, mmFromMap (Map.delete v m))]-         , mmFromList [(X, Map.findWithDefault 0 v m)]-         )---- | Multivariate division algorithm-polyMDivMod-  :: forall k v. (Eq k, Fractional k, Ord v)-  => MonomialOrder v -> Polynomial k v -> [Polynomial k v] -> ([Polynomial k v], Polynomial k v)-polyMDivMod cmp p fs = go IM.empty p-  where-    ls = [(leadingTerm cmp f, f) | f <- fs]--    go :: IM.IntMap (Polynomial k v) -> Polynomial k v -> ([Polynomial k v], Polynomial k v)-    go qs g =-      case xs of-        [] -> ([IM.findWithDefault 0 i qs | i <- [0 .. length fs - 1]], g)-        (i, b, g') : _ -> go (IM.insertWith (+) i b qs) g'-      where-        ms = sortBy (flip cmp `on` snd) (terms g)-        xs = do-          (i,(a,f)) <- zip [0..] ls-          h <- ms-          guard $ monomialDivisible h a-          let b = fromMonomial $ monomialDiv h a-          return (i, b, g - b * f)---- | Multivariate division algorithm-reduce-  :: (Eq k, Fractional k, Ord v)-  => MonomialOrder v -> Polynomial k v -> [Polynomial k v] -> Polynomial k v-reduce cmp p fs = go p-  where-    ls = [(leadingTerm cmp f, f) | f <- fs]-    go g = if null xs then g else go (head xs)-      where-        ms = sortBy (flip cmp `on` snd) (terms g)-        xs = do-          (a,f) <- ls-          h <- ms-          guard $ monomialDivisible h a-          return (g - fromMonomial (monomialDiv h a) * f)--{---------------------------------------------------------------------  Pretty printing---------------------------------------------------------------------}--data PrintOptions k v-  = PrintOptions-  { pOptPrintVar        :: PrettyLevel -> Rational -> v -> Doc-  , pOptPrintCoeff      :: PrettyLevel -> Rational -> k -> Doc-  , pOptIsNegativeCoeff :: k -> Bool-  , pOptMonomialOrder   :: MonomialOrder v-  }--defaultPrintOptions :: (PrettyCoeff k, PrettyVar v, Ord v) => PrintOptions k v-defaultPrintOptions-  = PrintOptions-  { pOptPrintVar        = pPrintVar-  , pOptPrintCoeff      = pPrintCoeff-  , pOptIsNegativeCoeff = isNegativeCoeff-  , pOptMonomialOrder   = grlex-  }--instance (Ord k, Num k, Ord v, PrettyCoeff k, PrettyVar v) => Pretty (Polynomial k v) where-  pPrintPrec = prettyPrint defaultPrintOptions--prettyPrint-  :: (Ord k, Num k, Ord v)-  => PrintOptions k v-  -> PrettyLevel -> Rational -> Polynomial k v -> Doc-prettyPrint opt lv prec p =-    case sortBy (flip (pOptMonomialOrder opt) `on` snd) $ terms p of-      [] -> PP.int 0-      [t] -> pLeadingTerm prec t-      t:ts ->-        prettyParen (prec > addPrec) $-          PP.hcat (pLeadingTerm addPrec t : map pTrailingTerm ts)-    where-      pLeadingTerm prec (c,xs) =-        if pOptIsNegativeCoeff opt c-        then prettyParen (prec > addPrec) $-               PP.char '-' <> prettyPrintMonomial opt lv (addPrec+1) (-c,xs)-        else prettyPrintMonomial opt lv prec (c,xs)--      pTrailingTerm (c,xs) =-        if pOptIsNegativeCoeff opt c-        then PP.space <> PP.char '-' <> PP.space <> prettyPrintMonomial opt lv (addPrec+1) (-c,xs)-        else PP.space <> PP.char '+' <> PP.space <> prettyPrintMonomial opt lv (addPrec+1) (c,xs)--prettyPrintMonomial-  :: (Ord k, Num k, Ord v)-  => PrintOptions k v-  -> PrettyLevel -> Rational -> Monomial k v -> Doc-prettyPrintMonomial opt lv prec (c,xs)-  | len == 0  = pOptPrintCoeff opt lv (appPrec+1) c-    -- intentionally specify (appPrec+1) to parenthesize any composite expression-  | len == 1 && c == 1 = pPow prec $ head (mmToList xs)-  | otherwise =-      prettyParen (prec > mulPrec) $-        PP.hcat $ intersperse (PP.char '*') fs-    where-      len = length $ mmToList xs-      fs  = [pOptPrintCoeff opt lv (appPrec+1) c | c /= 1] ++ [pPow (mulPrec+1) p | p <- mmToList xs]-      -- intentionally specify (appPrec+1) to parenthesize any composite expression--      pPow prec (x,1) = pOptPrintVar opt lv prec x-      pPow prec (x,n) =-        prettyParen (prec > expPrec) $-          pOptPrintVar opt lv (expPrec+1) x <> PP.char '^' <> PP.integer n--class PrettyCoeff a where-  pPrintCoeff :: PrettyLevel -> Rational -> a -> Doc-  isNegativeCoeff :: a -> Bool-  isNegativeCoeff _ = False--instance PrettyCoeff Integer where-  pPrintCoeff = pPrintPrec-  isNegativeCoeff = (0>)--instance (PrettyCoeff a, Integral a) => PrettyCoeff (Ratio a) where-  pPrintCoeff lv p r-    | denominator r == 1 = pPrintCoeff lv p (numerator r)-    | otherwise = -        prettyParen (p > ratPrec) $-          pPrintCoeff lv (ratPrec+1) (numerator r) <>-          PP.char '/' <>-          pPrintCoeff lv (ratPrec+1) (denominator r)-  isNegativeCoeff x = isNegativeCoeff (numerator x)--instance PrettyCoeff (FF.PrimeField a) where-  pPrintCoeff lv p a = pPrintCoeff lv p (FF.toInteger a)-  isNegativeCoeff _  = False--instance (Num c, Ord c, PrettyCoeff c, Ord v, PrettyVar v) => PrettyCoeff (Polynomial c v) where-  pPrintCoeff = pPrintPrec--class PrettyVar a where-  pPrintVar :: PrettyLevel -> Rational -> a -> Doc--instance PrettyVar Int where-  pPrintVar _ _ n = PP.char 'x' <> PP.int n--instance PrettyVar X where-  pPrintVar _ _ X = PP.char 'x'--addPrec, mulPrec, ratPrec, expPrec :: Rational-addPrec = 6 -- Precedence of '+'-mulPrec = 7 -- Precedence of '*'-ratPrec = 7 -- Precedence of '/'-expPrec = 8 -- Precedence of '^'-appPrec = 10 -- Precedence of function application--{---------------------------------------------------------------------  Univariate polynomials---------------------------------------------------------------------}---- | Univariate polynomials over commutative ring r-type UPolynomial r = Polynomial r X--data X = X-  deriving (Eq, Ord, Bounded, Enum, Show, Read, Typeable, Data)--instance NFData X---- | division of univariate polynomials-polyDiv :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k-polyDiv f1 f2 = fst (polyDivMod f1 f2)---- | division of univariate polynomials-polyMod :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k-polyMod f1 f2 = snd (polyDivMod f1 f2)---- | division of univariate polynomials-polyDivMod :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)-polyDivMod f g-  | isZero g  = error "polyDivMod: division by zero"-  | otherwise = go 0 f-  where-    lt_g = leadingTerm lex g-    go !q !r-      | deg r < deg g = (q,r)-      | otherwise     = go (q + t) (r - t * g)-        where-          lt_r = leadingTerm lex r-          t    = fromMonomial $ lt_r `monomialDiv` lt_g---- | GCD of univariate polynomials-polyGCD :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k-polyGCD f1 0  = associatedMonicPolynomial grlex f1-polyGCD f1 f2 = polyGCD f2 (f1 `polyMod` f2)---- | LCM of univariate polynomials-polyLCM :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k-polyLCM _ 0 = 0-polyLCM 0 _ = 0-polyLCM f1 f2 = associatedMonicPolynomial grlex $ (f1 `polyMod` (polyGCD f1 f2)) * f2---- | pseudo reminder-prem :: (Eq r, Integral r) => UPolynomial r -> UPolynomial r -> UPolynomial r-prem _ 0 = error "prem: division by 0"-prem f g-  | deg f < deg g = f-  | otherwise     = go (scale (lc_g ^ (deg f - deg g + 1)) f)-  where-    (lc_g, lm_g) = leadingTerm lex g-    deg_g    = deg g-    go !f1-      | deg_g > deg f1 = f1-      | otherwise =-          assert (lc_f1 `mod` lc_g == 0 && mmDivisible lm_f1 lm_g) $-            go (f1 - fromMonomial (lc_f1 `div` lc_g, lm_f1 `mmDiv` lm_g) * g)-          where-            (lc_f1, lm_f1) = leadingTerm lex f1---- | GCD of univariate polynomials-polyGCD' :: (Eq r, Integral r) => UPolynomial r -> UPolynomial r -> UPolynomial r-polyGCD' f1 0  = ppI f1-polyGCD' f1 f2 = polyGCD' f2 (f1 `prem` f2)--{---------------------------------------------------------------------  Monomial---------------------------------------------------------------------}--type Monomial k v = (k, MonicMonomial v)--monomialDegree :: Monomial k v -> Integer-monomialDegree (_,xs) = deg xs--monomialProd :: (Num k, Ord v) => Monomial k v -> Monomial k v -> Monomial k v-monomialProd (c1,xs1) (c2,xs2) = (c1*c2, xs1 `mmProd` xs2)--monomialDivisible :: (Fractional k, Ord v) => Monomial k v -> Monomial k v -> Bool-monomialDivisible (c1,xs1) (c2,xs2) = mmDivisible xs1 xs2--monomialDiv :: (Fractional k, Ord v) => Monomial k v -> Monomial k v -> Monomial k v-monomialDiv (c1,xs1) (c2,xs2) = (c1 / c2, xs1 `mmDiv` xs2)--monomialDeriv :: (Eq k, Num k, Ord v) => Monomial k v -> v -> Monomial k v-monomialDeriv (c,xs) x =-  case mmDeriv xs x of-    (s,ys) -> (c * fromIntegral s, ys)--monomialIntegral :: (Eq k, Fractional k, Ord v) => Monomial k v -> v -> Monomial k v-monomialIntegral (c,xs) x =-  case mmIntegral xs x of-    (s,ys) -> (c * fromRational s, ys)--{---------------------------------------------------------------------  Monic Monomial---------------------------------------------------------------------}---- 本当は変数の型に応じて type family で表現を変えたい---- | Monic monomials-newtype MonicMonomial v = MonicMonomial{ mmToMap :: Map.Map v Integer }-  deriving (Eq, Ord, Typeable)--instance (Ord v, Show v) => Show (MonicMonomial v) where-  showsPrec d m  = showParen (d > 10) $-    showString "mmFromList " . shows (mmToList m)--instance (NFData v) => NFData (MonicMonomial v) where-  rnf (MonicMonomial m) = rnf m--instance Degree (MonicMonomial v) where-  deg (MonicMonomial m) = sum $ Map.elems m--instance Variables MonicMonomial where-  var x        = MonicMonomial $ Map.singleton x 1-  variables mm = Map.keysSet (mmToMap mm)--mmNormalize :: Ord v => MonicMonomial v -> MonicMonomial v-mmNormalize (MonicMonomial m) = MonicMonomial $ Map.filter (>0) m--mmOne :: MonicMonomial v-mmOne = MonicMonomial $ Map.empty--mmFromList :: Ord v => [(v, Integer)] -> MonicMonomial v-mmFromList xs-  | any (\(x,e) -> 0>e) xs = error "mmFromList: negative exponent"-  | otherwise = MonicMonomial $ Map.fromListWith (+) [(x,e) | (x,e) <- xs, e > 0]--mmFromMap :: Ord v => Map.Map v Integer -> MonicMonomial v-mmFromMap m-  | any (\(x,e) -> 0>e) (Map.toList m) = error "mmFromFromMap: negative exponent"-  | otherwise = mmNormalize $ MonicMonomial m--mmFromIntMap :: IM.IntMap Integer -> MonicMonomial Int-mmFromIntMap = mmFromMap . Map.fromDistinctAscList . IM.toAscList--mmToList :: Ord v => MonicMonomial v -> [(v, Integer)]-mmToList (MonicMonomial m) = Map.toAscList m--mmToIntMap :: MonicMonomial Int -> IM.IntMap Integer-mmToIntMap (MonicMonomial m) = IM.fromDistinctAscList $ Map.toAscList m--mmProd :: Ord v => MonicMonomial v -> MonicMonomial v -> MonicMonomial v-mmProd (MonicMonomial xs1) (MonicMonomial xs2) = mmNormalize $ MonicMonomial $ Map.unionWith (+) xs1 xs2--mmDivisible :: Ord v => MonicMonomial v -> MonicMonomial v -> Bool-mmDivisible (MonicMonomial xs1) (MonicMonomial xs2) = Map.isSubmapOfBy (<=) xs2 xs1--mmDiv :: Ord v => MonicMonomial v -> MonicMonomial v -> MonicMonomial v-mmDiv (MonicMonomial xs1) (MonicMonomial xs2) = MonicMonomial $ Map.differenceWith f xs1 xs2-  where-    f m n-      | m <= n    = Nothing-      | otherwise = Just (m - n)--mmDeriv :: Ord v => MonicMonomial v -> v -> (Integer, MonicMonomial v)-mmDeriv (MonicMonomial xs) x-  | n==0      = (0, mmOne)-  | otherwise = (n, MonicMonomial $ Map.update f x xs)-  where-    n = Map.findWithDefault 0 x xs-    f m-      | m <= 1    = Nothing-      | otherwise = Just $! m - 1--mmIntegral :: Ord v => MonicMonomial v -> v -> (Rational, MonicMonomial v)-mmIntegral (MonicMonomial xs) x =-  (1 % fromIntegral (n + 1), MonicMonomial $ Map.insert x (n+1) xs)-  where-    n = Map.findWithDefault 0 x xs--mmLCM :: Ord v => MonicMonomial v -> MonicMonomial v -> MonicMonomial v-mmLCM (MonicMonomial m1) (MonicMonomial m2) = MonicMonomial $ Map.unionWith max m1 m2--mmGCD :: Ord v => MonicMonomial v -> MonicMonomial v -> MonicMonomial v-mmGCD (MonicMonomial m1) (MonicMonomial m2) = MonicMonomial $ Map.intersectionWith min m1 m2--mmMapVar :: (Ord v1, Ord v2) => (v1 -> v2) -> MonicMonomial v1 -> MonicMonomial v2-mmMapVar f (MonicMonomial m) = MonicMonomial $ Map.mapKeysWith (+) f m--{---------------------------------------------------------------------  Monomial Order---------------------------------------------------------------------}--type MonomialOrder v = MonicMonomial v -> MonicMonomial v -> Ordering---- | Lexicographic order-lex :: Ord v => MonomialOrder v-lex xs1 xs2 = go (mmToList xs1) (mmToList xs2)-  where-    go [] [] = EQ-    go [] _  = LT -- = cmpare 0 n2-    go _ []  = GT -- = cmpare n1 0-    go ((x1,n1):xs1) ((x2,n2):xs2) =-      case compare x1 x2 of-        LT -> GT -- = compare n1 0-        GT -> LT -- = compare 0 n2-        EQ -> compare n1 n2 `mappend` go xs1 xs2---- | Reverse lexicographic order--- Note that revlex is NOT a monomial order.-revlex :: Ord v => MonicMonomial v -> MonicMonomial v -> Ordering-revlex xs1 xs2 = go (reverse (mmToList xs1)) (reverse (mmToList xs2))-  where-    go [] [] = EQ-    go [] _  = GT -- = cmp 0 n2-    go _ []  = LT -- = cmp n1 0-    go ((x1,n1):xs1) ((x2,n2):xs2) =-      case compare x1 x2 of-        LT -> GT -- = cmp 0 n2-        GT -> LT -- = cmp n1 0-        EQ -> cmp n1 n2 `mappend` go xs1 xs2-    cmp n1 n2 = compare n2 n1---- | graded lexicographic order-grlex :: Ord v => MonomialOrder v-grlex = (compare `on` deg) `mappend` lex---- | graded reverse lexicographic order-grevlex :: Ord v => MonomialOrder v-grevlex = (compare `on` deg) `mappend` revlex+import Prelude hiding (lex, div, mod, divMod, gcd, lcm)+import Data.Polynomial.Base+import Data.Polynomial.Factorization.FiniteField ()+import Data.Polynomial.Factorization.Integer ()+import Data.Polynomial.Factorization.Rational ()

@@ -0,0 +1,802 @@+{-# LANGUAGE ScopedTypeVariables, FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies, TypeFamilies, BangPatterns, DeriveDataTypeable #-}+{-# OPTIONS_GHC -Wall #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Polynomial.Base+-- Copyright   :  (c) Masahiro Sakai 2012-2013+-- License     :  BSD-style+-- +-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable (ScopedTypeVariables, FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies, TypeFamilies, BangPatterns, DeriveDataTypeable)+--+-- Polynomials+--+-- References:+--+-- * Monomial order <http://en.wikipedia.org/wiki/Monomial_order>+--+-- * Polynomial class for Ruby <http://www.math.kobe-u.ac.jp/~kodama/tips-RubyPoly.html>+--+-- * constructive-algebra package <http://hackage.haskell.org/package/constructive-algebra>+-- +-----------------------------------------------------------------------------+module Data.Polynomial.Base+  (+  -- * Polynomial type+    Polynomial++  -- * Conversion+  , Var (..)+  , constant+  , terms+  , fromTerms+  , coeffMap+  , fromCoeffMap+  , fromTerm++  -- * Query+  , Degree (..)+  , Vars (..)+  , lt+  , lc+  , lm+  , coeff+  , lookupCoeff+  , isPrimitive+  , isRootOf++  -- * Operations+  , Factor (..)+  , SQFree (..)+  , ContPP (..)+  , deriv+  , integral+  , eval+  , subst+  , mapCoeff+  , toMonic+  , toUPolynomialOf+  , divModMP+  , reduce++  -- * Univariate polynomials+  , UPolynomial+  , X (..)+  , UTerm+  , UMonomial+  , div+  , mod+  , divMod+  , divides+  , gcd+  , lcm+  , exgcd+  , pdivMod+  , pdiv+  , pmod+  , gcd'+  , isSquareFree++  -- * Term+  , Term+  , tdeg+  , tmult+  , tdivides+  , tdiv+  , tderiv+  , tintegral++  -- * Monic monomial+  , Monomial+  , mone+  , mfromIndices+  , mfromIndicesMap+  , mindices+  , mindicesMap+  , mmult+  , mpow+  , mdivides+  , mdiv+  , mderiv+  , mintegral+  , mlcm+  , mgcd+  , mcoprime++  -- * Monomial order+  , MonomialOrder+  , lex+  , revlex+  , grlex+  , grevlex++  -- * Pretty Printing+  , PrintOptions (..)+  , defaultPrintOptions+  , prettyPrint+  , PrettyCoeff (..)+  , PrettyVar (..)+  ) where++import Prelude hiding (lex, div, mod, divMod, gcd, lcm)+import qualified Prelude+import Control.DeepSeq+import Control.Exception (assert)+import Control.Monad+import Data.Data+import qualified Data.FiniteField as FF+import Data.Function+import Data.List+import Data.Monoid+import Data.Ratio+import Data.Map (Map)+import qualified Data.Map as Map+import Data.Set (Set)+import qualified Data.Set as Set+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap+import Data.Typeable+import Data.VectorSpace+import qualified Text.PrettyPrint.HughesPJClass as PP+import Text.PrettyPrint.HughesPJClass (Doc, PrettyLevel, Pretty (..), prettyParen)++infixl 7  `div`, `mod`++{--------------------------------------------------------------------+  Classes+--------------------------------------------------------------------}++class Vars a v => Var a v | a -> v where+  var :: v -> a++class Ord v => Vars a v | a -> v where+  vars :: a -> Set v++-- | total degree of a given polynomial+class Degree t where+  deg :: t -> Integer++{--------------------------------------------------------------------+  Polynomial type+--------------------------------------------------------------------}++-- | Polynomial over commutative ring r+newtype Polynomial r v = Polynomial{ coeffMap :: Map (Monomial v) r }+  deriving (Eq, Ord, Typeable)++instance (Eq k, Num k, Ord v) => Num (Polynomial k v) where+  (+)      = plus+  (*)      = mult+  negate   = neg+  abs x    = x -- OK?+  signum _ = 1 -- OK?+  fromInteger x = constant (fromInteger x)++instance (Eq k, Num k, Ord v) => AdditiveGroup (Polynomial k v) where+  (^+^)   = plus+  zeroV   = zero+  negateV = neg++instance (Eq k, Num k, Ord v) => VectorSpace (Polynomial k v) where+  type Scalar (Polynomial k v) = k+  k *^ p = scale k p++instance (Show v, Ord v, Show k) => Show (Polynomial k v) where+  showsPrec d p  = showParen (d > 10) $+    showString "fromTerms " . shows (terms p)++instance (NFData k, NFData v) => NFData (Polynomial k v) where+  rnf (Polynomial m) = rnf m++instance (Eq k, Num k, Ord v) => Var (Polynomial k v) v where+  var x = fromTerm (1, var x)++instance (Eq k, Num k, Ord v) => Vars (Polynomial k v) v where+  vars p = Set.unions $ [vars mm | (_, mm) <- terms p]++instance Degree (Polynomial k v) where+  deg p+    | isZero p  = -1+    | otherwise = maximum [deg mm | (_,mm) <- terms p]++normalize :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v+normalize (Polynomial m) = Polynomial (Map.filter (0/=) m)++asConstant :: Num k => Polynomial k v -> Maybe k+asConstant p =+  case terms p of+    [] -> Just 0+    [(c,xs)] | Map.null (mindicesMap xs) -> Just c+    _ -> Nothing++scale :: (Eq k, Num k, Ord v) => k -> Polynomial k v -> Polynomial k v+scale 0 _ = zero+scale 1 p = p+scale a (Polynomial m) = normalize $ Polynomial (Map.map (a*) m)++zero :: (Eq k, Num k, Ord v) => Polynomial k v+zero = Polynomial $ Map.empty++plus :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v -> Polynomial k v+plus (Polynomial m1) (Polynomial m2) = normalize $ Polynomial $ Map.unionWith (+) m1 m2++neg :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v+neg (Polynomial m) = Polynomial $ Map.map negate m++mult :: (Eq k, Num k, Ord v) => Polynomial k v -> Polynomial k v -> Polynomial k v+mult a b+  | Just c <- asConstant a = scale c b+  | Just c <- asConstant b = scale c a+mult (Polynomial m1) (Polynomial m2) = normalize $ Polynomial $ Map.fromListWith (+)+      [ (xs1 `mmult` xs2, c1*c2)+      | (xs1,c1) <- Map.toList m1, (xs2,c2) <- Map.toList m2+      ]++isZero :: Polynomial k v -> Bool+isZero (Polynomial m) = Map.null m++-- | construct a polynomial from a constant+constant :: (Eq k, Num k, Ord v) => k -> Polynomial k v+constant c = fromTerm (c, mone)++-- | construct a polynomial from a list of monomials+fromTerms :: (Eq k, Num k, Ord v) => [Term k v] -> Polynomial k v+fromTerms = normalize . Polynomial . Map.fromListWith (+) . map (\(c,xs) -> (xs,c))++fromCoeffMap :: (Eq k, Num k, Ord v) => Map (Monomial v) k -> Polynomial k v+fromCoeffMap m = normalize $ Polynomial m++-- | construct a polynomial from a monomial+fromTerm :: (Eq k, Num k, Ord v) => Term k v -> Polynomial k v+fromTerm (c,xs) = normalize $ Polynomial $ Map.singleton xs c++-- | list of monomials+terms :: Polynomial k v -> [Term k v]+terms (Polynomial m) = [(c,xs) | (xs,c) <- Map.toList m]++-- | leading term with respect to a given monomial order+lt :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Term k v+lt cmp p =+  case terms p of+    [] -> (0, mone) -- should be error?+    ms -> maximumBy (cmp `on` snd) ms++-- | leading coefficient with respect to a given monomial order+lc :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> k+lc cmp = fst . lt cmp++-- | leading monomial with respect to a given monomial order+lm :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Monomial v+lm cmp = snd . lt cmp++coeff :: (Num k, Ord v) => Monomial v -> Polynomial k v -> k+coeff xs (Polynomial m) = Map.findWithDefault 0 xs m++lookupCoeff :: Ord v => Monomial v -> Polynomial k v -> Maybe k+lookupCoeff xs (Polynomial m) = Map.lookup xs m++contI :: (Integral r, Ord v) => Polynomial r v -> r+contI 0 = 1+contI p = foldl1' Prelude.gcd [abs c | (c,_) <- terms p]++ppI :: (Integral r, Ord v) => Polynomial r v -> Polynomial r v+ppI p = mapCoeff f p+  where+    c = contI p+    f x = assert (x `Prelude.mod` c == 0) $ x `Prelude.div` c++class ContPP k where+  -- | Content of a polynomial  +  cont :: (Ord v) => Polynomial k v -> k+  -- constructive-algebra-0.3.0 では cont 0 は error になる++  -- | Primitive part of a polynomial+  pp :: (Ord v) => Polynomial k v -> Polynomial k v++instance ContPP Integer where+  cont = contI+  pp   = ppI++instance Integral r => ContPP (Ratio r) where+  {-# SPECIALIZE instance ContPP (Ratio Integer) #-}++  cont 0 = 1+  cont p = foldl1' Prelude.gcd ns % foldl' Prelude.lcm 1 ds+    where+      ns = [abs (numerator c) | (c,_) <- terms p]+      ds = [denominator c     | (c,_) <- terms p]  ++  pp p = mapCoeff (/ c) p+    where+      c = cont p++isPrimitive :: (Eq k, Num k, ContPP k, Ord v) => Polynomial k v -> Bool+isPrimitive p = isZero p || cont p == 1++-- | Formal derivative of polynomials+deriv :: (Eq k, Num k, Ord v) => Polynomial k v -> v -> Polynomial k v+deriv p x = sumV [fromTerm (tderiv m x) | m <- terms p]++-- | Formal integral of polynomials+integral :: (Eq k, Fractional k, Ord v) => Polynomial k v -> v -> Polynomial k v+integral p x = sumV [fromTerm (tintegral m x) | m <- terms p]++-- | Evaluation+eval :: (Num k, Ord v) => (v -> k) -> Polynomial k v -> k+eval env p = sum [c * product [(env x) ^ e | (x,e) <- mindices xs] | (c,xs) <- terms p]++-- | Substitution or bind+subst+  :: (Eq k, Num k, Ord v1, Ord v2)+  => Polynomial k v1 -> (v1 -> Polynomial k v2) -> Polynomial k v2+subst p s =+  sumV [constant c * product [(s x)^e | (x,e) <- mindices xs] | (c, xs) <- terms p]++isRootOf :: (Eq k, Num k) => k -> UPolynomial k -> Bool+isRootOf x p = eval (\_ -> x) p == 0++isSquareFree :: (Eq k, Fractional k) => UPolynomial k -> Bool+isSquareFree p = gcd p (deriv p X) == 1++mapCoeff :: (Eq k1, Num k1, Ord v) => (k -> k1) -> Polynomial k v -> Polynomial k1 v+mapCoeff f (Polynomial m) = Polynomial $ Map.mapMaybe g m+  where+    g x = if y == 0 then Nothing else Just y+      where+        y = f x++toMonic :: (Eq r, Fractional r, Ord v) => MonomialOrder v -> Polynomial r v -> Polynomial r v+toMonic cmp p+  | c == 0 || c == 1 = p+  | otherwise = mapCoeff (/c) p+  where+    c = lc cmp p++toUPolynomialOf :: (Ord k, Num k, Ord v) => Polynomial k v -> v -> UPolynomial (Polynomial k v)+toUPolynomialOf p v = fromTerms $ do+  (c,mm) <- terms p+  let m = mindicesMap mm+  return ( fromTerms [(c, mfromIndicesMap (Map.delete v m))]+         , var X `mpow` Map.findWithDefault 0 v m+         )++-- | Multivariate division algorithm+divModMP+  :: forall k v. (Eq k, Fractional k, Ord v)+  => MonomialOrder v -> Polynomial k v -> [Polynomial k v] -> ([Polynomial k v], Polynomial k v)+divModMP cmp p fs = go IntMap.empty p+  where+    ls = [(lt cmp f, f) | f <- fs]++    go :: IntMap (Polynomial k v) -> Polynomial k v -> ([Polynomial k v], Polynomial k v)+    go qs g =+      case xs of+        [] -> ([IntMap.findWithDefault 0 i qs | i <- [0 .. length fs - 1]], g)+        (i, b, g') : _ -> go (IntMap.insertWith (+) i b qs) g'+      where+        ms = sortBy (flip cmp `on` snd) (terms g)+        xs = do+          (i,(a,f)) <- zip [0..] ls+          h <- ms+          guard $ a `tdivides` h+          let b = fromTerm $ tdiv h a+          return (i, b, g - b * f)++-- | Multivariate division algorithm+reduce+  :: (Eq k, Fractional k, Ord v)+  => MonomialOrder v -> Polynomial k v -> [Polynomial k v] -> Polynomial k v+reduce cmp p fs = go p+  where+    ls = [(lt cmp f, f) | f <- fs]+    go g = if null xs then g else go (head xs)+      where+        ms = sortBy (flip cmp `on` snd) (terms g)+        xs = do+          (a,f) <- ls+          h <- ms+          guard $ a `tdivides` h+          return (g - fromTerm (tdiv h a) * f)++-- | Factorization of polynomials+class Factor a where+  -- | factor a polynomial @p@ into @p1 ^ n1 + p2 ^ n2 + ..@ and+  -- return a list @[(p1,n1), (p2,n2), ..]@.+  factor :: a -> [(a, Integer)]++-- | Square-free factorization of polynomials+class SQFree a where+  -- | factor a polynomial @p@ into @p1 ^ n1 + p2 ^ n2 + ..@ and+  -- return a list @[(p1,n1), (p2,n2), ..]@.+  sqfree :: a -> [(a, Integer)]++{--------------------------------------------------------------------+  Pretty printing+--------------------------------------------------------------------}++data PrintOptions k v+  = PrintOptions+  { pOptPrintVar        :: PrettyLevel -> Rational -> v -> Doc+  , pOptPrintCoeff      :: PrettyLevel -> Rational -> k -> Doc+  , pOptIsNegativeCoeff :: k -> Bool+  , pOptMonomialOrder   :: MonomialOrder v+  }++defaultPrintOptions :: (PrettyCoeff k, PrettyVar v, Ord v) => PrintOptions k v+defaultPrintOptions+  = PrintOptions+  { pOptPrintVar        = pPrintVar+  , pOptPrintCoeff      = pPrintCoeff+  , pOptIsNegativeCoeff = isNegativeCoeff+  , pOptMonomialOrder   = grlex+  }++instance (Ord k, Num k, Ord v, PrettyCoeff k, PrettyVar v) => Pretty (Polynomial k v) where+  pPrintPrec = prettyPrint defaultPrintOptions++prettyPrint+  :: (Ord k, Num k, Ord v)+  => PrintOptions k v+  -> PrettyLevel -> Rational -> Polynomial k v -> Doc+prettyPrint opt lv prec p =+    case sortBy (flip (pOptMonomialOrder opt) `on` snd) $ terms p of+      [] -> PP.int 0+      [t] -> pLeadingTerm prec t+      t:ts ->+        prettyParen (prec > addPrec) $+          PP.hcat (pLeadingTerm addPrec t : map pTrailingTerm ts)+    where+      pLeadingTerm prec (c,xs) =+        if pOptIsNegativeCoeff opt c+        then prettyParen (prec > addPrec) $+               PP.char '-' <> prettyPrintTerm opt lv (addPrec+1) (-c,xs)+        else prettyPrintTerm opt lv prec (c,xs)++      pTrailingTerm (c,xs) =+        if pOptIsNegativeCoeff opt c+        then PP.space <> PP.char '-' <> PP.space <> prettyPrintTerm opt lv (addPrec+1) (-c,xs)+        else PP.space <> PP.char '+' <> PP.space <> prettyPrintTerm opt lv (addPrec+1) (c,xs)++prettyPrintTerm+  :: (Ord k, Num k, Ord v)+  => PrintOptions k v+  -> PrettyLevel -> Rational -> Term k v -> Doc+prettyPrintTerm opt lv prec (c,xs)+  | len == 0  = pOptPrintCoeff opt lv (appPrec+1) c+    -- intentionally specify (appPrec+1) to parenthesize any composite expression+  | len == 1 && c == 1 = pPow prec $ head (mindices xs)+  | otherwise =+      prettyParen (prec > mulPrec) $+        PP.hcat $ intersperse (PP.char '*') fs+    where+      len = Map.size $ mindicesMap xs+      fs  = [pOptPrintCoeff opt lv (appPrec+1) c | c /= 1] ++ [pPow (mulPrec+1) p | p <- mindices xs]+      -- intentionally specify (appPrec+1) to parenthesize any composite expression++      pPow prec (x,1) = pOptPrintVar opt lv prec x+      pPow prec (x,n) =+        prettyParen (prec > expPrec) $+          pOptPrintVar opt lv (expPrec+1) x <> PP.char '^' <> PP.integer n++class PrettyCoeff a where+  pPrintCoeff :: PrettyLevel -> Rational -> a -> Doc+  isNegativeCoeff :: a -> Bool+  isNegativeCoeff _ = False++instance PrettyCoeff Integer where+  pPrintCoeff = pPrintPrec+  isNegativeCoeff = (0>)++instance (PrettyCoeff a, Integral a) => PrettyCoeff (Ratio a) where+  pPrintCoeff lv p r+    | denominator r == 1 = pPrintCoeff lv p (numerator r)+    | otherwise = +        prettyParen (p > ratPrec) $+          pPrintCoeff lv (ratPrec+1) (numerator r) <>+          PP.char '/' <>+          pPrintCoeff lv (ratPrec+1) (denominator r)+  isNegativeCoeff x = isNegativeCoeff (numerator x)++instance PrettyCoeff (FF.PrimeField a) where+  pPrintCoeff lv p a = pPrintCoeff lv p (FF.toInteger a)+  isNegativeCoeff _  = False++instance (Num c, Ord c, PrettyCoeff c, Ord v, PrettyVar v) => PrettyCoeff (Polynomial c v) where+  pPrintCoeff = pPrintPrec++class PrettyVar a where+  pPrintVar :: PrettyLevel -> Rational -> a -> Doc++instance PrettyVar Int where+  pPrintVar _ _ n = PP.char 'x' <> PP.int n++instance PrettyVar X where+  pPrintVar _ _ X = PP.char 'x'++addPrec, mulPrec, ratPrec, expPrec, appPrec :: Rational+addPrec = 6 -- Precedence of '+'+mulPrec = 7 -- Precedence of '*'+ratPrec = 7 -- Precedence of '/'+expPrec = 8 -- Precedence of '^'+appPrec = 10 -- Precedence of function application++{--------------------------------------------------------------------+  Univariate polynomials+--------------------------------------------------------------------}++-- | Univariate polynomials over commutative ring r+type UPolynomial r = Polynomial r X++data X = X+  deriving (Eq, Ord, Bounded, Enum, Show, Read, Typeable, Data)++instance NFData X++ucmp :: MonomialOrder X+ucmp = grlex++-- | division of univariate polynomials+div :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k+div f1 f2 = fst (divMod f1 f2)++-- | division of univariate polynomials+mod :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k+mod f1 f2 = snd (divMod f1 f2)++-- | division of univariate polynomials+divMod :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, UPolynomial k)+divMod f g+  | isZero g  = error "divMod: division by zero"+  | otherwise = go 0 f+  where+    lt_g = lt ucmp g+    go !q !r+      | deg r < deg g = (q,r)+      | otherwise     = go (q + t) (r - t * g)+        where+          lt_r = lt ucmp r+          t    = fromTerm $ lt_r `tdiv` lt_g++divides :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> Bool+divides f1 f2 = f2 `mod` f1 == 0++-- | GCD of univariate polynomials+gcd :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k+gcd f1 0  = toMonic ucmp f1+gcd f1 f2 = gcd f2 (f1 `mod` f2)++-- | LCM of univariate polynomials+lcm :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> UPolynomial k+lcm _ 0 = 0+lcm 0 _ = 0+lcm f1 f2 = toMonic ucmp $ (f1 `mod` (gcd f1 f2)) * f2++-- | Extended GCD algorithm+exgcd+  :: (Eq k, Fractional k)+  => UPolynomial k+  -> UPolynomial k+  -> (UPolynomial k, UPolynomial k, UPolynomial k)+exgcd f1 f2 = f $ go f1 f2 1 0 0 1+  where+    go !r0 !r1 !s0 !s1 !t0 !t1+      | r1 == 0   = (r0, s0, t0)+      | otherwise = go r1 r2 s1 s2 t1 t2+      where+        (q, r2) = r0 `divMod` r1+        s2 = s0 - q*s1+        t2 = t0 - q*t1+    f (g,u,v)+      | lc_g == 0 = (g, u, v)+      | otherwise = (mapCoeff (/lc_g) g, mapCoeff (/lc_g) u, mapCoeff (/lc_g) v)+      where+        lc_g = lc ucmp g++-- | pseudo division+pdivMod :: (Eq r, Num r) => UPolynomial r -> UPolynomial r -> (r, UPolynomial r, UPolynomial r)+pdivMod _ 0 = error "pdivMod: division by 0"+pdivMod f g+  | deg f < deg g = (1, 0, f)+  | otherwise     = go (deg f - deg g + 1) f 0+  where+    (lc_g, lm_g) = lt ucmp g+    b = lc_g ^ (deg f - deg_g + 1)+    deg_g = deg g+    go !n !f1 !q+      | deg_g > deg f1 = (b, q, scale (lc_g ^ n) f1)+      | otherwise      = go (n - 1) (scale lc_g f1 - s * g) (q + scale (lc_g ^ (n-1)) s)+          where+            (lc_f1, lm_f1) = lt ucmp f1+            s = fromTerm (lc_f1, lm_f1 `mdiv` lm_g)++-- | pseudo quotient+pdiv :: (Eq r, Num r) => UPolynomial r -> UPolynomial r -> UPolynomial r+pdiv f g =+  case f `pdivMod` g of+    (_, q, _) -> q++-- | pseudo reminder+pmod :: (Eq r, Num r) => UPolynomial r -> UPolynomial r -> UPolynomial r+pmod _ 0 = error "pmod: division by 0"+pmod f g+  | deg f < deg g = f+  | otherwise     = go (deg f - deg g + 1) f+  where+    (lc_g, lm_g) = lt ucmp g+    deg_g = deg g+    go !n !f1+      | deg_g > deg f1 = scale (lc_g ^ n) f1+      | otherwise      = go (n - 1) (scale lc_g f1 - s * g)+          where+            (lc_f1, lm_f1) = lt ucmp f1+            s = fromTerm (lc_f1, lm_f1 `mdiv` lm_g)++-- | GCD of univariate polynomials+gcd' :: (Eq r, Integral r) => UPolynomial r -> UPolynomial r -> UPolynomial r+gcd' f1 0  = ppI f1+gcd' f1 f2 = gcd' f2 (f1 `pmod` f2)++{--------------------------------------------------------------------+  Term+--------------------------------------------------------------------}++type Term k v = (k, Monomial v)+type UTerm k = Term k X++tdeg :: Term k v -> Integer+tdeg (_,xs) = deg xs++tmult :: (Num k, Ord v) => Term k v -> Term k v -> Term k v+tmult (c1,xs1) (c2,xs2) = (c1*c2, xs1 `mmult` xs2)++tdivides :: (Fractional k, Ord v) => Term k v -> Term k v -> Bool+tdivides (_,xs1) (_,xs2) = xs1 `mdivides` xs2++tdiv :: (Fractional k, Ord v) => Term k v -> Term k v -> Term k v+tdiv (c1,xs1) (c2,xs2) = (c1 / c2, xs1 `mdiv` xs2)++tderiv :: (Eq k, Num k, Ord v) => Term k v -> v -> Term k v+tderiv (c,xs) x =+  case mderiv xs x of+    (s,ys) -> (c * fromIntegral s, ys)++tintegral :: (Eq k, Fractional k, Ord v) => Term k v -> v -> Term k v+tintegral (c,xs) x =+  case mintegral xs x of+    (s,ys) -> (c * fromRational s, ys)++{--------------------------------------------------------------------+  Monic Monomial+--------------------------------------------------------------------}++-- 本当は変数の型に応じて type family で表現を変えたい++-- | Monic monomials+newtype Monomial v = Monomial{ mindicesMap :: Map v Integer }+  deriving (Eq, Ord, Typeable)++type UMonomial = Monomial X++instance (Ord v, Show v) => Show (Monomial v) where+  showsPrec d m  = showParen (d > 10) $+    showString "mfromIndices " . shows (mindices m)++instance (NFData v) => NFData (Monomial v) where+  rnf (Monomial m) = rnf m++instance Degree (Monomial v) where+  deg (Monomial m) = sum $ Map.elems m++instance Ord v => Var (Monomial v) v where+  var x = Monomial $ Map.singleton x 1++instance Ord v => Vars (Monomial v) v where+  vars mm = Map.keysSet (mindicesMap mm)++mone :: Monomial v+mone = Monomial $ Map.empty++mfromIndices :: Ord v => [(v, Integer)] -> Monomial v+mfromIndices xs+  | any (\(_,e) -> 0>e) xs = error "mfromIndices: negative exponent"+  | otherwise = Monomial $ Map.fromListWith (+) [(x,e) | (x,e) <- xs, e > 0]++mfromIndicesMap :: Ord v => Map v Integer -> Monomial v+mfromIndicesMap m+  | any (\(_,e) -> 0>e) (Map.toList m) = error "mfromIndicesMap: negative exponent"+  | otherwise = mfromIndicesMap' m++mfromIndicesMap' :: Ord v => Map v Integer -> Monomial v+mfromIndicesMap' m = Monomial $ Map.filter (>0) m++mindices :: Ord v => Monomial v -> [(v, Integer)]+mindices = Map.toAscList . mindicesMap++mmult :: Ord v => Monomial v -> Monomial v -> Monomial v+mmult (Monomial xs1) (Monomial xs2) = mfromIndicesMap' $ Map.unionWith (+) xs1 xs2++mpow :: Ord v => Monomial v -> Integer -> Monomial v+mpow _ 0 = mone+mpow m 1 = m+mpow (Monomial xs) e+  | 0 > e     = error "mpow: negative exponent"+  | otherwise = Monomial $ Map.map (e*) xs++mdivides :: Ord v => Monomial v -> Monomial v -> Bool+mdivides (Monomial xs1) (Monomial xs2) = Map.isSubmapOfBy (<=) xs1 xs2++mdiv :: Ord v => Monomial v -> Monomial v -> Monomial v+mdiv (Monomial xs1) (Monomial xs2) = Monomial $ Map.differenceWith f xs1 xs2+  where+    f m n+      | m <= n    = Nothing+      | otherwise = Just (m - n)++mderiv :: Ord v => Monomial v -> v -> (Integer, Monomial v)+mderiv (Monomial xs) x+  | n==0      = (0, mone)+  | otherwise = (n, Monomial $ Map.update f x xs)+  where+    n = Map.findWithDefault 0 x xs+    f m+      | m <= 1    = Nothing+      | otherwise = Just $! m - 1++mintegral :: Ord v => Monomial v -> v -> (Rational, Monomial v)+mintegral (Monomial xs) x =+  (1 % fromIntegral (n + 1), Monomial $ Map.insert x (n+1) xs)+  where+    n = Map.findWithDefault 0 x xs++mlcm :: Ord v => Monomial v -> Monomial v -> Monomial v+mlcm (Monomial m1) (Monomial m2) = Monomial $ Map.unionWith max m1 m2++mgcd :: Ord v => Monomial v -> Monomial v -> Monomial v+mgcd (Monomial m1) (Monomial m2) = Monomial $ Map.intersectionWith min m1 m2++mcoprime :: Ord v => Monomial v -> Monomial v -> Bool+mcoprime m1 m2 = mgcd m1 m2 == mone++{--------------------------------------------------------------------+  Monomial Order+--------------------------------------------------------------------}++type MonomialOrder v = Monomial v -> Monomial v -> Ordering++-- | Lexicographic order+lex :: Ord v => MonomialOrder v+lex xs1 xs2 = go (mindices xs1) (mindices xs2)+  where+    go [] [] = EQ+    go [] _  = LT -- = compare 0 n2+    go _ []  = GT -- = compare n1 0+    go ((x1,n1):xs1) ((x2,n2):xs2) =+      case compare x1 x2 of+        LT -> GT -- = compare n1 0+        GT -> LT -- = compare 0 n2+        EQ -> compare n1 n2 `mappend` go xs1 xs2++-- | Reverse lexicographic order.+-- +-- Note that revlex is NOT a monomial order.+revlex :: Ord v => Monomial v -> Monomial v -> Ordering+revlex xs1 xs2 = go (Map.toDescList (mindicesMap xs1)) (Map.toDescList (mindicesMap xs2))+  where+    go [] [] = EQ+    go [] _  = GT -- = cmp 0 n2+    go _ []  = LT -- = cmp n1 0+    go ((x1,n1):xs1) ((x2,n2):xs2) =+      case compare x1 x2 of+        LT -> GT -- = cmp 0 n2+        GT -> LT -- = cmp n1 0+        EQ -> cmp n1 n2 `mappend` go xs1 xs2+    cmp n1 n2 = compare n2 n1++-- | Graded lexicographic order+grlex :: Ord v => MonomialOrder v+grlex = (compare `on` deg) `mappend` lex++-- | Graded reverse lexicographic order+grevlex :: Ord v => MonomialOrder v+grevlex = (compare `on` deg) `mappend` revlex

@@ -1,4 +1,4 @@-{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}+{-# LANGUAGE ScopedTypeVariables, BangPatterns, TypeSynonymInstances, FlexibleInstances #-} {-# OPTIONS_GHC -Wall #-} ----------------------------------------------------------------------------- -- |@@ -8,7 +8,7 @@ --  -- Maintainer  :  masahiro.sakai@gmail.com -- Stability   :  provisional--- Portability :  non-portable (ScopedTypeVariables, BangPatterns)+-- Portability :  non-portable (ScopedTypeVariables, BangPatterns, TypeSynonymInstances, FlexibleInstances) -- -- Factoriation of polynomial over a finite field. --@@ -34,13 +34,21 @@ import Data.List import Data.Set (Set) import qualified Data.Set as Set-import Data.Polynomial-import qualified Data.Polynomial.GBasis as GB+import Data.Polynomial.Base (Polynomial, UPolynomial, X (..), MonomialOrder)+import qualified Data.Polynomial.Base as P+import qualified Data.Polynomial.GroebnerBasis as GB+import qualified TypeLevel.Number.Nat as TL +instance TL.Nat p => P.Factor (UPolynomial (PrimeField p)) where+  factor = factor++instance TL.Nat p => P.SQFree (UPolynomial (PrimeField p)) where+  sqfree = sqfree+ factor :: forall k. (Ord k, FiniteField k) => UPolynomial k -> [(UPolynomial k, Integer)] factor f = do   (g,n) <- sqfree f-  if deg g > 0+  if P.deg g > 0     then do       h <- berlekamp g       return (h,n)@@ -51,9 +59,9 @@ sqfree :: forall k. (Eq k, FiniteField k) => UPolynomial k -> [(UPolynomial k, Integer)] sqfree f   | c == 1    = sqfree' f-  | otherwise = (constant c, 1) : sqfree' (mapCoeff (/c) f)+  | otherwise = (P.constant c, 1) : sqfree' (P.mapCoeff (/c) f)   where-    (c,_) = leadingTerm grlex f+    c = P.lc ucmp f  sqfree' :: forall k. (Eq k, FiniteField k) => UPolynomial k -> [(UPolynomial k, Integer)] sqfree' 0 = []@@ -62,9 +70,9 @@   | otherwise = go 1 c0 w0 []   where     p = char (undefined :: k)-    g = deriv f X-    c0 = polyGCD f g-    w0 = polyDiv f c0+    g = P.deriv f X+    c0 = P.gcd f g+    w0 = P.div f c0     go !i c w !result       | w == 1    =           if c == 1@@ -72,18 +80,21 @@           else result ++ [(h, n*p) | (h,n) <- sqfree' (polyPthRoot c)]       | otherwise = go (i+1) c' w' result'           where-            y  = polyGCD w c-            z  = w `polyDiv` y            -            c' = c `polyDiv` y+            y  = P.gcd w c+            z  = w `P.div` y            +            c' = c `P.div` y             w' = y             result' = [(z,i) | z /= 1] ++ result +ucmp :: MonomialOrder X+ucmp = P.grlex+ polyPthRoot :: forall k. (Eq k, FiniteField k) => UPolynomial k -> UPolynomial k-polyPthRoot f = assert (deriv f X == 0) $-  fromTerms [(pthRoot c, g mm) | (c,mm) <- terms f]+polyPthRoot f = assert (P.deriv f X == 0) $+  P.fromTerms [(pthRoot c, g mm) | (c,mm) <- P.terms f]   where     p = char (undefined :: k)-    g mm = mmFromList [(X, deg mm `div` p)]+    g mm = P.var X `P.mpow` (P.deg mm `div` p)  -- | Berlekamp algorithm for polynomial factorization. --@@ -99,38 +110,40 @@         where           func fi = Set.fromList $ hs2 ++ hs1             where-              hs1 = [h | k <- allValues, let h = polyGCD fi (b - constant k), deg h > 0]-              hs2 = if deg g > 0 then [g] else []+              hs1 = [h | k <- allValues, let h = P.gcd fi (b - P.constant k), P.deg h > 0]+              hs2 = if P.deg g > 0 then [g] else []                 where-                  g = fi `polyDiv` product hs1+                  g = fi `P.div` product hs1     basis = basisOfBerlekampSubalgebra f     r     = length basis  basisOfBerlekampSubalgebra :: forall k. (Ord k, FiniteField k) => UPolynomial k -> [UPolynomial k] basisOfBerlekampSubalgebra f =-  sortBy (flip compare `on` deg) $-    map (associatedMonicPolynomial grlex) $+  sortBy (flip compare `on` P.deg) $+    map (P.toMonic ucmp) $       basis   where     q    = order (undefined :: k)-    d    = deg f-    x    = var X+    d    = P.deg f+    x    = P.var X      qs :: [UPolynomial k]-    qs = [(x^(q*i)) `polyMod` f | i <- [0 .. d - 1]]+    qs = [(x^(q*i)) `P.mod` f | i <- [0 .. d - 1]] -    gb = GB.basis grlex [p3 | (p3,_) <- terms p2]+    gb :: [Polynomial k Int]+    gb = GB.basis P.grlex [p3 | (p3,_) <- P.terms p2]      p1 :: Polynomial k Int-    p1 = sum [var i * (subst qi (\X -> var (-1)) - (var (-1) ^ i)) | (i, qi) <- zip [0..] qs]+    p1 = sum [P.var i * (P.subst qi (\X -> P.var (-1)) - (P.var (-1) ^ i)) | (i, qi) <- zip [0..] qs]     p2 :: UPolynomial (Polynomial k Int)-    p2 = toUPolynomialOf p1 (-1)+    p2 = P.toUPolynomialOf p1 (-1) -    es  = [(i, reduce grlex (var i) gb) | i <- [0 .. fromIntegral d - 1]]-    vs1 = [i           | (i, gi_def) <- es, gi_def == var i]-    vs2 = [(i, gi_def) | (i, gi_def) <- es, gi_def /= var i]+    es  = [(i, P.reduce P.grlex (P.var i) gb) | i <- [0 .. fromIntegral d - 1]]+    vs1 = [i           | (i, gi_def) <- es, gi_def == P.var i]+    vs2 = [(i, gi_def) | (i, gi_def) <- es, gi_def /= P.var i] -    basis = [ x^i + sum [constant (eval (delta i) gj_def) * x^j | (j, gj_def) <- vs2] | i <- vs1 ]+    basis :: [UPolynomial k]+    basis = [ x^i + sum [P.constant (P.eval (delta i) gj_def) * x^j | (j, gj_def) <- vs2] | i <- vs1 ]       where         delta i k           | k==i      = 1

@@ -0,0 +1,147 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns, TemplateHaskell #-}+{-# OPTIONS_GHC -Wall #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Polynomial.Factorization.Hensel+-- Copyright   :  (c) Masahiro Sakai 2013+-- License     :  BSD-style+-- +-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable (ScopedTypeVariables, BangPatterns, TemplateHaskell)+--+-- References:+--+-- * <http://www.math.kobe-u.ac.jp/Asir/ca.pdf>+-- +-- * <http://www14.in.tum.de/konferenzen/Jass07/courses/1/Bulwahn/Buhlwahn_Paper.pdf>+--+-----------------------------------------------------------------------------+module Data.Polynomial.Factorization.Hensel+  ( hensel+  ) where++import Control.Exception (assert)+import Data.FiniteField+import Data.Polynomial.Base (UPolynomial, X (..))+import qualified Data.Polynomial.Base as P+import qualified TypeLevel.Number.Nat as TL++-- import Text.PrettyPrint.HughesPJClass++hensel :: forall p. TL.Nat p => UPolynomial Integer -> [UPolynomial (PrimeField p)] -> Integer -> [UPolynomial Integer]+hensel f fs1 k+  | k <= 0    = error "hensel; k <= 0"+  | otherwise = assert precondition $ go 1 (map (P.mapCoeff Data.FiniteField.toInteger) fs1)+  where+    precondition =+      P.mapCoeff fromInteger f == product fs1 && +      P.deg f == P.deg (product fs1)++    p :: Integer+    p = TL.toInt (undefined :: p)++    go :: Integer -> [UPolynomial Integer] -> [UPolynomial Integer]+    go !i fs+      | i==k      = assert (check i fs) $ fs+      | otherwise = assert (check i fs) $ go (i+1) (lift i fs)++    check :: Integer -> [UPolynomial Integer] -> Bool+    check k fs =+        and +        [ P.mapCoeff (`mod` pk) f == P.mapCoeff (`mod` pk) (product fs)+        , fs1 == map (P.mapCoeff fromInteger) fs+        , and [P.deg fi1 == P.deg fik | (fi1, fik) <- zip fs1 fs]+        ]+      where+        pk = p ^ k++    lift :: Integer -> [UPolynomial Integer] -> [UPolynomial Integer]+    lift k fs = fs'+      where+        pk  = p^k+        pk1 = p^(k+1)++        -- f ≡ product fs + p^k h  (mod p^(k+1))+        h :: UPolynomial Integer+        h = P.mapCoeff (\c -> (c `mod` pk1) `div` pk) (f - product fs)++        hs :: [UPolynomial (PrimeField p)]+        hs = prop_5_11 (map (P.mapCoeff fromInteger) fs) (P.mapCoeff fromInteger h)++        fs' :: [UPolynomial Integer]+        fs' = [ P.mapCoeff (`mod` pk1) (fi + P.constant pk * P.mapCoeff Data.FiniteField.toInteger hi)+              | (fi, hi) <- zip fs hs ]++-- http://www14.in.tum.de/konferenzen/Jass07/courses/1/Bulwahn/Buhlwahn_Paper.pdf+test_hensel :: Bool+test_hensel = and+  [ hensel f fs 2 == [x^(2::Int) + 5*x + 18, x + 5]+  , hensel f fs 3 == [x^(2::Int) + 105*x + 43, x + 30]+  , hensel f fs 4 == [x^(2::Int) + 605*x + 168, x + 30]+  ]+  where+    x :: forall k. (Eq k, Num k) => UPolynomial k+    x  = P.var X+    f :: UPolynomial Integer+    f  = x^(3::Int) + 10*x^(2::Int) - 432*x + 5040+    fs :: [UPolynomial $(primeField 5)]+    fs = [x^(2::Int)+3, x]++-- http://www.math.kobe-u.ac.jp/Asir/ca.pdf+prop_5_10 :: forall k. (Num k, Fractional k, Eq k) => [UPolynomial k] -> [UPolynomial k]+prop_5_10 fs = normalize (go fs)+  where+    check :: [UPolynomial k] -> [UPolynomial k] -> Bool+    check es fs = sum [ei * (product fs `P.div` fi) | (ei,fi) <- zip es fs] == 1++    go :: [UPolynomial k] -> [UPolynomial k]+    go [] = error "prop_5_10: empty list"+    go [fi] = assert (check [1] [fi]) [1]+    go fs@(fi : fs') = +      case P.exgcd (product fs') fi of+        (g,ei,v) ->+           assert (g == 1) $+             let es' = go fs'+                 es  = ei : map (v*) es'+             in assert (check es fs) es++    normalize :: [UPolynomial k] -> [UPolynomial k]+    normalize es = assert (check es2 fs) es2+      where+        es2 = zipWith P.mod es fs++test_prop_5_10 :: Bool+test_prop_5_10 = sum [ei * (product fs `P.div` fi) | (ei,fi) <- zip es fs] == 1+  where+    x :: UPolynomial Rational+    x = P.var P.X+    fs = [x, x+1, x+2]+    es = prop_5_10 fs++-- http://www.math.kobe-u.ac.jp/Asir/ca.pdf+prop_5_11 :: forall k. (Num k, Fractional k, Eq k, P.PrettyCoeff k, Ord k) => [UPolynomial k] -> UPolynomial k -> [UPolynomial k]+prop_5_11 fs g =+  assert (P.deg g <= P.deg (product fs)) $+  assert (P.deg c <= 0) $+  assert (check es2 fs g) $+    es2+  where+    es  = map (g*) $ prop_5_10 fs+    c   = sum [ei `P.div` fi | (ei,fi) <- zip es fs]+    es2 = case zipWith P.mod es fs of+            e2' : es2' -> e2' + c * head fs : es2'          ++    check :: [UPolynomial k] -> [UPolynomial k] -> UPolynomial k -> Bool+    check es fs g =+      sum [ei * (product fs `P.div` fi) | (ei,fi) <- zip es fs] == g &&+      and [P.deg ei <= P.deg fi | (ei,fi) <- zip es fs]++test_prop_5_11 :: Bool+test_prop_5_11 = sum [ei * (product fs `P.div` fi) | (ei,fi) <- zip es fs] == g+  where+    x :: UPolynomial Rational+    x = P.var P.X+    fs = [x, x+1, x+2]+    g  = x^(2::Int) + 1+    es = prop_5_11 fs g

@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Polynomial.Factorization.Integer@@ -7,127 +7,14 @@ --  -- Maintainer  :  masahiro.sakai@gmail.com -- Stability   :  provisional--- Portability :  non-portable (BangPatterns)------ Factoriation of integer-coefficient polynomial using Kronecker's method.------ References:------ * <http://en.wikipedia.org/wiki/Polynomial_factorization>+-- Portability :  non-portable (TypeSynonymInstances, FlexibleInstances) -- ------------------------------------------------------------------------------module Data.Polynomial.Factorization.Integer-  ( factor-  ) where--import Data.List-import Data.MultiSet (MultiSet)-import qualified Data.MultiSet as MultiSet-import Data.Numbers.Primes (primes)-import Data.Ratio-import Data.Polynomial-import qualified Data.Polynomial.Interpolation.Lagrange as Interpolation-import Util (isInteger)--factor :: UPolynomial Integer -> [(UPolynomial Integer, Integer)]-factor 0 = [(0,1)]-factor 1 = []-factor p | deg p == 0 = [(p,1)]-factor p = [(constant c, 1) | c /= 1] ++ [(q, fromIntegral m) | (q,m) <- MultiSet.toOccurList qs]-  where-    (c,qs) = normalize (cont p, factor' (pp p))--normalize :: (Integer, MultiSet (UPolynomial Integer)) -> (Integer, MultiSet (UPolynomial Integer))-normalize (c,ps) = go (MultiSet.toOccurList ps) c MultiSet.empty-  where-    go [] !c !qs = (c, qs)-    go ((p,m) : ps) !c !qs-      | deg p == 0 = go ps (c * (coeff (var X) p) ^ m) qs-      | fst (leadingTerm grlex p) < 0 = go ps (c * (-1)^m) (MultiSet.insertMany (-p) m qs)-      | otherwise = go ps c (MultiSet.insertMany p m qs)--factor' :: UPolynomial Integer -> MultiSet (UPolynomial Integer)-factor' p = go (MultiSet.singleton p) MultiSet.empty-  where-    go ps ret-      | MultiSet.null ps = ret-      | otherwise =-          case factor2 p of-            Nothing ->-              go ps' (MultiSet.insertMany p m ret)-            Just (q1,q2) ->-              go (MultiSet.insertMany q1 m $ MultiSet.insertMany q2 m ps') ret-          where-            p   = MultiSet.findMin ps-            m   = MultiSet.occur p ps-            ps' = MultiSet.deleteAll p ps--factor2 :: UPolynomial Integer -> Maybe (UPolynomial Integer, UPolynomial Integer)-factor2 p | p == var X = Nothing-factor2 p =-  case find (\(_,yi) -> yi==0) vs of-    Just (xi,_) ->-      let q1 = x - constant xi-          q2 = p' `polyDiv` mapCoeff fromInteger q1-      in Just (q1, toZ q2)-    Nothing ->-      let qs = map Interpolation.interpolate $-                  sequence [[(fromInteger xi, fromInteger z) | z <- factors yi] | (xi,yi) <- vs]-          zs = [ (q1,q2)-               | q1 <- qs, deg q1 > 0, isUPolyZ q1-               , let (q2,r) = p' `polyDivMod` q1-               , r == 0, deg q2 > 0, isUPolyZ q2-               ]-      in case zs of-           [] -> Nothing-           (q1,q2):_ -> Just (toZ q1, toZ q2)-  where-    n = (deg p `div` 2)-    xs = take (fromIntegral n + 1) xvalues-    vs = [(x, eval (\X -> x) p) | x <- xs]-    x = var X-    p' :: UPolynomial Rational-    p' = mapCoeff fromInteger p--isUPolyZ :: UPolynomial Rational -> Bool-isUPolyZ p = and [isInteger c | (c,_) <- terms p]--toZ :: Ord v => Polynomial Rational v -> Polynomial Integer v-toZ p = fromTerms [(numerator (c * fromInteger s), xs) | (c,xs) <- terms p]-  where-    s = foldl' lcm  1 [denominator c | (c,_) <- terms p]---- [0, 1, -1, 2, -2, 3, -3 ..]-xvalues :: [Integer]-xvalues = 0 : interleave [1,2..] [-1,-2..]--interleave :: [a] -> [a] -> [a]-interleave xs [] = xs-interleave [] ys     = ys-interleave (x:xs) ys = x : interleave ys xs--factors :: Integer -> [Integer]-factors 0 = []-factors x = xs ++ map negate xs-  where-    ps = primeFactors (abs x)-    xs = map product $ sequence [take (n+1) (iterate (p*) 1) | (p,n) <- MultiSet.toOccurList ps]+module Data.Polynomial.Factorization.Integer () where -primeFactors :: Integer -> MultiSet Integer-primeFactors 0 = MultiSet.empty-primeFactors n = go n primes MultiSet.empty-  where-    go :: Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer-    go 1 !_ !result = result-    go n (p:ps) !result-      | p*p > n   = MultiSet.insert n result-      | otherwise =-          case f p n of-            (m,n') -> go n' ps (MultiSet.insertMany p m result)+-- import Data.Polynomial.Factorization.Kronecker+import qualified Data.Polynomial.Base as P+import Data.Polynomial.Factorization.Zassenhaus -    f :: Integer -> Integer -> (Int, Integer)-    f p = go2 0-      where-        go2 !m !n-          | n `mod` p == 0 = go2 (m+1) (n `div` p)-          | otherwise = (m, n)+instance P.Factor (P.UPolynomial Integer) where+  factor = factor

@@ -0,0 +1,132 @@+{-# LANGUAGE BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Polynomial.Factorization.Kronecker+-- Copyright   :  (c) Masahiro Sakai 2012-2013+-- License     :  BSD-style+-- +-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable (BangPatterns)+--+-- Factoriation of integer-coefficient polynomial using Kronecker's method.+--+-- References:+--+-- * <http://en.wikipedia.org/wiki/Polynomial_factorization>+--+-----------------------------------------------------------------------------+module Data.Polynomial.Factorization.Kronecker+  ( factor+  ) where++import Data.List+import Data.MultiSet (MultiSet)+import qualified Data.MultiSet as MultiSet+import Data.Numbers.Primes (primes)+import Data.Ratio+import Data.Polynomial.Base (Polynomial, UPolynomial, X (..))+import qualified Data.Polynomial.Base as P+import qualified Data.Polynomial.Interpolation.Lagrange as Interpolation+import Util (isInteger)++factor :: UPolynomial Integer -> [(UPolynomial Integer, Integer)]+factor 0 = [(0,1)]+factor 1 = []+factor p | P.deg p == 0 = [(p,1)]+factor p = [(P.constant c, 1) | c /= 1] ++ [(q, fromIntegral m) | (q,m) <- MultiSet.toOccurList qs]+  where+    (c,qs) = normalize (P.cont p, factor' (P.pp p))++normalize :: (Integer, MultiSet (UPolynomial Integer)) -> (Integer, MultiSet (UPolynomial Integer))+normalize (c,ps) = go (MultiSet.toOccurList ps) c MultiSet.empty+  where+    go [] !c !qs = (c, qs)+    go ((p,m) : ps) !c !qs+      | P.deg p == 0 = go ps (c * (P.coeff (P.var X) p) ^ m) qs+      | P.lc P.grlex p < 0 = go ps (c * (-1)^m) (MultiSet.insertMany (-p) m qs)+      | otherwise = go ps c (MultiSet.insertMany p m qs)++factor' :: UPolynomial Integer -> MultiSet (UPolynomial Integer)+factor' p = go (MultiSet.singleton p) MultiSet.empty+  where+    go ps ret+      | MultiSet.null ps = ret+      | otherwise =+          case factor2 p of+            Nothing ->+              go ps' (MultiSet.insertMany p m ret)+            Just (q1,q2) ->+              go (MultiSet.insertMany q1 m $ MultiSet.insertMany q2 m ps') ret+          where+            p   = MultiSet.findMin ps+            m   = MultiSet.occur p ps+            ps' = MultiSet.deleteAll p ps++factor2 :: UPolynomial Integer -> Maybe (UPolynomial Integer, UPolynomial Integer)+factor2 p | p == P.var X = Nothing+factor2 p =+  case find (\(_,yi) -> yi==0) vs of+    Just (xi,_) ->+      let q1 = x - P.constant xi+          q2 = p' `P.div` P.mapCoeff fromInteger q1+      in Just (q1, toZ q2)+    Nothing ->+      let qs = map Interpolation.interpolate $+                  sequence [[(fromInteger xi, fromInteger z) | z <- factors yi] | (xi,yi) <- vs]+          zs = [ (q1,q2)+               | q1 <- qs, P.deg q1 > 0, isUPolyZ q1+               , let (q2,r) = p' `P.divMod` q1+               , r == 0, P.deg q2 > 0, isUPolyZ q2+               ]+      in case zs of+           [] -> Nothing+           (q1,q2):_ -> Just (toZ q1, toZ q2)+  where+    n = P.deg p `div` 2+    xs = take (fromIntegral n + 1) xvalues+    vs = [(x, P.eval (\X -> x) p) | x <- xs]+    x = P.var X+    p' :: UPolynomial Rational+    p' = P.mapCoeff fromInteger p++isUPolyZ :: UPolynomial Rational -> Bool+isUPolyZ p = and [isInteger c | (c,_) <- P.terms p]++toZ :: Ord v => Polynomial Rational v -> Polynomial Integer v+toZ = P.mapCoeff numerator . P.pp++-- [0, 1, -1, 2, -2, 3, -3 ..]+xvalues :: [Integer]+xvalues = 0 : interleave [1,2..] [-1,-2..]++interleave :: [a] -> [a] -> [a]+interleave xs [] = xs+interleave [] ys     = ys+interleave (x:xs) ys = x : interleave ys xs++factors :: Integer -> [Integer]+factors 0 = []+factors x = xs ++ map negate xs+  where+    ps = primeFactors (abs x)+    xs = map product $ sequence [take (n+1) (iterate (p*) 1) | (p,n) <- MultiSet.toOccurList ps]++primeFactors :: Integer -> MultiSet Integer+primeFactors 0 = MultiSet.empty+primeFactors n = go n primes MultiSet.empty+  where+    go :: Integer -> [Integer] -> MultiSet Integer -> MultiSet Integer+    go 1 !_ !result = result+    go n (p:ps) !result+      | p*p > n   = MultiSet.insert n result+      | otherwise =+          case f p n of+            (m,n') -> go n' ps (MultiSet.insertMany p m result)++    f :: Integer -> Integer -> (Int, Integer)+    f p = go2 0+      where+        go2 !m !n+          | n `mod` p == 0 = go2 (m+1) (n `div` p)+          | otherwise = (m, n)

@@ -1,18 +1,16 @@-module Data.Polynomial.Factorization.Rational-  ( factor-  ) where+{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}+module Data.Polynomial.Factorization.Rational () where  import Data.List (foldl')-import Data.Polynomial-import qualified Data.Polynomial.Factorization.Integer as FactorZ+import Data.Polynomial.Base (UPolynomial)+import qualified Data.Polynomial.Base as P+import Data.Polynomial.Factorization.Integer () import Data.Ratio -factor :: UPolynomial Rational -> [(UPolynomial Rational, Integer)]-factor 0 = [(0,1)]-factor p = [(constant c, 1) | c /= 1] ++ qs2-  where-    s   = foldl' lcm  1 [denominator c | (c,_) <- terms p]-    p'  = mapCoeff (\c -> numerator (c * fromInteger s)) p-    qs  = FactorZ.factor p'-    qs2 = [(mapCoeff fromInteger q, m) | (q,m) <- qs, deg q > 0]-    c   = toRational (product [(coeff mmOne q)^m | (q,m) <- qs, deg q == 0]) / toRational s+instance P.Factor (UPolynomial Rational) where+  factor 0 = [(0,1)]+  factor p = [(P.constant c, 1) | c /= 1] ++ qs2+    where+      qs  = P.factor $ P.mapCoeff numerator $ P.pp p+      qs2 = [(P.mapCoeff fromInteger q, m) | (q,m) <- qs, P.deg q > 0]+      c   = toRational (product [(P.coeff P.mone q)^m | (q,m) <- qs, P.deg q == 0]) * P.cont p

@@ -1,4 +1,4 @@-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BangPatterns, TypeSynonymInstances, FlexibleInstances #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Polynomial.Factorization.SquareFree@@ -7,9 +7,7 @@ --  -- Maintainer  :  masahiro.sakai@gmail.com -- Stability   :  provisional--- Portability :  non-portable (BangPatterns)------ Square-free decomposition of univariate polynomials over a field of characteristic 0.+-- Portability :  non-portable (BangPatterns, TypeSynonymInstances, FlexibleInstances) -- -- References: --@@ -17,24 +15,26 @@ -- ----------------------------------------------------------------------------- module Data.Polynomial.Factorization.SquareFree-  ( sqfree+  ( sqfreeChar0   ) where  import Control.Exception-import Data.Polynomial+import Data.Polynomial.Base (UPolynomial, X (..))+import qualified Data.Polynomial.Base as P+import Data.Ratio  -- | Square-free decomposition of univariate polynomials over a field of characteristic 0.-sqfree :: (Eq k, Fractional k) => UPolynomial k -> [(UPolynomial k, Integer)]-sqfree 0 = []-sqfree p = assert (product [q^m | (q,m) <- result] == p) $ result+sqfreeChar0 :: (Eq k, Fractional k) => UPolynomial k -> [(UPolynomial k, Integer)]+sqfreeChar0 0 = []+sqfreeChar0 p = assert (product [q^m | (q,m) <- result] == p) $ result   where-    result = go p (p `polyDiv` polyGCD p (deriv p X)) 0 []+    result = go p (p `P.div` P.gcd p (P.deriv p X)) 0 []     go p flat !m result-      | deg flat <= 0 = [(p,1) | p /= 1] ++ reverse result-      | otherwise     = go p' flat' m' ((flat `polyDiv` flat', m') : result)+      | P.deg flat <= 0 = [(p,1) | p /= 1] ++ reverse result+      | otherwise     = go p' flat' m' ((flat `P.div` flat', m') : result)           where             (p',n) = f p flat-            flat'  = polyGCD p' flat+            flat'  = P.gcd p' flat             m' = m + n  f :: (Eq k, Fractional k) => UPolynomial k -> UPolynomial k -> (UPolynomial k, Integer)@@ -42,6 +42,21 @@   where     result@(q, m) = go 0 p1     go !m p =-      case p `polyDivMod` p2 of+      case p `P.divMod` p2 of         (q, 0) -> go (m+1) q         _ -> (p, m)+++instance P.SQFree (UPolynomial Rational) where+  sqfree = sqfreeChar0++instance P.SQFree (UPolynomial Integer) where+  sqfree 0 = [(0,1)]+  sqfree f = go 1 [] (P.sqfree (P.mapCoeff fromIntegral f))+    where+      go !u ys [] =+        assert (denominator u == 1) $+          [(P.constant (numerator u), 1) | u /= 1] ++ ys+      go !u ys ((g,n):xs)+        | P.deg g <= 0 = go (u * P.coeff P.mone g) ys xs+        | otherwise    = go (u * (P.cont g)^n) ((P.mapCoeff numerator (P.pp g), n) : ys) xs

@@ -0,0 +1,171 @@+{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}+{-# OPTIONS_GHC -Wall #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Polynomial.Factorization.Zassenhaus+-- Copyright   :  (c) Masahiro Sakai 2012-2013+-- License     :  BSD-style+-- +-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable (BangPatterns, ScopedTypeVariables)+--+-- Factoriation of integer-coefficient polynomial using Zassenhaus algorithm.+--+-- References:+--+-- * <http://www.math.kobe-u.ac.jp/Asir/ca.pdf>+--+-----------------------------------------------------------------------------+module Data.Polynomial.Factorization.Zassenhaus+  ( factor+  ) where++import Control.Monad+import Control.Monad.ST+import Control.Exception (assert)+import Data.List+import Data.Maybe+import Data.Numbers.Primes (primes)+import Data.Ratio+import Data.STRef++import Data.Polynomial.Base (UPolynomial, X (..))+import qualified Data.Polynomial.Base as P+import Data.Polynomial.Factorization.FiniteField ()+import Data.Polynomial.Factorization.SquareFree ()+import qualified Data.Polynomial.Factorization.Hensel as Hensel++import qualified TypeLevel.Number.Nat as TL+import Data.FiniteField++-- import Text.PrettyPrint.HughesPJClass++factor :: UPolynomial Integer -> [(UPolynomial Integer, Integer)]+factor f = [(h,n) | (g,n) <- P.sqfree f, h <- if P.deg g > 0 then zassenhaus g else return g]++zassenhaus :: UPolynomial Integer -> [UPolynomial Integer]+zassenhaus f = fromJust $ msum [TL.withNat zassenhausWithP p | p <- primes]+  where+    zassenhausWithP :: forall p. TL.Nat p => p -> Maybe [UPolynomial Integer]+    zassenhausWithP _ = do+      let f_mod_p :: UPolynomial (PrimeField p)+          f_mod_p = P.mapCoeff fromInteger f+      guard $ P.deg f == P.deg f_mod_p -- 主係数を割り切らないことと同値+      guard $ P.isSquareFree f_mod_p+      let fs :: [UPolynomial (PrimeField p)]+          fs = [assert (n==1) fi | (fi,n) <- P.factor f_mod_p]+      return $ lift f fs++{-+Suppose @g@ is a factor of @f@.++From Landau-Mignotte inequality,+  @sum [abs c | (c,_) <- mapCoeff ((lc f / lc g) *) $ terms g] <= 2^(deg g) * norm2 f@ holds.++This together with @deg g <= deg f@ implies+  @all [- 2^(deg f) * norm2 f <= c <= 2^(deg f) * norm2 f | (c,_) <- terms ((lc f / lc g) * g)]@.++Choose smallest @k@ such that @p^k / 2 > 2^(deg f) * norm2 f@, so that+  @all [- (p^k)/2 < c < (p^k)/2 | (c,_) <- terms ((lc f / lc g) * g)]@.++Then it call @search@ to look for actual factorization.+-}+lift :: forall p. TL.Nat p => UPolynomial Integer -> [UPolynomial (PrimeField p)] -> [UPolynomial Integer]+lift f [_] = [f]+lift f fs  = search pk f (Hensel.hensel f fs k)+  where+    p = TL.toInt (undefined :: p)+    k, pk :: Integer+    (k,pk) = head [(k,pk) | k <- [1,2..], let pk = p^k, pk^(2::Int) > (2^(P.deg f + 1))^(2::Int) * norm2sq f]++search :: Integer -> UPolynomial Integer -> [UPolynomial Integer] -> [UPolynomial Integer]+search pk f0 fs0 = runST $ do+  let a = P.lc P.grlex f0+      m = length fs0++  fRef   <- newSTRef f0+  fsRef  <- newSTRef fs0+  retRef <- newSTRef []++  forM_ [1 .. m `div` 2] $ \l -> do+    fs <- readSTRef fsRef+    forM_ (comb fs l) $ \s -> do+      {-+          A factor @g@ of @f@ must satisfy @(lc f / lc g) * g ≡ product s (mod p^k)@ for some @s@.+          So we construct a candidate of @(lc f / lc g) * g@ from @product s@.+       -}+      let g0 = product s+          -- @g1@ is a candidate of @(lc f / lc g) * g@+          g1 :: UPolynomial Rational+          g1 = P.mapCoeff conv g0+          conv :: Integer -> Rational+          conv b = b3+            where+              b1  = (a % P.lc P.grlex g0) * fromIntegral b+              -- @b1 ≡ b2 (mod p^k)@ and @0 <= b2 < p^k@+              b2  = b1 - (fromIntegral (floor (b1 / pk') :: Integer) * pk')+              -- @b1 ≡ b2 ≡ b3 (mod p^k)@ and @-(p^k)/2 <= b3 <= (p^k)/2@+              b3  = if pk'/2 < b2 then b2 - pk' else b2+              pk' = fromIntegral pk++      f <- readSTRef fRef+      let f1 = P.mapCoeff fromInteger f++      when (P.deg g1 > 0 && g1 `P.divides` f1) $ do+        let g2 = P.mapCoeff numerator $ P.pp g1+            -- we choose leading coefficient to be positive.+            g :: UPolynomial Integer+            g = if P.lc P.grlex g2 < 0 then - g2 else g2+        writeSTRef fRef $! f `div'` g+        modifySTRef retRef (g :)+        modifySTRef fsRef (\\ s)++  f <- readSTRef fRef+  ret <- readSTRef retRef+  if f==1+    then return ret+    else return $ f : ret++-- |f|^2+norm2sq :: Num a => UPolynomial a -> a+norm2sq f = sum [c^(2::Int) | (c,_) <- P.terms f]++div' :: UPolynomial Integer -> UPolynomial Integer -> UPolynomial Integer+div' f1 f2 = assert (and [denominator c == 1 | (c,_) <- P.terms g3]) g4+  where+    g1, g2 :: UPolynomial Rational+    g1 = P.mapCoeff fromInteger f1+    g2 = P.mapCoeff fromInteger f2+    g3 = g1 `P.div` g2+    g4 = P.mapCoeff numerator g3++comb :: [a] -> Int -> [[a]]+comb _ 0      = [[]]+comb [] _     = []+comb (x:xs) n = [x:ys | ys <- comb xs (n-1)] ++ comb xs n++-- ---------------------------------------------------------------------------++test_zassenhaus :: [UPolynomial Integer]+test_zassenhaus = zassenhaus f+  where+    x = P.var X+    f = x^(4::Int) + 4++test_zassenhaus2 :: [UPolynomial Integer]+test_zassenhaus2 = zassenhaus f+  where+    x = P.var X+    f = x^(9::Int) - 15*x^(6::Int) - 87*x^(3::Int) - 125++test_foo :: [(UPolynomial Integer, Integer)]+test_foo = actual+  where+    x :: UPolynomial Integer+    x = P.var X   +    f = - (x^(5::Int) + x^(4::Int) + x^(2::Int) + x + 2)+    actual   = factor f+    expected = [(-1,1), (x^(2::Int)+x+1,1), (x^(3::Int)-x+2,1)]++-- ---------------------------------------------------------------------------

@@ -1,150 +0,0 @@-{-# LANGUAGE ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Data.Polynomial.GBasis--- Copyright   :  (c) Masahiro Sakai 2012-2013--- License     :  BSD-style--- --- Maintainer  :  masahiro.sakai@gmail.com--- Stability   :  provisional--- Portability :  non-portable (ScopedTypeVariables)--- --- Gröbner basis------ References:------ * Monomial order <http://en.wikipedia.org/wiki/Monomial_order>--- --- * Gröbner basis <http://en.wikipedia.org/wiki/Gr%C3%B6bner_basis>------ * グレブナー基底 <http://d.hatena.ne.jp/keyword/%A5%B0%A5%EC%A5%D6%A5%CA%A1%BC%B4%F0%C4%EC>------ * Gröbner Bases and Buchberger’s Algorithm <http://math.rice.edu/~cbruun/vigre/vigreHW6.pdf>------ * Docon <http://www.haskell.org/docon/>--- --------------------------------------------------------------------------------module Data.Polynomial.GBasis-  (-  -- * Options-    Options (..)-  , Strategy (..)-  , defaultOptions--  -- * Gröbner basis computation-  , basis-  , basis'-  , spolynomial-  , reduceGBasis-  ) where--import qualified Data.Set as Set-import qualified Data.Heap as H -- http://hackage.haskell.org/package/heaps-import Data.Polynomial--data Options-  = Options-  { optStrategy :: Strategy-  }--defaultOptions :: Options-defaultOptions =-  Options-  { optStrategy = NormalStrategy-  }--data Strategy-  = NormalStrategy-  | SugarStrategy  -- ^ sugar strategy (not implemented yet)-  deriving (Eq, Ord, Show, Read, Bounded, Enum)--spolynomial-  :: (Eq k, Fractional k, Ord v)-  => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Polynomial k v-spolynomial cmp f g =-      fromMonomial ((1,xs) `monomialDiv` (c1,xs1)) * f-    - fromMonomial ((1,xs) `monomialDiv` (c2,xs2)) * g-  where-    xs = mmLCM xs1 xs2-    (c1, xs1) = leadingTerm cmp f-    (c2, xs2) = leadingTerm cmp g--basis-  :: forall k v. (Eq k, Fractional k, Ord k, Ord v)-  => MonomialOrder v-  -> [Polynomial k v]-  -> [Polynomial k v]-basis = basis' defaultOptions--basis'-  :: forall k v. (Eq k, Fractional k, Ord k, Ord v)-  => Options-  -> MonomialOrder v-  -> [Polynomial k v]-  -> [Polynomial k v]-basis' opt cmp fs =-  reduceGBasis cmp $ go fs (H.fromList [item cmp fi fj | (fi,fj) <- pairs fs, checkGCD fi fj])-  where-    go :: [Polynomial k v] -> H.Heap (Item k v) -> [Polynomial k v]-    go gs h | H.null h = gs-    go gs h-      | r == 0    = go gs h'-      | otherwise = go (r:gs) (H.union h' (H.fromList [item cmp r g | g <- gs, checkGCD fi fj]))-      where-        Just (i, h') = H.viewMin h-        fi = iFst i-        fj = iSnd i-        spoly = spolynomial cmp fi fj-        r = reduce cmp spoly gs--    -- gcdが1となる組は選ばなくて良い-    checkGCD fi fj = mmGCD mm1 mm2 /= mmOne-      where-        (_, mm1) = leadingTerm cmp fi-        (_, mm2) = leadingTerm cmp fj--reduceGBasis-  :: forall k v. (Eq k, Ord k, Fractional k, Ord v)-  => MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]-reduceGBasis cmp ps = Set.toList $ Set.fromList $ go ps []-  where-    go [] qs = qs-    go (p:ps) qs-      | q == 0    = go ps qs-      | otherwise = go ps (constant (1/c) * q : qs)-      where-        q = reduce cmp p (ps++qs)-        (c,_) = leadingTerm cmp q--{---------------------------------------------------------------------  Item---------------------------------------------------------------------}--data Item k v-  = Item-  { iFst :: Polynomial k v-  , iSnd :: Polynomial k v-  , iCmp :: MonomialOrder v-  , iLCM :: MonicMonomial v-  }--item :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v-item cmp f g = Item f g cmp (mmLCM mm1 mm2)-  where-    (_, mm1) = leadingTerm cmp f-    (_, mm2) = leadingTerm cmp g--instance Ord v => Ord (Item k v) where-  a `compare` b = iCmp a (iLCM a) (iLCM b)--instance Ord v => Eq (Item k v) where-  a == b = compare a b == EQ--{---------------------------------------------------------------------  Utilities---------------------------------------------------------------------}--pairs :: [a] -> [(a,a)]-pairs [] = []-pairs (x:xs) = [(x,y) | y <- xs] ++ pairs xs

@@ -0,0 +1,144 @@+{-# LANGUAGE ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Polynomial.GroebnerBasis+-- Copyright   :  (c) Masahiro Sakai 2012-2013+-- License     :  BSD-style+-- +-- Maintainer  :  masahiro.sakai@gmail.com+-- Stability   :  provisional+-- Portability :  non-portable (ScopedTypeVariables)+-- +-- Gröbner basis+--+-- References:+--+-- * Monomial order <http://en.wikipedia.org/wiki/Monomial_order>+-- +-- * Gröbner basis <http://en.wikipedia.org/wiki/Gr%C3%B6bner_basis>+--+-- * グレブナー基底 <http://d.hatena.ne.jp/keyword/%A5%B0%A5%EC%A5%D6%A5%CA%A1%BC%B4%F0%C4%EC>+--+-- * Gröbner Bases and Buchberger’s Algorithm <http://math.rice.edu/~cbruun/vigre/vigreHW6.pdf>+--+-- * Docon <http://www.haskell.org/docon/>+-- +-----------------------------------------------------------------------------++module Data.Polynomial.GroebnerBasis+  (+  -- * Options+    Options (..)+  , Strategy (..)+  , defaultOptions++  -- * Gröbner basis computation+  , basis+  , basis'+  , spolynomial+  , reduceGBasis+  ) where++import qualified Data.Set as Set+import qualified Data.Heap as H -- http://hackage.haskell.org/package/heaps+import Data.Polynomial.Base (Polynomial, Monomial, MonomialOrder)+import qualified Data.Polynomial.Base as P++data Options+  = Options+  { optStrategy :: Strategy+  }++defaultOptions :: Options+defaultOptions =+  Options+  { optStrategy = NormalStrategy+  }++data Strategy+  = NormalStrategy+  | SugarStrategy  -- ^ sugar strategy (not implemented yet)+  deriving (Eq, Ord, Show, Read, Bounded, Enum)++spolynomial+  :: (Eq k, Fractional k, Ord v)+  => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Polynomial k v+spolynomial cmp f g =+      P.fromTerm ((1,xs) `P.tdiv` lt1) * f+    - P.fromTerm ((1,xs) `P.tdiv` lt2) * g+  where+    xs = P.mlcm xs1 xs2+    lt1@(c1, xs1) = P.lt cmp f+    lt2@(c2, xs2) = P.lt cmp g++basis+  :: forall k v. (Eq k, Fractional k, Ord k, Ord v)+  => MonomialOrder v+  -> [Polynomial k v]+  -> [Polynomial k v]+basis = basis' defaultOptions++basis'+  :: forall k v. (Eq k, Fractional k, Ord k, Ord v)+  => Options+  -> MonomialOrder v+  -> [Polynomial k v]+  -> [Polynomial k v]+basis' opt cmp fs =+  reduceGBasis cmp $ go fs (H.fromList [item cmp fi fj | (fi,fj) <- pairs fs, checkGCD fi fj])+  where+    go :: [Polynomial k v] -> H.Heap (Item k v) -> [Polynomial k v]+    go gs h | H.null h = gs+    go gs h+      | r == 0    = go gs h'+      | otherwise = go (r:gs) (H.union h' (H.fromList [item cmp r g | g <- gs, checkGCD  r g]))+      where+        Just (i, h') = H.viewMin h+        fi = iFst i+        fj = iSnd i+        spoly = spolynomial cmp fi fj+        r = P.reduce cmp spoly gs++    -- gcdが1となる組は選ばなくて良い+    checkGCD fi fj = not $ P.mcoprime (P.lm cmp fi) (P.lm cmp fj)++reduceGBasis+  :: forall k v. (Eq k, Ord k, Fractional k, Ord v)+  => MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]+reduceGBasis cmp ps = Set.toList $ Set.fromList $ go ps []+  where+    go [] qs = qs+    go (p:ps) qs+      | q == 0    = go ps qs+      | otherwise = go ps (P.toMonic cmp q : qs)+      where+        q = P.reduce cmp p (ps++qs)++{--------------------------------------------------------------------+  Item+--------------------------------------------------------------------}++data Item k v+  = Item+  { iFst :: Polynomial k v+  , iSnd :: Polynomial k v+  , iCmp :: MonomialOrder v+  , iLCM :: Monomial v+  }++item :: (Eq k, Num k, Ord v) => MonomialOrder v -> Polynomial k v -> Polynomial k v -> Item k v+item cmp f g = Item f g cmp (P.mlcm (P.lm cmp f) (P.lm cmp g))++instance Ord v => Ord (Item k v) where+  a `compare` b = iCmp a (iLCM a) (iLCM b)++instance Ord v => Eq (Item k v) where+  a == b = compare a b == EQ++{--------------------------------------------------------------------+  Utilities+--------------------------------------------------------------------}++pairs :: [a] -> [(a,a)]+pairs [] = []+pairs (x:xs) = [(x,y) | y <- xs] ++ pairs xs

@@ -3,11 +3,12 @@   ( interpolate   ) where -import Data.Polynomial+import Data.Polynomial (UPolynomial, X (..))+import qualified Data.Polynomial as P  interpolate :: (Eq k, Fractional k) => [(k,k)] -> UPolynomial k interpolate zs = sum $ do   (xj,yj) <- zs-  let lj x = product [constant (1 / (xj - xm)) * (x - constant xm) | (xm,_) <- zs, xj /= xm]-  let x = var X-  return $ constant yj * lj x+  let lj x = product [P.constant (1 / (xj - xm)) * (x - P.constant xm) | (xm,_) <- zs, xj /= xm]+  let x = P.var X+  return $ P.constant yj * lj x

@@ -25,7 +25,8 @@  import Control.Exception import qualified Data.IntMap as IM-import Data.Polynomial+import Data.Polynomial (UPolynomial, X (..))+import qualified Data.Polynomial as P  data NthRoot = NthRoot !Integer !Rational   deriving (Show)@@ -33,9 +34,9 @@ graeffesMethod :: UPolynomial Rational -> Int -> [NthRoot] graeffesMethod p v = xs !! (v - 1)   where-    xs = map (uncurry g) $ zip [1..] (tail $ iterate f $ associatedMonicPolynomial grlex p)+    xs = map (uncurry g) $ zip [1..] (tail $ iterate f $ P.toMonic P.grlex p) -    n = deg p+    n = P.deg p      g :: Int -> UPolynomial Rational -> [NthRoot]     g v p = do@@ -43,22 +44,22 @@       let yi = if i == 1 then - (b i) else - (b i / b (i-1))       return $ NthRoot (2 ^ fromIntegral v) yi       where-        bs = IM.fromList [(fromInteger i, b) | (b,ys) <- terms p, let i = n - deg ys, i /= 0]+        bs = IM.fromList [(fromInteger i, b) | (b,ys) <- P.terms p, let i = n - P.deg ys, i /= 0]         b i = IM.findWithDefault 0 i bs  f :: UPolynomial Rational -> UPolynomial Rational-f p = (-1) ^ (deg p) *-      fromTerms [ (c, mmFromList [assert (e `mod` 2 == 0) (x, e `div` 2) | (x,e) <- mmToList xs])-                | (c,xs) <- terms (p * subst p (\_ -> - var X)) ]+f p = (-1) ^ (P.deg p) *+      P.fromTerms [ (c, assert (P.deg xs `mod` 2 == 0) (P.var X `P.mpow` (P.deg xs `div` 2)))+                  | (c, xs) <- P.terms (p * P.subst p (\X -> - P.var X)) ]  f' :: UPolynomial Rational -> UPolynomial Rational-f' p = fromTerms [(b k, mmFromList [(X, n - k)]) | k <- [0..n]]+f' p = P.fromTerms [(b k, P.var X `P.mpow` (n - k)) | k <- [0..n]]   where-    n = deg p+    n = P.deg p      a :: Integer -> Rational     a k-      | n >= k    = coeff (mmFromList [(X, n - k)]) p+      | n >= k    = P.coeff (P.var X `P.mpow` (n - k)) p       | otherwise = 0      b :: Integer -> Rational@@ -66,10 +67,10 @@  test v = graeffesMethod p v   where-    x = var X+    x = P.var X     p = x^2 - 2  test2 v = graeffesMethod p v  where-    x = var X+    x = P.var X     p = x^5 - 3*x - 1

@@ -34,7 +34,8 @@   ) where  import Data.Maybe-import Data.Polynomial+import Data.Polynomial (UPolynomial, X (..))+import qualified Data.Polynomial as P import qualified Data.Interval as Interval import Data.Interval (Interval, EndPoint (..), (<..<=), (<=..<=)) @@ -46,10 +47,10 @@ sturmChain p = p0 : p1 : go p0 p1   where     p0 = p-    p1 = deriv p X+    p1 = P.deriv p P.X     go p q = if r==0 then [] else r : go q r       where-        r = - (p `polyMod` q)+        r = - (p `P.mod` q)  -- | The number of distinct real roots of @p@ in a given interval numRoots@@ -70,12 +71,12 @@       case (Interval.lowerBound ival2, Interval.upperBound ival2) of         (Finite lb, Finite ub) ->           (if lb==ub then 0 else (n lb - n ub)) +-          (if lb `Interval.member` ival2 && isRootOf lb p then 1 else 0) +-          (if ub `Interval.notMember` ival2 && isRootOf ub p then -1 else 0)+          (if lb `Interval.member` ival2 && lb `P.isRootOf` p then 1 else 0) ++          (if ub `Interval.notMember` ival2 && ub `P.isRootOf`  p then -1 else 0)         _ -> error "numRoots'': should not happen"   where     ival2 = boundInterval p ival-    n x = countSignChanges [eval (\X -> x) q | q <- chain]+    n x = countSignChanges [P.eval (\X -> x) q | q <- chain]  countSignChanges :: [Rational] -> Int countSignChanges rs = countChanges xs@@ -100,8 +101,8 @@   where     m = if p==0         then 0-        else max 1 (sum [abs (c/s) | (c,_) <- terms p] - 1)-    (s,_) = leadingTerm grlex p+        else max 1 (sum [abs (c/s) | (c,_) <- P.terms p] - 1)+    s = P.lc P.grlex p  boundInterval :: UPolynomial Rational -> Interval Rational -> Interval Rational boundInterval p ival = Interval.intersection ival (Finite lb <=..<= Finite ub)@@ -119,10 +120,10 @@ separate' :: SturmChain -> [Interval Rational] separate' chain@(p:_) = f (bounds p)   where-    n x = countSignChanges [eval (\X -> x) q | q <- chain]+    n x = countSignChanges [P.eval (\X -> x) q | q <- chain]      f (lb,ub) =-      if lb `isRootOf` p+      if lb `P.isRootOf` p       then Interval.singleton lb : g (lb,ub)       else g (lb,ub)     

@@ -1,4 +1,4 @@-{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleInstances, DeriveDataTypeable, CPP #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Sign@@ -7,7 +7,7 @@ --  -- Maintainer  :  masahiro.sakai@gmail.com -- Stability   :  provisional--- Portability :  non-portable (DeriveDataTypeable)+-- Portability :  non-portable (FlexibleInstances, DeriveDataTypeable, CPP) -- -- Algebra of Signs. --@@ -16,17 +16,22 @@   (   -- * Algebra of Sign     Sign (..)-  , signNegate-  , signMul-  , signRecip-  , signDiv-  , signPow+  , negate+  , mult+  , recip+  , div+  , pow   , signOf-  , showSign+  , symbol   ) where -import Algebra.Enumerable (Enumerable (..)) -- from lattices package+import Prelude hiding (negate, recip, div)+import Algebra.Enumerable (Enumerable (..), universeBounded) -- from lattices package+import qualified Algebra.Lattice as L -- from lattices package import Control.DeepSeq+import Data.Hashable+import Data.Set (Set)+import qualified Data.Set as Set import Data.Typeable import Data.Data import qualified Numeric.Algebra as Alg@@ -36,52 +41,54 @@  instance NFData Sign +instance Hashable Sign where hashWithSalt = hashUsing fromEnum+ instance Enumerable Sign where-  universe = [Neg .. Pos]+  universe = universeBounded  instance Alg.Multiplicative Sign where-  (*)   = signMul-  pow1p = signPow+  (*)   = mult+  pow1p s n = pow s (1+n)  instance Alg.Commutative Sign  instance Alg.Unital Sign where   one = Pos-  pow = signPow+  pow = pow  instance Alg.Division Sign where-  recip = signRecip-  (/)   = signDiv-  (\\)  = flip signDiv-  (^)   = signPow+  recip = recip+  (/)   = div+  (\\)  = flip div+  (^)   = pow -signNegate :: Sign -> Sign-signNegate Neg  = Pos-signNegate Zero = Zero-signNegate Pos  = Neg+negate :: Sign -> Sign+negate Neg  = Pos+negate Zero = Zero+negate Pos  = Neg -signMul :: Sign -> Sign -> Sign-signMul Pos s  = s-signMul s Pos  = s-signMul Neg s  = signNegate s-signMul s Neg  = signNegate s-signMul _ _    = Zero+mult :: Sign -> Sign -> Sign+mult Pos s  = s+mult s Pos  = s+mult Neg s  = negate s+mult s Neg  = negate s+mult _ _    = Zero -signRecip :: Sign -> Sign-signRecip Pos  = Pos-signRecip Zero = error "signRecip: division by Zero"-signRecip Neg  = Neg+recip :: Sign -> Sign+recip Pos  = Pos+recip Zero = error "Data.Sign.recip: division by Zero"+recip Neg  = Neg -signDiv :: Sign -> Sign -> Sign-signDiv s Pos  = s-signDiv _ Zero = error "signDiv: division by Zero"-signDiv s Neg  = signNegate s+div :: Sign -> Sign -> Sign+div s Pos  = s+div _ Zero = error "Data.Sign.div: division by Zero"+div s Neg  = negate s -signPow :: Integral x => Sign -> x -> Sign-signPow _ 0    = Pos-signPow Pos _  = Pos-signPow Zero _ = Zero-signPow Neg n  = if even n then Pos else Neg+pow :: Integral x => Sign -> x -> Sign+pow _ 0    = Pos+pow Pos _  = Pos+pow Zero _ = Zero+pow Neg n  = if even n then Pos else Neg  signOf :: Real a => a -> Sign signOf r =@@ -90,8 +97,45 @@     EQ -> Zero     GT -> Pos -showSign :: Sign -> String-showSign Pos  = "+"-showSign Neg  = "-"-showSign Zero = "0"+symbol :: Sign -> String+symbol Pos  = "+"+symbol Neg  = "-"+symbol Zero = "0"++instance L.MeetSemiLattice (Set Sign) where+  meet = Set.intersection++instance L.Lattice (Set Sign)++instance L.BoundedMeetSemiLattice (Set Sign) where+  top = Set.fromList universe++instance L.BoundedLattice (Set Sign)++#if !MIN_VERSION_hashable(1,2,0)+-- Copied from hashable-1.2.0.7:+-- Copyright   :  (c) Milan Straka 2010+--                (c) Johan Tibell 2011+--                (c) Bryan O'Sullivan 2011, 2012++-- | Transform a value into a 'Hashable' value, then hash the+-- transformed value using the given salt.+--+-- This is a useful shorthand in cases where a type can easily be+-- mapped to another type that is already an instance of 'Hashable'.+-- Example:+--+-- > data Foo = Foo | Bar+-- >          deriving (Enum)+-- >+-- > instance Hashable Foo where+-- >     hashWithSalt = hashUsing fromEnum+hashUsing :: (Hashable b) =>+             (a -> b)           -- ^ Transformation function.+          -> Int                -- ^ Salt.+          -> a                  -- ^ Value to transform.+          -> Int+hashUsing f salt x = hashWithSalt salt (f x)+{-# INLINE hashUsing #-}+#endif 

@@ -17,8 +17,9 @@   , Model   ) where -import qualified Data.IntMap as IM-import qualified Data.IntSet as IS+import Data.IntMap (IntMap)+import Data.IntSet (IntSet)+import qualified Data.IntSet as IntSet import Data.Ratio  -- ---------------------------------------------------------------------------@@ -27,17 +28,17 @@ type Var = Int  -- | Set of variables-type VarSet = IS.IntSet+type VarSet = IntSet  -- | Map from variables-type VarMap = IM.IntMap+type VarMap = IntMap  -- | collecting free variables class Variables a where   vars :: a -> VarSet  instance Variables a => Variables [a] where-  vars = IS.unions . map vars+  vars = IntSet.unions . map vars  -- | A @Model@ is a map from variables to values. type Model r = VarMap r

@@ -1,5 +1,10 @@ {-# OPTIONS_GHC -Wall -fno-warn-unused-do-bind #-}-{-# LANGUAGE BangPatterns, DoRec, ScopedTypeVariables, CPP, DeriveDataTypeable #-}+{-# LANGUAGE BangPatterns, ScopedTypeVariables, CPP, DeriveDataTypeable #-}+#if __GLASGOW_HASKELL__ < 706+{-# LANGUAGE DoRec #-}+#else+{-# LANGUAGE RecursiveDo #-}+#endif ----------------------------------------------------------------------------- -- | -- Module      :  SAT@@ -8,7 +13,7 @@ --  -- Maintainer  :  masahiro.sakai@gmail.com -- Stability   :  provisional--- Portability :  non-portable (BangPatterns, DoRec, ScopedTypeVariables, CPP, DeriveDataTypeable)+-- Portability :  non-portable (BangPatterns, RecursiveDo, ScopedTypeVariables, CPP, DeriveDataTypeable) -- -- A CDCL SAT solver. --

@@ -19,7 +19,8 @@   ) where  import Control.Monad-import qualified Data.IntMap as IM+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap import qualified SAT as SAT import qualified SAT.Types as SAT @@ -53,13 +54,13 @@       }  solveWBO :: SAT.Solver -> [(SAT.Lit, Integer)] -> Options -> IO (Maybe (SAT.Model, Integer))-solveWBO solver sels0 opt = loop 0 (IM.fromList sels0)+solveWBO solver sels0 opt = loop 0 (IntMap.fromList sels0)   where-    loop :: Integer -> IM.IntMap Integer -> IO (Maybe (SAT.Model, Integer))+    loop :: Integer -> IntMap Integer -> IO (Maybe (SAT.Model, Integer))     loop !lb sels = do       optUpdateLB opt lb -      ret <- SAT.solveWith solver (IM.keys sels)+      ret <- SAT.solveWith solver (IntMap.keys sels)       if ret       then do         m <- SAT.model solver@@ -71,7 +72,7 @@         case core of           [] -> return Nothing           _  -> do-            let !min_c = minimum [sels IM.! sel | sel <- core]+            let !min_c = minimum [sels IntMap.! sel | sel <- core]                 !lb' = lb + min_c              xs <- forM core $ \sel -> do@@ -80,13 +81,13 @@             SAT.addExactly solver (map snd xs) 1             SAT.addClause solver [-l | l <- core] -- optional constraint but sometimes useful -            ys <- liftM IM.unions $ forM xs $ \(sel, r) -> do+            ys <- liftM IntMap.unions $ forM xs $ \(sel, r) -> do               sel' <- SAT.newVar solver               SAT.addClause solver [-sel', r, sel]-              let c = sels IM.! sel+              let c = sels IntMap.! sel               if c > min_c-                then return $ IM.fromList [(sel', min_c), (sel, c - min_c)]-                else return $ IM.singleton sel' min_c-            let sels' = IM.union ys (IM.difference sels (IM.fromList [(sel, ()) | sel <- core]))+                then return $ IntMap.fromList [(sel', min_c), (sel, c - min_c)]+                else return $ IntMap.singleton sel' min_c+            let sels' = IntMap.union ys (IntMap.difference sels (IntMap.fromList [(sel, ()) | sel <- core]))              loop lb' sels'

@@ -50,8 +50,10 @@  import Control.Monad import Data.IORef+import Data.Map (Map) import qualified Data.Map as Map-import qualified Data.IntSet as IS+import Data.IntSet (IntSet)+import qualified Data.IntSet as IntSet import qualified SAT as SAT  -- | Arbitrary formula not restricted to CNF@@ -69,7 +71,7 @@   Encoder   { encSolver    :: SAT.Solver   , encUsePB     :: IORef Bool-  , encConjTable :: !(IORef (Map.Map IS.IntSet SAT.Lit))+  , encConjTable :: !(IORef (Map IntSet SAT.Lit))   }  -- | Create a @Encoder@ instance.@@ -161,7 +163,7 @@ encodeConj :: Encoder -> [SAT.Lit] -> IO SAT.Lit encodeConj _ [l] =  return l encodeConj encoder ls =  do-  let ls2 = IS.fromList ls+  let ls2 = IntSet.fromList ls   table <- readIORef (encConjTable encoder)   case Map.lookup ls2 table of     Just l -> return l

@@ -46,15 +46,16 @@ import Data.Array.Unboxed import Data.Ord import Data.List-import qualified Data.IntMap as IM-import qualified Data.IntSet as IS-import qualified Data.Set as Set+import Data.IntMap (IntMap)+import qualified Data.IntMap as IntMap+import Data.IntSet (IntSet)+import qualified Data.IntSet as IntSet  -- | Variable is represented as positive integers (DIMACS format). type Var = Int -type VarSet = IS.IntSet-type VarMap = IM.IntMap+type VarSet = IntSet+type VarMap = IntMap  {-# INLINE validVar #-} validVar :: Var -> Bool@@ -71,8 +72,8 @@ litUndef :: Lit litUndef = 0 -type LitSet = IS.IntSet-type LitMap = IM.IntMap+type LitSet = IntSet+type LitMap = IntMap  {-# INLINE validLit #-} validLit :: Lit -> Bool@@ -114,22 +115,22 @@ -- -- 'Nothing' if the clause is trivially true. normalizeClause :: Clause -> Maybe Clause-normalizeClause lits = assert (IS.size ys `mod` 2 == 0) $-  if IS.null ys-    then Just (IS.toList xs)+normalizeClause lits = assert (IntSet.size ys `mod` 2 == 0) $+  if IntSet.null ys+    then Just (IntSet.toList xs)     else Nothing   where-    xs = IS.fromList lits-    ys = xs `IS.intersection` (IS.map litNot xs)+    xs = IntSet.fromList lits+    ys = xs `IntSet.intersection` (IntSet.map litNot xs)  normalizeAtLeast :: ([Lit],Int) -> ([Lit],Int)-normalizeAtLeast (lits,n) = assert (IS.size ys `mod` 2 == 0) $-   (IS.toList lits', n')+normalizeAtLeast (lits,n) = assert (IntSet.size ys `mod` 2 == 0) $+   (IntSet.toList lits', n')    where-     xs = IS.fromList lits-     ys = xs `IS.intersection` (IS.map litNot xs)-     lits' = xs `IS.difference` ys-     n' = n - (IS.size ys `div` 2)+     xs = IntSet.fromList lits+     ys = xs `IntSet.intersection` (IntSet.map litNot xs)+     lits' = xs `IntSet.difference` ys+     n' = n - (IntSet.size ys `div` 2)  -- | normalizing PB term of the form /c1 x1 + c2 x2 ... cn xn + c/ into -- /d1 x1 + d2 x2 ... dm xm + d/ where d1,...,dm ≥ 1.@@ -139,15 +140,15 @@     -- 同じ変数が複数回現れないように、一度全部 @v@ に統一。     step1 :: ([(Integer,Lit)], Integer) -> ([(Integer,Lit)], Integer)     step1 (xs,n) =-      case loop (IM.empty,n) xs of-        (ys,n') -> ([(c,v) | (v,c) <- IM.toList ys], n')+      case loop (IntMap.empty,n) xs of+        (ys,n') -> ([(c,v) | (v,c) <- IntMap.toList ys], n')       where         loop :: (VarMap Integer, Integer) -> [(Integer,Lit)] -> (VarMap Integer, Integer)         loop (ys,m) [] = (ys,m)         loop (ys,m) ((c,l):zs) =           if litPolarity l-            then loop (IM.insertWith (+) l c ys, m) zs-            else loop (IM.insertWith (+) (litNot l) (negate c) ys, m+c) zs+            then loop (IntMap.insertWith (+) l c ys, m) zs+            else loop (IntMap.insertWith (+) (litNot l) (negate c) ys, m+c) zs      -- 係数が0のものも取り除き、係数が負のリテラルを反転することで、     -- 係数が正になるようにする。

@@ -3,10 +3,11 @@   , render   ) where +import Data.Map (Map) import qualified Data.Map as Map import Data.Ratio -type Model = Map.Map String Double+type Model = Map String Double  render :: Model -> Maybe Double -> String render m obj = unlines $ ls1 ++ ls2

@@ -56,7 +56,9 @@ import Data.List import Data.Maybe import Data.Ratio+import Data.Map (Map) import qualified Data.Map as Map+import Data.Set (Set) import qualified Data.Set as Set import Data.OptDir import Text.ParserCombinators.Parsec hiding (label)@@ -68,11 +70,11 @@ -- | Problem data LP   = LP-  { variables :: Set.Set Var+  { variables :: Set Var   , dir :: OptDir   , objectiveFunction :: ObjectiveFunction   , constraints :: [Constraint]-  , varInfo :: Map.Map Var VarInfo+  , varInfo :: Map Var VarInfo   , sos :: [SOS]   }   deriving (Show, Eq, Ord)@@ -155,7 +157,7 @@ type SOS = (Maybe Label, SOSType, [(Var, Rational)])  class Variables a where-  vars :: a -> Set.Set Var+  vars :: a -> Set Var  instance Variables a => Variables [a] where   vars = Set.unions . map vars@@ -199,12 +201,12 @@ intersectBounds :: Bounds -> Bounds -> Bounds intersectBounds (lb1,ub1) (lb2,ub2) = (max lb1 lb2, min ub1 ub2) -integerVariables :: LP -> Set.Set Var+integerVariables :: LP -> Set Var integerVariables lp = Map.keysSet $ Map.filter p (varInfo lp)   where     p VarInfo{ varType = vt } = vt == IntegerVariable -semiContinuousVariables :: LP -> Set.Set Var+semiContinuousVariables :: LP -> Set Var semiContinuousVariables lp = Map.keysSet $ Map.filter p (varInfo lp)   where     p VarInfo{ varType = vt } = vt == SemiContinuousVariable@@ -259,7 +261,7 @@   tok $ char ':'   return name -reserved :: Set.Set String+reserved :: Set String reserved = Set.fromList   [ "bound", "bounds"   , "gen", "general", "generals"@@ -400,7 +402,7 @@  type Bounds2 = (Maybe BoundExpr, Maybe BoundExpr) -boundsSection :: Parser (Map.Map Var Bounds)+boundsSection :: Parser (Map Var Bounds) boundsSection = do   tok $ string' "bound" >> optional (char' 's')   liftM (Map.map g . Map.fromListWith f) $ many (try bound)@@ -740,7 +742,7 @@ -- ---------------------------------------------------------------------------  {--compileExpr :: Expr -> Maybe (Map.Map Var Rational)+compileExpr :: Expr -> Maybe (Map Var Rational) compileExpr e = do   xs <- forM e $ \(Term c vs) ->     case vs of

@@ -29,7 +29,9 @@  import Control.Monad import Data.Maybe+import Data.Set (Set) import qualified Data.Set as Set+import Data.Map (Map) import qualified Data.Map as Map import Data.Ratio @@ -359,12 +361,12 @@   newline'   return (op, name) -colsSection :: Parser (Map.Map Column (Map.Map Row Rational), Set.Set Column)+colsSection :: Parser (Map Column (Map Row Rational), Set Column) colsSection = do   try $ stringLn "COLUMNS"   body False Map.empty Set.empty   where-    body :: Bool -> Map.Map Column (Map.Map Row Rational) -> Set.Set Column -> Parser (Map.Map Column (Map.Map Row Rational), Set.Set Column)+    body :: Bool -> Map Column (Map Row Rational) -> Set Column -> Parser (Map Column (Map Row Rational), Set Column)     body isInt rs ivs = msum       [ do isInt' <- try intMarker            body isInt' rs ivs@@ -386,7 +388,7 @@       newline'       return b -    entry :: Parser (Column, Map.Map Row Rational)+    entry :: Parser (Column, Map Row Rational)     entry = do       spaces1'       col <- ident@@ -397,13 +399,13 @@         Nothing -> return (col, rv1)         Just rv2 ->  return (col, Map.union rv1 rv2) -rowAndVal :: Parser (Map.Map Row Rational)+rowAndVal :: Parser (Map Row Rational) rowAndVal = do   row <- ident   val <- number   return $ Map.singleton row val -rhsSection :: Parser (Map.Map Row Rational)+rhsSection :: Parser (Map Row Rational) rhsSection = do   try $ stringLn "RHS"   liftM Map.unions $ many entry@@ -418,7 +420,7 @@         Nothing  -> return rv1         Just rv2 -> return $ Map.union rv1 rv2 -rangesSection :: Parser (Map.Map Row Rational)+rangesSection :: Parser (Map Row Rational) rangesSection = do   try $ stringLn "RANGES"   liftM Map.unions $ many entry@@ -519,7 +521,7 @@       newline'       return $ LPFile.Term val [col1, col2] -indicatorsSection :: Parser (Map.Map Row (Column, Rational))+indicatorsSection :: Parser (Map Row (Column, Rational)) indicatorsSection = do   try $ stringLn "INDICATORS"   liftM Map.fromList $ many entry