Source code for diofant.polys.rootisolation

"""Real and complex root isolation and refinement algorithms. """

import collections
import math

from ..core import I
from .densearith import (dmp_add, dmp_mul_ground, dmp_neg, dmp_pow, dmp_quo,
                         dmp_rem, dup_rshift)
from .densebasic import (dmp_convert, dmp_degree_in, dmp_LC, dmp_permute,
                         dmp_strip, dmp_TC, dmp_terms_gcd, dmp_to_tuple,
                         dup_reverse)
from .densetools import (dmp_clear_denoms, dmp_compose, dmp_diff_in,
                         dmp_eval_in, dmp_ground_primitive, dup_mirror,
                         dup_real_imag, dup_scale, dup_shift, dup_transform)
from .euclidtools import dmp_gcd, dmp_resultant
from .factortools import dmp_trial_division
from .polyerrors import DomainError, RefinementFailed
from .sqfreetools import dmp_sqf_list, dmp_sqf_part


[docs]def dup_sturm(f, K): """ Computes the Sturm sequence of ``f`` in ``F[x]``. Given a univariate, square-free polynomial ``f(x)`` returns the associated Sturm sequence (see e.g. :cite:`Davenport1988systems`) ``f_0(x), ..., f_n(x)`` defined by:: f_0(x), f_1(x) = f(x), f'(x) f_n = -rem(f_{n-2}(x), f_{n-1}(x)) Examples ======== >>> R, x = ring("x", QQ) >>> R.dup_sturm(x**3 - 2*x**2 + x - 3) [x**3 - 2*x**2 + x - 3, 3*x**2 - 4*x + 1, 2/9*x + 25/9, -2079/4] """ if not K.is_Field: raise DomainError("can't compute Sturm sequence over %s" % K) f = dmp_sqf_part(f, 0, K) sturm = [f, dmp_diff_in(f, 1, 0, 0, K)] while sturm[-1]: s = dmp_rem(sturm[-2], sturm[-1], 0, K) sturm.append(dmp_neg(s, 0, K)) return sturm[:-1]
[docs]def dup_sign_variations(f, K): """ Compute the number of sign variations of ``f`` in ``K[x]``. Examples ======== >>> R, x = ring("x", ZZ) >>> R.dup_sign_variations(x**4 - x**2 - x + 1) 2 """ prev, k = K.zero, 0 for coeff in f: if coeff*prev < 0: k += 1 if coeff: prev = coeff return k
[docs]def dup_root_upper_bound(f, K): """Compute the LMQ upper bound for the positive roots of `f`. LMQ (Local Max Quadratic) bound was developed by Akritas-Strzebonski-Vigklas :cite:`Alkiviadis2009bounds`. """ n, P = len(f), [] t = n * [1] if dmp_LC(f, K) < 0: f = dmp_neg(f, 0, K) f = list(reversed(f)) def ilog2(a): return int(math.log(a, 2)) for i in range(n): b = int(-f[i]) if b <= 0: continue a, QL = ilog2(b), [] for j in range(i + 1, n): b = int(f[j]) if b <= 0: continue q = t[j] + a - ilog2(b) QL.append([q // (j - i), j]) if not QL: continue q = min(QL) t[q[1]] = t[q[1]] + 1 P.append(q[0]) if P: return K(2)**int(max(P) + 1)
def _mobius_from_interval(I, field): """Convert an open interval to a Mobius transform.""" s, t = I a, c = map(field, (s.numerator, s.denominator)) b, d = map(field, (t.numerator, t.denominator)) return a, b, c, d def _mobius_to_interval(M): """Convert a Mobius transform to an open interval.""" a, b, c, d = M s, t = a/c, b/d return (s, t) if s <= t else (t, s)
[docs]def dup_step_refine_real_root(f, M, K): """One step of positive real root refinement algorithm.""" a, b, c, d = M if a == b and c == d: return f, (a, b, c, d) A = dup_root_upper_bound(dup_reverse(f), K) if A is not None: A = 1/K.convert(A) else: A = K.zero if A > 16: f = dup_scale(f, A, K) a, c, A = A*a, A*c, K.one if A >= 1: f = dup_shift(f, A, K) b, d = A*a + b, A*c + d assert dmp_TC(f, K) f, g = dup_shift(f, K.one, K), f a1, b1, c1, d1 = a, a + b, c, c + d if not dmp_eval_in(f, K.zero, 0, 0, K): return f, (b1, b1, d1, d1) k = dup_sign_variations(f, K) if k == 1: a, b, c, d = a1, b1, c1, d1 else: f = dup_shift(dup_reverse(g), K.one, K) assert dmp_TC(f, K) a, b, c, d = b, a + b, d, c + d return f, (a, b, c, d)
[docs]def dup_inner_refine_real_root(f, M, K, eps=None, steps=None, disjoint=None, mobius=False): """Refine a positive root of `f` given a Mobius transform or an interval.""" a, b, c, d = M while not c: f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K) if eps is not None and steps is not None: for i in range(steps): if abs(a/c - b/d) >= eps: f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K) else: break else: if eps is not None: while abs(a/c - b/d) >= eps: f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K) if steps is not None: for i in range(steps): f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K) if disjoint is not None: while True: u, v = _mobius_to_interval((a, b, c, d)) if u < disjoint < v: f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K) else: break return (f, (a, b, c, d)) if mobius else _mobius_to_interval((a, b, c, d))
[docs]def dup_outer_refine_real_root(f, s, t, K, eps=None, steps=None, disjoint=None): """Refine a positive root of `f` given an interval `(s, t)`.""" a, b, c, d = _mobius_from_interval((s, t), K) f = dup_transform(f, dmp_strip([a, b], 0), dmp_strip([c, d], 0), K) if dup_sign_variations(f, K) != 1: raise RefinementFailed("there should be exactly one root in (%s, %s) interval" % (s, t)) return dup_inner_refine_real_root(f, (a, b, c, d), K, eps=eps, steps=steps, disjoint=disjoint)
[docs]def dup_refine_real_root(f, s, t, K, eps=None, steps=None, disjoint=None): """Refine real root's approximating interval to the given precision.""" R, K = K, K.field f = dmp_convert(f, 0, R, K) f = dmp_clear_denoms(f, 0, K)[1] if not (K.is_RationalField or K.is_RealAlgebraicField): raise DomainError("real root refinement not supported over %s" % K) if s == t: return s, t if s > t: s, t = t, s negative = False if s < 0: if t <= 0: f, s, t, negative = dup_mirror(f, K), -t, -s, True else: raise ValueError("can't refine a real root in (%s, %s)" % (s, t)) if negative and disjoint is not None: if disjoint < 0: disjoint = -disjoint else: disjoint = None s, t = dup_outer_refine_real_root(f, s, t, K, eps=eps, steps=steps, disjoint=disjoint) return (-t, -s) if negative else (s, t)
[docs]def dup_inner_isolate_real_roots(f, K, eps=None): """Internal function for isolation positive roots up to given precision.""" a, b, c, d = K.one, K.zero, K.zero, K.one k = dup_sign_variations(f, K) roots, stack = [], [(a, b, c, d, f, k)] while stack: a, b, c, d, f, k = stack.pop() A = dup_root_upper_bound(dup_reverse(f), K) if A is not None: A = 1/K.convert(A) else: A = K.zero if A > 16: f = dup_scale(f, A, K) a, c, A = A*a, A*c, K.one if A >= 1: f = dup_shift(f, A, K) b, d = A*a + b, A*c + d assert dmp_TC(f, K) k = dup_sign_variations(f, K) if k == 0: continue if k == 1: roots.append(dup_inner_refine_real_root(f, (a, b, c, d), K, eps=eps, mobius=True)) continue f1 = dup_shift(f, K.one, K) a1, b1, c1, d1, r = a, a + b, c, c + d, 0 if not dmp_TC(f1, K): roots.append((f1, (b1, b1, d1, d1))) f1, r = dup_rshift(f1, 1, K), 1 k1 = dup_sign_variations(f1, K) k2 = k - k1 - r a2, b2, c2, d2 = b, a + b, d, c + d if k2 > 1: f2 = dup_shift(dup_reverse(f), K.one, K) if not dmp_TC(f2, K): f2 = dup_rshift(f2, 1, K) k2 = dup_sign_variations(f2, K) else: f2 = None if k1 < k2: a1, a2, b1, b2 = a2, a1, b2, b1 c1, c2, d1, d2 = c2, c1, d2, d1 f1, f2, k1, k2 = f2, f1, k2, k1 if not k1: continue if f1 is None: f1 = dup_shift(dup_reverse(f), K.one, K) if not dmp_TC(f1, K): f1 = dup_rshift(f1, 1, K) if k1 == 1: roots.append(dup_inner_refine_real_root(f1, (a1, b1, c1, d1), K, eps=eps, mobius=True)) else: stack.append((a1, b1, c1, d1, f1, k1)) if not k2: continue if f2 is None: f2 = dup_shift(dup_reverse(f), K.one, K) if not dmp_TC(f2, K): f2 = dup_rshift(f2, 1, K) if k2 == 1: roots.append(dup_inner_refine_real_root(f2, (a2, b2, c2, d2), K, eps=eps, mobius=True)) else: stack.append((a2, b2, c2, d2, f2, k2)) return roots
def _discard_if_outside_interval(f, M, inf, sup, K, negative, mobius): """Discard an isolating interval if outside ``(inf, sup)``.""" while True: u, v = _mobius_to_interval(M) if negative: u, v = -v, -u if (inf is None or u >= inf) and (sup is None or v <= sup): if not mobius: return u, v else: return f, M elif (sup is not None and u > sup) or (inf is not None and v < inf): return else: f, M = dup_step_refine_real_root(f, M, K)
[docs]def dup_inner_isolate_positive_roots(f, K, eps=None, inf=None, sup=None, mobius=False): """Iteratively compute disjoint positive root isolation intervals.""" if sup is not None and sup < 0: return [] roots = dup_inner_isolate_real_roots(f, K, eps=eps) results = [] if inf is not None or sup is not None: for f, M in roots: result = _discard_if_outside_interval(f, M, inf, sup, K, False, mobius) if result is not None: results.append(result) elif not mobius: for f, M in roots: u, v = _mobius_to_interval(M) results.append((u, v)) else: results = roots return results
[docs]def dup_inner_isolate_negative_roots(f, K, inf=None, sup=None, eps=None, mobius=False): """Iteratively compute disjoint negative root isolation intervals.""" if inf is not None and inf >= 0: return [] roots = dup_inner_isolate_real_roots(dup_mirror(f, K), K, eps=eps) results = [] if inf is not None or sup is not None: for f, M in roots: result = _discard_if_outside_interval(f, M, inf, sup, K, True, mobius) if result is not None: results.append(result) elif not mobius: for f, M in roots: u, v = _mobius_to_interval(M) results.append((-v, -u)) else: results = roots return results
def _isolate_zero(f, K, inf, sup, sqf=False): """Handle special case of CF algorithm when ``f`` is homogeneous.""" (j,), f = dmp_terms_gcd(f, 0, K) if j > 0: if (inf is None or inf <= 0) and (sup is None or 0 <= sup): if not sqf: return [((K.zero, K.zero), j)], f else: return [(K.zero, K.zero)], f return [], f
[docs]def dup_isolate_real_roots_sqf(f, K, eps=None, inf=None, sup=None, blackbox=False): """Isolate real roots of a square-free polynomial.""" if not K.is_Field: R, K = K, K.field f = dmp_convert(f, 0, R, K) if not (K.is_ComplexAlgebraicField or K.is_RationalField): raise DomainError("Can't isolate real roots in domain %s" % K) f = dmp_clear_denoms(f, 0, K)[1] if K.is_ComplexAlgebraicField and not K.is_RealAlgebraicField: A, K = K, K.domain polys = [dmp_eval_in(_, K.zero, 1, 1, K) for _ in dup_real_imag(f, A)] if not polys[1]: f = polys[0] else: roots = dup_isolate_real_roots_list(polys, K, eps=eps, inf=inf, sup=sup, strict=True) roots = [_[0] for _ in roots if _[1].keys() == {0, 1}] return [RealInterval((a, b), f, K) for (a, b) in roots] if blackbox else roots if dmp_degree_in(f, 0, 0) <= 0: return [] I_zero, f = _isolate_zero(f, K, inf, sup, sqf=True) I_neg = dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup) I_pos = dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup) roots = sorted(I_neg + I_zero + I_pos) return [RealInterval((a, b), f, K) for (a, b) in roots] if blackbox else roots
[docs]def dup_isolate_real_roots(f, K, eps=None, inf=None, sup=None): """Isolate real roots. Notes ===== Implemented algorithms use Vincent-Akritas-Strzebonski (VAS) continued fractions approach :cite:`Alkiviadis2005comp`, :cite:`Alkiviadis2008cf`. """ if not K.is_Field: R, K = K, K.field f = dmp_convert(f, 0, R, K) if not (K.is_ComplexAlgebraicField or K.is_RationalField): raise DomainError("isolation of real roots not supported over %s" % K) _, factors = dmp_sqf_list(f, 0, K) factors = [(dmp_clear_denoms(f, 0, K)[1], k) for f, k in factors] if len(factors) == 1: (f, k), = factors return [(r, k) for r in dup_isolate_real_roots_sqf(f, K, eps, inf, sup)] else: I_zero, f = _isolate_zero(f, K, inf, sup) I_neg, I_pos = _real_isolate_and_disjoin(factors, K, eps, inf, sup) return sorted(I_neg + I_zero + I_pos)
[docs]def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=False, basis=False): """Isolate real roots of a list of polynomials.""" R, K = K, K.field if not (K.is_RationalField or K.is_RealAlgebraicField): raise DomainError("isolation of real roots not supported over %s" % K) if (inf is None or inf <= 0) and (sup is None or 0 <= sup): zeros, zero_indices = True, {} else: zeros = False gcd = [] for i, p in enumerate(polys): p = dmp_convert(p, 0, R, K) p = dmp_clear_denoms(p, 0, K)[1] (j,), polys[i] = dmp_terms_gcd(p, 0, K) if zeros and j > 0: zero_indices[i] = j gcd = dmp_gcd(gcd, polys[i], 0, K) polys = [dmp_quo(p, gcd, 0, K) for p in polys] if len(polys) != 2: raise NotImplementedError factors = collections.defaultdict(dict) for i, p in enumerate(polys): for f, _ in dmp_sqf_list(dmp_gcd(p, gcd, 0, K), 0, K)[1]: k1 = dmp_trial_division(gcd, [f], 0, K)[0][1] k2 = dmp_trial_division(p, [f], 0, K)[0][1] factors[tuple(f)] = {i: k1 + k2, (i + 1) % 2: k1} gcd = dmp_quo(gcd, dmp_pow(f, k1, 0, K), 0, K) p = dmp_quo(p, dmp_pow(f, k2, 0, K), 0, K) for f, k in dmp_sqf_list(p, 0, K)[1]: factors[tuple(f)] = {i: k} for f, k in dmp_sqf_list(gcd, 0, K)[1]: factors[tuple(f)] = {0: k, 1: k} I_neg, I_pos = _real_isolate_and_disjoin(factors.items(), K, eps=eps, inf=inf, sup=sup, strict=strict, basis=basis) I_zero = [] if zeros and zero_indices: if not basis: I_zero = [((K.zero, K.zero), zero_indices)] else: I_zero = [((K.zero, K.zero), zero_indices, [K.one, K.zero])] return sorted(I_neg + I_zero + I_pos)
def _disjoint_p(M, N, strict=False): """Check if Mobius transforms define disjoint intervals.""" a1, b1, c1, d1 = M a2, b2, c2, d2 = N a1d1, b1c1 = a1*d1, b1*c1 a2d2, b2c2 = a2*d2, b2*c2 if a1d1 == b1c1 and a2d2 == b2c2: return True if a1d1 > b1c1: a1, c1, b1, d1 = b1, d1, a1, c1 if a2d2 > b2c2: a2, c2, b2, d2 = b2, d2, a2, c2 if not strict: return a2*d1 >= c2*b1 or b2*c1 <= d2*a1 else: return a2*d1 > c2*b1 or b2*c1 < d2*a1 def _real_isolate_and_disjoin(factors, K, eps=None, inf=None, sup=None, strict=False, basis=False): """Isolate real roots of a list of polynomials and disjoin intervals.""" I_pos, I_neg = [], [] for i, (f, k) in enumerate(factors): f = dmp_ground_primitive(f, 0, K)[1] for F, M in dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, mobius=True): I_pos.append((F, M, k, f)) for G, N in dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, mobius=True): I_neg.append((G, N, k, f)) for i, (f, M, k, F) in enumerate(I_pos): for j, (g, N, m, G) in enumerate(I_pos[i + 1:]): while not _disjoint_p(M, N, strict=strict): f, M = dup_inner_refine_real_root(f, M, K, steps=1, mobius=True) g, N = dup_inner_refine_real_root(g, N, K, steps=1, mobius=True) I_pos[i + j + 1] = g, N, m, G I_pos[i] = f, M, k, F for i, (f, M, k, F) in enumerate(I_neg): for j, (g, N, m, G) in enumerate(I_neg[i + 1:]): while not _disjoint_p(M, N, strict=strict): f, M = dup_inner_refine_real_root(f, M, K, steps=1, mobius=True) g, N = dup_inner_refine_real_root(g, N, K, steps=1, mobius=True) I_neg[i + j + 1] = g, N, m, G I_neg[i] = f, M, k, F if strict: for i, (f, M, k, F) in enumerate(I_neg): if not M[0]: while not M[0]: f, M = dup_inner_refine_real_root(f, M, K, steps=1, mobius=True) I_neg[i] = f, M, k, F break for j, (g, N, m, G) in enumerate(I_pos): if not N[0]: while not N[0]: g, N = dup_inner_refine_real_root(g, N, K, steps=1, mobius=True) I_pos[j] = g, N, m, G break I_neg = [(_mobius_to_interval(M), k, f) for (_, M, k, f) in I_neg] I_pos = [(_mobius_to_interval(M), k, f) for (_, M, k, f) in I_pos] if not basis: I_neg = [((-v, -u), k) for ((u, v), k, _) in I_neg] I_pos = [((+u, +v), k) for ((u, v), k, _) in I_pos] else: I_neg = [((-v, -u), k, f) for ((u, v), k, f) in I_neg] I_pos = [((+u, +v), k, f) for ((u, v), k, f) in I_pos] return I_neg, I_pos
[docs]def dup_count_real_roots(f, K, inf=None, sup=None): """Returns the number of distinct real roots of ``f`` in ``[inf, sup]``.""" if dmp_degree_in(f, 0, 0) <= 0: return 0 if not K.is_Field: R, K = K, K.field f = dmp_convert(f, 0, R, K) if not (K.is_ComplexAlgebraicField or K.is_RationalField): raise DomainError("Can't count real roots in domain %s" % K) if K.is_ComplexAlgebraicField and not K.is_RealAlgebraicField: return sum(k for *_, k in dup_isolate_real_roots(f, K, inf, sup)) sturm = dup_sturm(f, K) if inf is None: signs_inf = dup_sign_variations([dmp_LC(s, K)*(-1)**dmp_degree_in(s, 0, 0) for s in sturm], K) else: signs_inf = dup_sign_variations([dmp_eval_in(s, inf, 0, 0, K) for s in sturm], K) if sup is None: signs_sup = dup_sign_variations([dmp_LC(s, K) for s in sturm], K) else: signs_sup = dup_sign_variations([dmp_eval_in(s, sup, 0, 0, K) for s in sturm], K) count = abs(signs_inf - signs_sup) if inf is not None and not dmp_eval_in(f, inf, 0, 0, K): count += 1 return count
OO = 'OO' # Origin of (re, im) coordinate system Q1 = 'Q1' # Quadrant #1 (++): re > 0 and im > 0 Q2 = 'Q2' # Quadrant #2 (-+): re < 0 and im > 0 Q3 = 'Q3' # Quadrant #3 (--): re < 0 and im < 0 Q4 = 'Q4' # Quadrant #4 (+-): re > 0 and im < 0 A1 = 'A1' # Axis #1 (+0): re > 0 and im = 0 A2 = 'A2' # Axis #2 (0+): re = 0 and im > 0 A3 = 'A3' # Axis #3 (-0): re < 0 and im = 0 A4 = 'A4' # Axis #4 (0-): re = 0 and im < 0 _rules_simple = { # Q --> Q (same) => no change (Q1, Q1): 0, (Q2, Q2): 0, (Q3, Q3): 0, (Q4, Q4): 0, # A -- CCW --> Q => +1/4 (CCW) (A1, Q1): 1, (A2, Q2): 1, (A3, Q3): 1, (A4, Q4): 1, # A -- CW --> Q => -1/4 (CCW) (A1, Q4): 2, (A2, Q1): 2, (A3, Q2): 2, (A4, Q3): 2, # Q -- CCW --> A => +1/4 (CCW) (Q1, A2): 3, (Q2, A3): 3, (Q3, A4): 3, (Q4, A1): 3, # Q -- CW --> A => -1/4 (CCW) (Q1, A1): 4, (Q2, A2): 4, (Q3, A3): 4, (Q4, A4): 4, # Q -- CCW --> Q => +1/2 (CCW) (Q1, Q2): +5, (Q2, Q3): +5, (Q3, Q4): +5, (Q4, Q1): +5, # Q -- CW --> Q => -1/2 (CW) (Q1, Q4): -5, (Q2, Q1): -5, (Q3, Q2): -5, (Q4, Q3): -5, } _rules_ambiguous = { # A -- CCW --> Q => { +1/4 (CCW), -9/4 (CW) } (A1, OO, Q1): -1, (A2, OO, Q2): -1, (A3, OO, Q3): -1, (A4, OO, Q4): -1, # A -- CW --> Q => { -1/4 (CCW), +7/4 (CW) } (A1, OO, Q4): -2, (A2, OO, Q1): -2, (A3, OO, Q2): -2, (A4, OO, Q3): -2, # Q -- CCW --> A => { +1/4 (CCW), -9/4 (CW) } (Q1, OO, A2): -3, (Q2, OO, A3): -3, (Q3, OO, A4): -3, (Q4, OO, A1): -3, # Q -- CW --> A => { -1/4 (CCW), +7/4 (CW) } (Q1, OO, A1): -4, (Q2, OO, A2): -4, (Q3, OO, A3): -4, (Q4, OO, A4): -4, # A -- OO --> A => { +1 (CCW), -1 (CW) } (A1, A3): 7, (A2, A4): 7, (A3, A1): 7, (A4, A2): 7, (A1, OO, A3): 7, (A2, OO, A4): 7, (A3, OO, A1): 7, (A4, OO, A2): 7, # Q -- DIA --> Q => { +1 (CCW), -1 (CW) } (Q1, Q3): 8, (Q2, Q4): 8, (Q3, Q1): 8, (Q4, Q2): 8, (Q1, OO, Q3): 8, (Q2, OO, Q4): 8, (Q3, OO, Q1): 8, (Q4, OO, Q2): 8, # A --- R ---> A => { +1/2 (CCW), -3/2 (CW) } (A1, A2): 9, (A2, A3): 9, (A3, A4): 9, (A4, A1): 9, (A1, OO, A2): 9, (A2, OO, A3): 9, (A3, OO, A4): 9, (A4, OO, A1): 9, # A --- L ---> A => { +3/2 (CCW), -1/2 (CW) } (A1, A4): 10, (A2, A1): 10, (A3, A2): 10, (A4, A3): 10, (A1, OO, A4): 10, (A2, OO, A1): 10, (A3, OO, A2): 10, (A4, OO, A3): 10, # Q --- 1 ---> A => { +3/4 (CCW), -5/4 (CW) } (Q1, A3): 11, (Q2, A4): 11, (Q3, A1): 11, (Q4, A2): 11, (Q1, OO, A3): 11, (Q2, OO, A4): 11, (Q3, OO, A1): 11, (Q4, OO, A2): 11, # Q --- 2 ---> A => { +5/4 (CCW), -3/4 (CW) } (Q1, A4): 12, (Q2, A1): 12, (Q3, A2): 12, (Q4, A3): 12, (Q1, OO, A4): 12, (Q2, OO, A1): 12, (Q3, OO, A2): 12, (Q4, OO, A3): 12, # A --- 1 ---> Q => { +5/4 (CCW), -3/4 (CW) } (A1, Q3): 13, (A2, Q4): 13, (A3, Q1): 13, (A4, Q2): 13, (A1, OO, Q3): 13, (A2, OO, Q4): 13, (A3, OO, Q1): 13, (A4, OO, Q2): 13, # A --- 2 ---> Q => { +3/4 (CCW), -5/4 (CW) } (A1, Q2): 14, (A2, Q3): 14, (A3, Q4): 14, (A4, Q1): 14, (A1, OO, Q2): 14, (A2, OO, Q3): 14, (A3, OO, Q4): 14, (A4, OO, Q1): 14, # Q --> OO --> Q => { +1/2 (CCW), -3/2 (CW) } (Q1, OO, Q2): 15, (Q2, OO, Q3): 15, (Q3, OO, Q4): 15, (Q4, OO, Q1): 15, # Q --> OO --> Q => { +3/2 (CCW), -1/2 (CW) } (Q1, OO, Q4): 16, (Q2, OO, Q1): 16, (Q3, OO, Q2): 16, (Q4, OO, Q3): 16, # A --> OO --> A => { +2 (CCW), 0 (CW) } (A1, OO, A1): 17, (A2, OO, A2): 17, (A3, OO, A3): 17, (A4, OO, A4): 17, # Q --> OO --> Q => { +2 (CCW), 0 (CW) } (Q1, OO, Q1): 18, (Q2, OO, Q2): 18, (Q3, OO, Q3): 18, (Q4, OO, Q4): 18, } _values = { 0: [(+0, 1)], 1: [(+1, 4)], 2: [(-1, 4)], 3: [(+1, 4)], 4: [(-1, 4)], -1: [(+9, 4), (+1, 4)], -2: [(+7, 4), (-1, 4)], -3: [(+9, 4), (+1, 4)], -4: [(+7, 4), (-1, 4)], +5: [(+1, 2)], -5: [(-1, 2)], 7: [(+1, 1), (-1, 1)], 8: [(+1, 1), (-1, 1)], 9: [(+1, 2), (-3, 2)], 10: [(+3, 2), (-1, 2)], 11: [(+3, 4), (-5, 4)], 12: [(+5, 4), (-3, 4)], 13: [(+5, 4), (-3, 4)], 14: [(+3, 4), (-5, 4)], 15: [(+1, 2), (-3, 2)], 16: [(+3, 2), (-1, 2)], 17: [(+2, 1), (+0, 1)], 18: [(+2, 1), (+0, 1)], } def _classify_point(re, im): """Return the half-axis (or origin) on which (re, im) point is located.""" assert not re or not im if not re and not im: return OO if not re: if im > 0: return A2 else: return A4 else: if re > 0: return A1 else: return A3 def _intervals_to_quadrants(intervals, f1, f2, s, t, F): """Generate a sequence of extended quadrants from a list of critical points.""" if not intervals: return [] Q = [] if not f1: (a, b), _, _ = intervals[0] if a == b == s: if len(intervals) == 1: if dmp_eval_in(f2, t, 0, 0, F) > 0: return [OO, A2] else: return [OO, A4] else: (a, _), _, _ = intervals[1] if dmp_eval_in(f2, (s + a)/2, 0, 0, F) > 0: Q.extend([OO, A2]) f2_sgn = +1 else: Q.extend([OO, A4]) f2_sgn = -1 intervals = intervals[1:] else: if dmp_eval_in(f2, s, 0, 0, F) > 0: Q.append(A2) f2_sgn = +1 else: Q.append(A4) f2_sgn = -1 for (a, _), indices, _ in intervals: Q.append(OO) assert indices[1] % 2 == 1 f2_sgn = -f2_sgn if a != t: if f2_sgn > 0: Q.append(A2) else: Q.append(A4) return Q if not f2: (a, b), _, _ = intervals[0] if a == b == s: if len(intervals) == 1: if dmp_eval_in(f1, t, 0, 0, F) > 0: return [OO, A1] else: return [OO, A3] else: (a, _), _, _ = intervals[1] if dmp_eval_in(f1, (s + a)/2, 0, 0, F) > 0: Q.extend([OO, A1]) f1_sgn = +1 else: Q.extend([OO, A3]) f1_sgn = -1 intervals = intervals[1:] else: if dmp_eval_in(f1, s, 0, 0, F) > 0: Q.append(A1) f1_sgn = +1 else: Q.append(A3) f1_sgn = -1 for (a, _), indices, _ in intervals: Q.append(OO) assert indices[0] % 2 == 1 f1_sgn = -f1_sgn if a != t: if f1_sgn > 0: Q.append(A1) else: Q.append(A3) return Q re = dmp_eval_in(f1, s, 0, 0, F) im = dmp_eval_in(f2, s, 0, 0, F) if not re or not im: Q.append(_classify_point(re, im)) if len(intervals) == 1: re = dmp_eval_in(f1, t, 0, 0, F) im = dmp_eval_in(f2, t, 0, 0, F) else: (a, _), _, _ = intervals[1] re = dmp_eval_in(f1, (s + a)/2, 0, 0, F) im = dmp_eval_in(f2, (s + a)/2, 0, 0, F) intervals = intervals[1:] if re > 0: f1_sgn = +1 else: f1_sgn = -1 if im > 0: f2_sgn = +1 else: f2_sgn = -1 sgn = { (+1, +1): Q1, (-1, +1): Q2, (-1, -1): Q3, (+1, -1): Q4, } Q.append(sgn[(f1_sgn, f2_sgn)]) for (a, b), indices, _ in intervals: if a == b: re = dmp_eval_in(f1, a, 0, 0, F) im = dmp_eval_in(f2, a, 0, 0, F) Q.append(_classify_point(re, im)) if 0 in indices: if indices[0] % 2 == 1: f1_sgn = -f1_sgn if 1 in indices: if indices[1] % 2 == 1: f2_sgn = -f2_sgn if not (a == b and b == t): Q.append(sgn[(f1_sgn, f2_sgn)]) return Q def _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=None): """Transform sequences of quadrants to a sequence of rules.""" if exclude is True: edges = [1, 1, 0, 0] corners = { (0, 1): 1, (1, 2): 1, (2, 3): 0, (3, 0): 1, } else: edges = [0, 0, 0, 0] corners = { (0, 1): 0, (1, 2): 0, (2, 3): 0, (3, 0): 0, } if exclude is not None and exclude is not True: exclude = set(exclude) for i, edge in enumerate(['S', 'E', 'N', 'W']): if edge in exclude: edges[i] = 1 for i, corner in enumerate(['SW', 'SE', 'NE', 'NW']): if corner in exclude: corners[((i - 1) % 4, i)] = 1 QQ, rules = [Q_L1, Q_L2, Q_L3, Q_L4], [] for i, Q in enumerate(QQ): if not Q: continue if Q[-1] == OO: Q = Q[:-1] if Q[0] == OO: j, Q = (i - 1) % 4, Q[1:] qq = QQ[j][-2], OO, Q[0] rules.append((_rules_ambiguous[qq], corners[(j, i)])) q1, k = Q[0], 1 while k < len(Q): q2, k = Q[k], k + 1 if q2 != OO: qq = q1, q2 rules.append((_rules_simple[qq], 0)) else: qq, k = (q1, q2, Q[k]), k + 1 rules.append((_rules_ambiguous[qq], edges[i])) q1 = qq[-1] return rules def _reverse_intervals(intervals): """Reverse intervals for traversal from right to left and from top to bottom.""" return [((b, a), indices, f) for (a, b), indices, f in reversed(intervals)] def _winding_number(T, field): """Compute the winding number of the input polynomial, i.e. the number of roots.""" return int(sum((field(_values[t][i][0])/field(_values[t][i][1]) for t, i in T), field(0)) / field(2)) def _roots_bound(f, F): lc = F.to_expr(dmp_LC(f, F)) B = 2*max(abs(F.to_expr(c)/lc) for c in f) if not F.is_AlgebraicField: return F.convert(B) else: return F.domain(int(100*B) + 1)/F.domain(100) def _get_rectangle(f1, f2, F, inf, sup, exclude=None): (u, v), (s, t) = inf, sup f1L1 = dmp_eval_in(f1, v, 1, 1, F) f2L1 = dmp_eval_in(f2, v, 1, 1, F) f1L2 = dmp_eval_in(f1, s, 0, 1, F) f2L2 = dmp_eval_in(f2, s, 0, 1, F) f1L3 = dmp_eval_in(f1, t, 1, 1, F) f2L3 = dmp_eval_in(f2, t, 1, 1, F) f1L4 = dmp_eval_in(f1, u, 0, 1, F) f2L4 = dmp_eval_in(f2, u, 0, 1, F) S_L1 = [f1L1, f2L1] S_L2 = [f1L2, f2L2] S_L3 = [f1L3, f2L3] S_L4 = [f1L4, f2L4] I_L1 = dup_isolate_real_roots_list(S_L1, F, inf=u, sup=s, basis=True, strict=True) I_L2 = dup_isolate_real_roots_list(S_L2, F, inf=v, sup=t, basis=True, strict=True) I_L3 = dup_isolate_real_roots_list(S_L3, F, inf=u, sup=s, basis=True, strict=True) I_L4 = dup_isolate_real_roots_list(S_L4, F, inf=v, sup=t, basis=True, strict=True) I_L3 = _reverse_intervals(I_L3) I_L4 = _reverse_intervals(I_L4) Q_L1 = _intervals_to_quadrants(I_L1, f1L1, f2L1, u, s, F) Q_L2 = _intervals_to_quadrants(I_L2, f1L2, f2L2, v, t, F) Q_L3 = _intervals_to_quadrants(I_L3, f1L3, f2L3, s, u, F) Q_L4 = _intervals_to_quadrants(I_L4, f1L4, f2L4, t, v, F) T = _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=exclude) N = _winding_number(T, F) I_L = I_L1, I_L2, I_L3, I_L4 Q_L = Q_L1, Q_L2, Q_L3, Q_L4 F1 = f1L1, f1L2, f1L3, f1L4 F2 = f2L1, f2L2, f2L3, f2L4 return N, inf, sup, I_L, Q_L, F1, F2
[docs]def dup_count_complex_roots(f, K, inf=None, sup=None, exclude=None): """Count all roots in [u + v*I, s + t*I] rectangle using Collins-Krandick algorithm.""" if not K.is_Field: R, K = K, K.field f = dmp_convert(f, 0, R, K) if not (K.is_ComplexAlgebraicField or K.is_RationalField): raise DomainError("Can't count complex roots in domain %s" % K) if not all(isinstance(_, tuple) for _ in (inf, sup)): B = _roots_bound(f, K) if isinstance(inf, tuple): u, v = inf elif inf is not None: u, v = inf, -B else: u, v = -B, -B if isinstance(sup, tuple): s, t = sup elif sup is not None: s, t = sup, B else: s, t = B, B f1, f2 = dup_real_imag(f, K) if not (K.is_RationalField or K.is_RealAlgebraicField): K = K.domain return _get_rectangle(f1, f2, K, (u, v), (s, t), exclude)[0]
def _vertical_bisection(N, a, b, I, Q, F1, F2, f1, f2, F): """Vertical bisection step in Collins-Krandick root isolation algorithm.""" (u, v), (s, t) = a, b I_L1, I_L2, I_L3, I_L4 = I Q_L1, Q_L2, Q_L3, Q_L4 = Q f1L1F, f1L2F, f1L3F, f1L4F = F1 f2L1F, f2L2F, f2L3F, f2L4F = F2 x = (u + s) / F(2) f1V = dmp_eval_in(f1, x, 0, 1, F) f2V = dmp_eval_in(f2, x, 0, 1, F) I_V = dup_isolate_real_roots_list([f1V, f2V], F, inf=v, sup=t, strict=True, basis=True) I_L1_L, I_L1_R = [], [] I_L2_L, I_L2_R = I_V, I_L2 I_L3_L, I_L3_R = [], [] I_L4_L, I_L4_R = I_L4, _reverse_intervals(I_V) for I in I_L1: (a, b), indices, h = I if a == b: if a == x: I_L1_L.append(I) I_L1_R.append(I) elif a < x: I_L1_L.append(I) else: I_L1_R.append(I) else: if b <= x: I_L1_L.append(I) elif a >= x: I_L1_R.append(I) else: a, b = dup_refine_real_root(h, a, b, F, disjoint=x) if b <= x: I_L1_L.append(((a, b), indices, h)) if a >= x: I_L1_R.append(((a, b), indices, h)) for I in I_L3: (b, a), indices, h = I if a == b: if a == x: I_L3_L.append(I) I_L3_R.append(I) elif a < x: I_L3_L.append(I) else: I_L3_R.append(I) else: if b <= x: I_L3_L.append(I) elif a >= x: I_L3_R.append(I) else: a, b = dup_refine_real_root(h, a, b, F, disjoint=x) if b <= x: I_L3_L.append(((b, a), indices, h)) if a >= x: I_L3_R.append(((b, a), indices, h)) Q_L1_L = _intervals_to_quadrants(I_L1_L, f1L1F, f2L1F, u, x, F) Q_L2_L = _intervals_to_quadrants(I_L2_L, f1V, f2V, v, t, F) Q_L3_L = _intervals_to_quadrants(I_L3_L, f1L3F, f2L3F, x, u, F) Q_L4_L = Q_L4 Q_L1_R = _intervals_to_quadrants(I_L1_R, f1L1F, f2L1F, x, s, F) Q_L2_R = Q_L2 Q_L3_R = _intervals_to_quadrants(I_L3_R, f1L3F, f2L3F, s, x, F) Q_L4_R = _intervals_to_quadrants(I_L4_R, f1V, f2V, t, v, F) T_L = _traverse_quadrants(Q_L1_L, Q_L2_L, Q_L3_L, Q_L4_L, exclude=True) T_R = _traverse_quadrants(Q_L1_R, Q_L2_R, Q_L3_R, Q_L4_R, exclude=True) N_L = _winding_number(T_L, F) N_R = _winding_number(T_R, F) I_L = I_L1_L, I_L2_L, I_L3_L, I_L4_L Q_L = Q_L1_L, Q_L2_L, Q_L3_L, Q_L4_L I_R = I_L1_R, I_L2_R, I_L3_R, I_L4_R Q_R = Q_L1_R, Q_L2_R, Q_L3_R, Q_L4_R F1_L = f1L1F, f1V, f1L3F, f1L4F F2_L = f2L1F, f2V, f2L3F, f2L4F F1_R = f1L1F, f1L2F, f1L3F, f1V F2_R = f2L1F, f2L2F, f2L3F, f2V a, b = (u, v), (x, t) c, d = (x, v), (s, t) D_L = N_L, a, b, I_L, Q_L, F1_L, F2_L D_R = N_R, c, d, I_R, Q_R, F1_R, F2_R return D_L, D_R def _horizontal_bisection(N, a, b, I, Q, F1, F2, f1, f2, F): """Horizontal bisection step in Collins-Krandick root isolation algorithm.""" (u, v), (s, t) = a, b I_L1, I_L2, I_L3, I_L4 = I Q_L1, Q_L2, Q_L3, Q_L4 = Q f1L1F, f1L2F, f1L3F, f1L4F = F1 f2L1F, f2L2F, f2L3F, f2L4F = F2 y = (v + t) / F(2) f1H = dmp_eval_in(f1, y, 1, 1, F) f2H = dmp_eval_in(f2, y, 1, 1, F) I_H = dup_isolate_real_roots_list([f1H, f2H], F, inf=u, sup=s, strict=True, basis=True) I_L1_B, I_L1_U = I_L1, I_H I_L2_B, I_L2_U = [], [] I_L3_B, I_L3_U = _reverse_intervals(I_H), I_L3 I_L4_B, I_L4_U = [], [] for I in I_L2: (a, b), indices, h = I if a == b: if a == y: I_L2_B.append(I) I_L2_U.append(I) elif a < y: I_L2_B.append(I) else: I_L2_U.append(I) else: if b <= y: I_L2_B.append(I) elif a >= y: I_L2_U.append(I) else: a, b = dup_refine_real_root(h, a, b, F, disjoint=y) if b <= y: I_L2_B.append(((a, b), indices, h)) if a >= y: I_L2_U.append(((a, b), indices, h)) for I in I_L4: (b, a), indices, h = I if a == b: if a == y: I_L4_B.append(I) I_L4_U.append(I) elif a < y: I_L4_B.append(I) else: I_L4_U.append(I) else: if b <= y: I_L4_B.append(I) elif a >= y: I_L4_U.append(I) else: a, b = dup_refine_real_root(h, a, b, F, disjoint=y) if b <= y: I_L4_B.append(((b, a), indices, h)) if a >= y: I_L4_U.append(((b, a), indices, h)) Q_L1_B = Q_L1 Q_L2_B = _intervals_to_quadrants(I_L2_B, f1L2F, f2L2F, v, y, F) Q_L3_B = _intervals_to_quadrants(I_L3_B, f1H, f2H, s, u, F) Q_L4_B = _intervals_to_quadrants(I_L4_B, f1L4F, f2L4F, y, v, F) Q_L1_U = _intervals_to_quadrants(I_L1_U, f1H, f2H, u, s, F) Q_L2_U = _intervals_to_quadrants(I_L2_U, f1L2F, f2L2F, y, t, F) Q_L3_U = Q_L3 Q_L4_U = _intervals_to_quadrants(I_L4_U, f1L4F, f2L4F, t, y, F) T_B = _traverse_quadrants(Q_L1_B, Q_L2_B, Q_L3_B, Q_L4_B, exclude=True) T_U = _traverse_quadrants(Q_L1_U, Q_L2_U, Q_L3_U, Q_L4_U, exclude=True) N_B = _winding_number(T_B, F) N_U = _winding_number(T_U, F) I_B = I_L1_B, I_L2_B, I_L3_B, I_L4_B Q_B = Q_L1_B, Q_L2_B, Q_L3_B, Q_L4_B I_U = I_L1_U, I_L2_U, I_L3_U, I_L4_U Q_U = Q_L1_U, Q_L2_U, Q_L3_U, Q_L4_U F1_B = f1L1F, f1L2F, f1H, f1L4F F2_B = f2L1F, f2L2F, f2H, f2L4F F1_U = f1H, f1L2F, f1L3F, f1L4F F2_U = f2H, f2L2F, f2L3F, f2L4F a, b = (u, v), (s, y) c, d = (u, y), (s, t) D_B = N_B, a, b, I_B, Q_B, F1_B, F2_B D_U = N_U, c, d, I_U, Q_U, F1_U, F2_U return D_B, D_U def _depth_first_select(rectangles): """Find a rectangle of minimum area for bisection.""" min_area, j = None, None for i, (_, (u, v), (s, t), _, _, _, _) in enumerate(rectangles): area = (s - u)*(t - v) if min_area is None or area < min_area: min_area, j = area, i return rectangles.pop(j) def _rectangle_small_p(a, b, eps): """Return ``True`` if the given rectangle is small enough.""" (u, v), (s, t) = a, b if eps is not None: return s - u < eps and t - v < eps else: return True
[docs]def dup_isolate_complex_roots_sqf(f, K, eps=None, inf=None, sup=None, blackbox=False): """Isolate complex roots of a square-free polynomial using Collins-Krandick algorithm.""" if dmp_degree_in(f, 0, 0) <= 0: return [] if not K.is_Field: R, K = K, K.field f = dmp_convert(f, 0, R, K) if not (K.is_ComplexAlgebraicField or K.is_RationalField): raise DomainError("Can't isolate complex roots in domain %s" % K) if not all(isinstance(_, tuple) for _ in (inf, sup)): B = _roots_bound(f, K) if isinstance(inf, tuple): u, v = inf elif inf is not None: u, v = inf, -B else: u, v = -B, -B if isinstance(sup, tuple): s, t = sup elif sup is not None: s, t = sup, B else: s, t = B, B if t <= v or s <= u: raise ValueError("not a valid complex isolation rectangle") if v < 0 < t: roots = dup_isolate_complex_roots_sqf(f, K, eps=eps, inf=(u, 0), sup=(s, t), blackbox=True) if K.is_RationalField or K.is_RealAlgebraicField: _roots = [] for root in roots: croot = root.conjugate() if croot.ay >= v: _roots.append(croot) _roots.append(root) roots = _roots else: # Take conjugated polynomial to get solutions in the # bottom half-plane. f = [_.conjugate() for _ in f] roots += [_.conjugate() for _ in dup_isolate_complex_roots_sqf(f, K, eps=eps, inf=(u, 0), sup=(s, -v), blackbox=True)] roots = sorted(roots, key=lambda r: (r.ax, r.ay)) return roots if blackbox else [r.as_tuple() for r in roots] f1, f2 = dup_real_imag(f, K) if not (K.is_RationalField or K.is_RealAlgebraicField): K = K.domain N, a, b, I_L, Q_L, F1, F2 = _get_rectangle(f1, f2, K, (u, v), (s, t)) if not N: return [] rectangles, roots = [], [] if N == 1 and (v > 0 or t < 0): roots.append(ComplexInterval(a, b, I_L, Q_L, F1, F2, f1, f2, K)) else: rectangles.append((N, a, b, I_L, Q_L, F1, F2)) while rectangles: N, (u, v), (s, t), I_L, Q_L, F1, F2 = _depth_first_select(rectangles) if s - u > t - v: D_L, D_R = _vertical_bisection(N, (u, v), (s, t), I_L, Q_L, F1, F2, f1, f2, K) N_L, a, b, I_L, Q_L, F1_L, F2_L = D_L N_R, c, d, I_R, Q_R, F1_R, F2_R = D_R if N_L >= 1: if N_L == 1 and _rectangle_small_p(a, b, eps): roots.append(ComplexInterval(a, b, I_L, Q_L, F1_L, F2_L, f1, f2, K)) else: rectangles.append(D_L) if N_R >= 1: if N_R == 1 and _rectangle_small_p(c, d, eps): roots.append(ComplexInterval(c, d, I_R, Q_R, F1_R, F2_R, f1, f2, K)) else: rectangles.append(D_R) else: D_B, D_U = _horizontal_bisection(N, (u, v), (s, t), I_L, Q_L, F1, F2, f1, f2, K) N_B, a, b, I_B, Q_B, F1_B, F2_B = D_B N_U, c, d, I_U, Q_U, F1_U, F2_U = D_U if N_B >= 1: if N_B == 1 and _rectangle_small_p(a, b, eps): roots.append(ComplexInterval(a, b, I_B, Q_B, F1_B, F2_B, f1, f2, K)) else: rectangles.append(D_B) if N_U >= 1: if N_U == 1 and _rectangle_small_p(c, d, eps): roots.append(ComplexInterval(c, d, I_U, Q_U, F1_U, F2_U, f1, f2, K)) else: rectangles.append(D_U) roots = sorted(roots, key=lambda r: (r.ax, r.ay)) return roots if blackbox else [r.as_tuple() for r in roots]
[docs]def dup_isolate_all_roots_sqf(f, K, eps=None, inf=None, sup=None, blackbox=False): """Isolate real and complex roots of a square-free polynomial ``f``.""" return (dup_isolate_real_roots_sqf(f, K, eps=eps, inf=inf, sup=sup, blackbox=blackbox), dup_isolate_complex_roots_sqf(f, K, eps=eps, inf=inf, sup=sup, blackbox=blackbox))
[docs]def dup_isolate_all_roots(f, K, eps=None, inf=None, sup=None): """Isolate real and complex roots of a non-square-free polynomial ``f``.""" if not K.is_IntegerRing and not K.is_RationalField: raise DomainError("isolation of real and complex roots is not supported over %s" % K) _, factors = dmp_sqf_list(f, 0, K) if len(factors) == 1: (f, k), = factors real_part, complex_part = dup_isolate_all_roots_sqf(f, K, eps=eps, inf=inf, sup=sup) real_part = [((a, b), k) for (a, b) in real_part] complex_part = [((a, b), k) for (a, b) in complex_part] return real_part, complex_part else: raise NotImplementedError("only trivial square-free polynomials are supported")
[docs]class RealInterval: """A fully qualified representation of a real isolation interval.""" def __init__(self, data, f, dom): """Initialize new real interval with complete information.""" if len(data) == 2: s, t = data self.neg = False if s < 0: if t <= 0: f, s, t, self.neg = dup_mirror(f, dom), -t, -s, True else: raise ValueError("can't refine a real root in (%s, %s)" % (s, t)) a, b, c, d = _mobius_from_interval((s, t), dom.field) f = dup_transform(f, dmp_strip([a, b], 0), dmp_strip([c, d], 0), dom) self.mobius = a, b, c, d else: self.mobius = data[:-1] self.neg = data[-1] self.f, self.domain = f, dom @property def a(self): """Return the position of the left end.""" a, b, c, d = self.mobius if not self.neg: if a*d < b*c: return a/c return b/d else: if a*d > b*c: return -a/c return -b/d @property def b(self): """Return the position of the right end.""" was = self.neg self.neg = not was rv = -self.a self.neg = was return rv @property def center(self): """Return the center of the real isolating interval.""" return (self.a + self.b)/2
[docs] def as_tuple(self): """Return tuple representation of real isolating interval.""" return self.a, self.b
[docs] def is_disjoint(self, other): """Return ``True`` if two isolation intervals are disjoint.""" return self.b <= other.a or other.b <= self.a
[docs] def refine(self): """Perform one step of real root refinement algorithm.""" assert self.mobius is not None f, mobius = dup_inner_refine_real_root(self.f, self.mobius, self.domain, steps=1, mobius=True) return RealInterval(mobius + (self.neg,), f, self.domain)
[docs]class ComplexInterval: """A fully qualified representation of a complex isolation interval. The printed form is shown as (x1, y1) x (x2, y2): the southwest x northeast coordinates of the interval's rectangle. """ def __init__(self, a, b, I, Q, F1, F2, f1, f2, dom, conj=False): """Initialize new complex interval with complete information.""" self.a, self.b = a, b # the southwest and northeast corner: (x1, y1), (x2, y2) self.I, self.Q = I, Q self.f1, self.F1 = f1, F1 self.f2, self.F2 = f2, F2 self.domain = dom self.conj = conj @property def ax(self): """Return ``x`` coordinate of south-western corner.""" return self.a[0] @property def ay(self): """Return ``y`` coordinate of south-western corner.""" if not self.conj: return +self.a[1] else: return -self.b[1] @property def bx(self): """Return ``x`` coordinate of north-eastern corner.""" return self.b[0] @property def by(self): """Return ``y`` coordinate of north-eastern corner.""" if not self.conj: return +self.b[1] else: return -self.a[1] @property def center(self): """Return the center of the complex isolating interval.""" return (self.ax + self.bx)/2, (self.ay + self.by)/2
[docs] def as_tuple(self): """Return tuple representation of complex isolating interval.""" return (self.ax, self.ay), (self.bx, self.by)
[docs] def conjugate(self): """Return conjugated isolating interval.""" return ComplexInterval(self.a, self.b, self.I, self.Q, self.F1, self.F2, self.f1, self.f2, self.domain, not self.conj)
[docs] def is_disjoint(self, other, check_re_refinement=False, re_disjoint=False): """Return ``True`` if two isolation intervals are disjoint. Parameters ========== check_re_refinement : bool, optional If enabled, test that either real projections of isolation intervals are disjoint or roots share common real part. re_disjoint : bool, optional If enabled, return ``True`` only if real projections of isolation intervals are disjoint. """ test_re = self.bx <= other.ax or other.bx <= self.ax if test_re or re_disjoint: return test_re if not check_re_refinement: return self.by <= other.ay or other.by <= self.ay resultants = [] for i in (self, other): re = dmp_permute(i.f1, [1, 0], 1, i.domain) im = dmp_permute(i.f2, [1, 0], 1, i.domain) re, im = map(lambda x: dmp_to_tuple(x, 1), [re, im]) resultants.append(dmp_resultant(re, im, 1, i.domain)) dom = self.domain.unify(other.domain) gcd = dmp_gcd(*resultants, 0, dom) gcd_roots = dup_isolate_real_roots(gcd, dom, inf=max(self.ax, other.ax), sup=min(self.bx, other.bx)) if len(gcd_roots) != 1: return False l, r = gcd_roots[0][0] if l == r: # Can't use dup_count_complex_roots() here, as the new isolation # interval will be degenerate: a vertical line segment. Make # sure we only count roots on the northern/western edges and # on the north-western corner of the original isolation rectangle. for i in (self, other): dom = i.domain.algebraic_field(I) f1 = dmp_convert(dmp_eval_in(i.f1, l, 0, 1, i.domain), 0, i.domain, dom) f2 = dmp_convert(dmp_eval_in(i.f2, l, 0, 1, i.domain), 0, i.domain, dom) f = dmp_add(f1, dmp_mul_ground(f2, dom.unit, 0, dom), 0, dom) f = dmp_compose(f, [-dom.unit, 0], 0, dom) if i.conj: f = [_.conjugate() for _ in f] f = dmp_compose(f, [dom.unit, 0], 0, dom) r = dup_isolate_real_roots(f, dom, inf=i.ay, sup=i.by) if len([1 for _ in r if ((not i.conj and i.ay < _[0][0]) or (i.conj and _[0][1] < i.by)) and l < i.bx]) != 1: return False return True else: return all(_get_rectangle(i.f1, i.f2, i.domain, (l, i.ay), (r, i.by))[0] == 1 for i in (self, other))
[docs] def refine(self, vertical=False): """Perform one step of complex root refinement algorithm.""" a, b = (u, v), (s, t) = self.a, self.b I, Q = self.I, self.Q f1, F1 = self.f1, self.F1 f2, F2 = self.f2, self.F2 dom = self.domain if s - u > t - v or vertical: D_L, D_R = _vertical_bisection(1, a, b, I, Q, F1, F2, f1, f2, dom) if D_L[0] == 1: _, a, b, I, Q, F1, F2 = D_L else: _, a, b, I, Q, F1, F2 = D_R else: D_B, D_U = _horizontal_bisection(1, a, b, I, Q, F1, F2, f1, f2, dom) if D_B[0] == 1: _, a, b, I, Q, F1, F2 = D_B else: _, a, b, I, Q, F1, F2 = D_U return ComplexInterval(a, b, I, Q, F1, F2, f1, f2, dom, self.conj)