Welcome, guest | Sign In | My Account | Store | Cart
```"""\
This class implements polynomial functions over a single variable. It
represents the polynomial as a list of numbers and allows most
arithmetic operations, using conventional Python syntax. It does not
do symbolic manipulations. Instead, you can do things like this:

>>> x = SimplePolynomial()
'X**2 - 1'
15
>>> for i in range(4):
...     polynomial = (x+1)**i
...     i, str(polynomial), polynomial(1)
(0, '1', 1)
(1, 'X + 1', 2)
(2, 'X**2 + 2*X + 1', 4)
(3, 'X**3 + 3*X**2 + 3*X + 1', 8)
"""

from __future__ import division, generators

try:
from itertools import izip_longest
except ImportError:
# The izip_longest function was added in version 2.6
# If we can't find it, use an equivalent implementation.
from itertools import chain, repeat
class ZipExhausted(Exception):
pass
def izip_longest(*args, **kwds):
# izip_longest('ABCD', 'xy', fillvalue='-') --> Ax By C- D-
fillvalue = kwds.get('fillvalue')
counter = [len(args) - 1]
def sentinel():
if not counter[0]:
raise ZipExhausted
counter[0] -= 1
yield fillvalue
fillers = repeat(fillvalue)
iterators = [chain(it, sentinel(), fillers) for it in args]
try:
while iterators:
yield tuple(map(next, iterators))
except ZipExhausted:
pass

try:
from numbers import Number
except ImportError:
# The numbers module was added in version 2.6
# If we can't find it, use an equivalent implementation.
Number = (int, float, long, complex)

class SimplePolynomial(object):

def __init__(self, terms=[0,1]):
"""\
>>> str(SimplePolynomial())
'X'
"""
try:
while terms[-1] == 0:
del terms[-1]
except IndexError:
pass
self.terms = list(terms)

def __str__(self):
"""\
"""
l = len(self.terms)
if l == 0:
return '0'
if l == 1:
return str(self.terms[0])
result = []
for i, c in reversed(list(enumerate(self.terms))):
if c == 0:
continue
if c < 0:
result.append('-')
c = - c
else:
result.append('+')

if c == 1:
if i == 0:
result.append('1')
elif i == 1:
result.append('X')
else:
result.append('X**%g' % i)
else:
if i == 0:
result.append('%g' % c)
elif i == 1:
result.append('%g*X' % c)
else:
result.append('%g*X**%d' % (c, i))
if len(result) == 0:
return '0'
if result[0] == '-':
result[1] = '-'+result[1]
del result[0]
return ' '.join(result)

"""\
>>> str(SimplePolynomial() + 1)
'X + 1'
>>> str(1 + SimplePolynomial())
'X + 1'
>>> str(SimplePolynomial() + SimplePolynomial())
'2*X'
"""
if len(self.terms) == 0:
return other
if isinstance(other, Number):
terms = self.terms[:]
terms[0] += other
return SimplePolynomial(terms)
return SimplePolynomial([
for pair in izip_longest(self.terms, other.terms, fillvalue=0)
])

def __neg__(self):
"""\
>>> str(- SimplePolynomial())
'-X'
"""
return SimplePolynomial([-c for c in self.terms])

def __sub__(self, other):
"""\
>>> str(SimplePolynomial() - 1)
'X - 1'
>>> str(SimplePolynomial() - SimplePolynomial())
'0'
"""
return self + -other

def __rsub__(self, other):
"""\
>>> str(1 - SimplePolynomial())
'-X + 1'
"""
return -self + other

def __mul__(self, other):
"""\
>>> str(SimplePolynomial() * 2)
'2*X'
>>> str(2 * SimplePolynomial())
'2*X'
>>> str(SimplePolynomial() * SimplePolynomial())
'X**2'
"""
if isinstance(other, Number):
return SimplePolynomial([c * other for c in self.terms])
terms = [0]*(len(self.terms)+len(other.terms))
for i1, c1 in enumerate(self.terms):
for i2, c2 in enumerate(other.terms):
terms[i1+i2] += c1*c2
return SimplePolynomial(terms)

# Since multiplication is commutative, reuse __mul__
__rmul__ = __mul__

def __truediv__(self, other):
"""\
Implements some simple forms of division.  See
http://en.wikipedia.org/wiki/Synthetic_division
for details.

>>> str(SimplePolynomial() / 2)
'0.5*X'
>>> x = SimplePolynomial()
>>> quotient, remainder  = (x**3 - 12*x**2 - 42) / (x - 3)
>>> str(quotient), str(remainder)
('X**2 - 9*X - 27', '-123')
"""
if isinstance(other, Number):
return SimplePolynomial([c / other for c in self.terms])
if len(other.terms) == 1:
return self/other.terms[0]
assert len(other.terms) == 2
dividend = self.terms[:]
divisor = other.terms[:]
assert divisor[-1] == 1
xi = 0
result = []
for i in xrange(-1, -len(dividend) - 1, -1):
xi = dividend[i] - xi * divisor[0]
result.insert(0, xi)
return SimplePolynomial(result[1:]), result[0]
raise NotImplementedError('synthetic division')

__div__ = __truediv__

def __eq__(self, other):
"""\
>>> SimplePolynomial() == SimplePolynomial()
True
>>> SimplePolynomial() - SimplePolynomial() == 0
True
"""
if isinstance(other, SimplePolynomial):
return self.terms == other.terms
if isinstance(other, Number):
if len(self.terms) > 1:
return False
try:
self.terms[0] == other
except IndexError:
return other == 0
return False

def __ne__(self, other):
return not self == other

def copy(self):
return SimplePolynomial(self.terms[:])

def __pow__(self, exponent):
"""\
Uses the Russian Peasant Multiplication algorithm.

>>> str(SimplePolynomial() ** 2)
'X**2'
>>> str(SimplePolynomial() ** 0.5)
Traceback (most recent call last):
...
NotImplementedError: exponent is not an integer
>>> str(SimplePolynomial() ** -1)
Traceback (most recent call last):
...
NotImplementedError: exponent is less than zero
"""
if not isinstance(exponent, int):
raise NotImplementedError('exponent is not an integer')
if exponent < 0:
raise NotImplementedError('exponent is less than zero')
tmp = self.copy()
result = SimplePolynomial([1])
while exponent > 0:
if exponent & 1:
result *= tmp
tmp *= tmp
exponent >>= 1
return result

def __call__(self, x):
"""\
Evaluate the polynomial for the given value using the Horner scheme.

>>> SimplePolynomial()(1)
1
"""
result = 0
for c in reversed(self.terms):
result = result * x + c
return result

def derivative(self):
"""\
>>> str(SimplePolynomial().derivative())
'1'
"""
terms = [i*c for i, c in enumerate(self.terms)]
return SimplePolynomial(terms[1:])

def integrate(self, const=0):
"""\
>>> str(SimplePolynomial().integrate())
'0.5*X**2'
"""
terms = [const]
terms.extend([c/(i+1) for i, c in enumerate(self.terms)])
return SimplePolynomial(terms)

if __name__ == '__main__':
import doctest
doctest.testmod()
```

