Skip to content

Commit

Permalink
Add integer LLL
Browse files Browse the repository at this point in the history
  • Loading branch information
@privet-kitty
privet-kitty committed Nov 23, 2023
1 parent b13f23d commit cae19e4
Showing 1 changed file with 158 additions and 18 deletions.
176 changes: 158 additions & 18 deletions module/lll.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,28 +94,29 @@ NOTE: MAT should have full column rank. Otherwise the consequence is undefined."
(rotatef (aref mat i (- pos 1)) (aref mat i pos)))
(dotimes (j (- pos 1))
(rotatef (aref coefs (- pos 1) j) (aref coefs pos j)))
(let* ((mu-k-k-1 (aref coefs pos (- pos 1)))
(let* ((old-mu (aref coefs pos (- pos 1)))
(new-l2 (+ (aref l2s pos)
(* (expt mu-k-k-1 2) (aref l2s (- pos 1))))))
(declare (rational mu-k-k-1 new-l2))
(* (expt old-mu 2) (aref l2s (- pos 1))))))
(declare (rational old-mu new-l2))
(setf (aref coefs pos (- pos 1))
(/ (* mu-k-k-1 (aref l2s (- pos 1))) new-l2))
(/ (* old-mu (aref l2s (- pos 1))) new-l2))
(setf (aref l2s pos)
(/ (* (aref l2s (- pos 1)) (aref l2s pos)) new-l2))
(setf (aref l2s (- pos 1)) new-l2)
(dotimes (i m)
(let* ((old-ort (aref ort-mat i (- pos 1)))
(new-ort (+ (aref ort-mat i pos)
(* old-mu old-ort))))
(setf (aref ort-mat i (- pos 1)) new-ort)
(setf (aref ort-mat i pos)
(- old-ort (* (aref coefs pos (- pos 1)) new-ort)))))
(loop for j from (+ pos 1) to max-pos
for mu = (aref coefs j pos)
do (setf (aref coefs j pos)
(- (aref coefs j (- pos 1)) (* mu-k-k-1 mu)))
(- (aref coefs j (- pos 1)) (* old-mu mu)))
(setf (aref coefs j (- pos 1))
(+ mu (* (aref coefs pos (- pos 1))
(aref coefs j pos)))))
(dotimes (i m)
(let* ((old (aref ort-mat i (- pos 1)))
(new (+ (aref ort-mat i pos)
(* mu-k-k-1 old))))
(setf (aref ort-mat i (- pos 1)) new
(aref ort-mat i pos) (- old (* (aref coefs pos (- pos 1)) new)))))
(aref coefs j pos)))))
;; (%debug)
)))
(loop
Expand All @@ -131,12 +132,11 @@ NOTE: MAT should have full column rank. Otherwise the consequence is undefined."
(dotimes (i m)
(setf (aref ort-mat i pos) (aref mat i pos)))
(dotimes (j pos)
(setf (aref coefs pos j)
(/ (loop for i below m
sum (* (aref mat i pos) (aref ort-mat i j))
of-type rational)
(aref l2s j)))
(let ((mu (aref coefs pos j)))
(let ((mu (/ (loop for i below m
sum (* (aref mat i pos) (aref ort-mat i j))
of-type rational)
(aref l2s j))))
(setf (aref coefs pos j) mu)
(dotimes (i m)
(decf (aref ort-mat i pos) (* mu (aref ort-mat i j))))))
(setf (aref l2s pos)
Expand All @@ -159,3 +159,143 @@ NOTE: MAT should have full column rank. Otherwise the consequence is undefined."
(incf pos)
(return))))))
mat))))

;; FIXME: This is just for validating this implementation. I will discard it
;; when I believe this module is stable.
(declaim (inline %div))
(defun %div (x y)
(declare (integer x y))
(multiple-value-bind (quot rem) (floor x y)
(assert (zerop rem))
quot))

