Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/polynomials.py: 9%

1from ..libmp.backend import xrange 

2from .calculus import defun 

3 

4#----------------------------------------------------------------------------# 

5# Polynomials # 

6#----------------------------------------------------------------------------# 

7 

8# XXX: extra precision 

9@defun 

10def polyval(ctx, coeffs, x, derivative=False): 

11 r""" 

12 Given coefficients `[c_n, \ldots, c_2, c_1, c_0]` and a number `x`, 

13 :func:`~mpmath.polyval` evaluates the polynomial 

14 

15 .. math :: 

16 

17 P(x) = c_n x^n + \ldots + c_2 x^2 + c_1 x + c_0. 

18 

19 If *derivative=True* is set, :func:`~mpmath.polyval` simultaneously 

20 evaluates `P(x)` with the derivative, `P'(x)`, and returns the 

21 tuple `(P(x), P'(x))`. 

22 

23 >>> from mpmath import * 

24 >>> mp.pretty = True 

25 >>> polyval([3, 0, 2], 0.5) 

26 2.75 

27 >>> polyval([3, 0, 2], 0.5, derivative=True) 

28 (2.75, 3.0) 

29 

30 The coefficients and the evaluation point may be any combination 

31 of real or complex numbers. 

32 """ 

33 if not coeffs: 

34 return ctx.zero 

35 p = ctx.convert(coeffs[0]) 

36 q = ctx.zero 

37 for c in coeffs[1:]: 

38 if derivative: 

39 q = p + x*q 

40 p = c + x*p 

41 if derivative: 

42 return p, q 

43 else: 

44 return p 

45 

46@defun 

47def polyroots(ctx, coeffs, maxsteps=50, cleanup=True, extraprec=10, 

48 error=False, roots_init=None): 

49 """ 

50 Computes all roots (real or complex) of a given polynomial. 

51 

52 The roots are returned as a sorted list, where real roots appear first 

53 followed by complex conjugate roots as adjacent elements. The polynomial 

54 should be given as a list of coefficients, in the format used by 

55 :func:`~mpmath.polyval`. The leading coefficient must be nonzero. 

56 

57 With *error=True*, :func:`~mpmath.polyroots` returns a tuple *(roots, err)* 

58 where *err* is an estimate of the maximum error among the computed roots. 

59 

60 **Examples** 

61 

62 Finding the three real roots of `x^3 - x^2 - 14x + 24`:: 

63 

64 >>> from mpmath import * 

65 >>> mp.dps = 15; mp.pretty = True 

66 >>> nprint(polyroots([1,-1,-14,24]), 4) 

67 [-4.0, 2.0, 3.0] 

68 

69 Finding the two complex conjugate roots of `4x^2 + 3x + 2`, with an 

70 error estimate:: 

71 

72 >>> roots, err = polyroots([4,3,2], error=True) 

73 >>> for r in roots: 

74 ... print(r) 

75 ... 

76 (-0.375 + 0.59947894041409j) 

77 (-0.375 - 0.59947894041409j) 

78 >>> 

79 >>> err 

80 2.22044604925031e-16 

81 >>> 

82 >>> polyval([4,3,2], roots[0]) 

83 (2.22044604925031e-16 + 0.0j) 

84 >>> polyval([4,3,2], roots[1]) 

85 (2.22044604925031e-16 + 0.0j) 

86 

87 The following example computes all the 5th roots of unity; that is, 

88 the roots of `x^5 - 1`:: 

89 

90 >>> mp.dps = 20 

91 >>> for r in polyroots([1, 0, 0, 0, 0, -1]): 

92 ... print(r) 

93 ... 

94 1.0 

95 (-0.8090169943749474241 + 0.58778525229247312917j) 

96 (-0.8090169943749474241 - 0.58778525229247312917j) 

97 (0.3090169943749474241 + 0.95105651629515357212j) 

98 (0.3090169943749474241 - 0.95105651629515357212j) 

99 

100 **Precision and conditioning** 

101 

102 The roots are computed to the current working precision accuracy. If this 

103 accuracy cannot be achieved in ``maxsteps`` steps, then a 

104 ``NoConvergence`` exception is raised. The algorithm internally is using 

105 the current working precision extended by ``extraprec``. If 

106 ``NoConvergence`` was raised, that is caused either by not having enough 

107 extra precision to achieve convergence (in which case increasing 

108 ``extraprec`` should fix the problem) or too low ``maxsteps`` (in which 

109 case increasing ``maxsteps`` should fix the problem), or a combination of 

110 both. 

111 

112 The user should always do a convergence study with regards to 

113 ``extraprec`` to ensure accurate results. It is possible to get 

114 convergence to a wrong answer with too low ``extraprec``. 

115 

116 Provided there are no repeated roots, :func:`~mpmath.polyroots` can 

117 typically compute all roots of an arbitrary polynomial to high precision:: 

118 

119 >>> mp.dps = 60 

120 >>> for r in polyroots([1, 0, -10, 0, 1]): 

121 ... print(r) 

122 ... 

123 -3.14626436994197234232913506571557044551247712918732870123249 

124 -0.317837245195782244725757617296174288373133378433432554879127 

125 0.317837245195782244725757617296174288373133378433432554879127 

126 3.14626436994197234232913506571557044551247712918732870123249 

127 >>> 

128 >>> sqrt(3) + sqrt(2) 

129 3.14626436994197234232913506571557044551247712918732870123249 

130 >>> sqrt(3) - sqrt(2) 

131 0.317837245195782244725757617296174288373133378433432554879127 

132 

133 **Algorithm** 

134 

135 :func:`~mpmath.polyroots` implements the Durand-Kerner method [1], which 

136 uses complex arithmetic to locate all roots simultaneously. 

137 The Durand-Kerner method can be viewed as approximately performing 

138 simultaneous Newton iteration for all the roots. In particular, 

139 the convergence to simple roots is quadratic, just like Newton's 

140 method. 

141 

142 Although all roots are internally calculated using complex arithmetic, any 

143 root found to have an imaginary part smaller than the estimated numerical 

144 error is truncated to a real number (small real parts are also chopped). 

145 Real roots are placed first in the returned list, sorted by value. The 

146 remaining complex roots are sorted by their real parts so that conjugate 

147 roots end up next to each other. 

148 

149 **References** 

150 

151 1. http://en.wikipedia.org/wiki/Durand-Kerner_method 

152 

153 """ 

