Implement math.gcd in dart.

pkg/math/lib/math.dart

Index: pkg/math/lib/math.dart
diff --git a/pkg/math/lib/math.dart b/pkg/math/lib/math.dart
index 9958bec52ba1800eb71959658cd887967d62e4f5..d0cbeaec8f94ac0c84e477535db62fdfb57d2c5f 100644
--- a/pkg/math/lib/math.dart
+++ b/pkg/math/lib/math.dart
@@ -4,4 +4,143 @@
 
 library dart.pkg.math;
 
-// Placeholder library, reserved for future extension.
+/**
+ * Computes the greatest common divisor between [a] and [b].
+ *
+ * The result is always positive even if either `a` or `b` is negative.
+ */
+int gcd(int a, int b) {
+  if (a is! int) throw new ArgumentError(a);

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

Reduce this to just
  if (a == null) ...
the "int" part is checked by the parameters in checked mode and by the analyzer,
and this isn't essential system code, so a wrong argument can't break anything
deep.
Dittos below.

[srawlins, 2014/08/29 06:43:25]

Done.

+  if (b is! int) throw new ArgumentError(b);
+  a = a.abs();
+  b = b.abs();
+
+  // Iterative Binary GCD algorithm.
+  if (a == 0) return b;
+  if (b == 0) return a;
+  int powerOfTwo = 1;
+  int temp;
+  while (((a | b) & 1) == 0) {
+    powerOfTwo *= 2;
+    a ~/= 2;
+    b ~/= 2;
+  }
+
+  while (a.isEven) a ~/= 2;
+
+  do {
+    while (b.isEven) b ~/= 2;
+    if (a > b) {
+      temp = b;

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

declare temp here.

[srawlins, 2014/08/29 06:43:25]

Done.

+      b = a;
+      a = temp;
+    }
+    b -= a;
+  } while (b != 0);
+
+  return a * powerOfTwo;
+}
+
+/**
+ * Computes the greatest common divisor between [a] and [b], as well as [x] and
+ * [y] such that `ax+by == gcd(a,b)`.
+ *
+ * The return value is a List of three ints: the greatest common divisor, `x`,
+ * and `y`, in that order.
+ */
+List<int> gcdext(int a, int b) {
+  if (a is! int) throw new ArgumentError(a);
+  if (b is! int) throw new ArgumentError(b);
+
+  if (a < 0) {
+    List<int> results = gcdext(-a, b);

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

result[1] = -result[1];
return result;

No need to create a new result object.
Similar below.

[srawlins, 2014/08/29 06:43:25]

:P I should have caught that. Thanks.

+    return [results[0], -results[1], results[2]];
+  }
+  if (b < 0) {
+    List<int> results = gcdext(a, -b);
+    return [results[0], results[1], -results[2]];
+  }
+
+  int g, x, y;

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

None of g, x, y are used below?

[srawlins, 2014/08/29 06:43:25]

Done.

+  int r0 = a;
+  int r1 = b;
+  int x0, x1, y0, y1, q, tmp;
+  x0 = y1 = 1;
+  x1 = y0 = 0;
+
+  while (r1 != 0) {

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

Declare q and tmp inside the while.

[srawlins, 2014/08/29 06:43:25]

Done.

+    q = r0 ~/ r1;
+    tmp = r0;
+    r0 = r1;
+    r1 = tmp - q*r1;
+
+    tmp = x0;
+    x0 = x1;
+    x1 = tmp - q*x1;
+
+    tmp = y0;
+    y0 = y1;
+    y1 = tmp - q*y1;
+  }
+
+  return [r0, x0, y0];

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

Maybe use:
  new List<int>(3)..[0] = r0 ..[1] = x0 ..[2] = y0;
That makes it a fixed-length list, which has lower resource overhead.

Maybe even make it a class, like: GCDResults { final int gcd, factor1, factor2;
}, but it's probably not worth the effort.

[srawlins, 2014/08/29 06:43:25]

Done, but spread out onto multiple lines. Should I keep them in one line?

[Lasse Reichstein Nielsen, 2014/08/29 07:14:42]

Multiple lines are fine.

+}
+
+/**
+ * Computes the inverse of [a] modulo [m].
+ *
+ * Throws an IntegerDivisionByZeroException if `a` has no inverse modulo `m`:

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

[IntegerDivisionByZeroException]
(use square brackets the first time you refer to a type in the documentation,
then just use code formatting (backping) if it's mentioned again.

[srawlins, 2014/08/29 06:43:24]

Done.

+ *
+ *     invert(4, 7); // 2
+ *     invert(4, 10); // throws IntegerDivisionByZeroException
+ */
+int invert(int a, int m) {
+  List<int> results = gcdext(a, m);
+  int g = results[0];
+  int x = results[1];
+  if (g != 1) {
+    throw new IntegerDivisionByZeroException();
+  }
+  return x % m;
+}
+
+/**
+ * Computes the least common multiple between [a] and [b].
+ */
+int lcm(int a, int b) {
+  if (a is! int) throw new ArgumentError(a);
+  if (b is! int) throw new ArgumentError(b);
+  if (a == 0 && b == 0) return 0;
+
+  return a.abs() ~/ gcd(a, b) * b.abs();

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

Nice touch to avoid overflows :)

