1 (* Copyright (C) 1999-2006 Henry Cejtin, Matthew Fluet, Suresh
2 * Jagannathan, and Stephen Weeks.
4 * MLton is released under a BSD-style license.
5 * See the file MLton-LICENSE for details.
8 functor EuclideanRing(S: EUCLIDEAN_RING_STRUCTS)
9 :> EUCLIDEAN_RING where type t = S.t =
14 structure IntInf = Pervasive.IntInf
18 ("EuclideanRing.divMod",
19 Layout.tuple2(layout, layout),
20 Layout.tuple2(layout, layout),
21 fn (p, q) => (not(equals(q, zero)),
22 fn (d, m) => (equals(p, q * d + m)
23 andalso (equals(m, zero)
24 orelse IntInf.<(metric m, metric q)))))
27 fun p div q = #1(divMod(p, q))
29 fun p mod q = #2(divMod(p, q))
31 fun divides(d: t, x: t): bool = equals(x mod d, zero)
34 Trace.trace("EuclideanRing.divides", Layout.tuple2(layout, layout), Bool.layout) divides
36 (* Taken from page 812 of CLR. *)
37 fun extendedEuclidTerm(a: t, b: t, done: t * t -> bool, trace): t * t * t =
42 else let val (d, m) = divMod(a, b)
43 val (d', x', y') = loop(b, m)
44 in (d', y', x' - d * y')
49 fun makeTraceExtendedEuclid f =
51 ("EuclideanRing.extendedEuclid",
52 Layout.tuple2(layout, layout),
53 Layout.tuple3(layout, layout, layout),
54 fn (a, b) => (not(isZero a) andalso not(isZero b),
55 fn (d, x, y) => (f(d, x, y)
56 andalso equals(d, a * x + b * y))))
60 makeTraceExtendedEuclid
61 (fn (d, x, y) => divides(d, x) andalso divides(d, y))
63 (* Page 72 of Bach and Shallit. *)
64 (* Identical to algorithm on page 23 of Berlekamp. *)
65 (* This algorithm is slower (about 2x) than the recursive extendedEuclid
66 * given above, but stores only a constant number of ring elements.
67 * Thus, for now, it is overridden below.
69 fun extendedEuclid(u0: t, u1: t): t * t * t =
72 fn (r as {m11, m12, m21, m22, u, v, nEven}) =>
73 (Assert.assert("EuclideanRing.extendedEuclid", fn () =>
74 equals(u0, m11 * u + m12 * v)
75 andalso equals(u1, m21 * u + m22 * v)
76 andalso equals(if nEven then one else negOne,
77 m11 * m22 - m12 * m21))
81 let val (q, r) = divMod(u, v)
82 in loop{m11 = q * m11 + m12,
90 val {m12, m22, u, nEven, ...} =
91 loop{m11 = one, m12 = zero, m21 = zero, m22 = one,
92 u = u0, v = u1, nEven = true}
93 val (a, b) = if nEven then (m22, ~m12) else (~m22, m12)
97 val _ = extendedEuclid
99 fun extendedEuclid (a, b) =
100 extendedEuclidTerm (a, b, fn (_, b) => equals (b, zero), trace)
104 val trace = makeTraceExtendedEuclid(fn _ => true)
106 val extendedEuclidTerm =
107 fn (a, b, done) => extendedEuclidTerm(a, b, done, trace)
110 val lastPrime = ref one
112 fun gcd(a, b) = if isZero b then a else gcd(b, a mod b)
114 fun lcm(a, b) = (a * b) div gcd(a, b)
116 val primes: t Stream.t =
118 fun loop(s: t Stream.t) =
121 let val (p, s) = valOf(Stream.force s)
122 val _ = lastPrime := p
124 (p, loop(Stream.keep(s, fn x => not(divides(p, x)))))
133 val layout = Layout.str o toString
136 type factors = (t * Int.t) list
138 fun factor(n: t): factors =
140 fun loop(n: t, primes: t Stream.t, factors: factors) =
143 else let val (p, primes) = valOf(Stream.force primes)
147 let val (q, r) = divMod(n, p)
149 then loop(q, Int.+(k, 1))
157 else (p, k) :: factors)
159 in loop(n, primes, [])
164 ("EuclideanRing.factor", layout, List.layout (Layout.tuple2(layout, Int.layout)),
165 fn n => (not(isZero n), fn factors =>
166 equals(n, List.fold(factors, one, fn ((p, k), prod) =>
167 prod * pow (p, k)))))
170 fun existsPrimeOfSmallerMetric(m: IntInf.int, f: t -> bool): bool =
173 let val (p, primes) = valOf(Stream.force primes)
174 in IntInf.<(metric p, m)
175 andalso (f p orelse loop primes)
180 fun isPrime(r: t): bool =
181 let val r = unitEquivalent r
182 in existsPrimeOfSmallerMetric(IntInf.+ (metric r, 1),
183 fn p => equals(r, p))
186 fun isComposite(r: t): bool =
187 existsPrimeOfSmallerMetric(metric r, fn p => divides(p, r))