f(x) = a0 + a1x + anxn の零点を探索する。
その前に、補題を考える。
f(x) = anxn - (a0 + ... + an-1xn-1) n ≧ 1 ai ≧ 0 an > 0
は、x > 0 で f(x) > 0 なら x < x' で f(x') > 0
xf'(x) = nanxn - (a1 + ... + (n-1)an-1xn-1) > (n-1)a1 + ... + an-1xn-1 > 0
より、f(x) > 0 の x > 0 では常に f'(x) > 0。
一般に、
f(x) = a0 + a1x + anxn
に対して、
R > 0, |an|Rn > |a0| + ... + |an-1|Rn-1
なら、
| anRn | - ( | a0 | + ... + | an-1 | Rn-1) > 0 |
だから、
このようなRに対し、R < x で、|f(x)| > 0 となる。
だから、x = 0からxを増していって、f(x) > 0となるところで、Newton法でも用いてf(R) = 0となるRを求めれば、f(x)の零点は、すべて|x| ≦ Rとなる。
これを、
f(x) = 16-52x+48x2+68x3-30x4-39x5+7x6+4x7-10x8-2x9+6x10+x11-x12
対して行うと、
R = 3.33668330031
となった。
from poly import Polydef Newton(f, x):
eps = 1e-9;
g = f.differential();
while abs(f.value(x)) > eps:
x -= f.value(x) / g.value(x)
return xdef convert(f):
p = [ -abs(x) for x in f ]
p[len(p)-1] *= -1;
return Poly(p)f = Poly([ 16, -52, 48, 68, -30, -39, 7, 4, -10, -2, 6, 1, -1 ]);
h = convert(f)x = 1e-3
while h.value(x) < 0:
x += 1e-3
x = Newton(h, x)
print x, h.value(x)
'poly.py
class Poly(list):
def __add__(self, a):
res = Poly([])
m = len(self)
if isinstance(a, Poly):
n = len(a)
if m < n:
res[0:n] = a[0:n]
for i in range(m):
res[i] += self[i]
else:
res[0:m] = self[0:m]
for i in range(n):
res[i] += a[i]
else:
res[0:m] = self[0:m]
if len(res) == 0:
res.append(a)
else:
res[0] += a
res.normalize()
return res
def __radd__(self, a):
return self + a
def __sub__(self, a):
res = Poly([])
m = len(self)
if isinstance(a, Poly):
n = len(a)
if m < n:
res[0:n] = [ -x for x in a[0:n] ]
for i in range(m):
res[i] -= self[i]
else:
res[0:m] = self[0:m]
for i in range(n):
res[i] -= a[i]
else:
res[0:m] = self[0:m]
if len(res) == 0:
res.append(-a)
else:
res[0] -= a
res.normalize()
return res
def __rsub__(self, a):
return -self + a
def __mul__(self, a):
if isinstance(a, Poly):
m = len(self)
n = len(a)
res = Poly([ 0 for i in range(m + n - 1) ])
for i in range(m):
for j in range(n):
res[i+j] += self[i] * a[j]
return res
else:
return Poly([ x * a for x in self ])
def __rmul__(self, a):
return self * a
def __div__(self, a):
if isinstance(a, Poly):
m = len(self)
n = len(a)
res = Poly([ 0 for i in range(m - n + 1) ])
b = Poly([ x for x in self ])
for i in range(m - n, -1, -1):
d = b[i+n-1] / a[n-1]
for j in range(n):
b[i+j] -= d * a[j]
res[i] = d
else:
res = Poly([ x / a for x in self])
res.normalize()
return res
def __mod__(self, a):
m = len(self)
n = len(a)
res = Poly([ 0 for i in range(m - n + 1) ])
b = Poly([ x for x in self ])
for i in range(m - n, -1, -1):
d = b[i+n-1] / a[n-1]
for j in range(n):
b[i+j] -= d * a[j]
res[i] = d
b.normalize()
return b
def __neg__(self):
return Poly([ -x for x in self ])
def __eq__(self, a):
if isinstance(a, Poly):
if len(self) != len(a):
return False
for i in range(0, len(self)):
if self[i] != a[i]:
return False
return True
else:
if len(self) == 0:
return a == 0
elif len(self) == 1:
return self[0] == a
else:
return False
def __ne__(self, a):
return not self.__eq__(a)
def __gt__(self, a):
if isinstance(a, Poly):
if len(self) == len(a):
return abs(self[len(self)-1]) > abs(a[len(a)-1])
else:
return len(self) > len(a)
else:
return True
def normalize(self):
for i in range(len(self) - 1, -1, -1):
if self[i] == 0:
self[i:] = [ ]
else:
return
def value(self, x):
ret = 0
n = len(self)
for i in range(n - 1, -1, -1):
a = self[i]
ret *= x
ret += a
return ret
def differential(self):
n = len(self)
return Poly([ i * self[i] for i in range(1, n) ])
def factorize(self):
ret = [ ]
f = Poly([ x for x in self ])
a = 0
while abs(a) <= f[0]:
if f.value(a) == 0:
factor = Poly([a, -1])
f /= factor
ret.append(factor)
else:
if a > 0:
a = -a
else:
a = 1 - a
# divide by a - x^2
a = 2
while abs(a) <= f[0]:
factor = Poly([a, 0, -1])
mod = f % factor
if mod == 0:
f /= factor
ret.append(factor)
else:
if a > 0:
a = -a
else:
a = 1 - a
ret.append(f)
return ret
def __str__(self):
s = ""
n = len(self)
for i in range(0, n):
if self[i] == 0:
continue
if i == 0:
s = str(self[i])
else:
tmp = ""
if self[i] == 1:
pass
elif self[i] == -1:
tmp = "-"
else:
tmp = str(self[i])
if s != "" and self[i] > 0:
tmp = "+" + tmp
tmp += "x"
if i > 1:
tmp += "^" + str(i)
s += tmp
if s == "":
s = "0"
return s
