aboutsummaryrefslogtreecommitdiff
path: root/lib/poly.py
blob: a7c24eced38cb1b908b24366e238cecd8482d340 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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...