Initial revision
diff --git a/Lib/poly.py b/Lib/poly.py
new file mode 100644
index 0000000..abac4c8
--- /dev/null
+++ b/Lib/poly.py
@@ -0,0 +1,55 @@
+# module 'poly' -- Polynomials
+
+# A polynomial is represented by a list of coefficients, e.g.,
+# [1, 10, 5] represents 1*x**0 + 10*x**1 + 5*x**2 (or 1 + 10x + 5x**2).
+# There is no way to suppress internal zeros; trailing zeros are
+# taken out by normalize().
+
+def normalize(p): # Strip unnecessary zero coefficients
+ n = len(p)
+ while p:
+ if p[n-1]: return p[:n]
+ n = n-1
+ return []
+
+def plus(a, b):
+ if len(a) < len(b): a, b = b, a # make sure a is the longest
+ res = a[:] # make a copy
+ for i in range(len(b)):
+ res[i] = res[i] + b[i]
+ return normalize(res)
+
+def minus(a, b):
+ if len(a) < len(b): a, b = b, a # make sure a is the longest
+ res = a[:] # make a copy
+ for i in range(len(b)):
+ res[i] = res[i] - b[i]
+ return normalize(res)
+
+def one(power, coeff): # Representation of coeff * x**power
+ res = []
+ for i in range(power): res.append(0)
+ return res + [coeff]
+
+def times(a, b):
+ res = []
+ for i in range(len(a)):
+ for j in range(len(b)):
+ res = plus(res, one(i+j, a[i]*b[j]))
+ return res
+
+def power(a, n): # Raise polynomial a to the positive integral power n
+ if n = 0: return [1]
+ if n = 1: return a
+ if n/2*2 = n:
+ b = power(a, n/2)
+ return times(b, b)
+ return times(power(a, n-1), a)
+
+def der(a): # First derivative
+ res = a[1:]
+ for i in range(len(res)):
+ res[i] = res[i] * (i+1)
+ return res
+
+# Computing a primitive function would require rational arithmetic...