24 lines
1.0 KiB
Common Lisp
24 lines
1.0 KiB
Common Lisp
;;
|
|
;; Calculates the GCD of a and b based on the Extended Euclidean Algorithm. The function also returns
|
|
;; the Bézout coefficients s and t, such that gcd(a, b) = as + bt.
|
|
;;
|
|
;; The algorithm is described on page http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Iterative_method_2
|
|
;;
|
|
(defun egcd (a b)
|
|
(do ((r (cons b a) (cons (- (cdr r) (* (car r) q)) (car r))) ; (r+1 r) i.e. the latest is first.
|
|
(s (cons 0 1) (cons (- (cdr s) (* (car s) q)) (car s))) ; (s+1 s)
|
|
(u (cons 1 0) (cons (- (cdr u) (* (car u) q)) (car u))) ; (t+1 t)
|
|
(q nil))
|
|
((zerop (car r)) (values (cdr r) (cdr s) (cdr u))) ; exit when r+1 = 0 and return r s t
|
|
(setq q (floor (/ (cdr r) (car r)))))) ; inside loop; calculate the q
|
|
|
|
;;
|
|
;; Calculates the inverse module for a = 1 (mod m).
|
|
;;
|
|
;; Note: The inverse is only defined when a and m are coprimes, i.e. gcd(a, m) = 1.”
|
|
;;
|
|
(defun invmod (a m)
|
|
(multiple-value-bind (r s k) (egcd a m)
|
|
(unless (= 1 r) (error "invmod: Values ~a and ~a are not coprimes." a m))
|
|
s))
|