154 if len(coeffs) <= 1: 

155 if not coeffs or not coeffs[0]: 

156 raise ValueError("Input to polyroots must not be the zero polynomial") 

157 # Constant polynomial with no roots 

158 return [] 

159 

160 orig = ctx.prec 

161 tol = +ctx.eps 

162 with ctx.extraprec(extraprec): 

163 deg = len(coeffs) - 1 

164 # Must be monic 

165 lead = ctx.convert(coeffs[0]) 

166 if lead == 1: 

167 coeffs = [ctx.convert(c) for c in coeffs] 

168 else: 

169 coeffs = [c/lead for c in coeffs] 

170 f = lambda x: ctx.polyval(coeffs, x) 

171 if roots_init is None: 

172 roots = [ctx.mpc((0.4+0.9j)**n) for n in xrange(deg)] 

173 else: 

174 roots = [None]*deg; 

175 deg_init = min(deg, len(roots_init)) 

176 roots[:deg_init] = list(roots_init[:deg_init]) 

177 roots[deg_init:] = [ctx.mpc((0.4+0.9j)**n) for n 

178 in xrange(deg_init,deg)] 

179 err = [ctx.one for n in xrange(deg)] 

180 # Durand-Kerner iteration until convergence 

181 for step in xrange(maxsteps): 

182 if abs(max(err)) < tol: 

183 break 

184 for i in xrange(deg): 

185 p = roots[i] 

186 x = f(p) 

187 for j in range(deg): 

188 if i != j: 

189 try: 

190 x /= (p-roots[j]) 

191 except ZeroDivisionError: 

192 continue 

193 roots[i] = p - x 

194 err[i] = abs(x) 

195 if abs(max(err)) >= tol: 

196 raise ctx.NoConvergence("Didn't converge in maxsteps=%d steps." \ 

197 % maxsteps) 

198 # Remove small real or imaginary parts 

199 if cleanup: 

200 for i in xrange(deg): 

201 if abs(roots[i]) < tol: 

202 roots[i] = ctx.zero 

203 elif abs(ctx._im(roots[i])) < tol: 

204 roots[i] = roots[i].real 

205 elif abs(ctx._re(roots[i])) < tol: 

206 roots[i] = roots[i].imag * 1j 

207 roots.sort(key=lambda x: (abs(ctx._im(x)), ctx._re(x))) 

208 if error: 

209 err = max(err) 

210 err = max(err, ctx.ldexp(1, -orig+1)) 

211 return [+r for r in roots], +err 

212 else: 

213 return [+r for r in roots]