diff --git a/src/algorithm.py b/src/algorithm.py
index 43cc3e84508d1150d465dfc7fbb2219a6b45324a..846447b06a6a19464033788436af77b41f45591c 100644
--- a/src/algorithm.py
+++ b/src/algorithm.py
@@ -22,6 +22,12 @@ def integer_power(x, n, one=1):
   
   
   
+def gcd(q, r):
+    return extended_euclidean(q, r)[0]
+
+
+
+
 def extended_euclidean(q, r):
     """Return a tuple (p, a, b) such that p = aq + br, 
     where p is the greatest common divisor.
diff --git a/src/compiler.py b/src/compiler.py
index 0f424d2b9ed470fb2d9caaa4f4b5a0d94e3caf36..67318999a59b7de77cf550985f3f1d6a17e19f45 100644
--- a/src/compiler.py
+++ b/src/compiler.py
@@ -2,6 +2,41 @@ import math
 
 import pymbolic
 import pymbolic.mapper.dependency
+from pymbolic.mapper.stringifier import ReprMapper, PREC_NONE, PREC_SUM, PREC_POWER
+
+
+
+
+class CompileMapper(ReprMapper):
+    def map_polynomial(self, expr, enclosing_prec):
+        # Use Horner's scheme to evaluate the polynomial
+
+        sbase = self(expr.base, PREC_POWER)
+        def stringify_exp(exp):
+            if exp == 0:
+                return ""
+            elif exp == 1:
+                return "*%s" % sbase
+            else:
+                return "*%s**%s" % (sbase, exp)
+
+        result = ""
+        rev_data = expr.data[::-1]
+        for i, (exp, coeff) in enumerate(rev_data):
+            if i+1 < len(rev_data):
+                next_exp = rev_data[i+1][0]
+            else:
+                next_exp = 0
+            result = "(%s+%s)%s" % (result, self(coeff, PREC_SUM), 
+                    stringify_exp(exp-next_exp))
+        #print "A", result
+        #print "B", expr
+
+        if enclosing_prec > PREC_SUM and len(expr.data) > 1:
+            return "(%s)" % result
+        else:
+            return result
+
 
 
 
@@ -17,7 +52,6 @@ class CompiledExpression:
         self.__compile__()
 
     def __compile__(self):
-        from pymbolic.mapper.stringifier import ReprMapper, PREC_NONE
         ctx = self.context()
 
         used_variables = pymbolic.get_dependencies(self._Expression)
@@ -28,7 +62,7 @@ class CompiledExpression:
         all_variables = self._Variables + used_variables
 
         expr_s = "lambda %s:%s" % (",".join(str(v) for v in all_variables), 
-                                   ReprMapper()(self._Expression, PREC_NONE))
+                                   CompileMapper()(self._Expression, PREC_NONE))
         self.__call__ = eval(expr_s, ctx)
     
     def __getinitargs__(self):
diff --git a/src/mapper/evaluator.py b/src/mapper/evaluator.py
index dac5e051b9a23e0911883df7f4c924cd82a71edc..f205a0cdd9ea8d9a27b49a782fa745d2139ae320 100644
--- a/src/mapper/evaluator.py
+++ b/src/mapper/evaluator.py
@@ -54,9 +54,18 @@ class EvaluationMapper(RecursiveMapper):
 
     def map_polynomial(self, expr):
         import pymbolic
-        return pymbolic.sum(
-                self.rec(coeff)*self.rec(expr.base)**exp
-                for exp,coeff in expr.data)
+
+        # evaluate using Horner's scheme
+        result = 0
+        rev_data = expr.data[::-1]
+        for i, (exp, coeff) in enumerate(rev_data):
+            if i+1 < len(rev_data):
+                next_exp = rev_data[i+1][0]
+            else:
+                next_exp = 0
+            result = (result+coeff)*expr.base**(exp-next_exp)
+
+        return result
 
     def map_list(self, expr):
         return [self.rec(child) for child in expr.Children]
diff --git a/src/polynomial.py b/src/polynomial.py
index abffdf193a860bbffdac5c8bcf4701932818601e..fd73acbf9a836234e28620f8b95276a6bdb3352e 100644
--- a/src/polynomial.py
+++ b/src/polynomial.py
@@ -143,8 +143,11 @@ class Polynomial(Expression):
 
     def __mul__(self, other):
         if not isinstance(other, Polynomial):
-            return Polynomial(self.Base, [(exp, coeff * other)
-                                          for exp, coeff in self.Data])
+            if other == self.Base:
+                other = Polynomial(self.Base)
+            else:
+                return Polynomial(self.Base, [(exp, coeff * other)
+                                              for exp, coeff in self.Data])
 
         if other.Base != self.Base:
             assert self.VarLess == other.VarLess
@@ -342,14 +345,20 @@ if __name__ == "__main__":
     x = Polynomial(pymbolic.var("x"))
     y = Polynomial(pymbolic.var("y"))
 
-    # NOT WORKING INTRODUCE TESTS
-    u = (x+y)**5
-    v = x+y
-    #u = x+1
-    #v = 3*x+1
-    q, r = divmod(u, v)
-    print q, "R", r
-    print q*v 
-    print "REASSEMBLY:", q*v + r
+    u = (x+1)**5
+    v = pymbolic.evaluate_kw(u, x=x)
+    print u
+    print v
+
+    if False:
+        # NOT WORKING INTRODUCE TESTS
+        u = (x+y)**5
+        v = x+y
+        #u = x+1
+        #v = 3*x+1
+        q, r = divmod(u, v)
+        print q, "R", r
+        print q*v 
+        print "REASSEMBLY:", q*v + r