diff --git a/src/main/scala/leon/synthesis/GCD.scala b/src/main/scala/leon/synthesis/Algebra.scala
similarity index 63%
rename from src/main/scala/leon/synthesis/GCD.scala
rename to src/main/scala/leon/synthesis/Algebra.scala
index 47057dcedf7354b4f6074febde22fce84e288cbe..c92b08980c05a850926fa64878df353c8a07909e 100644
--- a/src/main/scala/leon/synthesis/GCD.scala
+++ b/src/main/scala/leon/synthesis/Algebra.scala
@@ -1,6 +1,14 @@
 package leon.synthesis
 
-object GCD {
+/*
+ * This provides some algebra/number theory functions, including operation such as true division,
+ * GCD and LCM as well as some matrix computation.
+ *
+ * Notice that all those functionalities are independent of the Leon language and
+ * are working for Integer (by opposition to real numbers).
+ */
+
+object Algebra {
 
   def remainder(x: Int, y: Int) = ((x % y) + y.abs) % y.abs
 
@@ -92,4 +100,32 @@ object GCD {
     else sys.error("shouldn't have forgot a case here")
   }
 
+
+  //val that the sol vector with the term in the equation
+  def eval(sol: Array[Int], equation: Array[Int]): Int = {
+    require(sol.size == equation.size)
+    sol.zip(equation).foldLeft(0)((acc, p) => acc + p._1 * p._2)
+  }
+
+  //multiply the matrix by the vector: [M1 M2 .. Mn] * [v1 .. vn] = v1*M1 + ... + vn*Mn]
+  def mult(matrix: Array[Array[Int]], vector: Array[Int]): Array[Int] = {
+    require(vector.size == matrix(0).size && vector.size > 0)
+    val tmat = matrix.transpose
+    var tmp: Array[Int] = null
+    tmp = mult(vector(0), tmat(0))
+    var i = 1
+    while(i < vector.size) {
+      tmp = add(tmp, mult(vector(i), tmat(i)))
+      i += 1
+    }
+    tmp
+  }
+
+  def mult(c: Int, v: Array[Int]): Array[Int] = v.map(_ * c)
+
+  def add(v1: Array[Int], v2: Array[Int]): Array[Int] = {
+    require(v1.size == v2.size)
+    v1.zip(v2).map(p => p._1 + p._2)
+  }
+
 }
diff --git a/src/main/scala/leon/synthesis/ArithmeticNormalization.scala b/src/main/scala/leon/synthesis/ArithmeticNormalization.scala
index f81311ffdc16388e64d5e7eb2468646f234df728..6c148e2ccfd5066d6d99a9035ad392aecbb475ef 100644
--- a/src/main/scala/leon/synthesis/ArithmeticNormalization.scala
+++ b/src/main/scala/leon/synthesis/ArithmeticNormalization.scala
@@ -4,10 +4,6 @@ import leon.purescala.Trees._
 import leon.purescala.TreeOps._
 import leon.purescala.Common._
 
-/*
- * TODO: move those functions to TreeOps
- */
-
 object ArithmeticNormalization {
 
   case class NonLinearExpressionException(msg: String) extends Exception
diff --git a/src/main/scala/leon/synthesis/LinearEquations.scala b/src/main/scala/leon/synthesis/LinearEquations.scala
index 3d9eb6230bd9ad18c1b41da40de1aa07d754c354..63e234320a3522d3e277947ee173b193580a09f1 100644
--- a/src/main/scala/leon/synthesis/LinearEquations.scala
+++ b/src/main/scala/leon/synthesis/LinearEquations.scala
@@ -4,6 +4,7 @@ import leon.purescala.Trees._
 import leon.purescala.TypeTrees._
 import leon.purescala.Common._
 import leon.Evaluator 
+import leon.synthesis.Algebra._
 
 object LinearEquations {
 
@@ -17,7 +18,7 @@ object LinearEquations {
     val orderedParams: Array[Identifier] = as.toArray
     val coefsParams0: List[Int] = ArithmeticNormalization(t, orderedParams).map{case IntLiteral(i) => i}.toList
     val coefsParams: List[Int] = if(coefsParams0.head == 0) coefsParams0.tail else coefsParams0
-    val d: Int = GCD.gcd((coefsParams ++ coefsVars).toSeq)
+    val d: Int = gcd((coefsParams ++ coefsVars).toSeq)
 
     if(coefsVars.size == 1) {
       val coef = coefsVars.head
@@ -45,33 +46,6 @@ object LinearEquations {
 
   }
 
-  //val that the sol vector with the term in the equation
-  def eval(sol: Array[Int], equation: Array[Int]): Int = {
-    require(sol.size == equation.size)
-    sol.zip(equation).foldLeft(0)((acc, p) => acc + p._1 * p._2)
-  }
-
-  //multiply the matrix by the vector: [M1 M2 .. Mn] * [v1 .. vn] = v1*M1 + ... + vn*Mn]
-  def mult(matrix: Array[Array[Int]], vector: Array[Int]): Array[Int] = {
-    require(vector.size == matrix(0).size && vector.size > 0)
-    val tmat = matrix.transpose
-    var tmp: Array[Int] = null
-    tmp = mult(vector(0), tmat(0))
-    var i = 1
-    while(i < vector.size) {
-      tmp = add(tmp, mult(vector(i), tmat(i)))
-      i += 1
-    }
-    tmp
-  }
-
-  def mult(c: Int, v: Array[Int]): Array[Int] = v.map(_ * c)
-
-  def add(v1: Array[Int], v2: Array[Int]): Array[Int] = {
-    require(v1.size == v2.size)
-    v1.zip(v2).map(p => p._1 + p._2)
-  }
-
   //compute a list of solution of the equation c1*x1 + ... + cn*xn where coef = [c1 ... cn]
   //return the solution in the form of a list of n-1 vectors that form a basis for the set
   //of solutions, that is res=(v1, ..., v{n-1}) and any solution x* to the original solution
@@ -87,7 +61,7 @@ object LinearEquations {
         if(i < j)
           K(i)(j) = 0
         else if(i == j) {
-          K(j)(j) = GCD.gcd(coef.drop(j+1))/GCD.gcd(coef.drop(j))
+          K(j)(j) = gcd(coef.drop(j+1))/gcd(coef.drop(j))
         }
       }
     }
@@ -117,8 +91,8 @@ object LinearEquations {
   //return a particular solution to t + c1x + c2y = 0, with (pre, (x0, y0))
   def particularSolution(as: Set[Identifier], t: Expr, c1: Expr, c2: Expr): (Expr, (Expr, Expr)) = {
     val (IntLiteral(i1), IntLiteral(i2)) = (c1, c2)
-    val (v1, v2) = GCD.extendedEuclid(i1, i2)
-    val d = GCD.gcd(i1, i2)
+    val (v1, v2) = extendedEuclid(i1, i2)
+    val d = gcd(i1, i2)
 
     val pre = Equals(Modulo(t, IntLiteral(d)), IntLiteral(0))
 
@@ -135,7 +109,7 @@ object LinearEquations {
     require(normalizedEquation.size >= 2)
     val t: Expr = normalizedEquation.head
     val coefs: List[Int] = normalizedEquation.tail.map{case IntLiteral(i) => i}
-    val d = GCD.gcd(coefs.toSeq)
+    val d = gcd(coefs.toSeq)
     val pre = Equals(Modulo(t, IntLiteral(d)), IntLiteral(0))
 
     if(normalizedEquation.size == 2) {
@@ -146,7 +120,7 @@ object LinearEquations {
     } else {
       val gamma1: Expr = normalizedEquation(1)
       val coefs: List[Int] = normalizedEquation.drop(2).map{case IntLiteral(i) => i}
-      val gamma2: Expr = IntLiteral(GCD.gcd(coefs.toSeq))
+      val gamma2: Expr = IntLiteral(gcd(coefs.toSeq))
       val (_, (w1, w)) = particularSolution(as, t, gamma1, gamma2)
       val (_, sols) = particularSolution(as, UMinus(Times(w, gamma2)) :: normalizedEquation.drop(2))
       (pre, w1 :: sols)
diff --git a/src/main/scala/leon/synthesis/rules/IntegerInequalities.scala b/src/main/scala/leon/synthesis/rules/IntegerInequalities.scala
index a4ff5fc7a8ddf1675b5d5bc9471a84ca48fae0de..d3fb488de3f2956e34fdd9f418d179c74d7e5002 100644
--- a/src/main/scala/leon/synthesis/rules/IntegerInequalities.scala
+++ b/src/main/scala/leon/synthesis/rules/IntegerInequalities.scala
@@ -9,6 +9,7 @@ import purescala.TreeOps._
 import purescala.TypeTrees._
 import purescala.Definitions._
 import LinearEquations.elimVariable
+import leon.synthesis.Algebra.lcm
 
 class IntegerInequalities(synth: Synthesizer) extends Rule("Integer Inequalities", synth, 300) {
   def attemptToApplyOn(problem: Problem): RuleResult = {
@@ -96,7 +97,7 @@ class IntegerInequalities(synth: Synthesizer) extends Rule("Integer Inequalities
             yield LessEquals(ceilingDiv(lb, IntLiteral(lc)), floorDiv(ub, IntLiteral(uc))))
         RuleFastSuccess(Solution(pre, Set(), witness))
       } else {
-        val L = GCD.lcm((upperBounds ::: lowerBounds).map(_._2))
+        val L = lcm((upperBounds ::: lowerBounds).map(_._2))
         val newUpperBounds: List[Expr] = upperBounds.map{case (bound, coef) => Times(IntLiteral(L/coef), bound)}
         val newLowerBounds: List[Expr] = lowerBounds.map{case (bound, coef) => Times(IntLiteral(L/coef), bound)}
 
diff --git a/src/test/scala/leon/test/synthesis/GCDSuite.scala b/src/test/scala/leon/test/synthesis/AlgebraSuite.scala
similarity index 83%
rename from src/test/scala/leon/test/synthesis/GCDSuite.scala
rename to src/test/scala/leon/test/synthesis/AlgebraSuite.scala
index eaaa20eb6ceccfd23e3075bed20176eaaafd123c..240b613dc9082efe7de72f8b3fdaf4db2f62a43c 100644
--- a/src/test/scala/leon/test/synthesis/GCDSuite.scala
+++ b/src/test/scala/leon/test/synthesis/AlgebraSuite.scala
@@ -2,9 +2,38 @@ package leon.test.synthesis
 
 import org.scalatest.FunSuite
 
-import leon.synthesis.GCD._
+import leon.synthesis.Algebra._
 
-class GCDSuite extends FunSuite {
+class AlgebraSuite extends FunSuite {
+
+  test("remainder") {
+    assert(remainder(1,1) === 0)
+    assert(remainder(2,2) === 0)
+    assert(remainder(2,1) === 0)
+    assert(remainder(0,1) === 0)
+    assert(remainder(0,4) === 0)
+    assert(remainder(1,3) === 1)
+    assert(remainder(1,8) === 1)
+    assert(remainder(4,2) === 0)
+    assert(remainder(17,3) === 2)
+    assert(remainder(25,5) === 0)
+    assert(remainder(26,5) === 1)
+    assert(remainder(29,5) === 4)
+
+    assert(remainder(1,-1) === 0)
+    assert(remainder(-1,1) === 0)
+    assert(remainder(2,-2) === 0)
+    assert(remainder(-2,-2) === 0)
+    assert(remainder(-2,1) === 0)
+    assert(remainder(0,-1) === 0)
+    assert(remainder(1,-2) === 1)
+    assert(remainder(-1,2) === 1)
+    assert(remainder(-1,3) === 2)
+    assert(remainder(-1,-3) === 2)
+    assert(remainder(1,-3) === 1)
+    assert(remainder(17,-3) === 2)
+    assert(remainder(-25,5) === 0)
+  }
 
   test("divide") {
     assert(divide(1,1) === (1, 0))
@@ -170,3 +199,5 @@ class GCDSuite extends FunSuite {
     checkExtendedEuclid(10, -15)
   }
 }
+
+// vim: set ts=4 sw=4 et: