Merge pull request #633 from non/topic/laws-field-gcdlcm

Adding GCDRing laws, fixing implementations

GUIDE.md

@@ -330,10 +330,26 @@ inheritance:
 
 Rings also provide a `pow` method (`**`) for doing repeated multiplication.
 
+#### GCDRings
+
+GCDRings are commutative rings (`CRing[A]`) with existence of a
+greatest-common-divisor and least-common-multiple.
+
+Spire's `GCDRing[A]` supports the following operations:
+
+ * `gcd` (`a gcd b`) find the greatest common divisor of `a` and `b`.
+ * `lcm` (`a lcm b`) find the lowest common multiple of `a` and `b`.
+
+Spire requires these operations to be commutative. Note that fields
+have leeway to define the GCD operation. In practice, instances of
+`Field[A]` provide either a trivial implementation `gcd(x != 0 , y != 0) == 1`
+or a definition that extends the one used for the integer ring
+(`gcd(a/b, c/d) == gcd(a, c)/lcm(b, d)`).
+
 #### EuclideanRings
 
 Spire supports euclidean domains (called `EuclideanRing[A]`). A
-euclidean domain is a commutative ring (`CRing[A]`) that also supports
+euclidean domain is a GCD ring (`GCDRing[A]`) that also supports
 euclidean division (e.g. floor division or integer division). This
 structure generalizes many useful properties of the integers (for
 instance, quotients and remainders, and greatest common divisors).
@@ -349,8 +365,6 @@ Spire's `EuclideanRing[A]` supports the following operations:
  * `quot` (`a /~ b`) finding the quotient (often integer division).
  * `mod` (`a % b`) the remainder from the quotient operation.
  * `quotmod` (`a /% b`) combines `quot` and `mod` into one operation.
- * `gcd` (`a gcd b`) find the greatest common divisor of `a` and `b`.
- * `lcm` (`a lcm b`) find the lowest common multiple of `a` and `b`.
 
 Spire requires that `b * (a /~ b) + (a % b)` is equivalent to `a`.
 

README.md

@@ -249,12 +249,13 @@ most users will expect.
    + zero: additive identity
    + one: multiplicative identity
  * *Ring* (Rng + Rig)
+ * *GCDRing*
+   + gcd: greatest-common-divisor
+   + lcm: least-common-multiple
  * *EuclideanRing*
-   + quot (`/~`): quotient (floor division)
+   + quot (`/~`): quotient (Euclidean division)
    + mod (`%`): remainder
    + quotmod (`/%`): quotient and mod
-   + gcd: greatest-common-divisor
-   + lcm: least-common-multiple
  * *Field*
    + reciprocal: multiplicative inverse
    + div (`/`): division

core/shared/src/main/scala/spire/algebra/EuclideanRing.scala

@@ -22,10 +22,6 @@ trait EuclideanRing[@sp(Int, Long, Float, Double) A] extends Any with GCDRing[A]
   def quot(a: A, b: A): A
   def mod(a: A, b: A): A
   def quotmod(a: A, b: A): (A, A) = (quot(a, b), mod(a, b))
-  def gcd(a: A, b: A)(implicit ev: Eq[A]): A =
-    EuclideanRing.euclid(a, b)(ev, EuclideanRing.this)
-  def lcm(a: A, b: A)(implicit ev: Eq[A]): A =
-    if (isZero(a) || isZero(b)) zero else times(quot(a, gcd(a, b)), b)
 }
 
 trait EuclideanRingFunctions[R[T] <: EuclideanRing[T]] extends GCDRingFunctions[R] {
@@ -43,4 +39,11 @@ object EuclideanRing extends EuclideanRingFunctions[EuclideanRing] {
     if (EuclideanRing[A].isZero(b)) a else euclid(b, EuclideanRing[A].mod(a, b))
   }
 
+  trait WithEuclideanAlgorithm[@sp(Int, Long, Float, Double) A] extends Any with EuclideanRing[A] { self =>
+    def gcd(a: A, b: A)(implicit ev: Eq[A]): A =
+      EuclideanRing.euclid(a, b)(ev, self)
+    def lcm(a: A, b: A)(implicit ev: Eq[A]): A =
+      if (isZero(a) || isZero(b)) zero else times(quot(a, gcd(a, b)), b)
+  }
+
 }

core/shared/src/main/scala/spire/algebra/Field.scala

@@ -1,21 +1,62 @@
 package spire
 package algebra
 