(defun lll (mat &optional (delta 3/4))
"Applies the basis reduction to the given column vectors.
NOTE: MAT should have full column rank. Otherwise the consequence is undefined."
(declare (optimize (speed 3))
((array * (* *)) mat)
((rational (1/4) 1) delta))
(destructuring-bind (m n) (array-dimensions mat)
(declare ((mod #.array-dimension-limit) m n))
(let ((mat
(let ((tmp (make-array (list m n) :element-type 'integer)))
(dotimes (i m tmp)
(dotimes (j n)
(setf (aref tmp i j) (aref mat i j))))))
(l2s (make-array n :element-type 'integer))
(det2s (make-array n :element-type 'integer))
(ort-mat (make-array (list m n) :element-type 'integer))
(coefs (make-array (list m n) :element-type 'integer))
(max-pos -1)
(pos 0))
(declare ((simple-array integer (*)) l2s det2s)
((simple-array integer (* *)) mat ort-mat coefs))
(setf (aref det2s 0) 1)
(labels ((%reduce (target-pos pivot-pos)
(let ((pivot-det2 (aref det2s (+ pivot-pos 1))))
(when (> (abs (* 2 (aref coefs target-pos pivot-pos))) pivot-det2)
(let ((rounded-mu (round (aref coefs target-pos pivot-pos) pivot-det2)))
(dotimes (i m)
(decf (aref mat i target-pos)
(* rounded-mu (aref mat i pivot-pos))))
(decf (aref coefs target-pos pivot-pos)
(* rounded-mu pivot-det2))
(dotimes (j pivot-pos)
(decf (aref coefs target-pos j)
(* rounded-mu (aref coefs pivot-pos j))))))))
(%swap (pos)
(dotimes (i m)
(rotatef (aref mat i (- pos 1)) (aref mat i pos)))
(dotimes (j (- pos 1))
(rotatef (aref coefs (- pos 1) j) (aref coefs pos j)))
(let* ((prec-det2 (aref det2s (- pos 1)))
(det2 (aref det2s pos))
(old-mu (aref coefs pos (- pos 1)))
(new-l2 (%div (+ (* (expt prec-det2 2)
(aref l2s pos))
(* (expt old-mu 2) (aref l2s (- pos 1))))
(expt det2 2)))
(new-det2 (%div (* det2 new-l2) (aref l2s (- pos 1)))))
(setf (aref coefs pos (- pos 1))
(%div (* old-mu (aref l2s (- pos 1))) new-l2))
(setf (aref l2s pos)
(%div (* (expt new-det2 2)
(aref l2s (- pos 1))
(aref l2s pos))
(* (expt det2 2) new-l2)))
(setf (aref l2s (- pos 1)) new-l2)
(dotimes (i m)
(let* ((old-ort (aref ort-mat i (- pos 1)))
(new-ort (%div (+ (* prec-det2 (aref ort-mat i pos))
(* old-mu old-ort))
det2)))
(setf (aref ort-mat i (- pos 1)) new-ort)
(setf (aref ort-mat i pos)
(%div (- (* new-det2 old-ort)
(* (aref coefs pos (- pos 1)) new-ort))
prec-det2))))
(loop with next-det2 = (aref det2s (+ pos 1))
for j from (+ pos 1) to max-pos
for mu = (aref coefs j pos)
do (setf (aref coefs j pos)
(%div (- (* next-det2 (aref coefs j (- pos 1)))
(* old-mu mu))
det2))
(setf (aref coefs j (- pos 1))
(%div (+ (* new-det2 mu)
(* (aref coefs pos (- pos 1))
(aref coefs j pos)))
next-det2)))
(setf (aref det2s pos) new-det2))))
(loop
(when (= pos n)
(return))
(when (> pos max-pos)
;; add new column
(assert (= (+ max-pos 1) pos))
(setq max-pos pos)
(dotimes (i m)
(setf (aref ort-mat i pos) (aref mat i pos)))
(dotimes (j pos)
(let* ((det2 (aref det2s j))
(next-det2 (aref det2s (+ j 1)))
(mu (%div (* (loop for i below m
sum (* (aref mat i pos) (aref ort-mat i j))
of-type rational)
det2
next-det2)
(aref l2s j))))
(setf (aref coefs pos j) mu)
(dotimes (i m)
(decf (aref ort-mat i pos)
(%div (- (* next-det2 (aref ort-mat i pos))
(* mu (aref ort-mat i j)))
det2)))))
(setf (aref l2s pos)
(loop for i below m
sum (expt (aref ort-mat i pos) 2)
of-type rational))
(setf (aref det2s (+ pos 1))
(%div (aref l2s pos) (aref det2s pos))))
(if (zerop pos)
(incf pos)
(loop
(%reduce pos (- pos 1))
(if (< (* (denominator delta)
(+ (* (expt (aref det2s (- pos 1)) 2)
(aref l2s pos))
(* (expt (aref coefs pos (- pos 1)) 2)
(aref l2s (- pos 1)))))
(* (numerator delta)
(expt (aref det2s pos) 2)
(aref l2s (- pos 1))))
(progn
(%swap pos)
(setq pos (max 1 (- pos 1))))
(progn
(loop for pivot-pos from (- pos 2) downto 0
do (%reduce pos pivot-pos))
(incf pos)
(return))))))
mat))))

0 comments on commit cae19e4

Please sign in to comment.