-
-
Notifications
You must be signed in to change notification settings - Fork 68
/
euclid.cljc
61 lines (55 loc) · 2.3 KB
/
euclid.cljc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
;;
;; Copyright © 2017 Colin Smith.
;; This work is based on the Scmutils system of MIT/GNU Scheme:
;; Copyright © 2002 Massachusetts Institute of Technology
;;
;; This is free software; you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation; either version 3 of the License, or (at
;; your option) any later version.
;;
;; This software is distributed in the hope that it will be useful, but
;; WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
;; General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this code; if not, see <http://www.gnu.org/licenses/>.
;;
(ns sicmutils.euclid
"Implementations of various [greatest common
divisor](https://en.wikipedia.org/wiki/Greatest_common_divisor) algorithms."
(:require [sicmutils.generic :as g]
[sicmutils.value :as v]))
(defn extended-gcd
"Returns a vector containing the [greatest common
divisor](https://en.wikipedia.org/wiki/Greatest_common_divisor) and
the [Bézout coefficients](https://en.wikipedia.org/wiki/Bézout%27s_identity)
corresponding to the inputs `a` and `b`.
For more info, see the Wikipedia article on the [Extended Euclidean
algorithm](http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm)."
[a b]
(cond (v/zero? a) [(g/abs b) 0 1]
(v/zero? b) [(g/abs a) 1 0]
:else (loop [s 0 s0 1 t 1 t0 0 r (g/abs b) r0 (g/abs a)]
(if (v/zero? r)
[r0 s0 t0]
(let [q (g/quotient r0 r)]
(recur (g/- s0 (g/* q s)) s
(g/- t0 (g/* q t)) t
(g/- r0 (g/* q r)) r))))))
(defn gcd
"Returns the [greatest common
divisor](https://en.wikipedia.org/wiki/Greatest_common_divisor) of the two
inputs `a` and `b`."
[a b]
(cond (v/zero? a) (g/abs b)
(v/zero? b) (g/abs a)
(not (and (v/integral? a) (v/integral? b))) 1
:else (loop [a (g/abs a) b (g/abs b)]
(if (v/zero? b)
a
(recur b (g/remainder a b))))))
;; multimethod implementation for basic numeric types.
(defmethod g/gcd :default [a b]
(gcd a b))