[srawlins, 2014/08/29 06:43:24]

Acknowledged.

+}
+
+/**
+ * Computes [base] raised to [exp] modulo [mod].
+ *
+ * The result has the same sign as `mod`.
+ *
+ * Throws an IntegerDivisionByZeroException if `exp` is negative and `base` has

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

[]-quote the exception type.

[srawlins, 2014/08/29 06:43:25]

Done.

+ * no inverse modulo `mod`.
+ */
+int powmod(int base, int exp, int mod) {

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

If this method is expected to be used often with a mod that is a power of 2, it
might be worth it to have a specialized version that uses bit-and instead of
modulo. Modulo is a slow operation.

[srawlins, 2014/08/29 06:43:25]

I don't think that is a common use case. One of the primary uses of this
function is in RSA encryption, where the modulus is prime, or at least coprime
with the base.

+  if (base is! int) throw new ArgumentError(base);
+  if (exp is! int) throw new ArgumentError(exp);
+  if (mod is! int) throw new ArgumentError(mod);
+
+  // Right-to-left binary method of modular exponentiation.
+  if (exp < 0) {
+    base = invert(base, mod);
+    exp = -exp;
+  }
+  int result = 1;
+  base = base % mod;
+  while (exp != 0) {
+    if (exp.isOdd) {
+      result = (result * base) % mod;
+    }
+    exp ~/= 2;
+    base = (base * base) % mod;

[Lasse Reichstein Nielsen, 2014/08/26 07:44:26]

This is executed one time too many. That's probably ok, if it requires too much
effort to avoid it.

[srawlins, 2014/08/29 06:43:25]

Good call. I'd not noticed that before. What about this new forever loop that
returns if exp == 0? Feels a little awkward...

[Lasse Reichstein Nielsen, 2014/08/29 07:14:42]

It does. And I'd speed-check it too - having the exit condition obvious in the
loop also helps compilers. 
Either is fine, so just pick the most efficient one :)

This is where I'd love to have composite do+while loops:
 
    do {
      if (exp.isOdd) result = (result * base) % mod;
      exp ~/= 2;
    } while (exp != 0) {
      base = (base * base) % mod;
    }

Alas, not in the language (yet!!!) :)

[srawlins, 2014/09/02 19:49:23]

Good call, I think the short circuit here helped more: benchmarking with
powmod(6714601027348841563, 1000000, 13464639003769154407), I found:
while (exp != 0) yielded 12.69 ops/ms
while (true) yielded 13.24 ops/ms

+  }
+  return result;
+}