-
Notifications
You must be signed in to change notification settings - Fork 0
/
roots.py
60 lines (54 loc) · 1.21 KB
/
roots.py
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
56
57
58
59
60
def differentiate(poly):
ret = []
n = len(poly) - 1
for i in xrange(n):
ret.append(poly[i] * (n - i))
return ret
def solveNaive(poly):
if len(poly) == 2:
return [ -poly[1]/poly[0] ]
a, b, c = poly
r = b**2 - 4*a*c
if r < 0:
return []
elif r == 0:
return [ -b/(2.0 * a) ]
else:
return [ (-b - r**0.5)/(2.0 * a), (-b + r**0.5)/(2.0 * a)]
def evaluate(poly, x0):
ret = 0
n = len(poly)
for i in xrange(n):
ret += poly[i] * (x0 ** (n-i-1))
return ret
def solve(poly):
n = len(poly) - 1
if n <= 2:
return solveNaive(poly)
derivative = differentiate(poly)
sols = solve(derivative)
sols = [-L-1] + sols + [L+1]
ret = []
for i in xrange (len(sols)-1):
x1, x2 = sols[i], sols[i+1]
y1, y2 = evaluate(poly, x1), evaluate(poly, x2)
if y1 * y2 > 0: continue
if y1 > y2:
y1, y2 = y2, y1
x1, x2 = x2, x1
for it in xrange (100):
mx = (x1 + x2) / 2
my = evaluate(poly, mx)
if y1 * my > 0:
x1, y1 = mx, my
else:
x2, y2 = mx, my
ret.append((x1 + x2) / 2)
ret.sort()
return ret
L = 25.0
for _ in xrange (input()):
n = input()
poly = map(float, raw_input().split())
ret = map(str, solve(poly))
print " ".join(ret)