1 files changed, 55 insertions, 0 deletions
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...