From 5f0c2508c53f1d90c2dc8e425825a7f76415847b Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <mathem@tiker.net>
Date: Thu, 30 Jun 2005 20:27:48 +0000
Subject: [PATCH] [pymbolic @
 Arch-1:inform@tiker.net--iam-2005%pymbolic--mainline--1.0--patch-15]
 Continuous work...

---
 LITERATURE        |   2 +
 src/algorithm.py  |  58 +++++++++++++++++
 src/integer.py    |  57 +++++++++++++++++
 src/operation.py  |  26 --------
 src/polynomial.py | 158 +++++++++++++++++++++++++++++++---------------
 src/primitives.py |  27 ++++----
 src/traits.py     |  64 +++++++++++++------
 7 files changed, 283 insertions(+), 109 deletions(-)
 create mode 100644 src/algorithm.py
 create mode 100644 src/integer.py
 delete mode 100644 src/operation.py

diff --git a/LITERATURE b/LITERATURE
index 9327269..f975c8b 100644
--- a/LITERATURE
+++ b/LITERATURE
@@ -1 +1,3 @@
 [Bosch]: S. Bosch, Algebra, 3rd Edition, Springer, 1999
+[Davenport]: J.H. Davenport, Y. Siret, E. Tournier, Computer Algebra,
+  Academic Press, 1988
diff --git a/src/algorithm.py b/src/algorithm.py
new file mode 100644
index 0000000..43cc3e8
--- /dev/null
+++ b/src/algorithm.py
@@ -0,0 +1,58 @@
+import traits
+  
+  
+  
+
+def integer_power(x, n, one=1):
+    # http://c2.com/cgi/wiki?IntegerPowerAlgorithm
+    assert isinstance(n, int)
+
+    if n < 0:
+        raise RuntimeError, "the integer power algorithm does not work for negative numbers"
+      
+    aux = one
+  
+    while n > 0:
+        if n & 1:
+            aux *= x
+            if n == 1:
+                return aux
+        x = x * x
+        n //= 2
+  
+  
+  
+def extended_euclidean(q, r):
+    """Return a tuple (p, a, b) such that p = aq + br, 
+    where p is the greatest common divisor.
+    """
+
+    # see [Davenport], Appendix, p. 214
+
+    t = traits.common_traits(q, r)
+    
+    if t.norm(q) < t.norm(r):
+        p, a, b = extended_euclidean(r, q)
+        return p, b, a
+  
+    Q = 1, 0
+    R = 0, 1
+  
+    while r:
+        quot, t = divmod(q, r)
+        T = Q[0] - quot*R[0], Q[1] - quot*R[1]
+        q, r = r, t
+        Q, R = R, T
+  
+    return q, Q[0], Q[1]
+  
+  
+  
+  
+if __name__ == "__main__":
+    import integer
+    q = integer.Integer(14)
+    r = integer.Integer(22)
+    gcd, a, b = extended_euclidean(q, r)
+    print gcd, "=", a, "*", q, "+", b, "*", r
+    print a*q + b*r
diff --git a/src/integer.py b/src/integer.py
new file mode 100644
index 0000000..a7fc79c
--- /dev/null
+++ b/src/integer.py
@@ -0,0 +1,57 @@
+import rational
+import traits
+
+
+
+
+def _int_or_ni(x):
+    if x is NotImplemented:
+        return x
+    else:
+        return Integer(x)
+
+
+
+
+class Integer(int):
+    def traits(self):
+        return IntegerTraits()
+
+    def __add__(self, other): return _int_or_ni(int.__add__(self, other))
+    def __sub__(self, other): return _int_or_ni(int.__sub__(self, other))
+    def __mul__(self, other): return _int_or_ni(int.__mul__(self, other))
+    def __floordiv__(self, other): return _int_or_ni(int.__floordiv__(self, other))
+    def __mod__(self, other): return _int_or_ni(int.__mod__(self, other))
+    def __pow__(self, other): return _int_or_ni(int.__pow__(self, other))
+
+    def __truediv__(self, other):
+        if self % other != 0:
+            return rational.Rational(self, other)
+        else:
+            return super.__truediv__(self, other)
+
+    __div__ = __truediv__
+
+    def __radd__(self, other): return _int_or_ni(int.__radd__(self, other))
+    def __rsub__(self, other): return _int_or_ni(int.__rsub__(self, other))
+    def __rmul__(self, other): return _int_or_ni(int.__rmul__(self, other))
+    def __rfloordiv__(self, other): return _int_or_ni(int.__rfloordiv__(self, other))
+    def __rmod__(self, other): return _int_or_ni(int.__rmod__(self, other))
+    def __rpow__(self, other): return _int_or_ni(int.__rpow__(self, other))
+
+    def __neg__(self, other): return Integer(int.__neg__(self))
+
+    def __rtruediv__(self, other):
+        if other % self != 0:
+            return rational.Rational(other, self)
+        else:
+            return int.__rtruediv__(self, other)
+
+    __rdiv__ = __rtruediv__
+
+
+
+
+class IntegerTraits(traits.EuclideanRingTraits):
+    def norm(self, x):
+        return abs(x)
diff --git a/src/operation.py b/src/operation.py
deleted file mode 100644
index ac96c1b..0000000
--- a/src/operation.py
+++ /dev/null
@@ -1,26 +0,0 @@
-def degree(x):
-    if isinstance(x, Polynomial):
-        return
-
-
-
-    
-def gcd_extended(x, y):
-    pass
-
-
-
-
-def gcd(x, y):
-    raise NotImplementedError
-
-
-
-
-def lcm_extended(x, y):
-    raise NotImplementedError
-
-
-
-def lcm(x, y):
-    raise NotImplementedError
diff --git a/src/polynomial.py b/src/polynomial.py
index 6565a5d..43fb7ac 100644
--- a/src/polynomial.py
+++ b/src/polynomial.py
@@ -1,10 +1,35 @@
-import tests
+from __future__ import division
+import algorithm
+import traits
+
+
+
+
+def _sort_uniq(data):
+    def sortkey((exp, coeff)): return exp
+    data.sort(key=sortkey)
+    
+    uniq_result = []
+    i = 0
+    last_exp = None
+    for exp, coeff in data:
+        if last_exp == exp:
+            newcoeff = uniq_result[-1][1]+coeff
+            if not newcoeff:
+                uniq_result.pop()
+            else:
+                uniq_result[-1] = last_exp, newcoeff
+
+        else:
+            uniq_result.append((exp, coeff))
+            last_exp = exp
+    return uniq_result
 
 
 
 
 class Polynomial(object):
