Coverage for /usr/lib/python3/dist-packages/sympy/polys/rootisolation.py: 10%
1092 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Real and complex root isolation and refinement algorithms. """
4from sympy.polys.densearith import (
5 dup_neg, dup_rshift, dup_rem,
6 dup_l2_norm_squared)
7from sympy.polys.densebasic import (
8 dup_LC, dup_TC, dup_degree,
9 dup_strip, dup_reverse,
10 dup_convert,
11 dup_terms_gcd)
12from sympy.polys.densetools import (
13 dup_clear_denoms,
14 dup_mirror, dup_scale, dup_shift,
15 dup_transform,
16 dup_diff,
17 dup_eval, dmp_eval_in,
18 dup_sign_variations,
19 dup_real_imag)
20from sympy.polys.euclidtools import (
21 dup_discriminant)
22from sympy.polys.factortools import (
23 dup_factor_list)
24from sympy.polys.polyerrors import (
25 RefinementFailed,
26 DomainError,
27 PolynomialError)
28from sympy.polys.sqfreetools import (
29 dup_sqf_part, dup_sqf_list)
32def dup_sturm(f, K):
33 """
34 Computes the Sturm sequence of ``f`` in ``F[x]``.
36 Given a univariate, square-free polynomial ``f(x)`` returns the
37 associated Sturm sequence ``f_0(x), ..., f_n(x)`` defined by::
39 f_0(x), f_1(x) = f(x), f'(x)
40 f_n = -rem(f_{n-2}(x), f_{n-1}(x))
42 Examples
43 ========
45 >>> from sympy.polys import ring, QQ
46 >>> R, x = ring("x", QQ)
48 >>> R.dup_sturm(x**3 - 2*x**2 + x - 3)
49 [x**3 - 2*x**2 + x - 3, 3*x**2 - 4*x + 1, 2/9*x + 25/9, -2079/4]
51 References
52 ==========
54 .. [1] [Davenport88]_
56 """
57 if not K.is_Field:
58 raise DomainError("Cannot compute Sturm sequence over %s" % K)
60 f = dup_sqf_part(f, K)
62 sturm = [f, dup_diff(f, 1, K)]
64 while sturm[-1]:
65 s = dup_rem(sturm[-2], sturm[-1], K)
66 sturm.append(dup_neg(s, K))
68 return sturm[:-1]
70def dup_root_upper_bound(f, K):
71 """Compute the LMQ upper bound for the positive roots of `f`;
72 LMQ (Local Max Quadratic) was developed by Akritas-Strzebonski-Vigklas.
74 References
75 ==========
76 .. [1] Alkiviadis G. Akritas: "Linear and Quadratic Complexity Bounds on the
77 Values of the Positive Roots of Polynomials"
78 Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009.
79 """
80 n, P = len(f), []
81 t = n * [K.one]
82 if dup_LC(f, K) < 0:
83 f = dup_neg(f, K)
84 f = list(reversed(f))
86 for i in range(0, n):
87 if f[i] >= 0:
88 continue
90 a, QL = K.log(-f[i], 2), []
92 for j in range(i + 1, n):
94 if f[j] <= 0:
95 continue
97 q = t[j] + a - K.log(f[j], 2)
98 QL.append([q // (j - i), j])
100 if not QL:
101 continue
103 q = min(QL)
105 t[q[1]] = t[q[1]] + 1
107 P.append(q[0])
109 if not P:
110 return None
111 else:
112 return K.get_field()(2)**(max(P) + 1)
114def dup_root_lower_bound(f, K):
115 """Compute the LMQ lower bound for the positive roots of `f`;
116 LMQ (Local Max Quadratic) was developed by Akritas-Strzebonski-Vigklas.
118 References
119 ==========
120 .. [1] Alkiviadis G. Akritas: "Linear and Quadratic Complexity Bounds on the
121 Values of the Positive Roots of Polynomials"
122 Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009.
123 """
124 bound = dup_root_upper_bound(dup_reverse(f), K)
126 if bound is not None:
127 return 1/bound
128 else:
129 return None
131def dup_cauchy_upper_bound(f, K):
132 """
133 Compute the Cauchy upper bound on the absolute value of all roots of f,
134 real or complex.
136 References
137 ==========
138 .. [1] https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Lagrange's_and_Cauchy's_bounds
139 """
140 n = dup_degree(f)
141 if n < 1:
142 raise PolynomialError('Polynomial has no roots.')
144 if K.is_ZZ:
145 L = K.get_field()
146 f, K = dup_convert(f, K, L), L
147 elif not K.is_QQ or K.is_RR or K.is_CC:
148 # We need to compute absolute value, and we are not supporting cases
149 # where this would take us outside the domain (or its quotient field).
150 raise DomainError('Cauchy bound not supported over %s' % K)
151 else:
152 f = f[:]
154 while K.is_zero(f[-1]):
155 f.pop()
156 if len(f) == 1:
157 # Monomial. All roots are zero.
158 return K.zero
160 lc = f[0]
161 return K.one + max(abs(n / lc) for n in f[1:])
163def dup_cauchy_lower_bound(f, K):
164 """Compute the Cauchy lower bound on the absolute value of all non-zero
165 roots of f, real or complex."""
166 g = dup_reverse(f)
167 if len(g) < 2:
168 raise PolynomialError('Polynomial has no non-zero roots.')
169 if K.is_ZZ:
170 K = K.get_field()
171 b = dup_cauchy_upper_bound(g, K)
172 return K.one / b
174def dup_mignotte_sep_bound_squared(f, K):
175 """
176 Return the square of the Mignotte lower bound on separation between
177 distinct roots of f. The square is returned so that the bound lies in
178 K or its quotient field.
180 References
181 ==========
183 .. [1] Mignotte, Maurice. "Some useful bounds." Computer algebra.
184 Springer, Vienna, 1982. 259-263.
185 https://people.dm.unipi.it/gianni/AC-EAG/Mignotte.pdf
186 """
187 n = dup_degree(f)
188 if n < 2:
189 raise PolynomialError('Polynomials of degree < 2 have no distinct roots.')
191 if K.is_ZZ:
192 L = K.get_field()
193 f, K = dup_convert(f, K, L), L
194 elif not K.is_QQ or K.is_RR or K.is_CC:
195 # We need to compute absolute value, and we are not supporting cases
196 # where this would take us outside the domain (or its quotient field).
197 raise DomainError('Mignotte bound not supported over %s' % K)
199 D = dup_discriminant(f, K)
200 l2sq = dup_l2_norm_squared(f, K)
201 return K(3)*K.abs(D) / ( K(n)**(n+1) * l2sq**(n-1) )
203def _mobius_from_interval(I, field):
204 """Convert an open interval to a Mobius transform. """
205 s, t = I
207 a, c = field.numer(s), field.denom(s)
208 b, d = field.numer(t), field.denom(t)
210 return a, b, c, d
212def _mobius_to_interval(M, field):
213 """Convert a Mobius transform to an open interval. """
214 a, b, c, d = M
216 s, t = field(a, c), field(b, d)
218 if s <= t:
219 return (s, t)
220 else:
221 return (t, s)
223def dup_step_refine_real_root(f, M, K, fast=False):
224 """One step of positive real root refinement algorithm. """
225 a, b, c, d = M
227 if a == b and c == d:
228 return f, (a, b, c, d)
230 A = dup_root_lower_bound(f, K)
232 if A is not None:
233 A = K(int(A))
234 else:
235 A = K.zero
237 if fast and A > 16:
238 f = dup_scale(f, A, K)
239 a, c, A = A*a, A*c, K.one
241 if A >= K.one:
242 f = dup_shift(f, A, K)
243 b, d = A*a + b, A*c + d
245 if not dup_eval(f, K.zero, K):
246 return f, (b, b, d, d)
248 f, g = dup_shift(f, K.one, K), f
250 a1, b1, c1, d1 = a, a + b, c, c + d
252 if not dup_eval(f, K.zero, K):
253 return f, (b1, b1, d1, d1)
255 k = dup_sign_variations(f, K)
257 if k == 1:
258 a, b, c, d = a1, b1, c1, d1
259 else:
260 f = dup_shift(dup_reverse(g), K.one, K)
262 if not dup_eval(f, K.zero, K):
263 f = dup_rshift(f, 1, K)
265 a, b, c, d = b, a + b, d, c + d
267 return f, (a, b, c, d)
269def dup_inner_refine_real_root(f, M, K, eps=None, steps=None, disjoint=None, fast=False, mobius=False):
270 """Refine a positive root of `f` given a Mobius transform or an interval. """
271 F = K.get_field()
273 if len(M) == 2:
274 a, b, c, d = _mobius_from_interval(M, F)
275 else:
276 a, b, c, d = M
278 while not c:
279 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c,
280 d), K, fast=fast)
282 if eps is not None and steps is not None:
283 for i in range(0, steps):
284 if abs(F(a, c) - F(b, d)) >= eps:
285 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast)
286 else:
287 break
288 else:
289 if eps is not None:
290 while abs(F(a, c) - F(b, d)) >= eps:
291 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast)
293 if steps is not None:
294 for i in range(0, steps):
295 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast)
297 if disjoint is not None:
298 while True:
299 u, v = _mobius_to_interval((a, b, c, d), F)
301 if v <= disjoint or disjoint <= u:
302 break
303 else:
304 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast)
306 if not mobius:
307 return _mobius_to_interval((a, b, c, d), F)
308 else:
309 return f, (a, b, c, d)
311def dup_outer_refine_real_root(f, s, t, K, eps=None, steps=None, disjoint=None, fast=False):
312 """Refine a positive root of `f` given an interval `(s, t)`. """
313 a, b, c, d = _mobius_from_interval((s, t), K.get_field())
315 f = dup_transform(f, dup_strip([a, b]),
316 dup_strip([c, d]), K)
318 if dup_sign_variations(f, K) != 1:
319 raise RefinementFailed("there should be exactly one root in (%s, %s) interval" % (s, t))
321 return dup_inner_refine_real_root(f, (a, b, c, d), K, eps=eps, steps=steps, disjoint=disjoint, fast=fast)
323def dup_refine_real_root(f, s, t, K, eps=None, steps=None, disjoint=None, fast=False):
324 """Refine real root's approximating interval to the given precision. """
325 if K.is_QQ:
326 (_, f), K = dup_clear_denoms(f, K, convert=True), K.get_ring()
327 elif not K.is_ZZ:
328 raise DomainError("real root refinement not supported over %s" % K)
330 if s == t:
331 return (s, t)
333 if s > t:
334 s, t = t, s
336 negative = False
338 if s < 0:
339 if t <= 0:
340 f, s, t, negative = dup_mirror(f, K), -t, -s, True
341 else:
342 raise ValueError("Cannot refine a real root in (%s, %s)" % (s, t))
344 if negative and disjoint is not None:
345 if disjoint < 0:
346 disjoint = -disjoint
347 else:
348 disjoint = None
350 s, t = dup_outer_refine_real_root(
351 f, s, t, K, eps=eps, steps=steps, disjoint=disjoint, fast=fast)
353 if negative:
354 return (-t, -s)
355 else:
356 return ( s, t)
358def dup_inner_isolate_real_roots(f, K, eps=None, fast=False):
359 """Internal function for isolation positive roots up to given precision.
361 References
362 ==========
363 1. Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative Study of Two Real Root
364 Isolation Methods . Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
365 2. Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the
366 Performance of the Continued Fractions Method Using new Bounds of Positive Roots. Nonlinear
367 Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
368 """
369 a, b, c, d = K.one, K.zero, K.zero, K.one
371 k = dup_sign_variations(f, K)
373 if k == 0:
374 return []
375 if k == 1:
376 roots = [dup_inner_refine_real_root(
377 f, (a, b, c, d), K, eps=eps, fast=fast, mobius=True)]
378 else:
379 roots, stack = [], [(a, b, c, d, f, k)]
381 while stack:
382 a, b, c, d, f, k = stack.pop()
384 A = dup_root_lower_bound(f, K)
386 if A is not None:
387 A = K(int(A))
388 else:
389 A = K.zero
391 if fast and A > 16:
392 f = dup_scale(f, A, K)
393 a, c, A = A*a, A*c, K.one
395 if A >= K.one:
396 f = dup_shift(f, A, K)
397 b, d = A*a + b, A*c + d
399 if not dup_TC(f, K):
400 roots.append((f, (b, b, d, d)))
401 f = dup_rshift(f, 1, K)
403 k = dup_sign_variations(f, K)
405 if k == 0:
406 continue
407 if k == 1:
408 roots.append(dup_inner_refine_real_root(
409 f, (a, b, c, d), K, eps=eps, fast=fast, mobius=True))
410 continue
412 f1 = dup_shift(f, K.one, K)
414 a1, b1, c1, d1, r = a, a + b, c, c + d, 0
416 if not dup_TC(f1, K):
417 roots.append((f1, (b1, b1, d1, d1)))
418 f1, r = dup_rshift(f1, 1, K), 1
420 k1 = dup_sign_variations(f1, K)
421 k2 = k - k1 - r
423 a2, b2, c2, d2 = b, a + b, d, c + d
425 if k2 > 1:
426 f2 = dup_shift(dup_reverse(f), K.one, K)
428 if not dup_TC(f2, K):
429 f2 = dup_rshift(f2, 1, K)
431 k2 = dup_sign_variations(f2, K)
432 else:
433 f2 = None
435 if k1 < k2:
436 a1, a2, b1, b2 = a2, a1, b2, b1
437 c1, c2, d1, d2 = c2, c1, d2, d1
438 f1, f2, k1, k2 = f2, f1, k2, k1
440 if not k1:
441 continue
443 if f1 is None:
444 f1 = dup_shift(dup_reverse(f), K.one, K)
446 if not dup_TC(f1, K):
447 f1 = dup_rshift(f1, 1, K)
449 if k1 == 1:
450 roots.append(dup_inner_refine_real_root(
451 f1, (a1, b1, c1, d1), K, eps=eps, fast=fast, mobius=True))
452 else:
453 stack.append((a1, b1, c1, d1, f1, k1))
455 if not k2:
456 continue
458 if f2 is None:
459 f2 = dup_shift(dup_reverse(f), K.one, K)
461 if not dup_TC(f2, K):
462 f2 = dup_rshift(f2, 1, K)
464 if k2 == 1:
465 roots.append(dup_inner_refine_real_root(
466 f2, (a2, b2, c2, d2), K, eps=eps, fast=fast, mobius=True))
467 else:
468 stack.append((a2, b2, c2, d2, f2, k2))
470 return roots
472def _discard_if_outside_interval(f, M, inf, sup, K, negative, fast, mobius):
473 """Discard an isolating interval if outside ``(inf, sup)``. """
474 F = K.get_field()
476 while True:
477 u, v = _mobius_to_interval(M, F)
479 if negative:
480 u, v = -v, -u
482 if (inf is None or u >= inf) and (sup is None or v <= sup):
483 if not mobius:
484 return u, v
485 else:
486 return f, M
487 elif (sup is not None and u > sup) or (inf is not None and v < inf):
488 return None
489 else:
490 f, M = dup_step_refine_real_root(f, M, K, fast=fast)
492def dup_inner_isolate_positive_roots(f, K, eps=None, inf=None, sup=None, fast=False, mobius=False):
493 """Iteratively compute disjoint positive root isolation intervals. """
494 if sup is not None and sup < 0:
495 return []
497 roots = dup_inner_isolate_real_roots(f, K, eps=eps, fast=fast)
499 F, results = K.get_field(), []
501 if inf is not None or sup is not None:
502 for f, M in roots:
503 result = _discard_if_outside_interval(f, M, inf, sup, K, False, fast, mobius)
505 if result is not None:
506 results.append(result)
507 elif not mobius:
508 for f, M in roots:
509 u, v = _mobius_to_interval(M, F)
510 results.append((u, v))
511 else:
512 results = roots
514 return results
516def dup_inner_isolate_negative_roots(f, K, inf=None, sup=None, eps=None, fast=False, mobius=False):
517 """Iteratively compute disjoint negative root isolation intervals. """
518 if inf is not None and inf >= 0:
519 return []
521 roots = dup_inner_isolate_real_roots(dup_mirror(f, K), K, eps=eps, fast=fast)
523 F, results = K.get_field(), []
525 if inf is not None or sup is not None:
526 for f, M in roots:
527 result = _discard_if_outside_interval(f, M, inf, sup, K, True, fast, mobius)
529 if result is not None:
530 results.append(result)
531 elif not mobius:
532 for f, M in roots:
533 u, v = _mobius_to_interval(M, F)
534 results.append((-v, -u))
535 else:
536 results = roots
538 return results
540def _isolate_zero(f, K, inf, sup, basis=False, sqf=False):
541 """Handle special case of CF algorithm when ``f`` is homogeneous. """
542 j, f = dup_terms_gcd(f, K)
544 if j > 0:
545 F = K.get_field()
547 if (inf is None or inf <= 0) and (sup is None or 0 <= sup):
548 if not sqf:
549 if not basis:
550 return [((F.zero, F.zero), j)], f
551 else:
552 return [((F.zero, F.zero), j, [K.one, K.zero])], f
553 else:
554 return [(F.zero, F.zero)], f
556 return [], f
558def dup_isolate_real_roots_sqf(f, K, eps=None, inf=None, sup=None, fast=False, blackbox=False):
559 """Isolate real roots of a square-free polynomial using the Vincent-Akritas-Strzebonski (VAS) CF approach.
561 References
562 ==========
563 .. [1] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative
564 Study of Two Real Root Isolation Methods. Nonlinear Analysis:
565 Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
566 .. [2] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S.
567 Vigklas: Improving the Performance of the Continued Fractions
568 Method Using New Bounds of Positive Roots. Nonlinear Analysis:
569 Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
571 """
572 if K.is_QQ:
573 (_, f), K = dup_clear_denoms(f, K, convert=True), K.get_ring()
574 elif not K.is_ZZ:
575 raise DomainError("isolation of real roots not supported over %s" % K)
577 if dup_degree(f) <= 0:
578 return []
580 I_zero, f = _isolate_zero(f, K, inf, sup, basis=False, sqf=True)
582 I_neg = dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast)
583 I_pos = dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast)
585 roots = sorted(I_neg + I_zero + I_pos)
587 if not blackbox:
588 return roots
589 else:
590 return [ RealInterval((a, b), f, K) for (a, b) in roots ]
592def dup_isolate_real_roots(f, K, eps=None, inf=None, sup=None, basis=False, fast=False):
593 """Isolate real roots using Vincent-Akritas-Strzebonski (VAS) continued fractions approach.
595 References
596 ==========
598 .. [1] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative
599 Study of Two Real Root Isolation Methods. Nonlinear Analysis:
600 Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
601 .. [2] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S.
602 Vigklas: Improving the Performance of the Continued Fractions
603 Method Using New Bounds of Positive Roots.
604 Nonlinear Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
606 """
607 if K.is_QQ:
608 (_, f), K = dup_clear_denoms(f, K, convert=True), K.get_ring()
609 elif not K.is_ZZ:
610 raise DomainError("isolation of real roots not supported over %s" % K)
612 if dup_degree(f) <= 0:
613 return []
615 I_zero, f = _isolate_zero(f, K, inf, sup, basis=basis, sqf=False)
617 _, factors = dup_sqf_list(f, K)
619 if len(factors) == 1:
620 ((f, k),) = factors
622 I_neg = dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast)
623 I_pos = dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast)
625 I_neg = [ ((u, v), k) for u, v in I_neg ]
626 I_pos = [ ((u, v), k) for u, v in I_pos ]
627 else:
628 I_neg, I_pos = _real_isolate_and_disjoin(factors, K,
629 eps=eps, inf=inf, sup=sup, basis=basis, fast=fast)
631 return sorted(I_neg + I_zero + I_pos)
633def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=False, basis=False, fast=False):
634 """Isolate real roots of a list of square-free polynomial using Vincent-Akritas-Strzebonski (VAS) CF approach.
636 References
637 ==========
639 .. [1] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative
640 Study of Two Real Root Isolation Methods. Nonlinear Analysis:
641 Modelling and Control, Vol. 10, No. 4, 297-304, 2005.
642 .. [2] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S.
643 Vigklas: Improving the Performance of the Continued Fractions
644 Method Using New Bounds of Positive Roots.
645 Nonlinear Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008.
647 """
648 if K.is_QQ:
649 K, F, polys = K.get_ring(), K, polys[:]
651 for i, p in enumerate(polys):
652 polys[i] = dup_clear_denoms(p, F, K, convert=True)[1]
653 elif not K.is_ZZ:
654 raise DomainError("isolation of real roots not supported over %s" % K)
656 zeros, factors_dict = False, {}
658 if (inf is None or inf <= 0) and (sup is None or 0 <= sup):
659 zeros, zero_indices = True, {}
661 for i, p in enumerate(polys):
662 j, p = dup_terms_gcd(p, K)
664 if zeros and j > 0:
665 zero_indices[i] = j
667 for f, k in dup_factor_list(p, K)[1]:
668 f = tuple(f)
670 if f not in factors_dict:
671 factors_dict[f] = {i: k}
672 else:
673 factors_dict[f][i] = k
675 factors_list = []
677 for f, indices in factors_dict.items():
678 factors_list.append((list(f), indices))
680 I_neg, I_pos = _real_isolate_and_disjoin(factors_list, K, eps=eps,
681 inf=inf, sup=sup, strict=strict, basis=basis, fast=fast)
683 F = K.get_field()
685 if not zeros or not zero_indices:
686 I_zero = []
687 else:
688 if not basis:
689 I_zero = [((F.zero, F.zero), zero_indices)]
690 else:
691 I_zero = [((F.zero, F.zero), zero_indices, [K.one, K.zero])]
693 return sorted(I_neg + I_zero + I_pos)
695def _disjoint_p(M, N, strict=False):
696 """Check if Mobius transforms define disjoint intervals. """
697 a1, b1, c1, d1 = M
698 a2, b2, c2, d2 = N
700 a1d1, b1c1 = a1*d1, b1*c1
701 a2d2, b2c2 = a2*d2, b2*c2
703 if a1d1 == b1c1 and a2d2 == b2c2:
704 return True
706 if a1d1 > b1c1:
707 a1, c1, b1, d1 = b1, d1, a1, c1
709 if a2d2 > b2c2:
710 a2, c2, b2, d2 = b2, d2, a2, c2
712 if not strict:
713 return a2*d1 >= c2*b1 or b2*c1 <= d2*a1
714 else:
715 return a2*d1 > c2*b1 or b2*c1 < d2*a1
717def _real_isolate_and_disjoin(factors, K, eps=None, inf=None, sup=None, strict=False, basis=False, fast=False):
718 """Isolate real roots of a list of polynomials and disjoin intervals. """
719 I_pos, I_neg = [], []
721 for i, (f, k) in enumerate(factors):
722 for F, M in dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast, mobius=True):
723 I_pos.append((F, M, k, f))
725 for G, N in dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast, mobius=True):
726 I_neg.append((G, N, k, f))
728 for i, (f, M, k, F) in enumerate(I_pos):
729 for j, (g, N, m, G) in enumerate(I_pos[i + 1:]):
730 while not _disjoint_p(M, N, strict=strict):
731 f, M = dup_inner_refine_real_root(f, M, K, steps=1, fast=fast, mobius=True)
732 g, N = dup_inner_refine_real_root(g, N, K, steps=1, fast=fast, mobius=True)
734 I_pos[i + j + 1] = (g, N, m, G)
736 I_pos[i] = (f, M, k, F)
738 for i, (f, M, k, F) in enumerate(I_neg):
739 for j, (g, N, m, G) in enumerate(I_neg[i + 1:]):
740 while not _disjoint_p(M, N, strict=strict):
741 f, M = dup_inner_refine_real_root(f, M, K, steps=1, fast=fast, mobius=True)
742 g, N = dup_inner_refine_real_root(g, N, K, steps=1, fast=fast, mobius=True)
744 I_neg[i + j + 1] = (g, N, m, G)
746 I_neg[i] = (f, M, k, F)
748 if strict:
749 for i, (f, M, k, F) in enumerate(I_neg):
750 if not M[0]:
751 while not M[0]:
752 f, M = dup_inner_refine_real_root(f, M, K, steps=1, fast=fast, mobius=True)
754 I_neg[i] = (f, M, k, F)
755 break
757 for j, (g, N, m, G) in enumerate(I_pos):
758 if not N[0]:
759 while not N[0]:
760 g, N = dup_inner_refine_real_root(g, N, K, steps=1, fast=fast, mobius=True)
762 I_pos[j] = (g, N, m, G)
763 break
765 field = K.get_field()
767 I_neg = [ (_mobius_to_interval(M, field), k, f) for (_, M, k, f) in I_neg ]
768 I_pos = [ (_mobius_to_interval(M, field), k, f) for (_, M, k, f) in I_pos ]
770 if not basis:
771 I_neg = [ ((-v, -u), k) for ((u, v), k, _) in I_neg ]
772 I_pos = [ (( u, v), k) for ((u, v), k, _) in I_pos ]
773 else:
774 I_neg = [ ((-v, -u), k, f) for ((u, v), k, f) in I_neg ]
775 I_pos = [ (( u, v), k, f) for ((u, v), k, f) in I_pos ]
777 return I_neg, I_pos
779def dup_count_real_roots(f, K, inf=None, sup=None):
780 """Returns the number of distinct real roots of ``f`` in ``[inf, sup]``. """
781 if dup_degree(f) <= 0:
782 return 0
784 if not K.is_Field:
785 R, K = K, K.get_field()
786 f = dup_convert(f, R, K)
788 sturm = dup_sturm(f, K)
790 if inf is None:
791 signs_inf = dup_sign_variations([ dup_LC(s, K)*(-1)**dup_degree(s) for s in sturm ], K)
792 else:
793 signs_inf = dup_sign_variations([ dup_eval(s, inf, K) for s in sturm ], K)
795 if sup is None:
796 signs_sup = dup_sign_variations([ dup_LC(s, K) for s in sturm ], K)
797 else:
798 signs_sup = dup_sign_variations([ dup_eval(s, sup, K) for s in sturm ], K)
800 count = abs(signs_inf - signs_sup)
802 if inf is not None and not dup_eval(f, inf, K):
803 count += 1
805 return count
807OO = 'OO' # Origin of (re, im) coordinate system
809Q1 = 'Q1' # Quadrant #1 (++): re > 0 and im > 0
810Q2 = 'Q2' # Quadrant #2 (-+): re < 0 and im > 0
811Q3 = 'Q3' # Quadrant #3 (--): re < 0 and im < 0
812Q4 = 'Q4' # Quadrant #4 (+-): re > 0 and im < 0
814A1 = 'A1' # Axis #1 (+0): re > 0 and im = 0
815A2 = 'A2' # Axis #2 (0+): re = 0 and im > 0
816A3 = 'A3' # Axis #3 (-0): re < 0 and im = 0
817A4 = 'A4' # Axis #4 (0-): re = 0 and im < 0
819_rules_simple = {
820 # Q --> Q (same) => no change
821 (Q1, Q1): 0,
822 (Q2, Q2): 0,
823 (Q3, Q3): 0,
824 (Q4, Q4): 0,
826 # A -- CCW --> Q => +1/4 (CCW)
827 (A1, Q1): 1,
828 (A2, Q2): 1,
829 (A3, Q3): 1,
830 (A4, Q4): 1,
832 # A -- CW --> Q => -1/4 (CCW)
833 (A1, Q4): 2,
834 (A2, Q1): 2,
835 (A3, Q2): 2,
836 (A4, Q3): 2,
838 # Q -- CCW --> A => +1/4 (CCW)
839 (Q1, A2): 3,
840 (Q2, A3): 3,
841 (Q3, A4): 3,
842 (Q4, A1): 3,
844 # Q -- CW --> A => -1/4 (CCW)
845 (Q1, A1): 4,
846 (Q2, A2): 4,
847 (Q3, A3): 4,
848 (Q4, A4): 4,
850 # Q -- CCW --> Q => +1/2 (CCW)
851 (Q1, Q2): +5,
852 (Q2, Q3): +5,
853 (Q3, Q4): +5,
854 (Q4, Q1): +5,
856 # Q -- CW --> Q => -1/2 (CW)
857 (Q1, Q4): -5,
858 (Q2, Q1): -5,
859 (Q3, Q2): -5,
860 (Q4, Q3): -5,
861}
863_rules_ambiguous = {
864 # A -- CCW --> Q => { +1/4 (CCW), -9/4 (CW) }
865 (A1, OO, Q1): -1,
866 (A2, OO, Q2): -1,
867 (A3, OO, Q3): -1,
868 (A4, OO, Q4): -1,
870 # A -- CW --> Q => { -1/4 (CCW), +7/4 (CW) }
871 (A1, OO, Q4): -2,
872 (A2, OO, Q1): -2,
873 (A3, OO, Q2): -2,
874 (A4, OO, Q3): -2,
876 # Q -- CCW --> A => { +1/4 (CCW), -9/4 (CW) }
877 (Q1, OO, A2): -3,
878 (Q2, OO, A3): -3,
879 (Q3, OO, A4): -3,
880 (Q4, OO, A1): -3,
882 # Q -- CW --> A => { -1/4 (CCW), +7/4 (CW) }
883 (Q1, OO, A1): -4,
884 (Q2, OO, A2): -4,
885 (Q3, OO, A3): -4,
886 (Q4, OO, A4): -4,
888 # A -- OO --> A => { +1 (CCW), -1 (CW) }
889 (A1, A3): 7,
890 (A2, A4): 7,
891 (A3, A1): 7,
892 (A4, A2): 7,
894 (A1, OO, A3): 7,
895 (A2, OO, A4): 7,
896 (A3, OO, A1): 7,
897 (A4, OO, A2): 7,
899 # Q -- DIA --> Q => { +1 (CCW), -1 (CW) }
900 (Q1, Q3): 8,
901 (Q2, Q4): 8,
902 (Q3, Q1): 8,
903 (Q4, Q2): 8,
905 (Q1, OO, Q3): 8,
906 (Q2, OO, Q4): 8,
907 (Q3, OO, Q1): 8,
908 (Q4, OO, Q2): 8,
910 # A --- R ---> A => { +1/2 (CCW), -3/2 (CW) }
911 (A1, A2): 9,
912 (A2, A3): 9,
913 (A3, A4): 9,
914 (A4, A1): 9,
916 (A1, OO, A2): 9,
917 (A2, OO, A3): 9,
918 (A3, OO, A4): 9,
919 (A4, OO, A1): 9,
921 # A --- L ---> A => { +3/2 (CCW), -1/2 (CW) }
922 (A1, A4): 10,
923 (A2, A1): 10,
924 (A3, A2): 10,
925 (A4, A3): 10,
927 (A1, OO, A4): 10,
928 (A2, OO, A1): 10,
929 (A3, OO, A2): 10,
930 (A4, OO, A3): 10,
932 # Q --- 1 ---> A => { +3/4 (CCW), -5/4 (CW) }
933 (Q1, A3): 11,
934 (Q2, A4): 11,
935 (Q3, A1): 11,
936 (Q4, A2): 11,
938 (Q1, OO, A3): 11,
939 (Q2, OO, A4): 11,
940 (Q3, OO, A1): 11,
941 (Q4, OO, A2): 11,
943 # Q --- 2 ---> A => { +5/4 (CCW), -3/4 (CW) }
944 (Q1, A4): 12,
945 (Q2, A1): 12,
946 (Q3, A2): 12,
947 (Q4, A3): 12,
949 (Q1, OO, A4): 12,
950 (Q2, OO, A1): 12,
951 (Q3, OO, A2): 12,
952 (Q4, OO, A3): 12,
954 # A --- 1 ---> Q => { +5/4 (CCW), -3/4 (CW) }
955 (A1, Q3): 13,
956 (A2, Q4): 13,
957 (A3, Q1): 13,
958 (A4, Q2): 13,
960 (A1, OO, Q3): 13,
961 (A2, OO, Q4): 13,
962 (A3, OO, Q1): 13,
963 (A4, OO, Q2): 13,
965 # A --- 2 ---> Q => { +3/4 (CCW), -5/4 (CW) }
966 (A1, Q2): 14,
967 (A2, Q3): 14,
968 (A3, Q4): 14,
969 (A4, Q1): 14,
971 (A1, OO, Q2): 14,
972 (A2, OO, Q3): 14,
973 (A3, OO, Q4): 14,
974 (A4, OO, Q1): 14,
976 # Q --> OO --> Q => { +1/2 (CCW), -3/2 (CW) }
977 (Q1, OO, Q2): 15,
978 (Q2, OO, Q3): 15,
979 (Q3, OO, Q4): 15,
980 (Q4, OO, Q1): 15,
982 # Q --> OO --> Q => { +3/2 (CCW), -1/2 (CW) }
983 (Q1, OO, Q4): 16,
984 (Q2, OO, Q1): 16,
985 (Q3, OO, Q2): 16,
986 (Q4, OO, Q3): 16,
988 # A --> OO --> A => { +2 (CCW), 0 (CW) }
989 (A1, OO, A1): 17,
990 (A2, OO, A2): 17,
991 (A3, OO, A3): 17,
992 (A4, OO, A4): 17,
994 # Q --> OO --> Q => { +2 (CCW), 0 (CW) }
995 (Q1, OO, Q1): 18,
996 (Q2, OO, Q2): 18,
997 (Q3, OO, Q3): 18,
998 (Q4, OO, Q4): 18,
999}
1001_values = {
1002 0: [( 0, 1)],
1003 1: [(+1, 4)],
1004 2: [(-1, 4)],
1005 3: [(+1, 4)],
1006 4: [(-1, 4)],
1007 -1: [(+9, 4), (+1, 4)],
1008 -2: [(+7, 4), (-1, 4)],
1009 -3: [(+9, 4), (+1, 4)],
1010 -4: [(+7, 4), (-1, 4)],
1011 +5: [(+1, 2)],
1012 -5: [(-1, 2)],
1013 7: [(+1, 1), (-1, 1)],
1014 8: [(+1, 1), (-1, 1)],
1015 9: [(+1, 2), (-3, 2)],
1016 10: [(+3, 2), (-1, 2)],
1017 11: [(+3, 4), (-5, 4)],
1018 12: [(+5, 4), (-3, 4)],
1019 13: [(+5, 4), (-3, 4)],
1020 14: [(+3, 4), (-5, 4)],
1021 15: [(+1, 2), (-3, 2)],
1022 16: [(+3, 2), (-1, 2)],
1023 17: [(+2, 1), ( 0, 1)],
1024 18: [(+2, 1), ( 0, 1)],
1025}
1027def _classify_point(re, im):
1028 """Return the half-axis (or origin) on which (re, im) point is located. """
1029 if not re and not im:
1030 return OO
1032 if not re:
1033 if im > 0:
1034 return A2
1035 else:
1036 return A4
1037 elif not im:
1038 if re > 0:
1039 return A1
1040 else:
1041 return A3
1043def _intervals_to_quadrants(intervals, f1, f2, s, t, F):
1044 """Generate a sequence of extended quadrants from a list of critical points. """
1045 if not intervals:
1046 return []
1048 Q = []
1050 if not f1:
1051 (a, b), _, _ = intervals[0]
1053 if a == b == s:
1054 if len(intervals) == 1:
1055 if dup_eval(f2, t, F) > 0:
1056 return [OO, A2]
1057 else:
1058 return [OO, A4]
1059 else:
1060 (a, _), _, _ = intervals[1]
1062 if dup_eval(f2, (s + a)/2, F) > 0:
1063 Q.extend([OO, A2])
1064 f2_sgn = +1
1065 else:
1066 Q.extend([OO, A4])
1067 f2_sgn = -1
1069 intervals = intervals[1:]
1070 else:
1071 if dup_eval(f2, s, F) > 0:
1072 Q.append(A2)
1073 f2_sgn = +1
1074 else:
1075 Q.append(A4)
1076 f2_sgn = -1
1078 for (a, _), indices, _ in intervals:
1079 Q.append(OO)
1081 if indices[1] % 2 == 1:
1082 f2_sgn = -f2_sgn
1084 if a != t:
1085 if f2_sgn > 0:
1086 Q.append(A2)
1087 else:
1088 Q.append(A4)
1090 return Q
1092 if not f2:
1093 (a, b), _, _ = intervals[0]
1095 if a == b == s:
1096 if len(intervals) == 1:
1097 if dup_eval(f1, t, F) > 0:
1098 return [OO, A1]
1099 else:
1100 return [OO, A3]
1101 else:
1102 (a, _), _, _ = intervals[1]
1104 if dup_eval(f1, (s + a)/2, F) > 0:
1105 Q.extend([OO, A1])
1106 f1_sgn = +1
1107 else:
1108 Q.extend([OO, A3])
1109 f1_sgn = -1
1111 intervals = intervals[1:]
1112 else:
1113 if dup_eval(f1, s, F) > 0:
1114 Q.append(A1)
1115 f1_sgn = +1
1116 else:
1117 Q.append(A3)
1118 f1_sgn = -1
1120 for (a, _), indices, _ in intervals:
1121 Q.append(OO)
1123 if indices[0] % 2 == 1:
1124 f1_sgn = -f1_sgn
1126 if a != t:
1127 if f1_sgn > 0:
1128 Q.append(A1)
1129 else:
1130 Q.append(A3)
1132 return Q
1134 re = dup_eval(f1, s, F)
1135 im = dup_eval(f2, s, F)
1137 if not re or not im:
1138 Q.append(_classify_point(re, im))
1140 if len(intervals) == 1:
1141 re = dup_eval(f1, t, F)
1142 im = dup_eval(f2, t, F)
1143 else:
1144 (a, _), _, _ = intervals[1]
1146 re = dup_eval(f1, (s + a)/2, F)
1147 im = dup_eval(f2, (s + a)/2, F)
1149 intervals = intervals[1:]
1151 if re > 0:
1152 f1_sgn = +1
1153 else:
1154 f1_sgn = -1
1156 if im > 0:
1157 f2_sgn = +1
1158 else:
1159 f2_sgn = -1
1161 sgn = {
1162 (+1, +1): Q1,
1163 (-1, +1): Q2,
1164 (-1, -1): Q3,
1165 (+1, -1): Q4,
1166 }
1168 Q.append(sgn[(f1_sgn, f2_sgn)])
1170 for (a, b), indices, _ in intervals:
1171 if a == b:
1172 re = dup_eval(f1, a, F)
1173 im = dup_eval(f2, a, F)
1175 cls = _classify_point(re, im)
1177 if cls is not None:
1178 Q.append(cls)
1180 if 0 in indices:
1181 if indices[0] % 2 == 1:
1182 f1_sgn = -f1_sgn
1184 if 1 in indices:
1185 if indices[1] % 2 == 1:
1186 f2_sgn = -f2_sgn
1188 if not (a == b and b == t):
1189 Q.append(sgn[(f1_sgn, f2_sgn)])
1191 return Q
1193def _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=None):
1194 """Transform sequences of quadrants to a sequence of rules. """
1195 if exclude is True:
1196 edges = [1, 1, 0, 0]
1198 corners = {
1199 (0, 1): 1,
1200 (1, 2): 1,
1201 (2, 3): 0,
1202 (3, 0): 1,
1203 }
1204 else:
1205 edges = [0, 0, 0, 0]
1207 corners = {
1208 (0, 1): 0,
1209 (1, 2): 0,
1210 (2, 3): 0,
1211 (3, 0): 0,
1212 }
1214 if exclude is not None and exclude is not True:
1215 exclude = set(exclude)
1217 for i, edge in enumerate(['S', 'E', 'N', 'W']):
1218 if edge in exclude:
1219 edges[i] = 1
1221 for i, corner in enumerate(['SW', 'SE', 'NE', 'NW']):
1222 if corner in exclude:
1223 corners[((i - 1) % 4, i)] = 1
1225 QQ, rules = [Q_L1, Q_L2, Q_L3, Q_L4], []
1227 for i, Q in enumerate(QQ):
1228 if not Q:
1229 continue
1231 if Q[-1] == OO:
1232 Q = Q[:-1]
1234 if Q[0] == OO:
1235 j, Q = (i - 1) % 4, Q[1:]
1236 qq = (QQ[j][-2], OO, Q[0])
1238 if qq in _rules_ambiguous:
1239 rules.append((_rules_ambiguous[qq], corners[(j, i)]))
1240 else:
1241 raise NotImplementedError("3 element rule (corner): " + str(qq))
1243 q1, k = Q[0], 1
1245 while k < len(Q):
1246 q2, k = Q[k], k + 1
1248 if q2 != OO:
1249 qq = (q1, q2)
1251 if qq in _rules_simple:
1252 rules.append((_rules_simple[qq], 0))
1253 elif qq in _rules_ambiguous:
1254 rules.append((_rules_ambiguous[qq], edges[i]))
1255 else:
1256 raise NotImplementedError("2 element rule (inside): " + str(qq))
1257 else:
1258 qq, k = (q1, q2, Q[k]), k + 1
1260 if qq in _rules_ambiguous:
1261 rules.append((_rules_ambiguous[qq], edges[i]))
1262 else:
1263 raise NotImplementedError("3 element rule (edge): " + str(qq))
1265 q1 = qq[-1]
1267 return rules
1269def _reverse_intervals(intervals):
1270 """Reverse intervals for traversal from right to left and from top to bottom. """
1271 return [ ((b, a), indices, f) for (a, b), indices, f in reversed(intervals) ]
1273def _winding_number(T, field):
1274 """Compute the winding number of the input polynomial, i.e. the number of roots. """
1275 return int(sum([ field(*_values[t][i]) for t, i in T ]) / field(2))
1277def dup_count_complex_roots(f, K, inf=None, sup=None, exclude=None):
1278 """Count all roots in [u + v*I, s + t*I] rectangle using Collins-Krandick algorithm. """
1279 if not K.is_ZZ and not K.is_QQ:
1280 raise DomainError("complex root counting is not supported over %s" % K)
1282 if K.is_ZZ:
1283 R, F = K, K.get_field()
1284 else:
1285 R, F = K.get_ring(), K
1287 f = dup_convert(f, K, F)
1289 if inf is None or sup is None:
1290 _, lc = dup_degree(f), abs(dup_LC(f, F))
1291 B = 2*max([ F.quo(abs(c), lc) for c in f ])
1293 if inf is None:
1294 (u, v) = (-B, -B)
1295 else:
1296 (u, v) = inf
1298 if sup is None:
1299 (s, t) = (+B, +B)
1300 else:
1301 (s, t) = sup
1303 f1, f2 = dup_real_imag(f, F)
1305 f1L1F = dmp_eval_in(f1, v, 1, 1, F)
1306 f2L1F = dmp_eval_in(f2, v, 1, 1, F)
1308 _, f1L1R = dup_clear_denoms(f1L1F, F, R, convert=True)
1309 _, f2L1R = dup_clear_denoms(f2L1F, F, R, convert=True)
1311 f1L2F = dmp_eval_in(f1, s, 0, 1, F)
1312 f2L2F = dmp_eval_in(f2, s, 0, 1, F)
1314 _, f1L2R = dup_clear_denoms(f1L2F, F, R, convert=True)
1315 _, f2L2R = dup_clear_denoms(f2L2F, F, R, convert=True)
1317 f1L3F = dmp_eval_in(f1, t, 1, 1, F)
1318 f2L3F = dmp_eval_in(f2, t, 1, 1, F)
1320 _, f1L3R = dup_clear_denoms(f1L3F, F, R, convert=True)
1321 _, f2L3R = dup_clear_denoms(f2L3F, F, R, convert=True)
1323 f1L4F = dmp_eval_in(f1, u, 0, 1, F)
1324 f2L4F = dmp_eval_in(f2, u, 0, 1, F)
1326 _, f1L4R = dup_clear_denoms(f1L4F, F, R, convert=True)
1327 _, f2L4R = dup_clear_denoms(f2L4F, F, R, convert=True)
1329 S_L1 = [f1L1R, f2L1R]
1330 S_L2 = [f1L2R, f2L2R]
1331 S_L3 = [f1L3R, f2L3R]
1332 S_L4 = [f1L4R, f2L4R]
1334 I_L1 = dup_isolate_real_roots_list(S_L1, R, inf=u, sup=s, fast=True, basis=True, strict=True)
1335 I_L2 = dup_isolate_real_roots_list(S_L2, R, inf=v, sup=t, fast=True, basis=True, strict=True)
1336 I_L3 = dup_isolate_real_roots_list(S_L3, R, inf=u, sup=s, fast=True, basis=True, strict=True)
1337 I_L4 = dup_isolate_real_roots_list(S_L4, R, inf=v, sup=t, fast=True, basis=True, strict=True)
1339 I_L3 = _reverse_intervals(I_L3)
1340 I_L4 = _reverse_intervals(I_L4)
1342 Q_L1 = _intervals_to_quadrants(I_L1, f1L1F, f2L1F, u, s, F)
1343 Q_L2 = _intervals_to_quadrants(I_L2, f1L2F, f2L2F, v, t, F)
1344 Q_L3 = _intervals_to_quadrants(I_L3, f1L3F, f2L3F, s, u, F)
1345 Q_L4 = _intervals_to_quadrants(I_L4, f1L4F, f2L4F, t, v, F)
1347 T = _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=exclude)
1349 return _winding_number(T, F)
1351def _vertical_bisection(N, a, b, I, Q, F1, F2, f1, f2, F):
1352 """Vertical bisection step in Collins-Krandick root isolation algorithm. """
1353 (u, v), (s, t) = a, b
1355 I_L1, I_L2, I_L3, I_L4 = I
1356 Q_L1, Q_L2, Q_L3, Q_L4 = Q
1358 f1L1F, f1L2F, f1L3F, f1L4F = F1
1359 f2L1F, f2L2F, f2L3F, f2L4F = F2
1361 x = (u + s) / 2
1363 f1V = dmp_eval_in(f1, x, 0, 1, F)
1364 f2V = dmp_eval_in(f2, x, 0, 1, F)
1366 I_V = dup_isolate_real_roots_list([f1V, f2V], F, inf=v, sup=t, fast=True, strict=True, basis=True)
1368 I_L1_L, I_L1_R = [], []
1369 I_L2_L, I_L2_R = I_V, I_L2
1370 I_L3_L, I_L3_R = [], []
1371 I_L4_L, I_L4_R = I_L4, _reverse_intervals(I_V)
1373 for I in I_L1:
1374 (a, b), indices, h = I
1376 if a == b:
1377 if a == x:
1378 I_L1_L.append(I)
1379 I_L1_R.append(I)
1380 elif a < x:
1381 I_L1_L.append(I)
1382 else:
1383 I_L1_R.append(I)
1384 else:
1385 if b <= x:
1386 I_L1_L.append(I)
1387 elif a >= x:
1388 I_L1_R.append(I)
1389 else:
1390 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=x, fast=True)
1392 if b <= x:
1393 I_L1_L.append(((a, b), indices, h))
1394 if a >= x:
1395 I_L1_R.append(((a, b), indices, h))
1397 for I in I_L3:
1398 (b, a), indices, h = I
1400 if a == b:
1401 if a == x:
1402 I_L3_L.append(I)
1403 I_L3_R.append(I)
1404 elif a < x:
1405 I_L3_L.append(I)
1406 else:
1407 I_L3_R.append(I)
1408 else:
1409 if b <= x:
1410 I_L3_L.append(I)
1411 elif a >= x:
1412 I_L3_R.append(I)
1413 else:
1414 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=x, fast=True)
1416 if b <= x:
1417 I_L3_L.append(((b, a), indices, h))
1418 if a >= x:
1419 I_L3_R.append(((b, a), indices, h))
1421 Q_L1_L = _intervals_to_quadrants(I_L1_L, f1L1F, f2L1F, u, x, F)
1422 Q_L2_L = _intervals_to_quadrants(I_L2_L, f1V, f2V, v, t, F)
1423 Q_L3_L = _intervals_to_quadrants(I_L3_L, f1L3F, f2L3F, x, u, F)
1424 Q_L4_L = Q_L4
1426 Q_L1_R = _intervals_to_quadrants(I_L1_R, f1L1F, f2L1F, x, s, F)
1427 Q_L2_R = Q_L2
1428 Q_L3_R = _intervals_to_quadrants(I_L3_R, f1L3F, f2L3F, s, x, F)
1429 Q_L4_R = _intervals_to_quadrants(I_L4_R, f1V, f2V, t, v, F)
1431 T_L = _traverse_quadrants(Q_L1_L, Q_L2_L, Q_L3_L, Q_L4_L, exclude=True)
1432 T_R = _traverse_quadrants(Q_L1_R, Q_L2_R, Q_L3_R, Q_L4_R, exclude=True)
1434 N_L = _winding_number(T_L, F)
1435 N_R = _winding_number(T_R, F)
1437 I_L = (I_L1_L, I_L2_L, I_L3_L, I_L4_L)
1438 Q_L = (Q_L1_L, Q_L2_L, Q_L3_L, Q_L4_L)
1440 I_R = (I_L1_R, I_L2_R, I_L3_R, I_L4_R)
1441 Q_R = (Q_L1_R, Q_L2_R, Q_L3_R, Q_L4_R)
1443 F1_L = (f1L1F, f1V, f1L3F, f1L4F)
1444 F2_L = (f2L1F, f2V, f2L3F, f2L4F)
1446 F1_R = (f1L1F, f1L2F, f1L3F, f1V)
1447 F2_R = (f2L1F, f2L2F, f2L3F, f2V)
1449 a, b = (u, v), (x, t)
1450 c, d = (x, v), (s, t)
1452 D_L = (N_L, a, b, I_L, Q_L, F1_L, F2_L)
1453 D_R = (N_R, c, d, I_R, Q_R, F1_R, F2_R)
1455 return D_L, D_R
1457def _horizontal_bisection(N, a, b, I, Q, F1, F2, f1, f2, F):
1458 """Horizontal bisection step in Collins-Krandick root isolation algorithm. """
1459 (u, v), (s, t) = a, b
1461 I_L1, I_L2, I_L3, I_L4 = I
1462 Q_L1, Q_L2, Q_L3, Q_L4 = Q
1464 f1L1F, f1L2F, f1L3F, f1L4F = F1
1465 f2L1F, f2L2F, f2L3F, f2L4F = F2
1467 y = (v + t) / 2
1469 f1H = dmp_eval_in(f1, y, 1, 1, F)
1470 f2H = dmp_eval_in(f2, y, 1, 1, F)
1472 I_H = dup_isolate_real_roots_list([f1H, f2H], F, inf=u, sup=s, fast=True, strict=True, basis=True)
1474 I_L1_B, I_L1_U = I_L1, I_H
1475 I_L2_B, I_L2_U = [], []
1476 I_L3_B, I_L3_U = _reverse_intervals(I_H), I_L3
1477 I_L4_B, I_L4_U = [], []
1479 for I in I_L2:
1480 (a, b), indices, h = I
1482 if a == b:
1483 if a == y:
1484 I_L2_B.append(I)
1485 I_L2_U.append(I)
1486 elif a < y:
1487 I_L2_B.append(I)
1488 else:
1489 I_L2_U.append(I)
1490 else:
1491 if b <= y:
1492 I_L2_B.append(I)
1493 elif a >= y:
1494 I_L2_U.append(I)
1495 else:
1496 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=y, fast=True)
1498 if b <= y:
1499 I_L2_B.append(((a, b), indices, h))
1500 if a >= y:
1501 I_L2_U.append(((a, b), indices, h))
1503 for I in I_L4:
1504 (b, a), indices, h = I
1506 if a == b:
1507 if a == y:
1508 I_L4_B.append(I)
1509 I_L4_U.append(I)
1510 elif a < y:
1511 I_L4_B.append(I)
1512 else:
1513 I_L4_U.append(I)
1514 else:
1515 if b <= y:
1516 I_L4_B.append(I)
1517 elif a >= y:
1518 I_L4_U.append(I)
1519 else:
1520 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=y, fast=True)
1522 if b <= y:
1523 I_L4_B.append(((b, a), indices, h))
1524 if a >= y:
1525 I_L4_U.append(((b, a), indices, h))
1527 Q_L1_B = Q_L1
1528 Q_L2_B = _intervals_to_quadrants(I_L2_B, f1L2F, f2L2F, v, y, F)
1529 Q_L3_B = _intervals_to_quadrants(I_L3_B, f1H, f2H, s, u, F)
1530 Q_L4_B = _intervals_to_quadrants(I_L4_B, f1L4F, f2L4F, y, v, F)
1532 Q_L1_U = _intervals_to_quadrants(I_L1_U, f1H, f2H, u, s, F)
1533 Q_L2_U = _intervals_to_quadrants(I_L2_U, f1L2F, f2L2F, y, t, F)
1534 Q_L3_U = Q_L3
1535 Q_L4_U = _intervals_to_quadrants(I_L4_U, f1L4F, f2L4F, t, y, F)
1537 T_B = _traverse_quadrants(Q_L1_B, Q_L2_B, Q_L3_B, Q_L4_B, exclude=True)
1538 T_U = _traverse_quadrants(Q_L1_U, Q_L2_U, Q_L3_U, Q_L4_U, exclude=True)
1540 N_B = _winding_number(T_B, F)
1541 N_U = _winding_number(T_U, F)
1543 I_B = (I_L1_B, I_L2_B, I_L3_B, I_L4_B)
1544 Q_B = (Q_L1_B, Q_L2_B, Q_L3_B, Q_L4_B)
1546 I_U = (I_L1_U, I_L2_U, I_L3_U, I_L4_U)
1547 Q_U = (Q_L1_U, Q_L2_U, Q_L3_U, Q_L4_U)
1549 F1_B = (f1L1F, f1L2F, f1H, f1L4F)
1550 F2_B = (f2L1F, f2L2F, f2H, f2L4F)
1552 F1_U = (f1H, f1L2F, f1L3F, f1L4F)
1553 F2_U = (f2H, f2L2F, f2L3F, f2L4F)
1555 a, b = (u, v), (s, y)
1556 c, d = (u, y), (s, t)
1558 D_B = (N_B, a, b, I_B, Q_B, F1_B, F2_B)
1559 D_U = (N_U, c, d, I_U, Q_U, F1_U, F2_U)
1561 return D_B, D_U
1563def _depth_first_select(rectangles):
1564 """Find a rectangle of minimum area for bisection. """
1565 min_area, j = None, None
1567 for i, (_, (u, v), (s, t), _, _, _, _) in enumerate(rectangles):
1568 area = (s - u)*(t - v)
1570 if min_area is None or area < min_area:
1571 min_area, j = area, i
1573 return rectangles.pop(j)
1575def _rectangle_small_p(a, b, eps):
1576 """Return ``True`` if the given rectangle is small enough. """
1577 (u, v), (s, t) = a, b
1579 if eps is not None:
1580 return s - u < eps and t - v < eps
1581 else:
1582 return True
1584def dup_isolate_complex_roots_sqf(f, K, eps=None, inf=None, sup=None, blackbox=False):
1585 """Isolate complex roots of a square-free polynomial using Collins-Krandick algorithm. """
1586 if not K.is_ZZ and not K.is_QQ:
1587 raise DomainError("isolation of complex roots is not supported over %s" % K)
1589 if dup_degree(f) <= 0:
1590 return []
1592 if K.is_ZZ:
1593 F = K.get_field()
1594 else:
1595 F = K
1597 f = dup_convert(f, K, F)
1599 lc = abs(dup_LC(f, F))
1600 B = 2*max([ F.quo(abs(c), lc) for c in f ])
1602 (u, v), (s, t) = (-B, F.zero), (B, B)
1604 if inf is not None:
1605 u = inf
1607 if sup is not None:
1608 s = sup
1610 if v < 0 or t <= v or s <= u:
1611 raise ValueError("not a valid complex isolation rectangle")
1613 f1, f2 = dup_real_imag(f, F)
1615 f1L1 = dmp_eval_in(f1, v, 1, 1, F)
1616 f2L1 = dmp_eval_in(f2, v, 1, 1, F)
1618 f1L2 = dmp_eval_in(f1, s, 0, 1, F)
1619 f2L2 = dmp_eval_in(f2, s, 0, 1, F)
1621 f1L3 = dmp_eval_in(f1, t, 1, 1, F)
1622 f2L3 = dmp_eval_in(f2, t, 1, 1, F)
1624 f1L4 = dmp_eval_in(f1, u, 0, 1, F)
1625 f2L4 = dmp_eval_in(f2, u, 0, 1, F)
1627 S_L1 = [f1L1, f2L1]
1628 S_L2 = [f1L2, f2L2]
1629 S_L3 = [f1L3, f2L3]
1630 S_L4 = [f1L4, f2L4]
1632 I_L1 = dup_isolate_real_roots_list(S_L1, F, inf=u, sup=s, fast=True, strict=True, basis=True)
1633 I_L2 = dup_isolate_real_roots_list(S_L2, F, inf=v, sup=t, fast=True, strict=True, basis=True)
1634 I_L3 = dup_isolate_real_roots_list(S_L3, F, inf=u, sup=s, fast=True, strict=True, basis=True)
1635 I_L4 = dup_isolate_real_roots_list(S_L4, F, inf=v, sup=t, fast=True, strict=True, basis=True)
1637 I_L3 = _reverse_intervals(I_L3)
1638 I_L4 = _reverse_intervals(I_L4)
1640 Q_L1 = _intervals_to_quadrants(I_L1, f1L1, f2L1, u, s, F)
1641 Q_L2 = _intervals_to_quadrants(I_L2, f1L2, f2L2, v, t, F)
1642 Q_L3 = _intervals_to_quadrants(I_L3, f1L3, f2L3, s, u, F)
1643 Q_L4 = _intervals_to_quadrants(I_L4, f1L4, f2L4, t, v, F)
1645 T = _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4)
1646 N = _winding_number(T, F)
1648 if not N:
1649 return []
1651 I = (I_L1, I_L2, I_L3, I_L4)
1652 Q = (Q_L1, Q_L2, Q_L3, Q_L4)
1654 F1 = (f1L1, f1L2, f1L3, f1L4)
1655 F2 = (f2L1, f2L2, f2L3, f2L4)
1657 rectangles, roots = [(N, (u, v), (s, t), I, Q, F1, F2)], []
1659 while rectangles:
1660 N, (u, v), (s, t), I, Q, F1, F2 = _depth_first_select(rectangles)
1662 if s - u > t - v:
1663 D_L, D_R = _vertical_bisection(N, (u, v), (s, t), I, Q, F1, F2, f1, f2, F)
1665 N_L, a, b, I_L, Q_L, F1_L, F2_L = D_L
1666 N_R, c, d, I_R, Q_R, F1_R, F2_R = D_R
1668 if N_L >= 1:
1669 if N_L == 1 and _rectangle_small_p(a, b, eps):
1670 roots.append(ComplexInterval(a, b, I_L, Q_L, F1_L, F2_L, f1, f2, F))
1671 else:
1672 rectangles.append(D_L)
1674 if N_R >= 1:
1675 if N_R == 1 and _rectangle_small_p(c, d, eps):
1676 roots.append(ComplexInterval(c, d, I_R, Q_R, F1_R, F2_R, f1, f2, F))
1677 else:
1678 rectangles.append(D_R)
1679 else:
1680 D_B, D_U = _horizontal_bisection(N, (u, v), (s, t), I, Q, F1, F2, f1, f2, F)
1682 N_B, a, b, I_B, Q_B, F1_B, F2_B = D_B
1683 N_U, c, d, I_U, Q_U, F1_U, F2_U = D_U
1685 if N_B >= 1:
1686 if N_B == 1 and _rectangle_small_p(a, b, eps):
1687 roots.append(ComplexInterval(
1688 a, b, I_B, Q_B, F1_B, F2_B, f1, f2, F))
1689 else:
1690 rectangles.append(D_B)
1692 if N_U >= 1:
1693 if N_U == 1 and _rectangle_small_p(c, d, eps):
1694 roots.append(ComplexInterval(
1695 c, d, I_U, Q_U, F1_U, F2_U, f1, f2, F))
1696 else:
1697 rectangles.append(D_U)
1699 _roots, roots = sorted(roots, key=lambda r: (r.ax, r.ay)), []
1701 for root in _roots:
1702 roots.extend([root.conjugate(), root])
1704 if blackbox:
1705 return roots
1706 else:
1707 return [ r.as_tuple() for r in roots ]
1709def dup_isolate_all_roots_sqf(f, K, eps=None, inf=None, sup=None, fast=False, blackbox=False):
1710 """Isolate real and complex roots of a square-free polynomial ``f``. """
1711 return (
1712 dup_isolate_real_roots_sqf( f, K, eps=eps, inf=inf, sup=sup, fast=fast, blackbox=blackbox),
1713 dup_isolate_complex_roots_sqf(f, K, eps=eps, inf=inf, sup=sup, blackbox=blackbox))
1715def dup_isolate_all_roots(f, K, eps=None, inf=None, sup=None, fast=False):
1716 """Isolate real and complex roots of a non-square-free polynomial ``f``. """
1717 if not K.is_ZZ and not K.is_QQ:
1718 raise DomainError("isolation of real and complex roots is not supported over %s" % K)
1720 _, factors = dup_sqf_list(f, K)
1722 if len(factors) == 1:
1723 ((f, k),) = factors
1725 real_part, complex_part = dup_isolate_all_roots_sqf(
1726 f, K, eps=eps, inf=inf, sup=sup, fast=fast)
1728 real_part = [ ((a, b), k) for (a, b) in real_part ]
1729 complex_part = [ ((a, b), k) for (a, b) in complex_part ]
1731 return real_part, complex_part
1732 else:
1733 raise NotImplementedError( "only trivial square-free polynomials are supported")
1735class RealInterval:
1736 """A fully qualified representation of a real isolation interval. """
1738 def __init__(self, data, f, dom):
1739 """Initialize new real interval with complete information. """
1740 if len(data) == 2:
1741 s, t = data
1743 self.neg = False
1745 if s < 0:
1746 if t <= 0:
1747 f, s, t, self.neg = dup_mirror(f, dom), -t, -s, True
1748 else:
1749 raise ValueError("Cannot refine a real root in (%s, %s)" % (s, t))
1751 a, b, c, d = _mobius_from_interval((s, t), dom.get_field())
1753 f = dup_transform(f, dup_strip([a, b]),
1754 dup_strip([c, d]), dom)
1756 self.mobius = a, b, c, d
1757 else:
1758 self.mobius = data[:-1]
1759 self.neg = data[-1]
1761 self.f, self.dom = f, dom
1763 @property
1764 def func(self):
1765 return RealInterval
1767 @property
1768 def args(self):
1769 i = self
1770 return (i.mobius + (i.neg,), i.f, i.dom)
1772 def __eq__(self, other):
1773 if type(other) is not type(self):
1774 return False
1775 return self.args == other.args
1777 @property
1778 def a(self):
1779 """Return the position of the left end. """
1780 field = self.dom.get_field()
1781 a, b, c, d = self.mobius
1783 if not self.neg:
1784 if a*d < b*c:
1785 return field(a, c)
1786 return field(b, d)
1787 else:
1788 if a*d > b*c:
1789 return -field(a, c)
1790 return -field(b, d)
1792 @property
1793 def b(self):
1794 """Return the position of the right end. """
1795 was = self.neg
1796 self.neg = not was
1797 rv = -self.a
1798 self.neg = was
1799 return rv
1801 @property
1802 def dx(self):
1803 """Return width of the real isolating interval. """
1804 return self.b - self.a
1806 @property
1807 def center(self):
1808 """Return the center of the real isolating interval. """
1809 return (self.a + self.b)/2
1811 @property
1812 def max_denom(self):
1813 """Return the largest denominator occurring in either endpoint. """
1814 return max(self.a.denominator, self.b.denominator)
1816 def as_tuple(self):
1817 """Return tuple representation of real isolating interval. """
1818 return (self.a, self.b)
1820 def __repr__(self):
1821 return "(%s, %s)" % (self.a, self.b)
1823 def __contains__(self, item):
1824 """
1825 Say whether a complex number belongs to this real interval.
1827 Parameters
1828 ==========
1830 item : pair (re, im) or number re
1831 Either a pair giving the real and imaginary parts of the number,
1832 or else a real number.
1834 """
1835 if isinstance(item, tuple):
1836 re, im = item
1837 else:
1838 re, im = item, 0
1839 return im == 0 and self.a <= re <= self.b
1841 def is_disjoint(self, other):
1842 """Return ``True`` if two isolation intervals are disjoint. """
1843 if isinstance(other, RealInterval):
1844 return (self.b < other.a or other.b < self.a)
1845 assert isinstance(other, ComplexInterval)
1846 return (self.b < other.ax or other.bx < self.a
1847 or other.ay*other.by > 0)
1849 def _inner_refine(self):
1850 """Internal one step real root refinement procedure. """
1851 if self.mobius is None:
1852 return self
1854 f, mobius = dup_inner_refine_real_root(
1855 self.f, self.mobius, self.dom, steps=1, mobius=True)
1857 return RealInterval(mobius + (self.neg,), f, self.dom)
1859 def refine_disjoint(self, other):
1860 """Refine an isolating interval until it is disjoint with another one. """
1861 expr = self
1862 while not expr.is_disjoint(other):
1863 expr, other = expr._inner_refine(), other._inner_refine()
1865 return expr, other
1867 def refine_size(self, dx):
1868 """Refine an isolating interval until it is of sufficiently small size. """
1869 expr = self
1870 while not (expr.dx < dx):
1871 expr = expr._inner_refine()
1873 return expr
1875 def refine_step(self, steps=1):
1876 """Perform several steps of real root refinement algorithm. """
1877 expr = self
1878 for _ in range(steps):
1879 expr = expr._inner_refine()
1881 return expr
1883 def refine(self):
1884 """Perform one step of real root refinement algorithm. """
1885 return self._inner_refine()
1888class ComplexInterval:
1889 """A fully qualified representation of a complex isolation interval.
1890 The printed form is shown as (ax, bx) x (ay, by) where (ax, ay)
1891 and (bx, by) are the coordinates of the southwest and northeast
1892 corners of the interval's rectangle, respectively.
1894 Examples
1895 ========
1897 >>> from sympy import CRootOf, S
1898 >>> from sympy.abc import x
1899 >>> CRootOf.clear_cache() # for doctest reproducibility
1900 >>> root = CRootOf(x**10 - 2*x + 3, 9)
1901 >>> i = root._get_interval(); i
1902 (3/64, 3/32) x (9/8, 75/64)
1904 The real part of the root lies within the range [0, 3/4] while
1905 the imaginary part lies within the range [9/8, 3/2]:
1907 >>> root.n(3)
1908 0.0766 + 1.14*I
1910 The width of the ranges in the x and y directions on the complex
1911 plane are:
1913 >>> i.dx, i.dy
1914 (3/64, 3/64)
1916 The center of the range is
1918 >>> i.center
1919 (9/128, 147/128)
1921 The northeast coordinate of the rectangle bounding the root in the
1922 complex plane is given by attribute b and the x and y components
1923 are accessed by bx and by:
1925 >>> i.b, i.bx, i.by
1926 ((3/32, 75/64), 3/32, 75/64)
1928 The southwest coordinate is similarly given by i.a
1930 >>> i.a, i.ax, i.ay
1931 ((3/64, 9/8), 3/64, 9/8)
1933 Although the interval prints to show only the real and imaginary
1934 range of the root, all the information of the underlying root
1935 is contained as properties of the interval.
1937 For example, an interval with a nonpositive imaginary range is
1938 considered to be the conjugate. Since the y values of y are in the
1939 range [0, 1/4] it is not the conjugate:
1941 >>> i.conj
1942 False
1944 The conjugate's interval is
1946 >>> ic = i.conjugate(); ic
1947 (3/64, 3/32) x (-75/64, -9/8)
1949 NOTE: the values printed still represent the x and y range
1950 in which the root -- conjugate, in this case -- is located,
1951 but the underlying a and b values of a root and its conjugate
1952 are the same:
1954 >>> assert i.a == ic.a and i.b == ic.b
1956 What changes are the reported coordinates of the bounding rectangle:
1958 >>> (i.ax, i.ay), (i.bx, i.by)
1959 ((3/64, 9/8), (3/32, 75/64))
1960 >>> (ic.ax, ic.ay), (ic.bx, ic.by)
1961 ((3/64, -75/64), (3/32, -9/8))
1963 The interval can be refined once:
1965 >>> i # for reference, this is the current interval
1966 (3/64, 3/32) x (9/8, 75/64)
1968 >>> i.refine()
1969 (3/64, 3/32) x (9/8, 147/128)
1971 Several refinement steps can be taken:
1973 >>> i.refine_step(2) # 2 steps
1974 (9/128, 3/32) x (9/8, 147/128)
1976 It is also possible to refine to a given tolerance:
1978 >>> tol = min(i.dx, i.dy)/2
1979 >>> i.refine_size(tol)
1980 (9/128, 21/256) x (9/8, 291/256)
1982 A disjoint interval is one whose bounding rectangle does not
1983 overlap with another. An interval, necessarily, is not disjoint with
1984 itself, but any interval is disjoint with a conjugate since the
1985 conjugate rectangle will always be in the lower half of the complex
1986 plane and the non-conjugate in the upper half:
1988 >>> i.is_disjoint(i), i.is_disjoint(i.conjugate())
1989 (False, True)
1991 The following interval j is not disjoint from i:
1993 >>> close = CRootOf(x**10 - 2*x + 300/S(101), 9)
1994 >>> j = close._get_interval(); j
1995 (75/1616, 75/808) x (225/202, 1875/1616)
1996 >>> i.is_disjoint(j)
1997 False
1999 The two can be made disjoint, however:
2001 >>> newi, newj = i.refine_disjoint(j)
2002 >>> newi
2003 (39/512, 159/2048) x (2325/2048, 4653/4096)
2004 >>> newj
2005 (3975/51712, 2025/25856) x (29325/25856, 117375/103424)
2007 Even though the real ranges overlap, the imaginary do not, so
2008 the roots have been resolved as distinct. Intervals are disjoint
2009 when either the real or imaginary component of the intervals is
2010 distinct. In the case above, the real components have not been
2011 resolved (so we do not know, yet, which root has the smaller real
2012 part) but the imaginary part of ``close`` is larger than ``root``:
2014 >>> close.n(3)
2015 0.0771 + 1.13*I
2016 >>> root.n(3)
2017 0.0766 + 1.14*I
2018 """
2020 def __init__(self, a, b, I, Q, F1, F2, f1, f2, dom, conj=False):
2021 """Initialize new complex interval with complete information. """
2022 # a and b are the SW and NE corner of the bounding interval,
2023 # (ax, ay) and (bx, by), respectively, for the NON-CONJUGATE
2024 # root (the one with the positive imaginary part); when working
2025 # with the conjugate, the a and b value are still non-negative
2026 # but the ay, by are reversed and have oppositite sign
2027 self.a, self.b = a, b
2028 self.I, self.Q = I, Q
2030 self.f1, self.F1 = f1, F1
2031 self.f2, self.F2 = f2, F2
2033 self.dom = dom
2034 self.conj = conj
2036 @property
2037 def func(self):
2038 return ComplexInterval
2040 @property
2041 def args(self):
2042 i = self
2043 return (i.a, i.b, i.I, i.Q, i.F1, i.F2, i.f1, i.f2, i.dom, i.conj)
2045 def __eq__(self, other):
2046 if type(other) is not type(self):
2047 return False
2048 return self.args == other.args
2050 @property
2051 def ax(self):
2052 """Return ``x`` coordinate of south-western corner. """
2053 return self.a[0]
2055 @property
2056 def ay(self):
2057 """Return ``y`` coordinate of south-western corner. """
2058 if not self.conj:
2059 return self.a[1]
2060 else:
2061 return -self.b[1]
2063 @property
2064 def bx(self):
2065 """Return ``x`` coordinate of north-eastern corner. """
2066 return self.b[0]
2068 @property
2069 def by(self):
2070 """Return ``y`` coordinate of north-eastern corner. """
2071 if not self.conj:
2072 return self.b[1]
2073 else:
2074 return -self.a[1]
2076 @property
2077 def dx(self):
2078 """Return width of the complex isolating interval. """
2079 return self.b[0] - self.a[0]
2081 @property
2082 def dy(self):
2083 """Return height of the complex isolating interval. """
2084 return self.b[1] - self.a[1]
2086 @property
2087 def center(self):
2088 """Return the center of the complex isolating interval. """
2089 return ((self.ax + self.bx)/2, (self.ay + self.by)/2)
2091 @property
2092 def max_denom(self):
2093 """Return the largest denominator occurring in either endpoint. """
2094 return max(self.ax.denominator, self.bx.denominator,
2095 self.ay.denominator, self.by.denominator)
2097 def as_tuple(self):
2098 """Return tuple representation of the complex isolating
2099 interval's SW and NE corners, respectively. """
2100 return ((self.ax, self.ay), (self.bx, self.by))
2102 def __repr__(self):
2103 return "(%s, %s) x (%s, %s)" % (self.ax, self.bx, self.ay, self.by)
2105 def conjugate(self):
2106 """This complex interval really is located in lower half-plane. """
2107 return ComplexInterval(self.a, self.b, self.I, self.Q,
2108 self.F1, self.F2, self.f1, self.f2, self.dom, conj=True)
2110 def __contains__(self, item):
2111 """
2112 Say whether a complex number belongs to this complex rectangular
2113 region.
2115 Parameters
2116 ==========
2118 item : pair (re, im) or number re
2119 Either a pair giving the real and imaginary parts of the number,
2120 or else a real number.
2122 """
2123 if isinstance(item, tuple):
2124 re, im = item
2125 else:
2126 re, im = item, 0
2127 return self.ax <= re <= self.bx and self.ay <= im <= self.by
2129 def is_disjoint(self, other):
2130 """Return ``True`` if two isolation intervals are disjoint. """
2131 if isinstance(other, RealInterval):
2132 return other.is_disjoint(self)
2133 if self.conj != other.conj: # above and below real axis
2134 return True
2135 re_distinct = (self.bx < other.ax or other.bx < self.ax)
2136 if re_distinct:
2137 return True
2138 im_distinct = (self.by < other.ay or other.by < self.ay)
2139 return im_distinct
2141 def _inner_refine(self):
2142 """Internal one step complex root refinement procedure. """
2143 (u, v), (s, t) = self.a, self.b
2145 I, Q = self.I, self.Q
2147 f1, F1 = self.f1, self.F1
2148 f2, F2 = self.f2, self.F2
2150 dom = self.dom
2152 if s - u > t - v:
2153 D_L, D_R = _vertical_bisection(1, (u, v), (s, t), I, Q, F1, F2, f1, f2, dom)
2155 if D_L[0] == 1:
2156 _, a, b, I, Q, F1, F2 = D_L
2157 else:
2158 _, a, b, I, Q, F1, F2 = D_R
2159 else:
2160 D_B, D_U = _horizontal_bisection(1, (u, v), (s, t), I, Q, F1, F2, f1, f2, dom)
2162 if D_B[0] == 1:
2163 _, a, b, I, Q, F1, F2 = D_B
2164 else:
2165 _, a, b, I, Q, F1, F2 = D_U
2167 return ComplexInterval(a, b, I, Q, F1, F2, f1, f2, dom, self.conj)
2169 def refine_disjoint(self, other):
2170 """Refine an isolating interval until it is disjoint with another one. """
2171 expr = self
2172 while not expr.is_disjoint(other):
2173 expr, other = expr._inner_refine(), other._inner_refine()
2175 return expr, other
2177 def refine_size(self, dx, dy=None):
2178 """Refine an isolating interval until it is of sufficiently small size. """
2179 if dy is None:
2180 dy = dx
2181 expr = self
2182 while not (expr.dx < dx and expr.dy < dy):
2183 expr = expr._inner_refine()
2185 return expr
2187 def refine_step(self, steps=1):
2188 """Perform several steps of complex root refinement algorithm. """
2189 expr = self
2190 for _ in range(steps):
2191 expr = expr._inner_refine()
2193 return expr
2195 def refine(self):
2196 """Perform one step of complex root refinement algorithm. """
2197 return self._inner_refine()