Commit | Line | Data |
---|---|---|
920dae64 AT |
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 n: | |
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): | |
23 | neg_b = map(lambda x: -x, b[:]) | |
24 | return plus(a, neg_b) | |
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 | |
39 | if n == 0: return [1] | |
40 | if n == 1: return a | |
41 | if n/2*2 == n: | |
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... |