1"""Real and complex root isolation and refinement algorithms. """ 2 3 4from sympy.polys.densearith import ( 5 dup_neg, dup_rshift, dup_rem) 6from sympy.polys.densebasic import ( 7 dup_LC, dup_TC, dup_degree, 8 dup_strip, dup_reverse, 9 dup_convert, 10 dup_terms_gcd) 11from sympy.polys.densetools import ( 12 dup_clear_denoms, 13 dup_mirror, dup_scale, dup_shift, 14 dup_transform, 15 dup_diff, 16 dup_eval, dmp_eval_in, 17 dup_sign_variations, 18 dup_real_imag) 19from sympy.polys.factortools import ( 20 dup_factor_list) 21from sympy.polys.polyerrors import ( 22 RefinementFailed, 23 DomainError) 24from sympy.polys.sqfreetools import ( 25 dup_sqf_part, dup_sqf_list) 26 27 28def dup_sturm(f, K): 29 """ 30 Computes the Sturm sequence of ``f`` in ``F[x]``. 31 32 Given a univariate, square-free polynomial ``f(x)`` returns the 33 associated Sturm sequence ``f_0(x), ..., f_n(x)`` defined by:: 34 35 f_0(x), f_1(x) = f(x), f'(x) 36 f_n = -rem(f_{n-2}(x), f_{n-1}(x)) 37 38 Examples 39 ======== 40 41 >>> from sympy.polys import ring, QQ 42 >>> R, x = ring("x", QQ) 43 44 >>> R.dup_sturm(x**3 - 2*x**2 + x - 3) 45 [x**3 - 2*x**2 + x - 3, 3*x**2 - 4*x + 1, 2/9*x + 25/9, -2079/4] 46 47 References 48 ========== 49 50 .. [1] [Davenport88]_ 51 52 """ 53 if not K.is_Field: 54 raise DomainError("can't compute Sturm sequence over %s" % K) 55 56 f = dup_sqf_part(f, K) 57 58 sturm = [f, dup_diff(f, 1, K)] 59 60 while sturm[-1]: 61 s = dup_rem(sturm[-2], sturm[-1], K) 62 sturm.append(dup_neg(s, K)) 63 64 return sturm[:-1] 65 66def dup_root_upper_bound(f, K): 67 """Compute the LMQ upper bound for the positive roots of `f`; 68 LMQ (Local Max Quadratic) was developed by Akritas-Strzebonski-Vigklas. 69 70 References 71 ========== 72 .. [1] Alkiviadis G. Akritas: "Linear and Quadratic Complexity Bounds on the 73 Values of the Positive Roots of Polynomials" 74 Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009. 75 """ 76 n, P = len(f), [] 77 t = n * [K.one] 78 if dup_LC(f, K) < 0: 79 f = dup_neg(f, K) 80 f = list(reversed(f)) 81 82 for i in range(0, n): 83 if f[i] >= 0: 84 continue 85 86 a, QL = K.log(-f[i], 2), [] 87 88 for j in range(i + 1, n): 89 90 if f[j] <= 0: 91 continue 92 93 q = t[j] + a - K.log(f[j], 2) 94 QL.append([q // (j - i) , j]) 95 96 if not QL: 97 continue 98 99 q = min(QL) 100 101 t[q[1]] = t[q[1]] + 1 102 103 P.append(q[0]) 104 105 if not P: 106 return None 107 else: 108 return K.get_field()(2)**(max(P) + 1) 109 110def dup_root_lower_bound(f, K): 111 """Compute the LMQ lower bound for the positive roots of `f`; 112 LMQ (Local Max Quadratic) was developed by Akritas-Strzebonski-Vigklas. 113 114 References 115 ========== 116 .. [1] Alkiviadis G. Akritas: "Linear and Quadratic Complexity Bounds on the 117 Values of the Positive Roots of Polynomials" 118 Journal of Universal Computer Science, Vol. 15, No. 3, 523-537, 2009. 119 """ 120 bound = dup_root_upper_bound(dup_reverse(f), K) 121 122 if bound is not None: 123 return 1/bound 124 else: 125 return None 126 127def _mobius_from_interval(I, field): 128 """Convert an open interval to a Mobius transform. """ 129 s, t = I 130 131 a, c = field.numer(s), field.denom(s) 132 b, d = field.numer(t), field.denom(t) 133 134 return a, b, c, d 135 136def _mobius_to_interval(M, field): 137 """Convert a Mobius transform to an open interval. """ 138 a, b, c, d = M 139 140 s, t = field(a, c), field(b, d) 141 142 if s <= t: 143 return (s, t) 144 else: 145 return (t, s) 146 147def dup_step_refine_real_root(f, M, K, fast=False): 148 """One step of positive real root refinement algorithm. """ 149 a, b, c, d = M 150 151 if a == b and c == d: 152 return f, (a, b, c, d) 153 154 A = dup_root_lower_bound(f, K) 155 156 if A is not None: 157 A = K(int(A)) 158 else: 159 A = K.zero 160 161 if fast and A > 16: 162 f = dup_scale(f, A, K) 163 a, c, A = A*a, A*c, K.one 164 165 if A >= K.one: 166 f = dup_shift(f, A, K) 167 b, d = A*a + b, A*c + d 168 169 if not dup_eval(f, K.zero, K): 170 return f, (b, b, d, d) 171 172 f, g = dup_shift(f, K.one, K), f 173 174 a1, b1, c1, d1 = a, a + b, c, c + d 175 176 if not dup_eval(f, K.zero, K): 177 return f, (b1, b1, d1, d1) 178 179 k = dup_sign_variations(f, K) 180 181 if k == 1: 182 a, b, c, d = a1, b1, c1, d1 183 else: 184 f = dup_shift(dup_reverse(g), K.one, K) 185 186 if not dup_eval(f, K.zero, K): 187 f = dup_rshift(f, 1, K) 188 189 a, b, c, d = b, a + b, d, c + d 190 191 return f, (a, b, c, d) 192 193def dup_inner_refine_real_root(f, M, K, eps=None, steps=None, disjoint=None, fast=False, mobius=False): 194 """Refine a positive root of `f` given a Mobius transform or an interval. """ 195 F = K.get_field() 196 197 if len(M) == 2: 198 a, b, c, d = _mobius_from_interval(M, F) 199 else: 200 a, b, c, d = M 201 202 while not c: 203 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, 204 d), K, fast=fast) 205 206 if eps is not None and steps is not None: 207 for i in range(0, steps): 208 if abs(F(a, c) - F(b, d)) >= eps: 209 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast) 210 else: 211 break 212 else: 213 if eps is not None: 214 while abs(F(a, c) - F(b, d)) >= eps: 215 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast) 216 217 if steps is not None: 218 for i in range(0, steps): 219 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast) 220 221 if disjoint is not None: 222 while True: 223 u, v = _mobius_to_interval((a, b, c, d), F) 224 225 if v <= disjoint or disjoint <= u: 226 break 227 else: 228 f, (a, b, c, d) = dup_step_refine_real_root(f, (a, b, c, d), K, fast=fast) 229 230 if not mobius: 231 return _mobius_to_interval((a, b, c, d), F) 232 else: 233 return f, (a, b, c, d) 234 235def dup_outer_refine_real_root(f, s, t, K, eps=None, steps=None, disjoint=None, fast=False): 236 """Refine a positive root of `f` given an interval `(s, t)`. """ 237 a, b, c, d = _mobius_from_interval((s, t), K.get_field()) 238 239 f = dup_transform(f, dup_strip([a, b]), 240 dup_strip([c, d]), K) 241 242 if dup_sign_variations(f, K) != 1: 243 raise RefinementFailed("there should be exactly one root in (%s, %s) interval" % (s, t)) 244 245 return dup_inner_refine_real_root(f, (a, b, c, d), K, eps=eps, steps=steps, disjoint=disjoint, fast=fast) 246 247def dup_refine_real_root(f, s, t, K, eps=None, steps=None, disjoint=None, fast=False): 248 """Refine real root's approximating interval to the given precision. """ 249 if K.is_QQ: 250 (_, f), K = dup_clear_denoms(f, K, convert=True), K.get_ring() 251 elif not K.is_ZZ: 252 raise DomainError("real root refinement not supported over %s" % K) 253 254 if s == t: 255 return (s, t) 256 257 if s > t: 258 s, t = t, s 259 260 negative = False 261 262 if s < 0: 263 if t <= 0: 264 f, s, t, negative = dup_mirror(f, K), -t, -s, True 265 else: 266 raise ValueError("can't refine a real root in (%s, %s)" % (s, t)) 267 268 if negative and disjoint is not None: 269 if disjoint < 0: 270 disjoint = -disjoint 271 else: 272 disjoint = None 273 274 s, t = dup_outer_refine_real_root( 275 f, s, t, K, eps=eps, steps=steps, disjoint=disjoint, fast=fast) 276 277 if negative: 278 return (-t, -s) 279 else: 280 return ( s, t) 281 282def dup_inner_isolate_real_roots(f, K, eps=None, fast=False): 283 """Internal function for isolation positive roots up to given precision. 284 285 References 286 ========== 287 1. Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative Study of Two Real Root 288 Isolation Methods . Nonlinear Analysis: Modelling and Control, Vol. 10, No. 4, 297-304, 2005. 289 2. Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. Vigklas: Improving the 290 Performance of the Continued Fractions Method Using new Bounds of Positive Roots. Nonlinear 291 Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008. 292 """ 293 a, b, c, d = K.one, K.zero, K.zero, K.one 294 295 k = dup_sign_variations(f, K) 296 297 if k == 0: 298 return [] 299 if k == 1: 300 roots = [dup_inner_refine_real_root( 301 f, (a, b, c, d), K, eps=eps, fast=fast, mobius=True)] 302 else: 303 roots, stack = [], [(a, b, c, d, f, k)] 304 305 while stack: 306 a, b, c, d, f, k = stack.pop() 307 308 A = dup_root_lower_bound(f, K) 309 310 if A is not None: 311 A = K(int(A)) 312 else: 313 A = K.zero 314 315 if fast and A > 16: 316 f = dup_scale(f, A, K) 317 a, c, A = A*a, A*c, K.one 318 319 if A >= K.one: 320 f = dup_shift(f, A, K) 321 b, d = A*a + b, A*c + d 322 323 if not dup_TC(f, K): 324 roots.append((f, (b, b, d, d))) 325 f = dup_rshift(f, 1, K) 326 327 k = dup_sign_variations(f, K) 328 329 if k == 0: 330 continue 331 if k == 1: 332 roots.append(dup_inner_refine_real_root( 333 f, (a, b, c, d), K, eps=eps, fast=fast, mobius=True)) 334 continue 335 336 f1 = dup_shift(f, K.one, K) 337 338 a1, b1, c1, d1, r = a, a + b, c, c + d, 0 339 340 if not dup_TC(f1, K): 341 roots.append((f1, (b1, b1, d1, d1))) 342 f1, r = dup_rshift(f1, 1, K), 1 343 344 k1 = dup_sign_variations(f1, K) 345 k2 = k - k1 - r 346 347 a2, b2, c2, d2 = b, a + b, d, c + d 348 349 if k2 > 1: 350 f2 = dup_shift(dup_reverse(f), K.one, K) 351 352 if not dup_TC(f2, K): 353 f2 = dup_rshift(f2, 1, K) 354 355 k2 = dup_sign_variations(f2, K) 356 else: 357 f2 = None 358 359 if k1 < k2: 360 a1, a2, b1, b2 = a2, a1, b2, b1 361 c1, c2, d1, d2 = c2, c1, d2, d1 362 f1, f2, k1, k2 = f2, f1, k2, k1 363 364 if not k1: 365 continue 366 367 if f1 is None: 368 f1 = dup_shift(dup_reverse(f), K.one, K) 369 370 if not dup_TC(f1, K): 371 f1 = dup_rshift(f1, 1, K) 372 373 if k1 == 1: 374 roots.append(dup_inner_refine_real_root( 375 f1, (a1, b1, c1, d1), K, eps=eps, fast=fast, mobius=True)) 376 else: 377 stack.append((a1, b1, c1, d1, f1, k1)) 378 379 if not k2: 380 continue 381 382 if f2 is None: 383 f2 = dup_shift(dup_reverse(f), K.one, K) 384 385 if not dup_TC(f2, K): 386 f2 = dup_rshift(f2, 1, K) 387 388 if k2 == 1: 389 roots.append(dup_inner_refine_real_root( 390 f2, (a2, b2, c2, d2), K, eps=eps, fast=fast, mobius=True)) 391 else: 392 stack.append((a2, b2, c2, d2, f2, k2)) 393 394 return roots 395 396def _discard_if_outside_interval(f, M, inf, sup, K, negative, fast, mobius): 397 """Discard an isolating interval if outside ``(inf, sup)``. """ 398 F = K.get_field() 399 400 while True: 401 u, v = _mobius_to_interval(M, F) 402 403 if negative: 404 u, v = -v, -u 405 406 if (inf is None or u >= inf) and (sup is None or v <= sup): 407 if not mobius: 408 return u, v 409 else: 410 return f, M 411 elif (sup is not None and u > sup) or (inf is not None and v < inf): 412 return None 413 else: 414 f, M = dup_step_refine_real_root(f, M, K, fast=fast) 415 416def dup_inner_isolate_positive_roots(f, K, eps=None, inf=None, sup=None, fast=False, mobius=False): 417 """Iteratively compute disjoint positive root isolation intervals. """ 418 if sup is not None and sup < 0: 419 return [] 420 421 roots = dup_inner_isolate_real_roots(f, K, eps=eps, fast=fast) 422 423 F, results = K.get_field(), [] 424 425 if inf is not None or sup is not None: 426 for f, M in roots: 427 result = _discard_if_outside_interval(f, M, inf, sup, K, False, fast, mobius) 428 429 if result is not None: 430 results.append(result) 431 elif not mobius: 432 for f, M in roots: 433 u, v = _mobius_to_interval(M, F) 434 results.append((u, v)) 435 else: 436 results = roots 437 438 return results 439 440def dup_inner_isolate_negative_roots(f, K, inf=None, sup=None, eps=None, fast=False, mobius=False): 441 """Iteratively compute disjoint negative root isolation intervals. """ 442 if inf is not None and inf >= 0: 443 return [] 444 445 roots = dup_inner_isolate_real_roots(dup_mirror(f, K), K, eps=eps, fast=fast) 446 447 F, results = K.get_field(), [] 448 449 if inf is not None or sup is not None: 450 for f, M in roots: 451 result = _discard_if_outside_interval(f, M, inf, sup, K, True, fast, mobius) 452 453 if result is not None: 454 results.append(result) 455 elif not mobius: 456 for f, M in roots: 457 u, v = _mobius_to_interval(M, F) 458 results.append((-v, -u)) 459 else: 460 results = roots 461 462 return results 463 464def _isolate_zero(f, K, inf, sup, basis=False, sqf=False): 465 """Handle special case of CF algorithm when ``f`` is homogeneous. """ 466 j, f = dup_terms_gcd(f, K) 467 468 if j > 0: 469 F = K.get_field() 470 471 if (inf is None or inf <= 0) and (sup is None or 0 <= sup): 472 if not sqf: 473 if not basis: 474 return [((F.zero, F.zero), j)], f 475 else: 476 return [((F.zero, F.zero), j, [K.one, K.zero])], f 477 else: 478 return [(F.zero, F.zero)], f 479 480 return [], f 481 482def dup_isolate_real_roots_sqf(f, K, eps=None, inf=None, sup=None, fast=False, blackbox=False): 483 """Isolate real roots of a square-free polynomial using the Vincent-Akritas-Strzebonski (VAS) CF approach. 484 485 References 486 ========== 487 .. [1] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative 488 Study of Two Real Root Isolation Methods. Nonlinear Analysis: 489 Modelling and Control, Vol. 10, No. 4, 297-304, 2005. 490 .. [2] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. 491 Vigklas: Improving the Performance of the Continued Fractions 492 Method Using New Bounds of Positive Roots. Nonlinear Analysis: 493 Modelling and Control, Vol. 13, No. 3, 265-279, 2008. 494 495 """ 496 if K.is_QQ: 497 (_, f), K = dup_clear_denoms(f, K, convert=True), K.get_ring() 498 elif not K.is_ZZ: 499 raise DomainError("isolation of real roots not supported over %s" % K) 500 501 if dup_degree(f) <= 0: 502 return [] 503 504 I_zero, f = _isolate_zero(f, K, inf, sup, basis=False, sqf=True) 505 506 I_neg = dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast) 507 I_pos = dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast) 508 509 roots = sorted(I_neg + I_zero + I_pos) 510 511 if not blackbox: 512 return roots 513 else: 514 return [ RealInterval((a, b), f, K) for (a, b) in roots ] 515 516def dup_isolate_real_roots(f, K, eps=None, inf=None, sup=None, basis=False, fast=False): 517 """Isolate real roots using Vincent-Akritas-Strzebonski (VAS) continued fractions approach. 518 519 References 520 ========== 521 522 .. [1] Alkiviadis G. Akritas and Adam W. Strzebonski: A Comparative 523 Study of Two Real Root Isolation Methods. Nonlinear Analysis: 524 Modelling and Control, Vol. 10, No. 4, 297-304, 2005. 525 .. [2] Alkiviadis G. Akritas, Adam W. Strzebonski and Panagiotis S. 526 Vigklas: Improving the Performance of the Continued Fractions 527 Method Using New Bounds of Positive Roots. 528 Nonlinear Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008. 529 530 """ 531 if K.is_QQ: 532 (_, f), K = dup_clear_denoms(f, K, convert=True), K.get_ring() 533 elif not K.is_ZZ: 534 raise DomainError("isolation of real roots not supported over %s" % K) 535 536 if dup_degree(f) <= 0: 537 return [] 538 539 I_zero, f = _isolate_zero(f, K, inf, sup, basis=basis, sqf=False) 540 541 _, factors = dup_sqf_list(f, K) 542 543 if len(factors) == 1: 544 ((f, k),) = factors 545 546 I_neg = dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast) 547 I_pos = dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast) 548 549 I_neg = [ ((u, v), k) for u, v in I_neg ] 550 I_pos = [ ((u, v), k) for u, v in I_pos ] 551 else: 552 I_neg, I_pos = _real_isolate_and_disjoin(factors, K, 553 eps=eps, inf=inf, sup=sup, basis=basis, fast=fast) 554 555 return sorted(I_neg + I_zero + I_pos) 556 557def dup_isolate_real_roots_list(polys, K, eps=None, inf=None, sup=None, strict=False, basis=False, fast=False): 558 """Isolate real roots of a list of square-free polynomial using Vincent-Akritas-Strzebonski (VAS) CF approach. 559 560 References 561 ========== 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. 569 Nonlinear Analysis: Modelling and Control, Vol. 13, No. 3, 265-279, 2008. 570 571 """ 572 if K.is_QQ: 573 K, F, polys = K.get_ring(), K, polys[:] 574 575 for i, p in enumerate(polys): 576 polys[i] = dup_clear_denoms(p, F, K, convert=True)[1] 577 elif not K.is_ZZ: 578 raise DomainError("isolation of real roots not supported over %s" % K) 579 580 zeros, factors_dict = False, {} 581 582 if (inf is None or inf <= 0) and (sup is None or 0 <= sup): 583 zeros, zero_indices = True, {} 584 585 for i, p in enumerate(polys): 586 j, p = dup_terms_gcd(p, K) 587 588 if zeros and j > 0: 589 zero_indices[i] = j 590 591 for f, k in dup_factor_list(p, K)[1]: 592 f = tuple(f) 593 594 if f not in factors_dict: 595 factors_dict[f] = {i: k} 596 else: 597 factors_dict[f][i] = k 598 599 factors_list = [] 600 601 for f, indices in factors_dict.items(): 602 factors_list.append((list(f), indices)) 603 604 I_neg, I_pos = _real_isolate_and_disjoin(factors_list, K, eps=eps, 605 inf=inf, sup=sup, strict=strict, basis=basis, fast=fast) 606 607 F = K.get_field() 608 609 if not zeros or not zero_indices: 610 I_zero = [] 611 else: 612 if not basis: 613 I_zero = [((F.zero, F.zero), zero_indices)] 614 else: 615 I_zero = [((F.zero, F.zero), zero_indices, [K.one, K.zero])] 616 617 return sorted(I_neg + I_zero + I_pos) 618 619def _disjoint_p(M, N, strict=False): 620 """Check if Mobius transforms define disjoint intervals. """ 621 a1, b1, c1, d1 = M 622 a2, b2, c2, d2 = N 623 624 a1d1, b1c1 = a1*d1, b1*c1 625 a2d2, b2c2 = a2*d2, b2*c2 626 627 if a1d1 == b1c1 and a2d2 == b2c2: 628 return True 629 630 if a1d1 > b1c1: 631 a1, c1, b1, d1 = b1, d1, a1, c1 632 633 if a2d2 > b2c2: 634 a2, c2, b2, d2 = b2, d2, a2, c2 635 636 if not strict: 637 return a2*d1 >= c2*b1 or b2*c1 <= d2*a1 638 else: 639 return a2*d1 > c2*b1 or b2*c1 < d2*a1 640 641def _real_isolate_and_disjoin(factors, K, eps=None, inf=None, sup=None, strict=False, basis=False, fast=False): 642 """Isolate real roots of a list of polynomials and disjoin intervals. """ 643 I_pos, I_neg = [], [] 644 645 for i, (f, k) in enumerate(factors): 646 for F, M in dup_inner_isolate_positive_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast, mobius=True): 647 I_pos.append((F, M, k, f)) 648 649 for G, N in dup_inner_isolate_negative_roots(f, K, eps=eps, inf=inf, sup=sup, fast=fast, mobius=True): 650 I_neg.append((G, N, k, f)) 651 652 for i, (f, M, k, F) in enumerate(I_pos): 653 for j, (g, N, m, G) in enumerate(I_pos[i + 1:]): 654 while not _disjoint_p(M, N, strict=strict): 655 f, M = dup_inner_refine_real_root(f, M, K, steps=1, fast=fast, mobius=True) 656 g, N = dup_inner_refine_real_root(g, N, K, steps=1, fast=fast, mobius=True) 657 658 I_pos[i + j + 1] = (g, N, m, G) 659 660 I_pos[i] = (f, M, k, F) 661 662 for i, (f, M, k, F) in enumerate(I_neg): 663 for j, (g, N, m, G) in enumerate(I_neg[i + 1:]): 664 while not _disjoint_p(M, N, strict=strict): 665 f, M = dup_inner_refine_real_root(f, M, K, steps=1, fast=fast, mobius=True) 666 g, N = dup_inner_refine_real_root(g, N, K, steps=1, fast=fast, mobius=True) 667 668 I_neg[i + j + 1] = (g, N, m, G) 669 670 I_neg[i] = (f, M, k, F) 671 672 if strict: 673 for i, (f, M, k, F) in enumerate(I_neg): 674 if not M[0]: 675 while not M[0]: 676 f, M = dup_inner_refine_real_root(f, M, K, steps=1, fast=fast, mobius=True) 677 678 I_neg[i] = (f, M, k, F) 679 break 680 681 for j, (g, N, m, G) in enumerate(I_pos): 682 if not N[0]: 683 while not N[0]: 684 g, N = dup_inner_refine_real_root(g, N, K, steps=1, fast=fast, mobius=True) 685 686 I_pos[j] = (g, N, m, G) 687 break 688 689 field = K.get_field() 690 691 I_neg = [ (_mobius_to_interval(M, field), k, f) for (_, M, k, f) in I_neg ] 692 I_pos = [ (_mobius_to_interval(M, field), k, f) for (_, M, k, f) in I_pos ] 693 694 if not basis: 695 I_neg = [ ((-v, -u), k) for ((u, v), k, _) in I_neg ] 696 I_pos = [ (( u, v), k) for ((u, v), k, _) in I_pos ] 697 else: 698 I_neg = [ ((-v, -u), k, f) for ((u, v), k, f) in I_neg ] 699 I_pos = [ (( u, v), k, f) for ((u, v), k, f) in I_pos ] 700 701 return I_neg, I_pos 702 703def dup_count_real_roots(f, K, inf=None, sup=None): 704 """Returns the number of distinct real roots of ``f`` in ``[inf, sup]``. """ 705 if dup_degree(f) <= 0: 706 return 0 707 708 if not K.is_Field: 709 R, K = K, K.get_field() 710 f = dup_convert(f, R, K) 711 712 sturm = dup_sturm(f, K) 713 714 if inf is None: 715 signs_inf = dup_sign_variations([ dup_LC(s, K)*(-1)**dup_degree(s) for s in sturm ], K) 716 else: 717 signs_inf = dup_sign_variations([ dup_eval(s, inf, K) for s in sturm ], K) 718 719 if sup is None: 720 signs_sup = dup_sign_variations([ dup_LC(s, K) for s in sturm ], K) 721 else: 722 signs_sup = dup_sign_variations([ dup_eval(s, sup, K) for s in sturm ], K) 723 724 count = abs(signs_inf - signs_sup) 725 726 if inf is not None and not dup_eval(f, inf, K): 727 count += 1 728 729 return count 730 731OO = 'OO' # Origin of (re, im) coordinate system 732 733Q1 = 'Q1' # Quadrant #1 (++): re > 0 and im > 0 734Q2 = 'Q2' # Quadrant #2 (-+): re < 0 and im > 0 735Q3 = 'Q3' # Quadrant #3 (--): re < 0 and im < 0 736Q4 = 'Q4' # Quadrant #4 (+-): re > 0 and im < 0 737 738A1 = 'A1' # Axis #1 (+0): re > 0 and im = 0 739A2 = 'A2' # Axis #2 (0+): re = 0 and im > 0 740A3 = 'A3' # Axis #3 (-0): re < 0 and im = 0 741A4 = 'A4' # Axis #4 (0-): re = 0 and im < 0 742 743_rules_simple = { 744 # Q --> Q (same) => no change 745 (Q1, Q1): 0, 746 (Q2, Q2): 0, 747 (Q3, Q3): 0, 748 (Q4, Q4): 0, 749 750 # A -- CCW --> Q => +1/4 (CCW) 751 (A1, Q1): 1, 752 (A2, Q2): 1, 753 (A3, Q3): 1, 754 (A4, Q4): 1, 755 756 # A -- CW --> Q => -1/4 (CCW) 757 (A1, Q4): 2, 758 (A2, Q1): 2, 759 (A3, Q2): 2, 760 (A4, Q3): 2, 761 762 # Q -- CCW --> A => +1/4 (CCW) 763 (Q1, A2): 3, 764 (Q2, A3): 3, 765 (Q3, A4): 3, 766 (Q4, A1): 3, 767 768 # Q -- CW --> A => -1/4 (CCW) 769 (Q1, A1): 4, 770 (Q2, A2): 4, 771 (Q3, A3): 4, 772 (Q4, A4): 4, 773 774 # Q -- CCW --> Q => +1/2 (CCW) 775 (Q1, Q2): +5, 776 (Q2, Q3): +5, 777 (Q3, Q4): +5, 778 (Q4, Q1): +5, 779 780 # Q -- CW --> Q => -1/2 (CW) 781 (Q1, Q4): -5, 782 (Q2, Q1): -5, 783 (Q3, Q2): -5, 784 (Q4, Q3): -5, 785} 786 787_rules_ambiguous = { 788 # A -- CCW --> Q => { +1/4 (CCW), -9/4 (CW) } 789 (A1, OO, Q1): -1, 790 (A2, OO, Q2): -1, 791 (A3, OO, Q3): -1, 792 (A4, OO, Q4): -1, 793 794 # A -- CW --> Q => { -1/4 (CCW), +7/4 (CW) } 795 (A1, OO, Q4): -2, 796 (A2, OO, Q1): -2, 797 (A3, OO, Q2): -2, 798 (A4, OO, Q3): -2, 799 800 # Q -- CCW --> A => { +1/4 (CCW), -9/4 (CW) } 801 (Q1, OO, A2): -3, 802 (Q2, OO, A3): -3, 803 (Q3, OO, A4): -3, 804 (Q4, OO, A1): -3, 805 806 # Q -- CW --> A => { -1/4 (CCW), +7/4 (CW) } 807 (Q1, OO, A1): -4, 808 (Q2, OO, A2): -4, 809 (Q3, OO, A3): -4, 810 (Q4, OO, A4): -4, 811 812 # A -- OO --> A => { +1 (CCW), -1 (CW) } 813 (A1, A3): 7, 814 (A2, A4): 7, 815 (A3, A1): 7, 816 (A4, A2): 7, 817 818 (A1, OO, A3): 7, 819 (A2, OO, A4): 7, 820 (A3, OO, A1): 7, 821 (A4, OO, A2): 7, 822 823 # Q -- DIA --> Q => { +1 (CCW), -1 (CW) } 824 (Q1, Q3): 8, 825 (Q2, Q4): 8, 826 (Q3, Q1): 8, 827 (Q4, Q2): 8, 828 829 (Q1, OO, Q3): 8, 830 (Q2, OO, Q4): 8, 831 (Q3, OO, Q1): 8, 832 (Q4, OO, Q2): 8, 833 834 # A --- R ---> A => { +1/2 (CCW), -3/2 (CW) } 835 (A1, A2): 9, 836 (A2, A3): 9, 837 (A3, A4): 9, 838 (A4, A1): 9, 839 840 (A1, OO, A2): 9, 841 (A2, OO, A3): 9, 842 (A3, OO, A4): 9, 843 (A4, OO, A1): 9, 844 845 # A --- L ---> A => { +3/2 (CCW), -1/2 (CW) } 846 (A1, A4): 10, 847 (A2, A1): 10, 848 (A3, A2): 10, 849 (A4, A3): 10, 850 851 (A1, OO, A4): 10, 852 (A2, OO, A1): 10, 853 (A3, OO, A2): 10, 854 (A4, OO, A3): 10, 855 856 # Q --- 1 ---> A => { +3/4 (CCW), -5/4 (CW) } 857 (Q1, A3): 11, 858 (Q2, A4): 11, 859 (Q3, A1): 11, 860 (Q4, A2): 11, 861 862 (Q1, OO, A3): 11, 863 (Q2, OO, A4): 11, 864 (Q3, OO, A1): 11, 865 (Q4, OO, A2): 11, 866 867 # Q --- 2 ---> A => { +5/4 (CCW), -3/4 (CW) } 868 (Q1, A4): 12, 869 (Q2, A1): 12, 870 (Q3, A2): 12, 871 (Q4, A3): 12, 872 873 (Q1, OO, A4): 12, 874 (Q2, OO, A1): 12, 875 (Q3, OO, A2): 12, 876 (Q4, OO, A3): 12, 877 878 # A --- 1 ---> Q => { +5/4 (CCW), -3/4 (CW) } 879 (A1, Q3): 13, 880 (A2, Q4): 13, 881 (A3, Q1): 13, 882 (A4, Q2): 13, 883 884 (A1, OO, Q3): 13, 885 (A2, OO, Q4): 13, 886 (A3, OO, Q1): 13, 887 (A4, OO, Q2): 13, 888 889 # A --- 2 ---> Q => { +3/4 (CCW), -5/4 (CW) } 890 (A1, Q2): 14, 891 (A2, Q3): 14, 892 (A3, Q4): 14, 893 (A4, Q1): 14, 894 895 (A1, OO, Q2): 14, 896 (A2, OO, Q3): 14, 897 (A3, OO, Q4): 14, 898 (A4, OO, Q1): 14, 899 900 # Q --> OO --> Q => { +1/2 (CCW), -3/2 (CW) } 901 (Q1, OO, Q2): 15, 902 (Q2, OO, Q3): 15, 903 (Q3, OO, Q4): 15, 904 (Q4, OO, Q1): 15, 905 906 # Q --> OO --> Q => { +3/2 (CCW), -1/2 (CW) } 907 (Q1, OO, Q4): 16, 908 (Q2, OO, Q1): 16, 909 (Q3, OO, Q2): 16, 910 (Q4, OO, Q3): 16, 911 912 # A --> OO --> A => { +2 (CCW), 0 (CW) } 913 (A1, OO, A1): 17, 914 (A2, OO, A2): 17, 915 (A3, OO, A3): 17, 916 (A4, OO, A4): 17, 917 918 # Q --> OO --> Q => { +2 (CCW), 0 (CW) } 919 (Q1, OO, Q1): 18, 920 (Q2, OO, Q2): 18, 921 (Q3, OO, Q3): 18, 922 (Q4, OO, Q4): 18, 923} 924 925_values = { 926 0: [( 0, 1)], 927 1: [(+1, 4)], 928 2: [(-1, 4)], 929 3: [(+1, 4)], 930 4: [(-1, 4)], 931 -1: [(+9, 4), (+1, 4)], 932 -2: [(+7, 4), (-1, 4)], 933 -3: [(+9, 4), (+1, 4)], 934 -4: [(+7, 4), (-1, 4)], 935 +5: [(+1, 2)], 936 -5: [(-1, 2)], 937 7: [(+1, 1), (-1, 1)], 938 8: [(+1, 1), (-1, 1)], 939 9: [(+1, 2), (-3, 2)], 940 10: [(+3, 2), (-1, 2)], 941 11: [(+3, 4), (-5, 4)], 942 12: [(+5, 4), (-3, 4)], 943 13: [(+5, 4), (-3, 4)], 944 14: [(+3, 4), (-5, 4)], 945 15: [(+1, 2), (-3, 2)], 946 16: [(+3, 2), (-1, 2)], 947 17: [(+2, 1), ( 0, 1)], 948 18: [(+2, 1), ( 0, 1)], 949} 950 951def _classify_point(re, im): 952 """Return the half-axis (or origin) on which (re, im) point is located. """ 953 if not re and not im: 954 return OO 955 956 if not re: 957 if im > 0: 958 return A2 959 else: 960 return A4 961 elif not im: 962 if re > 0: 963 return A1 964 else: 965 return A3 966 967def _intervals_to_quadrants(intervals, f1, f2, s, t, F): 968 """Generate a sequence of extended quadrants from a list of critical points. """ 969 if not intervals: 970 return [] 971 972 Q = [] 973 974 if not f1: 975 (a, b), _, _ = intervals[0] 976 977 if a == b == s: 978 if len(intervals) == 1: 979 if dup_eval(f2, t, F) > 0: 980 return [OO, A2] 981 else: 982 return [OO, A4] 983 else: 984 (a, _), _, _ = intervals[1] 985 986 if dup_eval(f2, (s + a)/2, F) > 0: 987 Q.extend([OO, A2]) 988 f2_sgn = +1 989 else: 990 Q.extend([OO, A4]) 991 f2_sgn = -1 992 993 intervals = intervals[1:] 994 else: 995 if dup_eval(f2, s, F) > 0: 996 Q.append(A2) 997 f2_sgn = +1 998 else: 999 Q.append(A4) 1000 f2_sgn = -1 1001 1002 for (a, _), indices, _ in intervals: 1003 Q.append(OO) 1004 1005 if indices[1] % 2 == 1: 1006 f2_sgn = -f2_sgn 1007 1008 if a != t: 1009 if f2_sgn > 0: 1010 Q.append(A2) 1011 else: 1012 Q.append(A4) 1013 1014 return Q 1015 1016 if not f2: 1017 (a, b), _, _ = intervals[0] 1018 1019 if a == b == s: 1020 if len(intervals) == 1: 1021 if dup_eval(f1, t, F) > 0: 1022 return [OO, A1] 1023 else: 1024 return [OO, A3] 1025 else: 1026 (a, _), _, _ = intervals[1] 1027 1028 if dup_eval(f1, (s + a)/2, F) > 0: 1029 Q.extend([OO, A1]) 1030 f1_sgn = +1 1031 else: 1032 Q.extend([OO, A3]) 1033 f1_sgn = -1 1034 1035 intervals = intervals[1:] 1036 else: 1037 if dup_eval(f1, s, F) > 0: 1038 Q.append(A1) 1039 f1_sgn = +1 1040 else: 1041 Q.append(A3) 1042 f1_sgn = -1 1043 1044 for (a, _), indices, _ in intervals: 1045 Q.append(OO) 1046 1047 if indices[0] % 2 == 1: 1048 f1_sgn = -f1_sgn 1049 1050 if a != t: 1051 if f1_sgn > 0: 1052 Q.append(A1) 1053 else: 1054 Q.append(A3) 1055 1056 return Q 1057 1058 re = dup_eval(f1, s, F) 1059 im = dup_eval(f2, s, F) 1060 1061 if not re or not im: 1062 Q.append(_classify_point(re, im)) 1063 1064 if len(intervals) == 1: 1065 re = dup_eval(f1, t, F) 1066 im = dup_eval(f2, t, F) 1067 else: 1068 (a, _), _, _ = intervals[1] 1069 1070 re = dup_eval(f1, (s + a)/2, F) 1071 im = dup_eval(f2, (s + a)/2, F) 1072 1073 intervals = intervals[1:] 1074 1075 if re > 0: 1076 f1_sgn = +1 1077 else: 1078 f1_sgn = -1 1079 1080 if im > 0: 1081 f2_sgn = +1 1082 else: 1083 f2_sgn = -1 1084 1085 sgn = { 1086 (+1, +1): Q1, 1087 (-1, +1): Q2, 1088 (-1, -1): Q3, 1089 (+1, -1): Q4, 1090 } 1091 1092 Q.append(sgn[(f1_sgn, f2_sgn)]) 1093 1094 for (a, b), indices, _ in intervals: 1095 if a == b: 1096 re = dup_eval(f1, a, F) 1097 im = dup_eval(f2, a, F) 1098 1099 cls = _classify_point(re, im) 1100 1101 if cls is not None: 1102 Q.append(cls) 1103 1104 if 0 in indices: 1105 if indices[0] % 2 == 1: 1106 f1_sgn = -f1_sgn 1107 1108 if 1 in indices: 1109 if indices[1] % 2 == 1: 1110 f2_sgn = -f2_sgn 1111 1112 if not (a == b and b == t): 1113 Q.append(sgn[(f1_sgn, f2_sgn)]) 1114 1115 return Q 1116 1117def _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=None): 1118 """Transform sequences of quadrants to a sequence of rules. """ 1119 if exclude is True: 1120 edges = [1, 1, 0, 0] 1121 1122 corners = { 1123 (0, 1): 1, 1124 (1, 2): 1, 1125 (2, 3): 0, 1126 (3, 0): 1, 1127 } 1128 else: 1129 edges = [0, 0, 0, 0] 1130 1131 corners = { 1132 (0, 1): 0, 1133 (1, 2): 0, 1134 (2, 3): 0, 1135 (3, 0): 0, 1136 } 1137 1138 if exclude is not None and exclude is not True: 1139 exclude = set(exclude) 1140 1141 for i, edge in enumerate(['S', 'E', 'N', 'W']): 1142 if edge in exclude: 1143 edges[i] = 1 1144 1145 for i, corner in enumerate(['SW', 'SE', 'NE', 'NW']): 1146 if corner in exclude: 1147 corners[((i - 1) % 4, i)] = 1 1148 1149 QQ, rules = [Q_L1, Q_L2, Q_L3, Q_L4], [] 1150 1151 for i, Q in enumerate(QQ): 1152 if not Q: 1153 continue 1154 1155 if Q[-1] == OO: 1156 Q = Q[:-1] 1157 1158 if Q[0] == OO: 1159 j, Q = (i - 1) % 4, Q[1:] 1160 qq = (QQ[j][-2], OO, Q[0]) 1161 1162 if qq in _rules_ambiguous: 1163 rules.append((_rules_ambiguous[qq], corners[(j, i)])) 1164 else: 1165 raise NotImplementedError("3 element rule (corner): " + str(qq)) 1166 1167 q1, k = Q[0], 1 1168 1169 while k < len(Q): 1170 q2, k = Q[k], k + 1 1171 1172 if q2 != OO: 1173 qq = (q1, q2) 1174 1175 if qq in _rules_simple: 1176 rules.append((_rules_simple[qq], 0)) 1177 elif qq in _rules_ambiguous: 1178 rules.append((_rules_ambiguous[qq], edges[i])) 1179 else: 1180 raise NotImplementedError("2 element rule (inside): " + str(qq)) 1181 else: 1182 qq, k = (q1, q2, Q[k]), k + 1 1183 1184 if qq in _rules_ambiguous: 1185 rules.append((_rules_ambiguous[qq], edges[i])) 1186 else: 1187 raise NotImplementedError("3 element rule (edge): " + str(qq)) 1188 1189 q1 = qq[-1] 1190 1191 return rules 1192 1193def _reverse_intervals(intervals): 1194 """Reverse intervals for traversal from right to left and from top to bottom. """ 1195 return [ ((b, a), indices, f) for (a, b), indices, f in reversed(intervals) ] 1196 1197def _winding_number(T, field): 1198 """Compute the winding number of the input polynomial, i.e. the number of roots. """ 1199 return int(sum([ field(*_values[t][i]) for t, i in T ]) / field(2)) 1200 1201def dup_count_complex_roots(f, K, inf=None, sup=None, exclude=None): 1202 """Count all roots in [u + v*I, s + t*I] rectangle using Collins-Krandick algorithm. """ 1203 if not K.is_ZZ and not K.is_QQ: 1204 raise DomainError("complex root counting is not supported over %s" % K) 1205 1206 if K.is_ZZ: 1207 R, F = K, K.get_field() 1208 else: 1209 R, F = K.get_ring(), K 1210 1211 f = dup_convert(f, K, F) 1212 1213 if inf is None or sup is None: 1214 _, lc = dup_degree(f), abs(dup_LC(f, F)) 1215 B = 2*max([ F.quo(abs(c), lc) for c in f ]) 1216 1217 if inf is None: 1218 (u, v) = (-B, -B) 1219 else: 1220 (u, v) = inf 1221 1222 if sup is None: 1223 (s, t) = (+B, +B) 1224 else: 1225 (s, t) = sup 1226 1227 f1, f2 = dup_real_imag(f, F) 1228 1229 f1L1F = dmp_eval_in(f1, v, 1, 1, F) 1230 f2L1F = dmp_eval_in(f2, v, 1, 1, F) 1231 1232 _, f1L1R = dup_clear_denoms(f1L1F, F, R, convert=True) 1233 _, f2L1R = dup_clear_denoms(f2L1F, F, R, convert=True) 1234 1235 f1L2F = dmp_eval_in(f1, s, 0, 1, F) 1236 f2L2F = dmp_eval_in(f2, s, 0, 1, F) 1237 1238 _, f1L2R = dup_clear_denoms(f1L2F, F, R, convert=True) 1239 _, f2L2R = dup_clear_denoms(f2L2F, F, R, convert=True) 1240 1241 f1L3F = dmp_eval_in(f1, t, 1, 1, F) 1242 f2L3F = dmp_eval_in(f2, t, 1, 1, F) 1243 1244 _, f1L3R = dup_clear_denoms(f1L3F, F, R, convert=True) 1245 _, f2L3R = dup_clear_denoms(f2L3F, F, R, convert=True) 1246 1247 f1L4F = dmp_eval_in(f1, u, 0, 1, F) 1248 f2L4F = dmp_eval_in(f2, u, 0, 1, F) 1249 1250 _, f1L4R = dup_clear_denoms(f1L4F, F, R, convert=True) 1251 _, f2L4R = dup_clear_denoms(f2L4F, F, R, convert=True) 1252 1253 S_L1 = [f1L1R, f2L1R] 1254 S_L2 = [f1L2R, f2L2R] 1255 S_L3 = [f1L3R, f2L3R] 1256 S_L4 = [f1L4R, f2L4R] 1257 1258 I_L1 = dup_isolate_real_roots_list(S_L1, R, inf=u, sup=s, fast=True, basis=True, strict=True) 1259 I_L2 = dup_isolate_real_roots_list(S_L2, R, inf=v, sup=t, fast=True, basis=True, strict=True) 1260 I_L3 = dup_isolate_real_roots_list(S_L3, R, inf=u, sup=s, fast=True, basis=True, strict=True) 1261 I_L4 = dup_isolate_real_roots_list(S_L4, R, inf=v, sup=t, fast=True, basis=True, strict=True) 1262 1263 I_L3 = _reverse_intervals(I_L3) 1264 I_L4 = _reverse_intervals(I_L4) 1265 1266 Q_L1 = _intervals_to_quadrants(I_L1, f1L1F, f2L1F, u, s, F) 1267 Q_L2 = _intervals_to_quadrants(I_L2, f1L2F, f2L2F, v, t, F) 1268 Q_L3 = _intervals_to_quadrants(I_L3, f1L3F, f2L3F, s, u, F) 1269 Q_L4 = _intervals_to_quadrants(I_L4, f1L4F, f2L4F, t, v, F) 1270 1271 T = _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=exclude) 1272 1273 return _winding_number(T, F) 1274 1275def _vertical_bisection(N, a, b, I, Q, F1, F2, f1, f2, F): 1276 """Vertical bisection step in Collins-Krandick root isolation algorithm. """ 1277 (u, v), (s, t) = a, b 1278 1279 I_L1, I_L2, I_L3, I_L4 = I 1280 Q_L1, Q_L2, Q_L3, Q_L4 = Q 1281 1282 f1L1F, f1L2F, f1L3F, f1L4F = F1 1283 f2L1F, f2L2F, f2L3F, f2L4F = F2 1284 1285 x = (u + s) / 2 1286 1287 f1V = dmp_eval_in(f1, x, 0, 1, F) 1288 f2V = dmp_eval_in(f2, x, 0, 1, F) 1289 1290 I_V = dup_isolate_real_roots_list([f1V, f2V], F, inf=v, sup=t, fast=True, strict=True, basis=True) 1291 1292 I_L1_L, I_L1_R = [], [] 1293 I_L2_L, I_L2_R = I_V, I_L2 1294 I_L3_L, I_L3_R = [], [] 1295 I_L4_L, I_L4_R = I_L4, _reverse_intervals(I_V) 1296 1297 for I in I_L1: 1298 (a, b), indices, h = I 1299 1300 if a == b: 1301 if a == x: 1302 I_L1_L.append(I) 1303 I_L1_R.append(I) 1304 elif a < x: 1305 I_L1_L.append(I) 1306 else: 1307 I_L1_R.append(I) 1308 else: 1309 if b <= x: 1310 I_L1_L.append(I) 1311 elif a >= x: 1312 I_L1_R.append(I) 1313 else: 1314 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=x, fast=True) 1315 1316 if b <= x: 1317 I_L1_L.append(((a, b), indices, h)) 1318 if a >= x: 1319 I_L1_R.append(((a, b), indices, h)) 1320 1321 for I in I_L3: 1322 (b, a), indices, h = I 1323 1324 if a == b: 1325 if a == x: 1326 I_L3_L.append(I) 1327 I_L3_R.append(I) 1328 elif a < x: 1329 I_L3_L.append(I) 1330 else: 1331 I_L3_R.append(I) 1332 else: 1333 if b <= x: 1334 I_L3_L.append(I) 1335 elif a >= x: 1336 I_L3_R.append(I) 1337 else: 1338 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=x, fast=True) 1339 1340 if b <= x: 1341 I_L3_L.append(((b, a), indices, h)) 1342 if a >= x: 1343 I_L3_R.append(((b, a), indices, h)) 1344 1345 Q_L1_L = _intervals_to_quadrants(I_L1_L, f1L1F, f2L1F, u, x, F) 1346 Q_L2_L = _intervals_to_quadrants(I_L2_L, f1V, f2V, v, t, F) 1347 Q_L3_L = _intervals_to_quadrants(I_L3_L, f1L3F, f2L3F, x, u, F) 1348 Q_L4_L = Q_L4 1349 1350 Q_L1_R = _intervals_to_quadrants(I_L1_R, f1L1F, f2L1F, x, s, F) 1351 Q_L2_R = Q_L2 1352 Q_L3_R = _intervals_to_quadrants(I_L3_R, f1L3F, f2L3F, s, x, F) 1353 Q_L4_R = _intervals_to_quadrants(I_L4_R, f1V, f2V, t, v, F) 1354 1355 T_L = _traverse_quadrants(Q_L1_L, Q_L2_L, Q_L3_L, Q_L4_L, exclude=True) 1356 T_R = _traverse_quadrants(Q_L1_R, Q_L2_R, Q_L3_R, Q_L4_R, exclude=True) 1357 1358 N_L = _winding_number(T_L, F) 1359 N_R = _winding_number(T_R, F) 1360 1361 I_L = (I_L1_L, I_L2_L, I_L3_L, I_L4_L) 1362 Q_L = (Q_L1_L, Q_L2_L, Q_L3_L, Q_L4_L) 1363 1364 I_R = (I_L1_R, I_L2_R, I_L3_R, I_L4_R) 1365 Q_R = (Q_L1_R, Q_L2_R, Q_L3_R, Q_L4_R) 1366 1367 F1_L = (f1L1F, f1V, f1L3F, f1L4F) 1368 F2_L = (f2L1F, f2V, f2L3F, f2L4F) 1369 1370 F1_R = (f1L1F, f1L2F, f1L3F, f1V) 1371 F2_R = (f2L1F, f2L2F, f2L3F, f2V) 1372 1373 a, b = (u, v), (x, t) 1374 c, d = (x, v), (s, t) 1375 1376 D_L = (N_L, a, b, I_L, Q_L, F1_L, F2_L) 1377 D_R = (N_R, c, d, I_R, Q_R, F1_R, F2_R) 1378 1379 return D_L, D_R 1380 1381def _horizontal_bisection(N, a, b, I, Q, F1, F2, f1, f2, F): 1382 """Horizontal bisection step in Collins-Krandick root isolation algorithm. """ 1383 (u, v), (s, t) = a, b 1384 1385 I_L1, I_L2, I_L3, I_L4 = I 1386 Q_L1, Q_L2, Q_L3, Q_L4 = Q 1387 1388 f1L1F, f1L2F, f1L3F, f1L4F = F1 1389 f2L1F, f2L2F, f2L3F, f2L4F = F2 1390 1391 y = (v + t) / 2 1392 1393 f1H = dmp_eval_in(f1, y, 1, 1, F) 1394 f2H = dmp_eval_in(f2, y, 1, 1, F) 1395 1396 I_H = dup_isolate_real_roots_list([f1H, f2H], F, inf=u, sup=s, fast=True, strict=True, basis=True) 1397 1398 I_L1_B, I_L1_U = I_L1, I_H 1399 I_L2_B, I_L2_U = [], [] 1400 I_L3_B, I_L3_U = _reverse_intervals(I_H), I_L3 1401 I_L4_B, I_L4_U = [], [] 1402 1403 for I in I_L2: 1404 (a, b), indices, h = I 1405 1406 if a == b: 1407 if a == y: 1408 I_L2_B.append(I) 1409 I_L2_U.append(I) 1410 elif a < y: 1411 I_L2_B.append(I) 1412 else: 1413 I_L2_U.append(I) 1414 else: 1415 if b <= y: 1416 I_L2_B.append(I) 1417 elif a >= y: 1418 I_L2_U.append(I) 1419 else: 1420 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=y, fast=True) 1421 1422 if b <= y: 1423 I_L2_B.append(((a, b), indices, h)) 1424 if a >= y: 1425 I_L2_U.append(((a, b), indices, h)) 1426 1427 for I in I_L4: 1428 (b, a), indices, h = I 1429 1430 if a == b: 1431 if a == y: 1432 I_L4_B.append(I) 1433 I_L4_U.append(I) 1434 elif a < y: 1435 I_L4_B.append(I) 1436 else: 1437 I_L4_U.append(I) 1438 else: 1439 if b <= y: 1440 I_L4_B.append(I) 1441 elif a >= y: 1442 I_L4_U.append(I) 1443 else: 1444 a, b = dup_refine_real_root(h, a, b, F.get_ring(), disjoint=y, fast=True) 1445 1446 if b <= y: 1447 I_L4_B.append(((b, a), indices, h)) 1448 if a >= y: 1449 I_L4_U.append(((b, a), indices, h)) 1450 1451 Q_L1_B = Q_L1 1452 Q_L2_B = _intervals_to_quadrants(I_L2_B, f1L2F, f2L2F, v, y, F) 1453 Q_L3_B = _intervals_to_quadrants(I_L3_B, f1H, f2H, s, u, F) 1454 Q_L4_B = _intervals_to_quadrants(I_L4_B, f1L4F, f2L4F, y, v, F) 1455 1456 Q_L1_U = _intervals_to_quadrants(I_L1_U, f1H, f2H, u, s, F) 1457 Q_L2_U = _intervals_to_quadrants(I_L2_U, f1L2F, f2L2F, y, t, F) 1458 Q_L3_U = Q_L3 1459 Q_L4_U = _intervals_to_quadrants(I_L4_U, f1L4F, f2L4F, t, y, F) 1460 1461 T_B = _traverse_quadrants(Q_L1_B, Q_L2_B, Q_L3_B, Q_L4_B, exclude=True) 1462 T_U = _traverse_quadrants(Q_L1_U, Q_L2_U, Q_L3_U, Q_L4_U, exclude=True) 1463 1464 N_B = _winding_number(T_B, F) 1465 N_U = _winding_number(T_U, F) 1466 1467 I_B = (I_L1_B, I_L2_B, I_L3_B, I_L4_B) 1468 Q_B = (Q_L1_B, Q_L2_B, Q_L3_B, Q_L4_B) 1469 1470 I_U = (I_L1_U, I_L2_U, I_L3_U, I_L4_U) 1471 Q_U = (Q_L1_U, Q_L2_U, Q_L3_U, Q_L4_U) 1472 1473 F1_B = (f1L1F, f1L2F, f1H, f1L4F) 1474 F2_B = (f2L1F, f2L2F, f2H, f2L4F) 1475 1476 F1_U = (f1H, f1L2F, f1L3F, f1L4F) 1477 F2_U = (f2H, f2L2F, f2L3F, f2L4F) 1478 1479 a, b = (u, v), (s, y) 1480 c, d = (u, y), (s, t) 1481 1482 D_B = (N_B, a, b, I_B, Q_B, F1_B, F2_B) 1483 D_U = (N_U, c, d, I_U, Q_U, F1_U, F2_U) 1484 1485 return D_B, D_U 1486 1487def _depth_first_select(rectangles): 1488 """Find a rectangle of minimum area for bisection. """ 1489 min_area, j = None, None 1490 1491 for i, (_, (u, v), (s, t), _, _, _, _) in enumerate(rectangles): 1492 area = (s - u)*(t - v) 1493 1494 if min_area is None or area < min_area: 1495 min_area, j = area, i 1496 1497 return rectangles.pop(j) 1498 1499def _rectangle_small_p(a, b, eps): 1500 """Return ``True`` if the given rectangle is small enough. """ 1501 (u, v), (s, t) = a, b 1502 1503 if eps is not None: 1504 return s - u < eps and t - v < eps 1505 else: 1506 return True 1507 1508def dup_isolate_complex_roots_sqf(f, K, eps=None, inf=None, sup=None, blackbox=False): 1509 """Isolate complex roots of a square-free polynomial using Collins-Krandick algorithm. """ 1510 if not K.is_ZZ and not K.is_QQ: 1511 raise DomainError("isolation of complex roots is not supported over %s" % K) 1512 1513 if dup_degree(f) <= 0: 1514 return [] 1515 1516 if K.is_ZZ: 1517 F = K.get_field() 1518 else: 1519 F = K 1520 1521 f = dup_convert(f, K, F) 1522 1523 lc = abs(dup_LC(f, F)) 1524 B = 2*max([ F.quo(abs(c), lc) for c in f ]) 1525 1526 (u, v), (s, t) = (-B, F.zero), (B, B) 1527 1528 if inf is not None: 1529 u = inf 1530 1531 if sup is not None: 1532 s = sup 1533 1534 if v < 0 or t <= v or s <= u: 1535 raise ValueError("not a valid complex isolation rectangle") 1536 1537 f1, f2 = dup_real_imag(f, F) 1538 1539 f1L1 = dmp_eval_in(f1, v, 1, 1, F) 1540 f2L1 = dmp_eval_in(f2, v, 1, 1, F) 1541 1542 f1L2 = dmp_eval_in(f1, s, 0, 1, F) 1543 f2L2 = dmp_eval_in(f2, s, 0, 1, F) 1544 1545 f1L3 = dmp_eval_in(f1, t, 1, 1, F) 1546 f2L3 = dmp_eval_in(f2, t, 1, 1, F) 1547 1548 f1L4 = dmp_eval_in(f1, u, 0, 1, F) 1549 f2L4 = dmp_eval_in(f2, u, 0, 1, F) 1550 1551 S_L1 = [f1L1, f2L1] 1552 S_L2 = [f1L2, f2L2] 1553 S_L3 = [f1L3, f2L3] 1554 S_L4 = [f1L4, f2L4] 1555 1556 I_L1 = dup_isolate_real_roots_list(S_L1, F, inf=u, sup=s, fast=True, strict=True, basis=True) 1557 I_L2 = dup_isolate_real_roots_list(S_L2, F, inf=v, sup=t, fast=True, strict=True, basis=True) 1558 I_L3 = dup_isolate_real_roots_list(S_L3, F, inf=u, sup=s, fast=True, strict=True, basis=True) 1559 I_L4 = dup_isolate_real_roots_list(S_L4, F, inf=v, sup=t, fast=True, strict=True, basis=True) 1560 1561 I_L3 = _reverse_intervals(I_L3) 1562 I_L4 = _reverse_intervals(I_L4) 1563 1564 Q_L1 = _intervals_to_quadrants(I_L1, f1L1, f2L1, u, s, F) 1565 Q_L2 = _intervals_to_quadrants(I_L2, f1L2, f2L2, v, t, F) 1566 Q_L3 = _intervals_to_quadrants(I_L3, f1L3, f2L3, s, u, F) 1567 Q_L4 = _intervals_to_quadrants(I_L4, f1L4, f2L4, t, v, F) 1568 1569 T = _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4) 1570 N = _winding_number(T, F) 1571 1572 if not N: 1573 return [] 1574 1575 I = (I_L1, I_L2, I_L3, I_L4) 1576 Q = (Q_L1, Q_L2, Q_L3, Q_L4) 1577 1578 F1 = (f1L1, f1L2, f1L3, f1L4) 1579 F2 = (f2L1, f2L2, f2L3, f2L4) 1580 1581 rectangles, roots = [(N, (u, v), (s, t), I, Q, F1, F2)], [] 1582 1583 while rectangles: 1584 N, (u, v), (s, t), I, Q, F1, F2 = _depth_first_select(rectangles) 1585 1586 if s - u > t - v: 1587 D_L, D_R = _vertical_bisection(N, (u, v), (s, t), I, Q, F1, F2, f1, f2, F) 1588 1589 N_L, a, b, I_L, Q_L, F1_L, F2_L = D_L 1590 N_R, c, d, I_R, Q_R, F1_R, F2_R = D_R 1591 1592 if N_L >= 1: 1593 if N_L == 1 and _rectangle_small_p(a, b, eps): 1594 roots.append(ComplexInterval(a, b, I_L, Q_L, F1_L, F2_L, f1, f2, F)) 1595 else: 1596 rectangles.append(D_L) 1597 1598 if N_R >= 1: 1599 if N_R == 1 and _rectangle_small_p(c, d, eps): 1600 roots.append(ComplexInterval(c, d, I_R, Q_R, F1_R, F2_R, f1, f2, F)) 1601 else: 1602 rectangles.append(D_R) 1603 else: 1604 D_B, D_U = _horizontal_bisection(N, (u, v), (s, t), I, Q, F1, F2, f1, f2, F) 1605 1606 N_B, a, b, I_B, Q_B, F1_B, F2_B = D_B 1607 N_U, c, d, I_U, Q_U, F1_U, F2_U = D_U 1608 1609 if N_B >= 1: 1610 if N_B == 1 and _rectangle_small_p(a, b, eps): 1611 roots.append(ComplexInterval( 1612 a, b, I_B, Q_B, F1_B, F2_B, f1, f2, F)) 1613 else: 1614 rectangles.append(D_B) 1615 1616 if N_U >= 1: 1617 if N_U == 1 and _rectangle_small_p(c, d, eps): 1618 roots.append(ComplexInterval( 1619 c, d, I_U, Q_U, F1_U, F2_U, f1, f2, F)) 1620 else: 1621 rectangles.append(D_U) 1622 1623 _roots, roots = sorted(roots, key=lambda r: (r.ax, r.ay)), [] 1624 1625 for root in _roots: 1626 roots.extend([root.conjugate(), root]) 1627 1628 if blackbox: 1629 return roots 1630 else: 1631 return [ r.as_tuple() for r in roots ] 1632 1633def dup_isolate_all_roots_sqf(f, K, eps=None, inf=None, sup=None, fast=False, blackbox=False): 1634 """Isolate real and complex roots of a square-free polynomial ``f``. """ 1635 return ( 1636 dup_isolate_real_roots_sqf( f, K, eps=eps, inf=inf, sup=sup, fast=fast, blackbox=blackbox), 1637 dup_isolate_complex_roots_sqf(f, K, eps=eps, inf=inf, sup=sup, blackbox=blackbox)) 1638 1639def dup_isolate_all_roots(f, K, eps=None, inf=None, sup=None, fast=False): 1640 """Isolate real and complex roots of a non-square-free polynomial ``f``. """ 1641 if not K.is_ZZ and not K.is_QQ: 1642 raise DomainError("isolation of real and complex roots is not supported over %s" % K) 1643 1644 _, factors = dup_sqf_list(f, K) 1645 1646 if len(factors) == 1: 1647 ((f, k),) = factors 1648 1649 real_part, complex_part = dup_isolate_all_roots_sqf( 1650 f, K, eps=eps, inf=inf, sup=sup, fast=fast) 1651 1652 real_part = [ ((a, b), k) for (a, b) in real_part ] 1653 complex_part = [ ((a, b), k) for (a, b) in complex_part ] 1654 1655 return real_part, complex_part 1656 else: 1657 raise NotImplementedError( "only trivial square-free polynomials are supported") 1658 1659class RealInterval: 1660 """A fully qualified representation of a real isolation interval. """ 1661 1662 def __init__(self, data, f, dom): 1663 """Initialize new real interval with complete information. """ 1664 if len(data) == 2: 1665 s, t = data 1666 1667 self.neg = False 1668 1669 if s < 0: 1670 if t <= 0: 1671 f, s, t, self.neg = dup_mirror(f, dom), -t, -s, True 1672 else: 1673 raise ValueError("can't refine a real root in (%s, %s)" % (s, t)) 1674 1675 a, b, c, d = _mobius_from_interval((s, t), dom.get_field()) 1676 1677 f = dup_transform(f, dup_strip([a, b]), 1678 dup_strip([c, d]), dom) 1679 1680 self.mobius = a, b, c, d 1681 else: 1682 self.mobius = data[:-1] 1683 self.neg = data[-1] 1684 1685 self.f, self.dom = f, dom 1686 1687 @property 1688 def func(self): 1689 return RealInterval 1690 1691 @property 1692 def args(self): 1693 i = self 1694 return (i.mobius + (i.neg,), i.f, i.dom) 1695 1696 def __eq__(self, other): 1697 if type(other) != type(self): 1698 return False 1699 return self.args == other.args 1700 1701 @property 1702 def a(self): 1703 """Return the position of the left end. """ 1704 field = self.dom.get_field() 1705 a, b, c, d = self.mobius 1706 1707 if not self.neg: 1708 if a*d < b*c: 1709 return field(a, c) 1710 return field(b, d) 1711 else: 1712 if a*d > b*c: 1713 return -field(a, c) 1714 return -field(b, d) 1715 1716 @property 1717 def b(self): 1718 """Return the position of the right end. """ 1719 was = self.neg 1720 self.neg = not was 1721 rv = -self.a 1722 self.neg = was 1723 return rv 1724 1725 @property 1726 def dx(self): 1727 """Return width of the real isolating interval. """ 1728 return self.b - self.a 1729 1730 @property 1731 def center(self): 1732 """Return the center of the real isolating interval. """ 1733 return (self.a + self.b)/2 1734 1735 def as_tuple(self): 1736 """Return tuple representation of real isolating interval. """ 1737 return (self.a, self.b) 1738 1739 def __repr__(self): 1740 return "(%s, %s)" % (self.a, self.b) 1741 1742 def is_disjoint(self, other): 1743 """Return ``True`` if two isolation intervals are disjoint. """ 1744 if isinstance(other, RealInterval): 1745 return (self.b < other.a or other.b < self.a) 1746 assert isinstance(other, ComplexInterval) 1747 return (self.b < other.ax or other.bx < self.a 1748 or other.ay*other.by > 0) 1749 1750 def _inner_refine(self): 1751 """Internal one step real root refinement procedure. """ 1752 if self.mobius is None: 1753 return self 1754 1755 f, mobius = dup_inner_refine_real_root( 1756 self.f, self.mobius, self.dom, steps=1, mobius=True) 1757 1758 return RealInterval(mobius + (self.neg,), f, self.dom) 1759 1760 def refine_disjoint(self, other): 1761 """Refine an isolating interval until it is disjoint with another one. """ 1762 expr = self 1763 while not expr.is_disjoint(other): 1764 expr, other = expr._inner_refine(), other._inner_refine() 1765 1766 return expr, other 1767 1768 def refine_size(self, dx): 1769 """Refine an isolating interval until it is of sufficiently small size. """ 1770 expr = self 1771 while not (expr.dx < dx): 1772 expr = expr._inner_refine() 1773 1774 return expr 1775 1776 def refine_step(self, steps=1): 1777 """Perform several steps of real root refinement algorithm. """ 1778 expr = self 1779 for _ in range(steps): 1780 expr = expr._inner_refine() 1781 1782 return expr 1783 1784 def refine(self): 1785 """Perform one step of real root refinement algorithm. """ 1786 return self._inner_refine() 1787 1788 1789class ComplexInterval: 1790 """A fully qualified representation of a complex isolation interval. 1791 The printed form is shown as (ax, bx) x (ay, by) where (ax, ay) 1792 and (bx, by) are the coordinates of the southwest and northeast 1793 corners of the interval's rectangle, respectively. 1794 1795 Examples 1796 ======== 1797 1798 >>> from sympy import CRootOf, S 1799 >>> from sympy.abc import x 1800 >>> CRootOf.clear_cache() # for doctest reproducibility 1801 >>> root = CRootOf(x**10 - 2*x + 3, 9) 1802 >>> i = root._get_interval(); i 1803 (3/64, 3/32) x (9/8, 75/64) 1804 1805 The real part of the root lies within the range [0, 3/4] while 1806 the imaginary part lies within the range [9/8, 3/2]: 1807 1808 >>> root.n(3) 1809 0.0766 + 1.14*I 1810 1811 The width of the ranges in the x and y directions on the complex 1812 plane are: 1813 1814 >>> i.dx, i.dy 1815 (3/64, 3/64) 1816 1817 The center of the range is 1818 1819 >>> i.center 1820 (9/128, 147/128) 1821 1822 The northeast coordinate of the rectangle bounding the root in the 1823 complex plane is given by attribute b and the x and y components 1824 are accessed by bx and by: 1825 1826 >>> i.b, i.bx, i.by 1827 ((3/32, 75/64), 3/32, 75/64) 1828 1829 The southwest coordinate is similarly given by i.a 1830 1831 >>> i.a, i.ax, i.ay 1832 ((3/64, 9/8), 3/64, 9/8) 1833 1834 Although the interval prints to show only the real and imaginary 1835 range of the root, all the information of the underlying root 1836 is contained as properties of the interval. 1837 1838 For example, an interval with a nonpositive imaginary range is 1839 considered to be the conjugate. Since the y values of y are in the 1840 range [0, 1/4] it is not the conjugate: 1841 1842 >>> i.conj 1843 False 1844 1845 The conjugate's interval is 1846 1847 >>> ic = i.conjugate(); ic 1848 (3/64, 3/32) x (-75/64, -9/8) 1849 1850 NOTE: the values printed still represent the x and y range 1851 in which the root -- conjugate, in this case -- is located, 1852 but the underlying a and b values of a root and its conjugate 1853 are the same: 1854 1855 >>> assert i.a == ic.a and i.b == ic.b 1856 1857 What changes are the reported coordinates of the bounding rectangle: 1858 1859 >>> (i.ax, i.ay), (i.bx, i.by) 1860 ((3/64, 9/8), (3/32, 75/64)) 1861 >>> (ic.ax, ic.ay), (ic.bx, ic.by) 1862 ((3/64, -75/64), (3/32, -9/8)) 1863 1864 The interval can be refined once: 1865 1866 >>> i # for reference, this is the current interval 1867 (3/64, 3/32) x (9/8, 75/64) 1868 1869 >>> i.refine() 1870 (3/64, 3/32) x (9/8, 147/128) 1871 1872 Several refinement steps can be taken: 1873 1874 >>> i.refine_step(2) # 2 steps 1875 (9/128, 3/32) x (9/8, 147/128) 1876 1877 It is also possible to refine to a given tolerance: 1878 1879 >>> tol = min(i.dx, i.dy)/2 1880 >>> i.refine_size(tol) 1881 (9/128, 21/256) x (9/8, 291/256) 1882 1883 A disjoint interval is one whose bounding rectangle does not 1884 overlap with another. An interval, necessarily, is not disjoint with 1885 itself, but any interval is disjoint with a conjugate since the 1886 conjugate rectangle will always be in the lower half of the complex 1887 plane and the non-conjugate in the upper half: 1888 1889 >>> i.is_disjoint(i), i.is_disjoint(i.conjugate()) 1890 (False, True) 1891 1892 The following interval j is not disjoint from i: 1893 1894 >>> close = CRootOf(x**10 - 2*x + 300/S(101), 9) 1895 >>> j = close._get_interval(); j 1896 (75/1616, 75/808) x (225/202, 1875/1616) 1897 >>> i.is_disjoint(j) 1898 False 1899 1900 The two can be made disjoint, however: 1901 1902 >>> newi, newj = i.refine_disjoint(j) 1903 >>> newi 1904 (39/512, 159/2048) x (2325/2048, 4653/4096) 1905 >>> newj 1906 (3975/51712, 2025/25856) x (29325/25856, 117375/103424) 1907 1908 Even though the real ranges overlap, the imaginary do not, so 1909 the roots have been resolved as distinct. Intervals are disjoint 1910 when either the real or imaginary component of the intervals is 1911 distinct. In the case above, the real components have not been 1912 resolved (so we don't know, yet, which root has the smaller real 1913 part) but the imaginary part of ``close`` is larger than ``root``: 1914 1915 >>> close.n(3) 1916 0.0771 + 1.13*I 1917 >>> root.n(3) 1918 0.0766 + 1.14*I 1919 """ 1920 1921 def __init__(self, a, b, I, Q, F1, F2, f1, f2, dom, conj=False): 1922 """Initialize new complex interval with complete information. """ 1923 # a and b are the SW and NE corner of the bounding interval, 1924 # (ax, ay) and (bx, by), respectively, for the NON-CONJUGATE 1925 # root (the one with the positive imaginary part); when working 1926 # with the conjugate, the a and b value are still non-negative 1927 # but the ay, by are reversed and have oppositite sign 1928 self.a, self.b = a, b 1929 self.I, self.Q = I, Q 1930 1931 self.f1, self.F1 = f1, F1 1932 self.f2, self.F2 = f2, F2 1933 1934 self.dom = dom 1935 self.conj = conj 1936 1937 @property 1938 def func(self): 1939 return ComplexInterval 1940 1941 @property 1942 def args(self): 1943 i = self 1944 return (i.a, i.b, i.I, i.Q, i.F1, i.F2, i.f1, i.f2, i.dom, i.conj) 1945 1946 def __eq__(self, other): 1947 if type(other) != type(self): 1948 return False 1949 return self.args == other.args 1950 1951 @property 1952 def ax(self): 1953 """Return ``x`` coordinate of south-western corner. """ 1954 return self.a[0] 1955 1956 @property 1957 def ay(self): 1958 """Return ``y`` coordinate of south-western corner. """ 1959 if not self.conj: 1960 return self.a[1] 1961 else: 1962 return -self.b[1] 1963 1964 @property 1965 def bx(self): 1966 """Return ``x`` coordinate of north-eastern corner. """ 1967 return self.b[0] 1968 1969 @property 1970 def by(self): 1971 """Return ``y`` coordinate of north-eastern corner. """ 1972 if not self.conj: 1973 return self.b[1] 1974 else: 1975 return -self.a[1] 1976 1977 @property 1978 def dx(self): 1979 """Return width of the complex isolating interval. """ 1980 return self.b[0] - self.a[0] 1981 1982 @property 1983 def dy(self): 1984 """Return height of the complex isolating interval. """ 1985 return self.b[1] - self.a[1] 1986 1987 @property 1988 def center(self): 1989 """Return the center of the complex isolating interval. """ 1990 return ((self.ax + self.bx)/2, (self.ay + self.by)/2) 1991 1992 def as_tuple(self): 1993 """Return tuple representation of the complex isolating 1994 interval's SW and NE corners, respectively. """ 1995 return ((self.ax, self.ay), (self.bx, self.by)) 1996 1997 def __repr__(self): 1998 return "(%s, %s) x (%s, %s)" % (self.ax, self.bx, self.ay, self.by) 1999 2000 def conjugate(self): 2001 """This complex interval really is located in lower half-plane. """ 2002 return ComplexInterval(self.a, self.b, self.I, self.Q, 2003 self.F1, self.F2, self.f1, self.f2, self.dom, conj=True) 2004 2005 def is_disjoint(self, other): 2006 """Return ``True`` if two isolation intervals are disjoint. """ 2007 if isinstance(other, RealInterval): 2008 return other.is_disjoint(self) 2009 if self.conj != other.conj: # above and below real axis 2010 return True 2011 re_distinct = (self.bx < other.ax or other.bx < self.ax) 2012 if re_distinct: 2013 return True 2014 im_distinct = (self.by < other.ay or other.by < self.ay) 2015 return im_distinct 2016 2017 def _inner_refine(self): 2018 """Internal one step complex root refinement procedure. """ 2019 (u, v), (s, t) = self.a, self.b 2020 2021 I, Q = self.I, self.Q 2022 2023 f1, F1 = self.f1, self.F1 2024 f2, F2 = self.f2, self.F2 2025 2026 dom = self.dom 2027 2028 if s - u > t - v: 2029 D_L, D_R = _vertical_bisection(1, (u, v), (s, t), I, Q, F1, F2, f1, f2, dom) 2030 2031 if D_L[0] == 1: 2032 _, a, b, I, Q, F1, F2 = D_L 2033 else: 2034 _, a, b, I, Q, F1, F2 = D_R 2035 else: 2036 D_B, D_U = _horizontal_bisection(1, (u, v), (s, t), I, Q, F1, F2, f1, f2, dom) 2037 2038 if D_B[0] == 1: 2039 _, a, b, I, Q, F1, F2 = D_B 2040 else: 2041 _, a, b, I, Q, F1, F2 = D_U 2042 2043 return ComplexInterval(a, b, I, Q, F1, F2, f1, f2, dom, self.conj) 2044 2045 def refine_disjoint(self, other): 2046 """Refine an isolating interval until it is disjoint with another one. """ 2047 expr = self 2048 while not expr.is_disjoint(other): 2049 expr, other = expr._inner_refine(), other._inner_refine() 2050 2051 return expr, other 2052 2053 def refine_size(self, dx, dy=None): 2054 """Refine an isolating interval until it is of sufficiently small size. """ 2055 if dy is None: 2056 dy = dx 2057 expr = self 2058 while not (expr.dx < dx and expr.dy < dy): 2059 expr = expr._inner_refine() 2060 2061 return expr 2062 2063 def refine_step(self, steps=1): 2064 """Perform several steps of complex root refinement algorithm. """ 2065 expr = self 2066 for _ in range(steps): 2067 expr = expr._inner_refine() 2068 2069 return expr 2070 2071 def refine(self): 2072 """Perform one step of complex root refinement algorithm. """ 2073 return self._inner_refine() 2074