hermite-normal-form: Fix bug; add tests

privet-kitty · Nov 3, 2023 · 7b22847 · 7b22847
1 parent dfd347a
commit 7b22847
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 0 deletions.
diff --git a/module/hermite-normal-form.lisp b/module/hermite-normal-form.lisp
@@ -393,6 +393,7 @@ returns (VALUES NIL NIL HNF)."
            ((array * (*)) b))
   (destructuring-bind (m n) (array-dimensions a)
     (declare ((mod #.array-dimension-limit) m n))
+    (assert (= m (length b)))
     (let* ((hnf (hnf a t))
            (h (hnf-matrix hnf)) ; AU = H
            (u (hnf-unimodular-matrix hnf))
@@ -405,6 +406,9 @@ returns (VALUES NIL NIL HNF)."
                 do (incf row))
           (setf (aref pivot-rows col) row)))
       ;; solve Hy = b
+      (dotimes (row (if (zerop rank) m (aref pivot-rows 0)))
+        (unless (zerop (the integer (aref b row)))
+          (return-from solve-integer-linear-system (values nil nil hnf))))
       (dotimes (col rank)
         (let ((pivot-row (aref pivot-rows col)))
           (multiple-value-bind (quot rem)

diff --git a/module/test/hermite-normal-form.lisp b/module/test/hermite-normal-form.lisp
@@ -302,6 +302,54 @@
                                           'simple-vector)))
                             (is (equalp restored-vector target-vector)))))))))))
 
+(test solve-integer-linear-system/hand
+  ;; zero-size case
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a() #())
+    (is (equalp intercepts #()))
+    (is (equalp coefs #2a())))
+  (multiple-value-bind (intercepts coefs)
+      (solve-integer-linear-system (make-array '(0 3)) #())
+    (is (equalp intercepts #(0 0 0)))
+    (is (equalp coefs #2a((1 0 0) (0 1 0) (0 0 1)))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a(() ()) #(0 0))
+    (is (equalp intercepts #()))
+    (is (equalp coefs #2a())))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a(() ()) #(0 1))
+    (is (null intercepts))
+    (is (null coefs)))  
+  ;; one-size case
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2)) #(3))
+    (is (null intercepts))
+    (is (null coefs)))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2)) #(4))
+    (is (equalp intercepts #(2)))
+    (is (equalp coefs #2a(()))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2 3)) #(-1))
+    (is (equalp intercepts #(1 -1)))
+    (is (equalp coefs #2a((-3) (2)))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2) (3)) #(-2 -3))
+    (is (equalp intercepts #(-1)))
+    (is (equalp coefs #2a(()))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((0 0 0)) #(0))
+    (is (equalp intercepts #(0 0 0)))
+    (is (equalp coefs #2a((1 0 0) (0 1 0) (0 0 1)))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((0 0 0)) #(1))
+    (is (null intercepts))
+    (is (null coefs)))
+  ;; more size
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2 3) (0 0)) #(-1 0))
+    (is (equalp intercepts #(1 -1)))
+    (is (equalp coefs #2a((-3) (2)))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2 3) (4 6)) #(-1 -2))
+    (is (equalp intercepts #(1 -1)))
+    (is (equalp coefs #2a((-3) (2)))))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2 3) (4 6)) #(-1 -3))
+    (is (null intercepts))
+    (is (null coefs)))
+  (multiple-value-bind (intercepts coefs) (solve-integer-linear-system #2a((2 3) (1 -1)) #(-1 -3))
+    (is (equalp intercepts #(-2 1)))
+    (is (equalp coefs #2a(() ())))))
+
 (test hnf-p
   (declare (notinline hnf-p))
   (is (eql 0 (hnf-p #2a())))