(ns emmy.modint
"This namespace contains an implementation of a [[ModInt]] datatype and various
operations for creating and working with [[ModInt]] instances. See [\"Modular
Arithmetic\"](https://en.wikipedia.org/wiki/Modular_arithmetic) on Wikipedia
for more details about modular arithmetic.
[[emmy.modint]] also extends many Emmy generic operations
to the [[ModInt]] datatype."
(:require #?(:cljs [emmy.euclid :as e])
[emmy.generic :as g]
[emmy.util :as u]
[emmy.value :as v]))
(declare mod:=)
#object[emmy.modint$mod_COLON__EQ_ 0x55247b23 "
emmy.modint$mod_COLON__EQ_@55247b23"
]
(deftype ModInt [i m]
v/IKind
(kind [_] ::modint)
#?@(:clj
[Object
(equals [this that] (mod:= this that))
(toString [_] (str "[" i " mod " m "]"))]
:cljs
[IEquiv
(-equiv [this that] (mod:= this that))
Object
(toString [_] (str "[" i " mod " m "]"))
IPrintWithWriter
(-pr-writer [x writer _]
(write-all writer
"#object[emmy.modint.ModInt \""
(.toString x)
"\"]"))]))
#object[emmy.modint$eval86784$__GT_ModInt__86786 0x2801a738 "
emmy.modint$eval86784$__GT_ModInt__86786@2801a738"
]
(defn modint?
"Returns true if `x` is an instance of [[ModInt]], false otherwise."
[x]
(instance? ModInt x))
#object[emmy.modint$modint_QMARK_ 0x785cd93f "
emmy.modint$modint_QMARK_@785cd93f"
]
(defn residue [x]
(.-i ^ModInt x))
#object[emmy.modint$residue 0x305d1af8 "
emmy.modint$residue@305d1af8"
]
(defn modulus [x]
(.-m ^ModInt x))
#object[emmy.modint$modulus 0x465682b3 "
emmy.modint$modulus@465682b3"
]
(defn- mod:= [this that]
(cond (modint? that)
(and (= (modulus this)
(modulus that))
(v/= (residue this)
(residue that)))
(v/number? that)
(v/= (residue this)
(g/modulo that (modulus this)))
:else false))
#object[emmy.modint$mod_COLON__EQ_ 0x55247b23 "
emmy.modint$mod_COLON__EQ_@55247b23"
]
(defn make
"Returns an instance of [[ModInt]] that represents integer `i` with integral
modulus `m`."
[i m]
{:pre [(v/integral? i)
(v/integral? m)]}
(->ModInt (g/modulo i m) m))
#object[emmy.modint$make 0x63cbc526 "
emmy.modint$make@63cbc526"
]
(defn- modular-binop [op]
(fn [a b]
(if-not (= (modulus a) (modulus b))
(u/arithmetic-ex "unequal moduli")
(make (op (residue a) (residue b)) (modulus a)))))
#object[emmy.modint$modular_binop 0x7d407d41 "
emmy.modint$modular_binop@7d407d41"
]
(defn- invert
"Modular inverse. JVM implementation uses the native BigInt implementation."
([m] (invert (residue m) (modulus m)))
([i modulus]
#?(:clj
(try (-> (biginteger i)
(.modInverse (biginteger modulus))
(int)
(->ModInt modulus))
(catch ArithmeticException _
(u/arithmetic-ex (str i " is not invertible mod " modulus))))
:cljs
(let [[g a _] (e/extended-gcd i modulus)]
(if (< g 2)
(make a modulus)
(u/arithmetic-ex (str i " is not invertible mod " modulus)))))))
#object[emmy.modint$invert 0x251222cb "
emmy.modint$invert@251222cb"
]
(defn- mod-expt
"Modular exponentiation, more efficient on the JVM."
[base pow modulus]
#?(:clj (let [base (if (neg? pow)
(residue (invert base modulus))
base)]
(-> (.modPow (biginteger base)
(.abs (biginteger pow))
(biginteger modulus))
(int)
(->ModInt modulus)))
:cljs (-> (g/expt (u/bigint base)
(u/bigint pow))
(g/modulo modulus)
(js/Number)
(->ModInt modulus))))
#object[emmy.modint$mod_expt 0x57236fc0 "
emmy.modint$mod_expt@57236fc0"
]
(defn chinese-remainder
"[Chinese Remainder Algorithm](https://en.wikipedia.org/wiki/Chinese_remainder_theorem).
Accepts a sequence of [[ModInt]] instances (where the `modulus` of
all [[ModInt]] instances are relatively prime), and returns a [[ModInt]] `x`
such that `(residue input) == (mod x (modulus input))`.
For example:
```clojure
(let [a1 (m/make 2 5)
a2 (m/make 3 13)]
[(= 42 (chinese-remainder a1 a2))
(= (residue a1) (mod cr (modulus a1)))
(= (residue a2) (mod cr (modulus a2)))])
;;=> [true true true]
```"
[& modints]
(let [prod (transduce (map modulus) g/* modints)
xform (map (fn [mi]
(let [i (residue mi)
m (modulus mi)
c (g/quotient prod m)]
(g/* i c (residue (invert c m))))))]
(-> (transduce xform g/+ modints)
(g/modulo prod))))
#object[emmy.modint$chinese_remainder 0x6afc62f7 "
emmy.modint$chinese_remainder@6afc62f7"
]
(def ^:private add (modular-binop g/add))
#object[emmy.modint$modular_binop$fn__86795 0x1428a09c "
emmy.modint$modular_binop$fn__86795@1428a09c"
]
(def ^:private sub (modular-binop g/sub))
#object[emmy.modint$modular_binop$fn__86795 0x2018cca8 "
emmy.modint$modular_binop$fn__86795@2018cca8"
]
(def ^:private mul (modular-binop g/mul))
#object[emmy.modint$modular_binop$fn__86795 0x6c89ea04 "
emmy.modint$modular_binop$fn__86795@6c89ea04"
]
(def ^:private remainder (modular-binop g/remainder))
#object[emmy.modint$modular_binop$fn__86795 0x656bd6b8 "
emmy.modint$modular_binop$fn__86795@656bd6b8"
]
(def ^:private modulo (modular-binop g/modulo))
#object[emmy.modint$modular_binop$fn__86795 0x3496b6c9 "
emmy.modint$modular_binop$fn__86795@3496b6c9"
]
(defn- div [a b]
(mul a (invert b)))
#object[emmy.modint$div 0x528d517f "
emmy.modint$div@528d517f"
]
(defmethod v/= [::v/number ::modint] [l r] (mod:= r l))
#object[clojure.lang.MultiFn 0x744f207b "
clojure.lang.MultiFn@744f207b"
]
(defmethod v/= [::modint ::v/number] [l r] (mod:= l r))
#object[clojure.lang.MultiFn 0x744f207b "
clojure.lang.MultiFn@744f207b"
]
(defmethod g/zero? [::modint] [^ModInt a] (g/zero? (.-i a)))
#object[clojure.lang.MultiFn 0x78e2923f "
clojure.lang.MultiFn@78e2923f"
]
(defmethod g/one? [::modint] [^ModInt a] (g/one? (.-i a)))
#object[clojure.lang.MultiFn 0x710fa3bd "
clojure.lang.MultiFn@710fa3bd"
]
(defmethod g/identity? [::modint] [^ModInt a] (g/one? (.-i a)))
#object[clojure.lang.MultiFn 0x3e853a24 "
clojure.lang.MultiFn@3e853a24"
]
(defmethod g/zero-like [::modint] [^ModInt a] (ModInt. (g/zero-like (.-i a)) (.-m a)))
#object[clojure.lang.MultiFn 0x79337717 "
clojure.lang.MultiFn@79337717"
]
(defmethod g/one-like [::modint] [^ModInt a] (ModInt. (g/one-like (.-i a)) (.-m a)))
#object[clojure.lang.MultiFn 0x46a24adf "
clojure.lang.MultiFn@46a24adf"
]
(defmethod g/identity-like [::modint] [^ModInt a] (ModInt. (g/one-like (.-i a)) (.-m a)))
#object[clojure.lang.MultiFn 0x787144b0 "
clojure.lang.MultiFn@787144b0"
]
(defmethod g/freeze [::modint] [^ModInt a] (list 'modint (.-i a) (.-m a)))
#object[clojure.lang.MultiFn 0x7a0db79a "
clojure.lang.MultiFn@7a0db79a"
]
(defmethod g/exact? [::modint] [_] true)
#object[clojure.lang.MultiFn 0x1864e04e "
clojure.lang.MultiFn@1864e04e"
]
(defmethod g/integer-part [::modint] [^ModInt a] (.-i a))
#object[clojure.lang.MultiFn 0x12b60d1 "
clojure.lang.MultiFn@12b60d1"
]
(defmethod g/fractional-part [::modint] [_] 0)
#object[clojure.lang.MultiFn 0x4430be99 "
clojure.lang.MultiFn@4430be99"
]
(defmethod g/floor [::modint] [a] a)
#object[clojure.lang.MultiFn 0x6b6c4b9d "
clojure.lang.MultiFn@6b6c4b9d"
]
(defmethod g/ceiling [::modint] [a] a)
#object[clojure.lang.MultiFn 0x775a6a7c "
clojure.lang.MultiFn@775a6a7c"
]
(defmethod g/add [::modint ::modint] [a b] (add a b))
#object[clojure.lang.MultiFn 0x7e1fb1be "
clojure.lang.MultiFn@7e1fb1be"
]
(defmethod g/mul [::modint ::modint] [a b] (mul a b))
#object[clojure.lang.MultiFn 0x4485ec13 "
clojure.lang.MultiFn@4485ec13"
]
(defmethod g/div [::modint ::modint] [a b] (div a b))
#object[clojure.lang.MultiFn 0x233dab3a "
clojure.lang.MultiFn@233dab3a"
]
(defmethod g/sub [::modint ::modint] [a b] (sub a b))
#object[clojure.lang.MultiFn 0x3b212ce4 "
clojure.lang.MultiFn@3b212ce4"
]
(defmethod g/negate [::modint] [a] (make (g/negate (residue a)) (modulus a)))
#object[clojure.lang.MultiFn 0x1b9c4031 "
clojure.lang.MultiFn@1b9c4031"
]
(defmethod g/invert [::modint] [a] (invert a))
#object[clojure.lang.MultiFn 0x165df7c1 "
clojure.lang.MultiFn@165df7c1"
]
(defmethod g/magnitude [::modint] [a]
(g/modulo (residue a)
(modulus a)))
#object[clojure.lang.MultiFn 0x2df77656 "
clojure.lang.MultiFn@2df77656"
]
(defmethod g/abs [::modint] [a]
(let [i (residue a)]
(if (g/negative? i)
(make i (modulus a))
a)))
#object[clojure.lang.MultiFn 0x721e14f0 "
clojure.lang.MultiFn@721e14f0"
]
(defmethod g/quotient [::modint ::modint] [a b] (mul a (invert b)))
#object[clojure.lang.MultiFn 0x7a1d117a "
clojure.lang.MultiFn@7a1d117a"
]
(defmethod g/remainder [::modint ::modint] [a b] (remainder a b))
#object[clojure.lang.MultiFn 0x1d8867c7 "
clojure.lang.MultiFn@1d8867c7"
]
(defmethod g/modulo [::modint ::modint] [a b] (modulo a b))
#object[clojure.lang.MultiFn 0xf575ff8 "
clojure.lang.MultiFn@f575ff8"
]
(defmethod g/exact-divide [::modint ::modint] [a b] (mul a (invert b)))
#object[clojure.lang.MultiFn 0x198b6b28 "
clojure.lang.MultiFn@198b6b28"
]
(defmethod g/negative? [::modint] [a] (g/negative? (residue a)))
#object[clojure.lang.MultiFn 0x6e310c69 "
clojure.lang.MultiFn@6e310c69"
]

A more efficient exponent implementation is available on the JVM.

(defmethod g/expt [::v/integral ::modint] [a b] (mod-expt a (residue b) (modulus b)))
#object[clojure.lang.MultiFn 0x310104f9 "
clojure.lang.MultiFn@310104f9"
]
(defmethod g/expt [::modint ::v/integral] [a b] (mod-expt (residue a) b (modulus a)))
#object[clojure.lang.MultiFn 0x310104f9 "
clojure.lang.MultiFn@310104f9"
]
(defmethod g/solve-linear [::modint ::modint] [a b] (div b a))
#object[clojure.lang.MultiFn 0x48e9c356 "
clojure.lang.MultiFn@48e9c356"
]
(defmethod g/solve-linear-right [::modint ::modint] [a b] (div a b))
#object[clojure.lang.MultiFn 0x738d50cf "
clojure.lang.MultiFn@738d50cf"
]

Methods that allow interaction with other integral types. The first block is perhaps slightly more efficient:

(doseq [op [g/add g/mul g/sub]]
(defmethod op [::v/integral ::modint] [a b] (make (op a (residue b)) (modulus b)))
(defmethod op [::modint ::v/integral] [a b] (make (op (residue a) b) (modulus a))))
nil

The second block promotes any integral type to a ModInt before operating.

(doseq [op [g/div g/solve-linear g/solve-linear-right
g/quotient g/remainder g/exact-divide]]
(defmethod op [::v/integral ::modint] [a b] (op (make a (modulus b)) b))
(defmethod op [::modint ::v/integral] [a b] (op a (make b (modulus a)))))
nil