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