+/** Field type class. While algebra already provides one, we provide one in Spire
+  * that integrates with the commutative ring hierarchy, in particular `GCDRing`
+  * and `EuclideanRing`.
+  * 
+  * On a field, all nonzero elements are invertible, so the remainder of the
+  * division is always 0. The Euclidean function can take an arbitrary value on
+  * nonzero elements (it is undefined for zero); for compatibility with the degree
+  * of polynomials, we use the constant 0.
+  *
+  * The GCD and LCM are defined up to a unit; on a field, it means that either the GCD or LCM
+  * can be fixed arbitrarily. Some conventions with consistent defaults are provided in the 
+  * spire.algebra.Field companion object.
+  */
 trait Field[@sp(Int, Long, Float, Double) A] extends Any with AlgebraField[A] with EuclideanRing[A] {
-  /* On a field, all nonzero elements are invertible, so the remainder of the
-   division is always 0. The Euclidean function can take an arbitrary value on
-   nonzero elements (it is undefined for zero); for compatibility with the degree
-   of polynomials, we use the constant 0.
-   */
   def euclideanFunction(a: A): BigInt = BigInt(0)
   def quot(a: A, b: A): A = div(a, b)
   def mod(a: A, b: A): A = zero
   override def quotmod(a: A, b: A): (A, A) = (div(a, b), zero)
-  override def gcd(a: A, b: A)(implicit eqA: Eq[A]): A =
-    if (isZero(a) && isZero(b)) zero else one
-  override def lcm(a: A, b: A)(implicit eqA: Eq[A]): A = times(a, b)
 }
 
+
 object Field extends _root_.algebra.ring.FieldFunctions[Field] with EuclideanRingFunctions[Field] {
+
   @inline def apply[A](implicit ev: Field[A]): Field[A] = ev
+
+  /** Field with simple default GCD/LCM implementations:
+    * gcd(a, b) = 1 (except gcd(0, 0) = 0) while lcm(a, b) = a * b. */
+  trait WithDefaultGCD[@sp(Int, Long, Float, Double) A] extends Any with Field[A] {
+    override def gcd(a: A, b: A)(implicit eqA: Eq[A]): A =
+      if (isZero(a) && isZero(b)) zero else one
+    override def lcm(a: A, b: A)(implicit eqA: Eq[A]): A = times(a, b)
+  }
+
+  /** Field defined as a field of fractions with a default implementation of GCD/LCM such that
+    * - gcd(a/b, c/d) = gcd(a, c) / lcm(b, d)
+    * - lcm(a/b, c/d) = lcm(a, c) / gcd(b, d)
+    * which corresponds to the convention of the GCD domains of SageMath; on rational numbers, it
+    * "yields the unique extension of gcd from integers to rationals presuming the natural extension
+    * of the divisibility relation from integers to rationals", see http://math.stackexchange.com/a/151431
+    */
+  trait FieldOfFractionsGCD[A, R] extends Any with Field[A] {
+    implicit def ringR: GCDRing[R]
+    implicit def eqR: Eq[R]
+    def numerator(a: A): R
+    def denominator(a: A): R
+    def fraction(num: R, den: R): A
+    override def gcd(x: A, y: A)(implicit ev: Eq[A]): A = {
+      val num = ringR.gcd(numerator(x), numerator(y))
+      val den = ringR.lcm(denominator(x), denominator(y))
+      fraction(num, den)
+    }
+    override def lcm(x: A, y: A)(implicit ev: Eq[A]): A = {
+      val num = ringR.lcm(numerator(x), numerator(y))
+      val den = ringR.gcd(denominator(x), denominator(y))
+      fraction(num, den)
+    }
+  }
+
 }

core/shared/src/main/scala/spire/macros/Auto.scala

