(ns emmy.euclid
"Implementations of various [greatest common
divisor](https://en.wikipedia.org/wiki/Greatest_common_divisor) algorithms."
(:require [emmy.generic :as g]
[emmy.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 (g/zero? a) [(g/abs b) 0 1]
(g/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 (g/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))))))
#object[emmy.euclid$extended_gcd 0x26095ca6 "
emmy.euclid$extended_gcd@26095ca6"
]
(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 (g/zero? a) (g/abs b)
(g/zero? b) (g/abs a)
(or (v/= a b) (v/= a (g/negate b)))
(g/abs a)
(not (and (v/integral? a) (v/integral? b))) 1
:else (loop [a (g/abs a) b (g/abs b)]
(if (g/zero? b)
a
(recur b (g/remainder a b))))))
#object[emmy.euclid$gcd 0x5d452c7 "
emmy.euclid$gcd@5d452c7"
]

multimethod implementation for basic numeric types.

(defmethod g/gcd :default [a b]
(gcd a b))
#object[clojure.lang.MultiFn 0x724dd6e2 "
clojure.lang.MultiFn@724dd6e2"
]