|
5
|
1 # module 'poly' -- Polynomials
|
|
|
2
|
|
|
3 # A polynomial is represented by a list of coefficients, e.g.,
|
|
|
4 # [1, 10, 5] represents 1*x**0 + 10*x**1 + 5*x**2 (or 1 + 10x + 5x**2).
|
|
|
5 # There is no way to suppress internal zeros; trailing zeros are
|
|
|
6 # taken out by normalize().
|
|
|
7
|
|
|
8 def normalize(p): # Strip unnecessary zero coefficients
|
|
|
9 n = len(p)
|
|
|
10 while p:
|
|
|
11 if p[n-1]: return p[:n]
|
|
|
12 n = n-1
|
|
|
13 return []
|
|
|
14
|
|
|
15 def plus(a, b):
|
|
|
16 if len(a) < len(b): a, b = b, a # make sure a is the longest
|
|
|
17 res = a[:] # make a copy
|
|
|
18 for i in range(len(b)):
|
|
|
19 res[i] = res[i] + b[i]
|
|
|
20 return normalize(res)
|
|
|
21
|
|
|
22 def minus(a, b):
|
|
1989
|
23 neg_b = map(lambda x: -x, b[:])
|
|
|
24 return plus(a, neg_b)
|
|
5
|
25
|
|
|
26 def one(power, coeff): # Representation of coeff * x**power
|
|
|
27 res = []
|
|
|
28 for i in range(power): res.append(0)
|
|
|
29 return res + [coeff]
|
|
|
30
|
|
|
31 def times(a, b):
|
|
|
32 res = []
|
|
|
33 for i in range(len(a)):
|
|
|
34 for j in range(len(b)):
|
|
|
35 res = plus(res, one(i+j, a[i]*b[j]))
|
|
|
36 return res
|
|
|
37
|
|
|
38 def power(a, n): # Raise polynomial a to the positive integral power n
|
|
559
|
39 if n == 0: return [1]
|
|
|
40 if n == 1: return a
|
|
|
41 if n/2*2 == n:
|
|
5
|
42 b = power(a, n/2)
|
|
|
43 return times(b, b)
|
|
|
44 return times(power(a, n-1), a)
|
|
|
45
|
|
|
46 def der(a): # First derivative
|
|
|
47 res = a[1:]
|
|
|
48 for i in range(len(res)):
|
|
|
49 res[i] = res[i] * (i+1)
|
|
|
50 return res
|
|
|
51
|
|
|
52 # Computing a primitive function would require rational arithmetic...
|
