"""\
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()
>>> quadratic = (x+1)*(x-1)
>>> str(quadratic)
'X**2 - 1'
>>> quadratic(4)
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

from operator import add

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):
        """\
Needs some work, but adequate.
"""
        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)

    def __add__(self, other):
        """\
>>> 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([
            add(*pair)
            for pair in izip_longest(self.terms, other.terms, fillvalue=0)
            ])

    # Since addition is commutative, reuse __add__
    __radd__ = __add__

    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()