Needed, symbolic math in Python

[Dan Bishop]

Paul Rubin <http://phr.cx@NOSPAM.invalid> wrote in message news:<7x659zp1o1.fsf_-_ at ruckus.brouhaha.com>...
> I wonder if anyone can recommend any simple symbolic algebra modules
> in Python.  I just need to do some straightforward polynomial
> arithmetic, no fancy calculus or special functions anything like that.
> But I do need arbitrary precision rational coefficients in the
> polynomials.  Thanks.

Here's a polynomial class I wrote a couple of years ago:

from __future__ import division

def _reverse(seq):       # returns seq's elements in reverse order
   return [seq[len(seq)-i-1] for i in range(len(seq))]

def _strterm(coef, power):
   # returns the simplest string representation for coef*x^power
   if coef == 0:
      return ""
   if power == 0:
      return str(coef)
   if power == 1:
      if coef == 1:
         return "x"
      return "%gx" % coef
   if coef == 1:
      return "x^%d" % power
   return "%gx^%d" % (coef, power)

class Polynomial(object):
   def __init__(self, coefs):
      # coefs[i] == coefficient for x**i term
      self.__coefs = coefs[:]
      while len(self.__coefs) > 0 and self.__coefs[-1] == 0:
         self.__coefs.pop()
   def __call__(self, x):
      # Use Horner's algorithm to evaluate the polynomial
      result = 0.
      for c in _reverse(self.__coefs):
         result = result * x + c
      return result
   def __repr__(self):
      return "Polynomial(%s)" % self.__coefs
   def __str__(self):
      # string representation as an expression in terms of x
      return "+".join([_strterm(self.__coefs[i], i) \
                      for i in range(len(self.__coefs)) \
                      if self.__coefs[i]]).replace("+-", "-")
   def __nonzero__(self):
      return len(self.__coefs) > 0
   def __pos__(self):
      return self
   def __neg__(self):
      return Polynomial([-c for c in self.__coefs])
   def __add__(self, other):
      # ensure operands are of the same type
      if not isinstance(other, Polynomial):
         other = Polynomial([other])
      # ensure coefficient lists are the same length
      coefs1 = self.__coefs[:]
      coefs2 = other.__coefs[:]
      if len(coefs1) < len(coefs2):
         coefs1, coefs2 = coefs2, coefs1
      for i in range(len(coefs1) - len(coefs2)):
         coefs2.append(0)
      # finally, perform the addition
      return Polynomial([coefs1[i] + coefs2[i] for i in range(len(coefs1))])
   __radd__ = __add__
   def __sub__(self, other):
      return self + -other
   def __mul__(self, other):
      if isinstance(other, Polynomial):
         coefs = [0] * (len(self.__coefs) + len(other.__coefs) - 1)
         for i in range(len(self.__coefs)):
            for j in range(len(other.__coefs)):
               coefs[i+j] += self.__coefs[i] * other.__coefs[j]
         return Polynomial(coefs)
      return Polynomial([other * c for c in self.__coefs])
   __rmul__ = __mul__
   def __div__(self, other):
      try:
         other(0)
         return lambda x: self(x) / other(x)
      except:
         return Polynomial([c / other for c in self.__coefs])