diff --git a/LITERATURE b/LITERATURE
index 9327269fec600ddf57e06b08f873282a08fd362a..f975c8b56497071f2b54ead1b271ab779932e817 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 0000000000000000000000000000000000000000..43cc3e84508d1150d465dfc7fbb2219a6b45324a
--- /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 0000000000000000000000000000000000000000..a7fc79c2c2c37433922068ac7a869364d8c7cffb
--- /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 ac96c1b117d8feb77bc8318aac4a07c2cb76e347..0000000000000000000000000000000000000000
--- 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 6565a5dba32d01a81cf7cac633650fb5733043d5..43fb7ac8ea3e9fa983c9230679e3eac9ed132f56 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 21c31e5d0e4d645129aed96dde03eb27a42722cf..b03380c35866c3fb45e9732497317e5531391de9 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 9fe8d0ddae0c0315cd2a63e110c5b75844a31d93..7855c8980c023abbaba66274c9f12ca45c39af50 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