#### Diff to Previous Revision

```--- revision 6 2010-01-24 05:22:57
+++ revision 7 2014-02-13 20:55:55
@@ -1,65 +1,79 @@
-import operator
-
-Number = (int, float, long, complex)
-#sum_v = lambda *a: sum(a)
+"""\
+This class implements polynomial functions over a single variable. It
+represents the polynomial as a list of numbers and allows most
+arithmetic operations, using conventional Python syntax. It does not
+do symbolic manipulations. Instead, you can do things like this:
+
+>>> x = SimplePolynomial()
+'X**2 - 1'
+15
+>>> for i in range(4):
+...     polynomial = (x+1)**i
+...     i, str(polynomial), polynomial(1)
+(0, '1', 1)
+(1, 'X + 1', 2)
+(2, 'X**2 + 2*X + 1', 4)
+(3, 'X**3 + 3*X**2 + 3*X + 1', 8)
+"""
+
+from __future__ import division, generators
+
+
+try:
+    from itertools import izip_longest
+except ImportError:
+    # The izip_longest function was added in version 2.6
+    # If we can't find it, use an equivalent implementation.
+    from itertools import chain, repeat
+    class ZipExhausted(Exception):
+        pass
+    def izip_longest(*args, **kwds):
+        # izip_longest('ABCD', 'xy', fillvalue='-') --> Ax By C- D-
+        fillvalue = kwds.get('fillvalue')
+        counter = [len(args) - 1]
+        def sentinel():
+            if not counter[0]:
+                raise ZipExhausted
+            counter[0] -= 1
+            yield fillvalue
+        fillers = repeat(fillvalue)
+        iterators = [chain(it, sentinel(), fillers) for it in args]
+        try:
+            while iterators:
+                yield tuple(map(next, iterators))
+        except ZipExhausted:
+            pass
+
+try:
+    from numbers import Number
+except ImportError:
+    # The numbers module was added in version 2.6
+    # If we can't find it, use an equivalent implementation.
+    Number = (int, float, long, complex)
+

class SimplePolynomial(object):
+
def __init__(self, terms=[0,1]):
+        """\
+>>> str(SimplePolynomial())
+'X'
+"""
try:
while terms[-1] == 0:
del terms[-1]
except IndexError:
pass
-        self.terms = terms
-        if len(self.terms) == 0:
-            return other
-        if isinstance(other, Number):
-            terms = self.terms[:]
-            terms[0] += other
-        elif len(self.terms) <= len(other.terms):
-            terms = self.terms[:]
-            terms.extend([0]*(len(other.terms) - len(self.terms)))
-            terms = map(operator.add, other.terms, terms)
-        else:
-            return other + self
-        return SimplePolynomial(terms)
-        "addition is commutative, so we can just switch the arguments"
-        return self + other
-    def __neg__(self):
-        return SimplePolynomial([-c for c in self.terms])
-    def __sub__(self, other):
-        return self + -other
-    def __rsub__(self, other):
-        return -self + other
-    def __mul__(self, other):
-        if isinstance(other, Number):
-            return SimplePolynomial([other*c for c in self.terms])
-        terms = [0]*(len(self.terms)+len(other.terms))
-        for i1, c1 in enumerate(self.terms):
-            for i2, c2 in enumerate(other.terms):
-                terms[i1+i2] += c1*c2
-        return SimplePolynomial(terms)
-    def __rmul__(self, other):
-        "multiplication is commutative, so we can just switch the arguments"
-        return self * other
-    def copy(self):
-        return SimplePolynomial(self.terms[:])
-    def __pow__(self, exponent):
-        "uses the Russian Peasant Multiplication algorithm"
-        if not isinstance(exponent, int):
-            raise NotImplementedError
-        tmp = self.copy()
-        result = SimplePolynomial([1])
-        while exponent > 0:
-            if exponent & 1:
-                result *= tmp
-            tmp *= tmp
-            exponent >>= 1
-        return result
+        self.terms = list(terms)
+
def __str__(self):
-        "needs some work, but adequate"
+        """\
+"""
l = len(self.terms)
if l == 0:
return '0'
@@ -95,25 +109,186 @@
result[1] = '-'+result[1]
del result[0]
return ' '.join(result)
+
+        """\
+>>> str(SimplePolynomial() + 1)
+'X + 1'
+>>> str(1 + SimplePolynomial())
+'X + 1'
+>>> str(SimplePolynomial() + SimplePolynomial())
+'2*X'
+"""
+        if len(self.terms) == 0:
+            return other
+        if isinstance(other, Number):
+            terms = self.terms[:]
+            terms[0] += other
+            return SimplePolynomial(terms)
+        return SimplePolynomial([
+            for pair in izip_longest(self.terms, other.terms, fillvalue=0)
+            ])
+
+
+    def __neg__(self):
+        """\
+>>> str(- SimplePolynomial())
+'-X'
+"""
+        return SimplePolynomial([-c for c in self.terms])
+
+    def __sub__(self, other):
+        """\
+>>> str(SimplePolynomial() - 1)
+'X - 1'
+>>> str(SimplePolynomial() - SimplePolynomial())
+'0'
+"""
+        return self + -other
+
+    def __rsub__(self, other):
+        """\
+>>> str(1 - SimplePolynomial())
+'-X + 1'
+"""
+        return -self + other
+
+    def __mul__(self, other):
+        """\
+>>> str(SimplePolynomial() * 2)
+'2*X'
+>>> str(2 * SimplePolynomial())
+'2*X'
+>>> str(SimplePolynomial() * SimplePolynomial())
+'X**2'
+"""
+        if isinstance(other, Number):
+            return SimplePolynomial([c * other for c in self.terms])
+        terms = [0]*(len(self.terms)+len(other.terms))
+        for i1, c1 in enumerate(self.terms):
+            for i2, c2 in enumerate(other.terms):
+                terms[i1+i2] += c1*c2
+        return SimplePolynomial(terms)
+
+    # Since multiplication is commutative, reuse __mul__
+    __rmul__ = __mul__
+
+    def __truediv__(self, other):
+        """\
+Implements some simple forms of division.  See
+http://en.wikipedia.org/wiki/Synthetic_division
+for details.
+
+>>> str(SimplePolynomial() / 2)
+'0.5*X'
+>>> x = SimplePolynomial()
+>>> quotient, remainder  = (x**3 - 12*x**2 - 42) / (x - 3)
+>>> str(quotient), str(remainder)
+('X**2 - 9*X - 27', '-123')
+"""
+        if isinstance(other, Number):
+            return SimplePolynomial([c / other for c in self.terms])
+        if len(other.terms) == 1:
+            return self/other.terms[0]
+        assert len(other.terms) == 2
+        dividend = self.terms[:]
+        divisor = other.terms[:]
+        assert divisor[-1] == 1
+        xi = 0
+        result = []
+        for i in xrange(-1, -len(dividend) - 1, -1):
+            xi = dividend[i] - xi * divisor[0]
+            result.insert(0, xi)
+        return SimplePolynomial(result[1:]), result[0]
+        raise NotImplementedError('synthetic division')
+
+    __div__ = __truediv__
+
+    def __eq__(self, other):
+        """\
+>>> SimplePolynomial() == SimplePolynomial()
+True
+>>> SimplePolynomial() - SimplePolynomial() == 0
+True
+"""
+        if isinstance(other, SimplePolynomial):
+            return self.terms == other.terms
+        if isinstance(other, Number):
+            if len(self.terms) > 1:
+                return False
+            try:
+                self.terms[0] == other
+            except IndexError:
+                return other == 0
+        return False
+
+    def __ne__(self, other):
+        return not self == other
+
+    def copy(self):
+        return SimplePolynomial(self.terms[:])
+
+    def __pow__(self, exponent):
+        """\
+Uses the Russian Peasant Multiplication algorithm.
+
+>>> str(SimplePolynomial() ** 2)
+'X**2'
+>>> str(SimplePolynomial() ** 0.5)
+Traceback (most recent call last):
+...
+NotImplementedError: exponent is not an integer
+>>> str(SimplePolynomial() ** -1)
+Traceback (most recent call last):
+...
+NotImplementedError: exponent is less than zero
+"""
+        if not isinstance(exponent, int):
+            raise NotImplementedError('exponent is not an integer')
+        if exponent < 0:
+            raise NotImplementedError('exponent is less than zero')
+        tmp = self.copy()
+        result = SimplePolynomial([1])
+        while exponent > 0:
+            if exponent & 1:
+                result *= tmp
+            tmp *= tmp
+            exponent >>= 1
+        return result
+
def __call__(self, x):
-        "evaluate the polynomial for the given value using the Horner scheme"
+        """\
+Evaluate the polynomial for the given value using the Horner scheme.
+
+>>> SimplePolynomial()(1)
+1
+"""
result = 0
for c in reversed(self.terms):
-            result = result*x + c
+            result = result * x + c
return result
+
def derivative(self):
+        """\
+>>> str(SimplePolynomial().derivative())
+'1'
+"""
terms = [i*c for i, c in enumerate(self.terms)]
return SimplePolynomial(terms[1:])
+
def integrate(self, const=0):
+        """\
+>>> str(SimplePolynomial().integrate())
+'0.5*X**2'
+"""
terms = [const]
-        terms.extend([c/float(i+1) for i, c in enumerate(self.terms)])
+        terms.extend([c/(i+1) for i, c in enumerate(self.terms)])
return SimplePolynomial(terms)

+
if __name__ == '__main__':
-    print 'testing...'
-    x = SimplePolynomial()
-    for expression in 'x', 'x+1', '1+x', 'x+x', '-x', 'x-1', '1-x', '2*x', 'x*2', 'x*x', '(x+1)*(x-1)':
-        print expression, '\t=>', eval(expression)
-    print 'testing exponentiation...'
-    for i in range(4):
-        print '(x+1)**%d' % i, '\t=>', (x+1)**i
+    import doctest
+    doctest.testmod()
```