-    def __init__(self, base, data):
+    def __init__(self, base, data = ((1,1),)):
         self.Base = base
 
         # list of (exponent, coefficient tuples)
@@ -14,6 +39,21 @@ class Polynomial(object):
         
         # Remember the Zen, Luke: Sparse is better than dense.
 
+    def coefficients(self):
+        return [coeff for (exp, coeff) in self.Data]
+
+    def degree(self):
+        try:
+            return self.Data[-1][0]
+        except IndexError:
+            return -1
+
+    def traits(self):
+        return PolynomialTraits()
+
+    def __nonzero__(self):
+        return len(self.Data) != 0
+
     def __neg__(self):
         return Polynomial(self.Base,
                           [(exp, -coeff)
@@ -21,7 +61,7 @@ class Polynomial(object):
 
     def __add__(self, other):
         if not isinstance(other, Polynomial) or other.Base != self.Base:
-            other = Polynomial(other, ((1,1),))
+            other = Polynomial(self.Base, ((0, other),))
 
         iself = 0
         iother = 0
@@ -32,7 +72,7 @@ class Polynomial(object):
             exp_other = other.Data[iother][0]
             if exp_self == exp_other:
                 coeff = self.Data[iself][1] + other.Data[iother][1]
-                if not coeff:
+                if coeff:
                     result.append((exp_self, coeff))
                 iself += 1
                 iother += 1
@@ -76,64 +116,64 @@ class Polynomial(object):
             for o_exp, o_coeff in other.Data:
                 result.append((s_exp+o_exp, s_coeff*o_coeff))
 
-        def sortkey((exp, coeff)): return exp
-        result.sort(key=sortkey)
-
-        uniq_result = []
-        i = 0
-        last_exp = None
-        for exp, coeff in result:
-            if last_exp == exp:
-                newcoeff = uniq_result[-1][1]+coeff
-                if newcoeff:
-                    uniq_result.pop()
-                else:
-                    uniq_result[-1] = last_exp, newcoeff
-
-            else:
-                uniq_result.append((exp, coeff))
-                last_exp = exp
-
-        return Polynomial(self.Base, tuple(uniq_result))
+        return Polynomial(self.Base, tuple(_sort_uniq(result)))
 
     def __rmul__(self, other):
         return Polynomial(self.Base, [(exp, other * coeff)
                                       for exp, coeff in self.Data])
 
     def __pow__(self, other):
-        n = int(other)
-        if n < 0:
-            raise RuntimeError, "negative powers of polynomials not defined"
-        
-        aux = Polynomial(self.Base, ((0, 1),))
-        x = self
-
-        # http://c2.com/cgi/wiki?IntegerPowerAlgorithm
-        while n > 0:
-            if n & 1:
-                aux *= x
-                if n == 1:
-                    return aux
-            x = x * x
-            n //= 2
-    
-        return aux
+        return algorithm.integer_power(self, int(other),
+                                       Polynomial(self.Base, ((0, 1),)))
 
     def __divmod__(self, other):
-        pass
+        if not isinstance(other, Polynomial) or other.Base != self.Base:
+            dm_list = [(exp, divmod(coeff * other)) for exp, coeff in self.Data]
+            return Polynomial(self.Base, [(exp, quot) for (exp, (quot, rem)) in dm_list]),\
+                   Polynomial(self.Base, [(exp, rem) for (exp, (quot, rem)) in dm_list])
+
+        if other.degree() == -1:
+            raise DivisionByZeroError
+
+        quotient = Polynomial(self.Base, ())
+        remainder = self
+        other_lead_coeff = other.Data[-1][1]
+        other_lead_exp = other.Data[-1][0]
+        while remainder.degree() >= other.degree():
+            coeff_factor = remainder.Data[-1][1] / other_lead_coeff
+            deg_diff = remainder.Data[-1][0] - other_lead_exp
+
+            this_fac = Polynomial(self.Base, ((deg_diff, coeff_factor),))
+            quotient += this_fac
+            remainder -= this_fac * other
+        return quotient, remainder
+
+    def __div__(self):
+        q, r = self.__divmod__(self, other)
+        if r.degree() != -1:
+            raise ValueError, "division yielded a remainder"
+        return q
+
+    __truediv__ = __div__
+
+    def __floordiv__(self):
+        return self.__divmod__(self, other)[0]
+
+    def __mod__(self):
+        return self.__divmod__(self, other)[1]
 
     def __str__(self):
         sbase = str(self.Base)
         def stringify_expcoeff((exp, coeff)):
             if exp == 1:
-                if tests.is_one(coeff):
+                if not (coeff-1):
                     return sbase
                 else:
                     return "%s*%s" % (str(coeff), sbase)
             elif exp == 0:
                 return str(coeff)
             else:
-                if tests.is_one(coeff):
+                if not (coeff-1):
                     return "%s**%s" % (sbase, exp) 
                 else:
                     return "%s*%s**%s" % (str(coeff), sbase, exp) 
@@ -141,7 +181,7 @@ class Polynomial(object):
         return "(%s)" % " + ".join(stringify_expcoeff(i) for i in self.Data[::-1])
 
     def __repr__(self):
-        return "%s(%s, %s)" % (self.__class__.__name,
+        return "%s(%s, %s)" % (self.__class__.__name__,
                                repr(self.Base), 
                                repr(self.Data))
         
@@ -151,10 +191,28 @@ class Polynomial(object):
 
 
 
+class PolynomialTraits(traits.EuclideanRingTraits):
+    @staticmethod
+    def norm(x):
+        return x.degree()
+    
+
+
+   
 if __name__ == "__main__":
-    xpoly = Polynomial("x", ((0,1), (2,1)))
-    ypoly = Polynomial("y", ((0,xpoly), (1,2), (2,-xpoly)))
-    print xpoly**3
-    print xpoly + 1
-    print xpoly**3-xpoly
-    print (xpoly*ypoly)**7
+    x = Polynomial("x")
+    y = Polynomial("y")
+    xpoly = x**2 + 1
+    ypoly = -y**2*xpoly + xpoly
+    print xpoly
+    print ypoly
+    u = xpoly*ypoly
+    print u
+    print u**4
+    print
+
+    print 3*xpoly**3 + 1
+    print xpoly 
+    q,r = divmod(3*xpoly**3 + 1, xpoly)
+    print q, r
+
diff --git a/src/primitives.py b/src/primitives.py
index 21c31e5..b03380c 100644
--- a/src/primitives.py
+++ b/src/primitives.py
@@ -1,4 +1,3 @@
-import tests
 import mapper.stringifier
 
 
@@ -8,14 +7,14 @@ class Expression(object):
     def __add__(self, other):
         if not isinstance(other, Expression):
             other = Constant(other)
-        if tests.is_zero(other):
+        if not other:
             return self
         return Sum((self, other))
 
     def __radd__(self, other):
         assert not isinstance(other, Expression)
 
-        if other == 0:
+        if not other:
             return self
         else:
             other = Constant(other)
@@ -24,14 +23,14 @@ class Expression(object):
     def __sub__(self, other):
         if not isinstance(other, Expression):
                 other = Constant(other)
-        if tests.is_zero(other):
+        if not other:
             return self
         return Sum((self, -other))
 
     def __rsub__(self, other):
         assert not isinstance(other, Expression)
 
-        if other == 0:
+        if not other:
             return Negation(self)
         else:
             other = Constant(other)
@@ -40,18 +39,18 @@ class Expression(object):
     def __mul__(self, other):
         if not isinstance(other, Expression):
             other = Constant(other)
-        if tests.is_one(other):
+        if not (other - 1):
             return self
-        if tests.is_zero(other):
+        if not other:
             return Constant(0)
         return Product((self, other))
 
     def __rmul__(self, other):
         assert not isinstance(other, Expression)
 
-        if other == 1:
+        if not (other-1):
             return self
-        elif other == 0:
+        elif not other:
             return Constant(0)
         else:
             return Product((other, self))
@@ -59,7 +58,7 @@ class Expression(object):
     def __div__(self, other):
         if not isinstance(other, Expression):
             other = Constant(other)
-        if isinstance(other, Constant) and other.Value == 1:
+        if not (other-1):
             return self
         return RationalExpression(self, other)
 
@@ -70,17 +69,17 @@ class Expression(object):
     def __pow__(self, other):
         if not isinstance(other, Expression):
             other = Constant(other)
-        if tests.is_zero(other): # exponent zero
+        if not other: # exponent zero
             return Constant(1)
-        elif tests.is_one(other): # exponent one
+        elif not (other-1): # exponent one
             return self
         return Power(self, other)
 
     def __rpow__(self, other):
         assert not isinstance(other, Expression)
-        if tests.is_zero(other): # base zero
+        if not other: # base zero
             return Constant(0)
-        elif tests.is_one(other): # base one
+        elif not (other-1): # base one
             return Constant(1)
         return Power(other, self)
 
diff --git a/src/traits.py b/src/traits.py
index 9fe8d0d..7855c89 100644
--- a/src/traits.py
+++ b/src/traits.py
@@ -1,24 +1,26 @@
 def traits(x):
     try:
-        return x.traits
+        return x.traits()
     except AttributeError:
         if isinstance(x, complex): return ComplexTraits
-        if isinstance(x, float): return FloatTraits
-        if isinstance(x, int): return IntegerTraits
-        return NotImplementedError
+        if isinstance(x, float, int): return RealTraits
+        raise NotImplementedError
+
 
 
 
+def common_traits(*args):
+    def common_traits_two(t_x, t_y):
+        if isinstance(t_y, t_x.__class__):
+            return t_y
+        elif isinstance(t_x, t_y.__class__):
+            return t_x
+        else:
+            raise RuntimeError, "No common traits type between '%s' and '%s'" % \
+                  (t_x.__class__.__name__, t_y.__class__.__name__)
 
-def most_special_traits(x, y):
-    t_x = traits(x)
-    t_y = traits(y)
+    return reduce(common_traits_two, (traits(arg) for arg in args))
 
-    if isinstance(t_y, t_x.__class__):
-        return t_y
-    else:
-        return t_x
-    
 
 
 
@@ -26,17 +28,41 @@ class Traits:
     def one(self): raise NotImplementedError
     def zero(self): raise NotImplementedError
 
-    def degree(self, x):
-        """Returns the degree of the element x.
 
-        "Degree" is used as in the definition of a Euclidean ring,
+
+
+class EuclideanRingTraits(Traits):
+    @staticmethod
+    def norm(x):
+        """Returns the algebraic norm of the element x.
+ 
+        "Norm" is used as in the definition of a Euclidean ring,
         see [Bosch], p. 42
+        """"
         """
         raise NotImplementedError
 
+    @staticmethod
+    def gcd_extended(q, r): 
+        """Returns a tuple """
+        raise NotImplementedError
+ 
     def gcd_extended(self, q, r): 
-        """Returns a tuple ..."""
+        """Returns a tuple """
+
+    @staticmethod
+    def gcd(self): 
+        raise NotImplementedError
+
+    @staticmethod
+    def lcm(a, b): raise NotImplementedError
+
+    @staticmethod
+    def lcm_extended(a, b): raise NotImplementedError
+
+
+
+
+class FieldTraits(Traits):
+    pass
 
-    def gcd(self): raise NotImplementedError
-    def lcm(self): raise NotImplementedError
-    def lcm_extended(self): raise NotImplementedError
-- 
GitLab