@@ -152,7 +152,15 @@ abstract class AutoAlgebra extends AutoOps { ops =>
   def EuclideanRing[A: c.WeakTypeTag](z: c.Expr[A], o: c.Expr[A])
       (ev: c.Expr[Eq[A]]): c.Expr[EuclideanRing[A]] = {
     c.universe.reify {
-      new EuclideanRing[A] {
+      new EuclideanRing[A] { self =>
+        // default implementations from EuclideanRing.WithEuclideanAlgorithm
+        @tailrec final def euclid(a: A, b: A)(implicit ev: Eq[A]): A =
+          if (isZero(b)) a else euclid(b, mod(a, b))
+        def gcd(a: A, b: A)(implicit ev: Eq[A]): A =
+          euclid(a, b)(ev)
+        def lcm(a: A, b: A)(implicit ev: Eq[A]): A =
+          if (isZero(a) || isZero(b)) zero else times(quot(a, gcd(a, b)), b)
+
         def zero: A = z.splice
         def one: A = o.splice
         def plus(x: A, y: A): A = ops.plus[A].splice
@@ -170,6 +178,11 @@ abstract class AutoAlgebra extends AutoOps { ops =>
       (z: c.Expr[A], o: c.Expr[A])(ev: c.Expr[Eq[A]]): c.Expr[Field[A]] = {
     c.universe.reify {
       new Field[A] {
+        // default implementations from Field.WithDefaultGCD
+        override def gcd(a: A, b: A)(implicit eqA: Eq[A]): A =
+              if (isZero(a) && isZero(b)) zero else one
+        override def lcm(a: A, b: A)(implicit eqA: Eq[A]): A = times(a, b)
+
         def zero: A = z.splice
         def one: A = o.splice
         def plus(x: A, y: A): A = ops.plus[A].splice

core/shared/src/main/scala/spire/math/Algebraic.scala

@@ -1527,7 +1527,7 @@ trait AlgebraicInstances {
   implicit final val AlgebraicTag = new LargeTag[Algebraic](Exact, Algebraic(0))
 }
 
-private[math] trait AlgebraicIsField extends Field[Algebraic] {
+private[math] trait AlgebraicIsField extends Field.WithDefaultGCD[Algebraic] {
   def zero: Algebraic = Algebraic.Zero
   def one: Algebraic = Algebraic.One
   def plus(a: Algebraic, b: Algebraic): Algebraic = a + b

core/shared/src/main/scala/spire/math/Complex.scala

@@ -600,7 +600,7 @@ private[math] trait ComplexIsRing[@sp(Float, Double) A] extends Ring[Complex[A]]
   override def fromInt(n: Int): Complex[A] = Complex.fromInt[A](n)
 }
 
-private[math] trait ComplexIsField[@sp(Float,Double) A] extends ComplexIsRing[A] with Field[Complex[A]] {
+private[math] trait ComplexIsField[@sp(Float,Double) A] extends ComplexIsRing[A] with Field.WithDefaultGCD[Complex[A]] {
 
   implicit def algebra: Field[A]
 

core/shared/src/main/scala/spire/math/Jet.scala

@@ -553,7 +553,7 @@ private[math] trait JetIsEuclideanRing[@sp(Float,Double) T]
 
 /* TODO: Jet[T] is probably not a genuine Field */
 private[math] trait JetIsField[@sp(Float,Double) T]
-    extends JetIsEuclideanRing[T] with Field[Jet[T]] {
+    extends JetIsEuclideanRing[T] with Field.WithDefaultGCD[Jet[T]] {
   /* TODO: what are exactly the laws of Jet with respect to EuclideanRing ? */
   // duplicating methods because super[..].call() does not work on 2.10 and 2.11
   override def fromDouble(n: Double): Jet[T] = Jet(f.fromDouble(n))

core/shared/src/main/scala/spire/math/Number.scala

@@ -647,7 +647,7 @@ private[math] trait NumberIsEuclideanRing extends EuclideanRing[Number] with Num
 }
  */
 
-private[math] trait NumberIsField extends Field[Number] with NumberIsCRing {
+private[math] trait NumberIsField extends Field.WithDefaultGCD[Number] with NumberIsCRing {
   def div(a:Number, b:Number): Number = a / b
   override def fromDouble(a: Double): Number = Number(a)
 }

core/shared/src/main/scala/spire/math/Polynomial.scala

@@ -518,7 +518,7 @@ with Ring[Polynomial[C]] {
 }
 
 trait PolynomialOverField[@sp(Double) C] extends PolynomialOverRing[C]
-    with EuclideanRing[Polynomial[C]] with VectorSpace[Polynomial[C], C] { self =>
+    with EuclideanRing.WithEuclideanAlgorithm[Polynomial[C]] with VectorSpace[Polynomial[C], C] { self =>
   implicit override val scalar: Field[C]
 
   override def divr(x: Polynomial[C], k: C): Polynomial[C] = x :/ k
@@ -539,18 +539,6 @@ trait PolynomialOverField[@sp(Double) C] extends PolynomialOverRing[C]
     }
   }
 
-  // TODO: why final?
-  final override def gcd(x: Polynomial[C], y: Polynomial[C])(implicit ev: Eq[Polynomial[C]]): Polynomial[C] = {
-    val result = EuclideanRing.euclid(x, y)(ev, self)
-    if (result.degree > 0) {
-      result
-    } else {
-      // TODO: I don't like so much this special case, can we fold it in the general result,
-      // TODO: as the gcd is defined up to a unit anyway?
-      // return the gcd of all coefficients when there is no higher degree divisor
-      Polynomial.constant(spire.math.gcd(x.coeffsArray ++ y.coeffsArray))
-    }
-  }
 }
 
 trait PolynomialEq[@sp(Double) C] extends Eq[Polynomial[C]] {